# Overview: N-body

We provide some methods for performing some very simple fixed grid N-body simulations (Particle Mesh). The domain decomposition we use in this library is not suitable for huge high-resolution simuations so we only provide these methods and there are some really good and efficient N-body codes like RAMSES and GADGET that one should use for proper N-body anyway.

# Methods

In a N-body simulation we solve the equations $$\frac{d^2 x}{d\tau^2} = - \nabla\Phi$$ $$\nabla^2\Phi = C \delta$$ where, for cosmological simulations, $d\tau = H_0 dt / a^2$ and $C = \frac{3}{2}\Omega_M a$ and for Newtonian ones $C = 4\pi G \rho_{\rm mean} / B_{\rm Boxsize}^2$. This can be solved with the Leapfrog method $$x_{n+1} = x_n + v_{n+1/2}\Delta\tau$$ $$v_{n+1/2} = v_{n-1/2} + (\nabla\Phi)_n\Delta\tau$$ The first step is moving the partices (Drift) and the last step is updating the velocities (Kick). The force is computed by performing ${\rm NDIM}+1$ FFTs: $$\nabla\Phi = C F^{-1}[ -\frac{i\vec{k}}{k^2}F[\delta]]$$ The methods below can perform these steps:

// This is solving the Poisson equation D^2phi = norm_poisson_equation * delta and returns DPhi
template <int N, class T>
void compute_force_from_density_real(
const FFTWGrid<N> & density_grid_real,
std::array<FFTWGrid<N>, N> & force_real,
double norm_poisson_equation = 1.0);

// Same as above, but where the density field is already in fourier space
template <int N, class T>
void compute_force_from_density_fourier(
const FFTWGrid<N> & density_grid_fourier,
std::array<FFTWGrid<N>, N> & force_real,
double norm_poisson_equation = 1.0);

// Advance positions from velocity. Takes care of communication in case particles move past domain boundary
template <int N, class T>
void DriftParticles(
FML::PARTICLE::MPIParticles <T> & part,
double delta_time,
bool periodic_box = true);

// Does *not* takes care of communication in case particles move past domain boundary
template <int N, class T>
void DriftParticles(
T * p,
size_t NumPart,
double delta_time,
bool periodic_box = true);

template <int N, class T>
void KickParticles(std::array<FFTWGrid<N>, N> & force_grid,
FML::PARTICLE::MPIParticles<T> & part,
double delta_time,
std::string interpolation_method);

template <int N, class T>
void KickParticles(std::array<FFTWGrid<N>, N> & force_grid,
T * p,
size_t NumPart,
double delta_time,
std::string interpolation_method);

// Do one kick-drift-kick 2nd order leapfrog N-body step
template <int N, class T>
void KickDriftKickNBodyStep(
int Nmesh,
FML::PARTICLE::MPIParticles<T> & part,
double delta_time,
std::string density_assignment_method,
double norm_poisson_equation)

// Do one 4th order leapfrog N-body step
template <int N, class T>
void YoshidaNBodyStep(
int Nmesh,
FML::PARTICLE::MPIParticles<T> & part,
double delta_time,
std::string density_assignment_method,
double norm_poisson_equation);

// From a given density field delta(k,zini) generate N-body IC
// using 1, 2 or 3 LPT.
template <int N, class T>
void NBodyInitialConditions(FML::PARTICLE::MPIParticles<T> & part,
int Npart_1D,
double buffer_factor,

FFTWGrid<N> & delta_fourier,
int LPT_order,

double box,
double zini,
std::function<double(double)> & H_over_H0_of_loga,
std::vector<std::function<double(double)>> & growth_rate_f_of_loga);

// Generate a random field (gaussian or non-local non-gaussian) of a given type
// and then call the method above
template <int N, class T>
void NBodyInitialConditions(
FML::PARTICLE::MPIParticles<T> & part,
int Npart_1D,
double buffer_factor,

int Nmesh,
bool fix_amplitude,
FML::RANDOM::RandomGenerator * rng,
std::function<double(double)> Pofk_of_kBox_over_Pofk_primordal,
std::function<double(double)> Pofk_of_kBox_over_volume_primordial,
int LPT_order,
std::string type_of_random_field,
double fNL,

double box,
double zini,
std::function<double(double)> & H_over_H0_of_loga,
std::vector<std::function<double(double)>> & growth_rate_f_of_loga);



# COLA simulations

In COLA simulations we solve for displacements on top of that predicted by LPT. With this we can do simulations much faster than usual N-body (by sacrificing accuracy on small scales). Instead of propagating $x$ and $v$ via $$\frac{dx}{d\tau} = v$$ $$\frac{d v}{d\tau} = - \nabla\Phi$$ we instead define $x_{\rm COLA} = x - x_{\rm LPT}$ so that the equation becomes $$\frac{dx}{d\tau} = v_{\rm COLA} + \frac{dx_{\rm LPT}}{d\tau}$$ $$\frac{d v_{\rm COLA}}{d\tau} = - \nabla\Phi - \frac{d^2x_{\rm LPT}}{d\tau^2}$$ and we propagate $x$ and $v_{\rm COLA}$ instead. Note that in this frame the initial velocity is $v_{\rm COLA} = 0$ and this helps us keep in line with linear theory on the largest scales even with large time-steps. Also note that $x_{\rm LPT} \sim D_1(\tau)\Psi_1(q,\tau_{\rm ini}) + D_2(\tau)\Psi_2(q,\tau_{\rm ini}) + \ldots$ so the derivatives of $x_{\rm LPT}$ above just becomes some growth-factor times the initial displacement field at particle positions when we start the simulation. The library is set up such that when we create initial conditions for N-body then we store displacement fields at the particle positions if they have the related get_D_1LPT (and 2LPT, 3LPTa, 3LPTb if you want to use higher order) methods. With this its easy to turn a normal N-body into a COLA simulations - all we need is to compute some ODEs to get the growth-factors and add in the additional terms in the equations of motion (see the N-body example for this).

# Modified gravity simulations

The library contains some methods for computing the additional fifth-force in modified gravity models. We can compute the fifth-force in any linear model exactly, we can estimate the fifth-force in models with a screening mechanism using an approximate method and if you want you can always use the multigrid solver to solve the exact equation.

The methods below compute the fifth-force in terms of a fifth-force "density contrast" that is such that the total force (in fourier space) is just $$F[\nabla \Phi_{\rm tot}] = F[\nabla \Phi + \nabla \Phi_{\rm fifth-force}] \propto \frac{i\vec{k}}{k^2}(\delta(k) + \delta_{\rm fifth-force})$$ i.e. after computing this we only need to add this to the density contrast we would use to compute the usual gravitational force in order to get the total-force.

// For models where the fifth-force is given as F = Dphi where
// F[D^2 phi] = m^2 F[phi] + f(k)F[delta]
// Here the coupling_factor is f(k)k^2/(k^2 + m^2)
// e.g. for linearized f(R) the coupling factor is 1/3 * k^2/(k^2 + m^2)
// and for linearized DGP the coupling factor is 1/(3betaDGP)
template <int N>
void compute_delta_fifth_force(
const FFTWGrid<N> & density_fourier,
FFTWGrid<N> & density_mg_fourier,
std::function<double(double)> coupling_factor_of_kBox);

// Approximate method of Winther and Ferreira 2015
// For MG models that has a potential dependent screening (f(R) and the alike)
// D^2 phi = m^2 phi + f(Phi)delta
// e.g. for f(R) the screening factor is f(Phi) = 1.5 |f_R / Phi|
// and the coupling factor is 1/3 * k^2 / (k^2 + m^2)
// poisson_norm is used to convert delta to Phi and is defined via D^2 Phi = poisson_norm * delta
template <int N>
void compute_delta_fifth_force_potential_screening(
const FFTWGrid<N> & density_fourier,
FFTWGrid<N> & density_mg_fourier,
std::function<double(double)> coupling_factor_of_kBox);
std::function<double(double)> screening_factor_of_newtonian_potential,
double poisson_norm);

// Approximate method of Winther and Ferreira 2015
// For MG models that has a density dependent screening (DGP and the alike)
// D^2 phi = m^2 phi + f(rho)delta
// e.g. for DGP m=0 and the screening factor f(rho) = 2(sqrt(1 + C) - 1)/C where C = const * rho