Accounting for expansion in $N$-Body gravitational simulations

Question

I have an $N$-Body simulation for the evolution of the universe. It uses a hybrid Tree Particle-Mesh scheme to compute the forces. The tree code used Newtons law of gravity:

$$\textbf{F}_i = -G \sum_{j\neq i} \frac{m_i m_j}{|\textbf{x}_i - \textbf{x}_j|^3} (\textbf{x}_i - \textbf{x}_j). $$

The particle mesh method computes force by solving the poisson equation for potential in through a Fourier transform:

$$\nabla^2 \phi = 4 \pi G \rho.$$

The particle mesh method accounts for periodic boundary conditions and the tree code is used only for dense regions. Currently the simulation does not account for the expansion of the universe. Is there any way of adapting the equations or the methods to take expansion into account?

Answer 1

If you treat your $\mathbf{x}_i$ as physical coordinates, then at the level of Newtonian physics, you can match the result of a spacetime expansion in general relativity by giving the particles an initial outward velocity. That is, if the Hubble constant is $H$, and the initial proper velocities were $\mathbf{v}_i$, then the right initial velocities to account for expansion are $$H \mathbf{x}_i + \mathbf{v}_i.$$ You don't have to change anything else.

However, if you're doing the computation using periodic boundary conditions in a box of fixed size, this won't work. In this case, you want your $\mathbf{x}_i$ to represent comoving coordinates (i.e. physical distances are comoving distances multiplied by the scale factor $a$), in which case one can show that the equations of motion are $$\mathbf{a} = - H \mathbf{v} - \frac{1}{a} \nabla \Phi, \quad \frac{1}{a^2}\nabla^2 \Phi = 4 \pi G \rho.$$ In other words, we throw in a factor of $1/a$ every time there's a comoving space derivative, and there's an additional "Hubble friction" term $-H \mathbf{v}$ representing cosmological redshift. The time-dependence of $a$ and $H$ can be computed from the Friedmann equation for a flat, matter-dominated universe, $$\frac{\ddot{a}}{a} = - \frac{4 \pi}{3} G \bar{\rho}.$$ I believe this is the standard way Newtonian cosmological simulations are set up.