Lagrangian cosmological perturbation theory at shell-crossing
Perturbation theory is often used to describe gravitational dynamics of the large scale structure of the Universe in the weakly nonlinear regime and is valid in the single flow regime. However little is known on how the perturbative series expansion performs when the system reaches a significant level of nonlinearity, in particular when approaching shell crossing. To study this issue, we consider the growth of primordial dark matter halos seeded by three crossed initial sine waves of various amplitudes. Using a Lagrangian treatment of cosmological gravitational dynamics, we examine the convergence properties of a high-order perturbative expansion in the vicinity of shell-crossing by comparing the analytical results with state-of-the-art Vlasov-Poisson simulations. In agreement with intuition, we show that convergence of the perturbative series is the best in quasi one-dimensional shell-crossing (one sine-wave dominating) and slows down drastically when approaching triaxial symmetry (three sine-waves with same amplitude). In all cases, however, the perturbative series exhibit a generic behavior as a function of the order of the expansion, which can be fitted with a simple analytic form, allowing us to formally extrapolate the expansion to infinite order. The results of such an extrapolation agree remarkably well with the simulations, even at shell-crossing.
Introduction.—The observed large-scale structures of the Universe are thought to be mainly the result of the evolution through gravitational instability of small initial fluctuations in the matter distribution, with a dominant component given by collisionless dark matter. In the concordance model, Cold Dark Matter (CDM) Peebles (1982, 1984); Blumenthal et al. (1984), the dark matter particles form, at the macroscopic level, a smooth distribution with a virtually null initial local velocity dispersion. This means that dark matter dynamics follow Vlasov-Poisson equations and that the dark matter distribution can be represented as a three-dimensional sheet evolving in six-dimensional phase-space.
In the standard CDM picture, the phase-space sheet is initially perturbed by small Gaussian random fluctuations in the density distribution, with a cut-off scale related to the mass of the dark matter particle Peebles (1982); Berezinsky et al. (2003). The first nonlinear structures form at shell crossing, i.e. in regions of space where the phase-space sheet first self-intersects. Primordial dark matter halos with a power-law density profile of logarithmic slope around Diemand et al. (2005); Ishiyama (2014); Angulo et al. (2017); Ogiya and Hahn (2018) grow around these initial singularities through violent relaxation and then merge together hierarchically to form larger halos with universal but different properties Navarro et al. (1996, 1997, 2010); Diemand and Moore (2011); Fukushige and Makino (1997). Due to the complexity of post-collapse dynamics, the actual processes leading to the shape of a primordial or an evolved halo remain a subject of debate. Indeed, there is no exact analytical theory to predict the results of numerical simulations, although many approaches have been proposed to tackle the problem, relying e.g. on self-similarity Fillmore and Goldreich (1984); Bertschinger (1985); Alard (2013) or on entropy maximization Lynden-Bell (1967); Hjorth and Williams (2010); Carron and Szapudi (2013); Pontzen and Governato (2013). In this debate, the detailed knowledge of the structure of primordial dark matter halos prior to shell-crossing seems essential and even this remains a challenge in the general case.
To describe gravitational dynamics before shell-crossing, it is however possible to employ perturbation theory (PT) as long as fluctuations in the density field remain small Bernardeau et al. (2002). With Lagrangian perturbation theory (LPT) Zel’dovich (1970); Shandarin and Zeldovich (1989); Bouchet et al. (1992); Buchert (1992); Buchert and Ehlers (1993); Bouchet et al. (1995), which uses the displacement field as a small parameter in the expansion of the equations of motion, one can follow in a rather realistic way the evolution of a system in the nonlinear regime up to shell-crossing or even slightly beyond. Zel’dovich approximation Zel’dovich (1970), which corresponds to linear order, has been widely used in the literature. It already gives the exact solution until shell-crossing in the one-dimensional case Novikov (1969) and also provides, in higher dimension, a firm framework to study the families of singularities that first form at shell crossing Arnold et al. (1982); Hidding et al. (2014); Feldbrugge et al. (2017). In general, higher order PT is required to have a sufficiently accurate description of pre-collapse dynamics Moutarde et al. (1991); Buchert et al. (1997) and this obviously depends on the nature of initial conditions. While convergence of the perturbative series has been proven for LPT Zheligovsky and Frisch (2014); Rampf et al. (2015), knowing the speed of such convergence with the order of the expansion remains an important question.
In this Letter, we study, in the framework of LPT, the phase-space structure at shell-crossing of primordial halos initially seeded by three crossed sine waves of various amplitudes. We perform LPT expansion up to tenth order, examine its convergence properties and perform extrapolation to infinite order. The analytical results are validated with state-of-the-art Vlasov-Poisson simulations.
Setup.—In the presence of gravity, the Lagrangian equation of motion of a matter element in the expanding Universe is given by
where is the comoving position, the Newton potential obtained by solving Poisson equation,
and where the quantities , , , and correspond to the scale factor of the Universe, Hubble parameter, background mass density and matter density contrast, respectively. The symbol corresponds to Lagrangian time derivative. In this framework, the velocity of each mass element is given by .
Here, we use a Lagrangian approach i.e. we follow motion as a function of initial position, the Lagrangian coordinate . The subsequent evolution of the position is expressed as , where represents the displacement from initial position. Then, the velocity field is expressed as , and, in the single flow regime (prior to shell crossing), mass conservation implies , which leads to .
We consider the dynamics of mass flow in a finite box with periodic boundaries in the interval and start with initial conditions defined by the following three crossed sine waves:
where is the linear growth factor. The initial time, , and parameters, , , are chosen such that and a small density peak appears at or . This density peak will grow with time due to gravitational instability and the way this happens depends on the relative amplitudes of parameters . Here, we are going to study analytically and numerically the evolution of such a system for three configurations: quasi one dimensional, denoted by Q1D-S, with , triaxial asymmetric with , denoted by ASY-S and triaxial symmetric, i.e. , denoted by SYM-S. When passing from the quasi one-dimensional to the triaxial symmetric set-ups, one expects increasingly contrasted dynamical behavior of the system.
Lagrangian Perturbation Theory.—In LPT, the displacement field is the fundamental building block which is considered as a small quantity. As long as function is a single-valued function of , there is no singularity in the density field and a systematic perturbative expansion is possible, . The evolution equation of the displacement field at each order is obtained from Eq. (1). Taking the divergence and the curl of this equation in Eulerian coordinates space and rewriting the expressions in terms of the Lagrangian coordinates, a set of recurrence relations is obtained by substituting the perturbative expansion into Eq. (1) Rampf (2012); Zheligovsky and Frisch (2014); Rampf et al. (2015); Matsubara (2015):
for the longitudinal part, and
for the transverse part. Here, is the Levi-Civita symbol, and stands for a differential operator, , where . Thus, according to the initial conditions given by Eq. (3), one builds up from Eqs. (4) and (5) the perturbative solutions for the two kinds of derivatives and . Then, the displacement field can be consistently reconstructed from these derivatives. This last step is non-trivial and generally involves intricate calculations. However, thanks to the sinusoidal nature of the initial conditions we consider, the final form of is given as a series of harmonic functions which can be obtained through simple algebraic manipulations.
Since the initial set-up assumes a small amplitude of the fluctuations, , and we consider time sufficiently close to collapse, , taking only the fastest growing-mode provides an accurate description of the dynamics. With this additional simplification in mind, we perform the perturbative calculation in two different ways, detailed as follows. The first approach, quite standard, consists in using Eq. (3) as the first-order solution, i.e., , and develop higher-order calculations (LPT) up to tenth order. The second approach, that we denote by Q1D, assumes quasi one dimensional dynamics, following the footsteps of Rampf and Frisch (2017), i.e. . In this case, one takes the exact solution of the 1D problem along -axis given by Zel’dovich approximation as the unperturbed zeroth-order state: , , with being the Kronecker delta function. Transverse fluctuations are considered as small first order perturbations to this set-up, . The perturbative expansion is then performed by keeping terms proportional to and up to second order, (so in this sense we go one order beyond Rampf and Frisch (2017)), while keeping terms proportional to up to tenth order, .
Vlasov simulations and phase-space structure.—To quantitatively investigate the dynamics of our system up to shell-crossing, we perform high resolution simulations with the state-of-the-art Vlasov-Poisson code ColDICE Sousbie and Colombi (2016). This code uses a tessellation, i.e. tetrahedra, to represent the phase-space sheet. The vertices of the tessellation form initially a homogeneous mesh of resolution (which corresponds to simplices) and then follow Lagrangian equations of motion. Poisson equation is solved by Fast-Fourier-Transform on a regular cartesian grid of resolution , after exact deposit at linear order of the phase-space sheet on the grid. A number of simulations, as listed in Table 1, were performed assuming Einstein-de Sitter cosmology.
Fig. 1 shows representative examples of phase-space portraits at shell-crossing time. As the ratios increase, Zel’dovich approximation, which is exact in the strictly 1D case , starts to deviate significantly from the simulation, as expected. The Q1D prediction provides a substantial improvement, with an excellent match of the simulation measurements for (top panel). Still, it cannot catch up the shell-crossing structure when ratios approach unity (middle and bottom panels), because convergence of LPT gets slower. However, with a systematic calculation of all the contributions up to tenth-order, LPT prediction improves considerably and accurately reproduces the shell-crossing structure seen in the simulations for most of the parameter space, except in the vicinity of (bottom panel).
In the case, the phase-space structure is highly stretched along velocity axis, similarly as in spherical symmetry in the absence of angular momentum Fillmore and Goldreich (1984); Bertschinger (1985). This means that the inward mass flow along each direction is accelerated near , which is a natural consequence of the mutual interactions resulting from the simultaneous collapses taking place along the three axes. As a result, convergence of the LPT expansion is much slower near and even the tenth-order calculation is insufficient (see also Rampf (2017) for a discussion about the spherically symmetric case).
However, by studying the sequence of LPT predictions as a function of order up to , it is possible to extrapolate the asymptotic convergence at . Indeed, we find that the position of matter elements at collapse time calculated at order with LPT, , is well described by the following fitting form:
with and where , , , and are fitting parameters which depend both on initial conditions and Lagrangian position . Taking the limit gives the extrapolated result at infinite order, . The same treatment can be applied just before and after shell-crossing: by taking the central finite time difference, the extrapolated velocity, , is also obtained.
Examination of Fig. 1 shows that the result of this procedure (cyan curves) reproduces very well simulation measurements, even the spiky structure in the case. Disagreement is at worse a few percent when . Even if there are still some small mismatches, partly attributable to a small desynchronization due to the finite time step in the simulations, the extrapolation based on Eq. (6) is unquestionably successful and can be used to perform quantitative predictions over the entire parameter space covered by . Note that the same fitting technique has been employed to predict the shell-crossing time, , used as the output time for the simulations shown in Fig. 1.
Analysis of the shell-crossing structure.—Making use of the generic convergence behavior described in Eq. (6), we can explore the global parameter dependence of the shell-crossing structure, by studying the properties of the shape of the 2D section shown on Fig. 1. Fig. 2 summarizes the results of a parameter space survey, focusing on the maximum of the normalized velocity (top panel) and the corresponding Lagrangian position (bottom panel). Again, as expected, for the sub-space of parameters probed by our runs, the theoretical predictions given by the black dots are found to be in good agreement with the simulations.
What is more striking about Fig. 2 is the sudden transition in the vicinity of the point , where the maximum velocity suddenly augments when these ratios increase towards unity (upper right part of top panel) while the corresponding Lagrangian coordinate gets very close to zero, suggesting the convergence to some singular behavior (upper right part of bottom panel). This sudden change is in fact also associated to a drastic variation in the shape of the structure of the phase-space sheet at collapse, as illustrated by Fig. 3 in sub-space, where we consider the case and values of increasing from to unity. As seen in this figure, the cross-section of the phase-space sheet changes drastically from a smooth “S” shape, which is the normal behavior for most of values of the ratios , to a spiky structure when both these ratios approach unity, . Again, the case is special because of the (close to) simultaneous collapses along three main directions of motion and one expects a different behavior of post-collapse dynamics compared to the typical cases. This might correspond in the CDM paradigm to a population of rare halos or subhalos with particular properties and is worth investigating more in a future work.
Conclusion and outlook.—With Lagrangian perturbation theory (LPT) extrapolated to infinite order, we found a way to describe accurately the phase-space structure of proto-halos growing from three initial sine waves of various amplitudes, , and , until collapse time. To validate the theory, we used the state of the art Vlasov code of Sousbie and Colombi (2016). During our investigations, we found, as expected, that convergence of the LPT series expansion gets worse when going from quasi-one dimensional to tri-axial symmetric initial conditions. An exploration of parameter space shows that a spiky structure in phase-space appears when approaching tri-axial symmetry, . To which extent such a spiky structure can have consequences on the subsequent evolution of the system is an interesting question and will be the object of future investigations.
Given the high level of symmetry of the initial conditions we considered, some additional investigations are required to address the general case. However, we are confident that our results are quite generic. They set up the framework for accurate theoretical investigations beyond collapse time. Indeed, except for rare cases such as the triaxial symmetric configuration, collapse is generally expected to produce a planar singularity (as illustrated by two top panels of Fig. 1), which means that near collapse time, dynamics is quasi unidimensional. Starting from the state of the system at crossing time, that we can now describe accurately from the analytical point of view, it is certainly possible to develop a post-collapse Lagrangian perturbation theory formalism, generalizing the investigations performed in one dimension by Colombi (2015); Taruya and Colombi (2017) to the fully three-dimensional case and using the fact that the actual dynamical structure remains close to one-dimensional when considering epochs close to collapse time. We leave this project for future work.
Acknowledgements.—This work was supported in part by JSPS Research Fellow Grant-in-Aid Number 17J10553 (SS) and MEXT/JSPS KAKENHI Grants, Numbers JP15H05889 and JP16H03977 (AT), as well as ANR grant ANR-13-MONU-0003 (SC). Numerical computation was carried out using the HPC resources of CINES (Occigen supercomputer) under the GENCI allocation 2017-A0040407568. It has also made use of the Yukawa Institute Computer Facility.
- Peebles (1982) P. J. E. Peebles, ApJ 263, L1 (1982).
- Peebles (1984) P. J. E. Peebles, ApJ 277, 470 (1984).
- Blumenthal et al. (1984) G. R. Blumenthal, S. M. Faber, J. R. Primack, and M. J. Rees, Nature 311, 517 (1984).
- Berezinsky et al. (2003) V. Berezinsky, V. Dokuchaev, and Y. Eroshenko, Phys. Rev. D 68, 103003 (2003), astro-ph/0301551 .
- Diemand et al. (2005) J. Diemand, B. Moore, and J. Stadel, Nature 433, 389 (2005), astro-ph/0501589 .
- Ishiyama (2014) T. Ishiyama, ApJ 788, 27 (2014), arXiv:1404.1650 .
- Angulo et al. (2017) R. E. Angulo, O. Hahn, A. D. Ludlow, and S. Bonoli, MNRAS 471, 4687 (2017), arXiv:1604.03131 .
- Ogiya and Hahn (2018) G. Ogiya and O. Hahn, MNRAS 473, 4339 (2018), arXiv:1707.07693 .
- Navarro et al. (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 462, 563 (1996), astro-ph/9508025 .
- Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490, 493 (1997), astro-ph/9611107 .
- Navarro et al. (2010) J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, and A. Helmi, MNRAS 402, 21 (2010), arXiv:0810.1522 .
- Diemand and Moore (2011) J. Diemand and B. Moore, Advanced Science Letters 4, 297 (2011), arXiv:0906.4340 .
- Fukushige and Makino (1997) T. Fukushige and J. Makino, ApJ 477, L9 (1997), astro-ph/9610005 .
- Fillmore and Goldreich (1984) J. A. Fillmore and P. Goldreich, ApJ 281, 1 (1984).
- Bertschinger (1985) E. Bertschinger, ApJS 58, 39 (1985).
- Alard (2013) C. Alard, MNRAS 428, 340 (2013), arXiv:1206.6548 [astro-ph.CO] .
- Lynden-Bell (1967) D. Lynden-Bell, MNRAS 136, 101 (1967).
- Hjorth and Williams (2010) J. Hjorth and L. L. R. Williams, ApJ 722, 851 (2010), arXiv:1010.0265 [astro-ph.CO] .
- Carron and Szapudi (2013) J. Carron and I. Szapudi, MNRAS 432, 3161 (2013), arXiv:1301.6760 .
- Pontzen and Governato (2013) A. Pontzen and F. Governato, MNRAS 430, 121 (2013), arXiv:1210.1849 .
- Bernardeau et al. (2002) F. Bernardeau, S. Colombi, E. Gaztañaga, and R. Scoccimarro, Phys. Rep. 367, 1 (2002), astro-ph/0112551 .
- Zel’dovich (1970) Y. B. Zel’dovich, A&A 5, 84 (1970).
- Shandarin and Zeldovich (1989) S. F. Shandarin and Y. B. Zeldovich, Reviews of Modern Physics 61, 185 (1989).
- Bouchet et al. (1992) F. R. Bouchet, R. Juszkiewicz, S. Colombi, and R. Pellat, ApJ 394, L5 (1992).
- Buchert (1992) T. Buchert, MNRAS 254, 729 (1992).
- Buchert and Ehlers (1993) T. Buchert and J. Ehlers, MNRAS 264, 375 (1993).
- Bouchet et al. (1995) F. R. Bouchet, S. Colombi, E. Hivon, and R. Juszkiewicz, A&A 296, 575 (1995), astro-ph/9406013 .
- Novikov (1969) E. A. Novikov, Soviet Journal of Experimental and Theoretical Physics 30, 512 (1969).
- Arnold et al. (1982) V. I. Arnold, S. F. Shandarin, and I. B. Zeldovich, Geophysical and Astrophysical Fluid Dynamics 20, 111 (1982).
- Hidding et al. (2014) J. Hidding, S. F. Shandarin, and R. van de Weygaert, MNRAS 437, 3442 (2014), arXiv:1311.7134 [astro-ph.CO] .
- Feldbrugge et al. (2017) J. Feldbrugge, R. van de Weygaert, J. Hidding, and J. Feldbrugge, ArXiv e-prints (2017), arXiv:1703.09598 .
- Moutarde et al. (1991) F. Moutarde, J.-M. Alimi, F. R. Bouchet, R. Pellat, and A. Ramani, ApJ 382, 377 (1991).
- Buchert et al. (1997) T. Buchert, G. Karakatsanis, R. Klaffl, and P. Schiller, A&A 318, 1 (1997), astro-ph/9510133 .
- Zheligovsky and Frisch (2014) V. Zheligovsky and U. Frisch, Journal of Fluid Mechanics 749, 404 (2014), arXiv:1312.6320 [math.AP] .
- Rampf et al. (2015) C. Rampf, B. Villone, and U. Frisch, MNRAS 452, 1421 (2015), arXiv:1504.00032 .
- Rampf (2012) C. Rampf, J. Cosmology Astropart. Phys 12, 004 (2012), arXiv:1205.5274 .
- Matsubara (2015) T. Matsubara, Phys. Rev. D 92, 023534 (2015), arXiv:1505.01481 .
- Rampf and Frisch (2017) C. Rampf and U. Frisch, MNRAS 471, 671 (2017), arXiv:1705.08456 .
- Sousbie and Colombi (2016) T. Sousbie and S. Colombi, Journal of Computational Physics 321, 644 (2016), arXiv:1509.07720 [physics.comp-ph] .
- Rampf (2017) C. Rampf, ArXiv e-prints (2017), arXiv:1712.01878 .
- Colombi (2015) S. Colombi, MNRAS 446, 2902 (2015), arXiv:1411.4165 .
- Taruya and Colombi (2017) A. Taruya and S. Colombi, MNRAS 470, 4858 (2017), arXiv:1701.09088 .