Modelling non-linear large scale structure using Lagrangian perturbation theory re-expansions

Lagrangian perturbation theory (LPT) has been widely used to model the nonlinear growth of large scale structure analytically. However, the LPT series converges only for a finite time. In recent work, the authors have examined this issue in great detail, and proposed a new algorithm (called LPT re-expansions), to extend the domain of validity of the series. This article outlines the main ideas of the algorithm, first developed for the spherical top-hat system, and later generalized to deal with inhomogeneous initial conditions arising from random fields.


Introduction
We have presented a new algorithm to model large scale structure in the mildly non-linear regime. It is based on Lagrangian perturbation theory (LPT), a technique introduced by Zeldovich in the 1970s, and developed further by various authors. In the Lagrangian formalism, the position of a particle is the fundamental variable, and is expressed as a function of the initial position and time. The density is then reconstructed from the position by imposing mass conservation. This definition involves the Jacobian of the transformation, relating the initial (Lagrangian) coordinate space and the final (Eulerian) coordinate space, and hence, is intrinsically nonlinear. A fundamental assumption in this framework is that the matter is 'cold'. Once particle trajectories cross the initial space to final space, mapping becomes many-to-one and the density formally diverges.
However, LPT has another major drawback. The series is not guaranteed to converge at all times before orbit crossing. In a recent paper [1], this issue has been investigated in detail for the spherical top-hat system. Our analysis has estimated the exact 'time of validity' of the series, and we have proposed an algorithm (dubbed LPT re-expansions) to extend the domain of validity of the series beyond the initial time interval. We have shown that this technique has worked, and is computationally feasible for the top-hat system. Convergence could be achieved by increasing the re-expansion frequency and/or the Lagrangian order.
In a subsequent paper [2], this idea of the algorithm has been applied to inhomogeneous initial conditions arising from random fields. One important additional feature emerges. The solution outlined by the traditional framework in [3] is unique only up to a purely time-dependent transformation. This time dependent function needs to be added to get a complete solution that has the correct convergence properties. This paper summarizes the main results of our work.

Spherical top-hat system
The spherical top-hat perturbation consists of a inner over-dense or under-dense sphere surrounded by a mass compensating region. The perturbation is described by two parameters: fractional over-density δ and fractional Hubble parameter δ v . In a δ − δ v phase-space, the perturbation can be represented as a point with polar coordinates (Δ, θ). Δ is taken as the expansion parameter for the series. The scale factors a(t) and b(t) describe the time evolution of the background and perturbed spheres respectively. The exact evolution for b(t) is known as a function of (Δ, θ). By treating the variable Δ to be complex, one can estimate when the solution b(t) is singular and the roots to the singularity condition give the radius of convergence R(Δ, t). This is represented on a 'root plot'. It is generally different for different values of θ. From the root plot, the time of validity can be read off for the Δ corresponding to the initial value. Figure 1 shows a specific example of the problem.
To remedy this problem, we have proposed the algorithm of LPT re-expansions. The main idea is somewhat analogous to that of 'analytic continuation'. The series is expanded up to some final time, less than the time of validity given by the initial parameters. At that time, the system is re-initialized by defining a new Δ and θ, and the series is propagated forward until a new final time, less than the new time of validity. This is done successively until the desired end point is reached. A detailed examination of the dynamics in the δ − δ v phase-space proves that this procedure, indeed works, and implementation is realistically feasible. Figure 2 shows the implementation of this algorithm on the same problem considered in Figure 1. The upper panel is b(t) vs. a(t) and the lower panels show the log of the error between the various components. It is clear that taking more steps decreases the error. Comparing the second and third plots, it is also seen that higher orders converge faster. Compare, for example, the red line (order 1) with the black line (order 9), i.e., the convergence rates improve as the order of the approximation improves.  Figure 2. The LPT re-expansion algorithm applied to the example in Figure 1. The three plots correspond to one step (left), three steps (middle), five steps (right).

General initial conditions
In the general case, the equation of motion for the particle position r is given by the geodesic equation:r −ä a r = −∇ r ψ(r, t), and (1) The formal expansion of Ehlers and Buchert [3] uses a derivative form: where ψ(r, t) is the gravitational potential, ρ m (r, t) and δρ m (r, t) are the total and perturbed matter densities. Generally, one substitutes in Eqs. (3) and (4) the ansatz r = aX + p(X, t), where X is the Lagrangian coordinate, and solves for the displacement vector p(X, t). The solutions to Eqs. (1) and (2) differ from those to Eqs. (3) and (4) by a time dependent transformation (called frame shift). Satisfying the geodesic equation implies that the total momentum in the box decays as 1/a(t) for the expanding universe. Ignoring the shifts violate this law and break convergence. Figure 3 shows the effect on convergence with (left panels) and without (right panels) the shifts for generic initial conditions (where the velocity can have a non-zero transverse part).

Conclusion
We be fine tuned by changing the order of the series or the frequency of re-expansion. Higher orders converge faster. The algorithm is then generalized to inhomogeneous initial conditions. In this case, the series solution to the differential form of the geodesic equation must be supplemented by time dependent gauge functions called frame shifts. Omission of these functions violates momentum conservation, and breaks the expected convergence.
The new technique has three major advantages: (i) The displacement and hence, the density and velocity are treated as smooth fields, making the results free from particle induced shotnoise, (ii) The errors can be precisely controlled by changing the Lagrangian order and/or the number of steps, and (iii) It is capable of handling arbitrary velocity data including cases with a initial rotational component.
Nothing in this work deals with what happens once orbits cross. Some sort of effective pressure needs to be introduced to handle the onset of orbit crossing.