Back to the main page Back to the documentation 1. Overview: Lagrangian Perturbation Theory (1LPT+2LPT+3LPT) 2. Methods we provide 3. Other LPT related stuff: Reconstruction

This contains methods and algorithms for working with LPT. Can be used for initial conditions generation, COLA simulations, reconstruction and many more things. But if you just want to generate initial conditions for cosmological N-body simulation then you probably want to check out Oliver Hahn's MUSIC code.

The position of particles $x$ are related to the Lagrangian positions $q$ via $$x = q + \Psi(q,t)$$ where $\Psi = \Psi^{\rm 1LPT} + \Psi^{\rm 2LPT} + \ldots$ is the displacement field. We provide functions to computing displacement fields up to third order. Below all positions (and gradients) are in code units, i.e. in $[0,1)$.

Given a density field $\delta$ we compute the 1LPT potential $\phi^{\rm 1LPT}$ defined via $\Psi^{\rm 1LPT} = D_1 \nabla\phi^{\rm 1LPT}$ by solving the Poisson equation $$\nabla^2\phi^{\rm 1LPT} = -\delta$$ Note that we use the sane sign convention $\Psi^{\rm iLPT} = +\nabla\phi^{\rm iLPT}$ here and not the silly one often found in the litterature where $\Psi^{\rm 1LPT} = -\nabla\phi^{\rm 1LPT}$ and $\Psi^{\rm 2LPT} = +\nabla\phi^{\rm 2LPT}$!

Likewise the 2LPT potential is computed by solving $$\nabla^2\phi^{\rm 2LPT} = -\frac{3}{7}\sum_{i>j} [\nabla_{ii}\phi^{\rm 1LPT}\nabla_{jj}\phi^{\rm 1LPT} - (\nabla_{ij}\phi^{\rm 1LPT})^2]$$ where all terms are evaluated via Fourier transforms (where the factor $-3/7 = D_2/D_1^2$).

The full displacement field can then be computed as $$\vec{\Psi} = \frac{D_1}{D_1^{\rm ini}}\nabla\phi^{\rm 1LPT} + \frac{D_2}{D_2^{\rm ini}}\nabla\phi^{\rm 2LPT}$$ and the velocity can be found as $$\frac{d\vec{x}}{d\tau} = BH\left[\frac{d\log D_1}{d\tau}\frac{D_1}{D_1^{\rm ini}}\nabla\phi^{\rm 1LPT} + \frac{d\log D_2}{d\tau}\frac{D_2}{D_2^{\rm ini}}\nabla\phi^{\rm 2LPT}\right]$$ where $B$ is the physial size of the box ($BH$ has units of velocity). Note that $D_1=D_2=1$ at the initial time by definition, however in Einstein deSitter - matter dominated $\Lambda$CDM - the physically relevant solution has $D_2 = -\frac{3}{7} D_1^2$ so we need to add an additional factor of $-\frac{3}{7}$ infront of the 2LPT terms when assigning positions and velocities. We include this factor by default.

Its only in simple models that the growth factor is scale-independent. In general the displacement field do not factor in time-space, but rather in time-fourier-mode and the growth factor will depend on $k$. The methods here can handle this case. For scale-independent growth we only need to compute the displacement field once and then automatically have it for all times. For the case where we have scaledependent growth we need to provide the function(s) $D_i(k,t)/D_i(k,t_{\rm ini})$ and then evaluate it from $$\phi^{\rm iLPT}(k,t) = \frac{D_i(k,t)}{D_i(k,t_{\rm ini})} \phi^{\rm iLPT}(k,t_{\rm ini})$$ The velocity in this case can similary by computed from $$\vec{v} = HB \cdot F^{-1}\left[ \frac{d\log D_1(k,t)}{d\tau} \frac{D_1(k,t)}{D_1(k,t_{\rm ini})}F[\nabla\phi^{\rm 1LPT}] + \frac{d\log D_2(k,t)}{d\tau} \frac{D_2(k,t)}{D_2(k,t_{\rm ini})}F[\nabla\phi^{\rm 2LPT}]\right]$$

The equation for 3LPT is given by $$\Psi^{\rm 3LPT} = \left(\frac{D_1}{D_1^{\rm ini}}\right)^3\left(\nabla \phi^{3a} + \nabla \phi^{3b} + \nabla\times A\right)$$ where $$\nabla^2 \phi^{3a} = -\frac{1}{3}\det[\phi^{(1)}_{ij}]$$ $$\nabla^2 \phi^{3b} = -\frac{5}{9}[\nabla^2\phi^{(1)}\nabla^2\phi^{(2)} - \sum_{i,j} \phi^{(1)}_{ij} \phi^{(2)}_{ij} ]$$ $$\nabla^2 A = \frac{1}{3}\nabla\phi^{(2)}_i\times \nabla\phi^{(1)}_i$$ We normalize the potentials we provide such that $$\Psi^{\rm 3LPT} = \left(\frac{D_1}{D_1^{\rm ini}}\right)^3\left(\nabla \phi^{3a} + \nabla \phi^{3b} + \nabla\times A\right)$$ The main difference wrt 1LPT and 2LPT is that the displacement field is not a pure potential flow (it has non-zero curl - well its still curl free in Eulerian coordinates, but the transformation to Lagrangian coordinates introduces a curl. This curl part does not contribute to the initial power-spectrum of $\Psi$) so the only method we provide for 3LPT is to compute the full displacment vector $\Psi$ with all 3 orders included. The only thing missing from our implementation is dealiasing of the 3LPT term (see Michaux et al. 2020 for why its needed) and proper *testing* which I will do one day.

See LPT/example for an example on how to compute displacement fields from a given power-spectrum and then make particles from it (IC generation).

// From a given phi^iLPT(k, tini) compute the vector Dphi^iLPT(x,t)
// The vector is allocated in the method
// The factor DoverDini = D_nLPT(t) / D_nLPT(tini) where tini is the redshift we have phi_fourier at
template<int N>
  void from_LPT_potential_to_displacement_vector(
      const FFTWGrid<N> & phi_fourier,
      std::vector<FFTWGrid<N> > & psi_real,
      double DoverDini = 1.0);

// Compute phi^1LPT(k,tini) in Fourier space
// Where tini is the time we have delta at
// The phi grid is allocated in the method
template<int N>
  void compute_1LPT_potential_fourier(
      const FFTWGrid<N> & delta_fourier,
      FFTWGrid<N> & phi_1LPT_fourier);

// Compute phi^2LPT(k,tini) in Fourier space
// Where tini is the time we have delta at
// The phi grid is allocated in the method
template<int N>
  void compute_2LPT_potential_fourier(
      const FFTWGrid<N> & delta_fourier,
      FFTWGrid<N> & phi_2LPT_fourier);

// From phi^iLPT(k,t_ini) to Psi(x,t) using growth-factors: Di(k,t)/Di(k, t_ini)
// The vector is allocated in the method
template<int N>
  void from_LPT_potential_to_displacement_vector_scaledependent(
      const FFTWGrid<N> & phi_fourier,
      std::vector<FFTWGrid<N> > & psi_real,
      std::function<double(double)> & growth_function_ratio);

// Compute 1LPT, 2LPT and 3LPT potentials at (k,tini) at the same time
// Where tini is the time we have delta at
// With ignore_curl_term we don't compute the vector potential A
template<int N>
  void compute_3LPT_potential_fourier(
      const FFTWGrid<N> & delta_fourier,
      FFTWGrid<N> & phi_1LPT_fourier,
      FFTWGrid<N> & phi_2LPT_fourier,
      FFTWGrid<N> & phi_3LPT_a_fourier,
      FFTWGrid<N> & phi_3LPT_b_fourier,
      std::vector<FFTWGrid<N>> & phi_3LPT_Avec_fourier,
      bool ignore_curl_term);

// This method is not really needed as it follows from the above, but anyway
// Compute the displacement field Psi(k,tini) and its derivative up to 3LPT
// Where tini is the time we have delta at
// dlogDdt is the derivative of the (log) growth factor at the initial time
// With ignore_curl_term we don't compute the vector potential A
template<int N>
  void compute_3LPT_displacement_field(
      const FFTWGrid<N> & delta_fourier,
      std::vector<FFTWGrid<N>> & Psi,
      std::vector<FFTWGrid<N>> & dPsidt,
      double dlogDdt,
      bool ignore_curl_term)

// Compute augmented LPT (ALPT) Psi(x,t) potential (Kitaura and Hess 2013)
// DoverDini_1LPT = D(t)/D(tini) where tini is the time we have delta at
template<int N>
  void compute_ALPT_potential_fourier(
      const FFTWGrid<N> & delta_fourier,
      const FFTWGrid<N> & phi_2LPT_fourier,
      FFTWGrid<N> & phi_ALPT_fourier,
      double smoothing_scale,
      std::string smoothing_method,
      double DoverDini_1LPT = 1.0,
      double DoverDini_2LPT = 1.0);

These methods are templated on the dimension $N$ because why not (and N=2 is useful for testing).

Growth factors are easily evaluated using the ODESolver and Spline and for $\Lambda$CDM one can call the method below to generate splines of the growth-factors up to 3LPT (not included the vector potential at 3LPT yet). The growth-factors are normalized to unity at redshift $z_{\rm ini}$. The IC as set as to match the EdS solution so $D_1=1$, $D_2 = -3/7$ and so on.

void compute_LPT_growth_factors_LCDM(
      double OmegaM,
      double zini,
      std::function<double(double)> HoverH0_of_a,
      FML::INTERPOLATION::SPLINE::Spline & D_1LPT_of_loga,
      FML::INTERPOLATION::SPLINE::Spline & D_2LPT_of_loga,
      FML::INTERPOLATION::SPLINE::Spline & D_3LPTa_of_loga,
      FML::INTERPOLATION::SPLINE::Spline & D_3LPTb_of_loga);

Defining $\Theta$ via $q = \nabla_x \Theta(x)$ then $$\det[\Theta_{,ij}] = 1 + \delta(x)$$ this equation reduces to a qubic equation in the Laplacian $$(\nabla^2\Theta)^3 + p \nabla^2\Theta = q$$ where $p,q$ depends on cross-derivatives and $\delta$. Thus without expanding in a perturbative series we can solve this equation exactly using the multigrid solver, see examples for how to do this. This can again be used to find the displacement field $\Psi(q) = q-x$ by the following method: Compute $\Theta(x)$ in a $x$-box and create 'particles' that have both $x$ and $q$ coordinates. Swap $q$ and $x$ so that we now have a $q$-box and communicate particles so that they end up at their $q$ positions. The particles are now in a $q$-box and we can bin $q-x = \Psi(q)$ to a grid to get the displacement field. This has many applications like for example initial density reconstruction. Though it does not give us time-derivatives.

We have methods for solving the usual 1LPT reconstruction equation (i.e. given observed galaxies try to move them back to the initial density field or to remove RSD). From observed $x$ compute (and smooth) $\delta_{\rm obs}$ and solve the LPT equation $$\nabla\cdot \Psi + \frac{f}{b}\nabla[(\Psi\cdot r)r] = -\frac{\delta_{\rm obs}}{b}$$ where $r$ is the unit norm line of sight vector, $f$ is the growth rate and $b$ the galaxy bias. Given $\Psi$ subtract $\Psi$ from $x$ to get the $q$'s, or subtract $-\frac{f}{b}(\Psi\cdot r)r$ to remove the RSD contribution. This can be solved with the multigrid solver (see examples) or with Fourier methods (see Reconstruction.h; thought this is just approximate as the term $(\Psi\cdot r)r$ is not curl-free, but can be eliviated by doing it iteratively)

template <int N, class T>
  void RSDReconstructionFourierMethod(
      MPIParticles<T> & part,
      std::string density_assignment_method,
      std::vector<double> los_direction,
      int Nmesh,
      int niterations,
      double beta,
      std::pair<std::string, double> smoothing_options,
      bool survey_data);