Issue 
A&A
Volume 664, August 2022



Article Number  A3  
Number of page(s)  31  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202142756  
Published online  03 August 2022 
Cold dark matter protohalo structure around collapse: Lagrangian cosmological perturbation theory versus Vlasov simulations
^{1}
Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris, CNRS,
92190
Meudon, France
email: shohei.saga@obspm.fr
^{2}
Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kyoto,
6068502
Japan
^{3}
Kavli Institute for the Physics and Mathematics of the Universe (WPI), Todai institute for Advanced Study, University of Tokyo,
Kashiwa,
Chiba,
2778568
Japan
^{4}
Sorbonne Universtié, CNRS, UMR7095, Institut d’Astrophysique de Paris,
98bis boulevard Arago,
75014
Paris, France
Received:
26
November
2021
Accepted:
8
May
2022
We explore the structure around the shellcrossing time of cold dark matter protohaloes seeded by two or three crossed sine waves of various relative initial amplitudes, by comparing Lagrangian perturbation theory (LPT) up to the tenth order with highresolution cosmological simulations performed with the public Vlasov code ColDICE. Accurate analyses of the density, the velocity, and related quantities such as the vorticity are performed by exploiting the fact that ColDICE can follow the phasespace sheet locally at the quadratic level. To test LPT predictions beyond the shellcrossing, we employ a ballistic approximation, which assumes that the velocity field is frozen just after the shellcrossing. In the generic case, where the amplitudes of the sine waves are all different, highorder LPT predictions match the exact solution very well, even beyond collapse. As expected, convergence slows down when going from quasi1D dynamics, where one wave dominates over the two others, to the axialsymmetric configuration, where all the amplitudes of the waves are equal. We also notice that LPT convergence is slower when considering velocityrelated quantities. Additionally, the structure of the system at and beyond collapse given by LPT and the simulations agrees very well with singularity theory predictions, in particular with respect to the caustic and vorticity patterns that develop beyond collapse. Again, this does not apply to axialsymmetric configurations, which are still correct from the qualitative point of view, but rather when multiple foldings of the phasespace sheet produce very high density contrasts and hence a strong backreaction of the gravitational force.
Key words: gravitation / largescale structure of Universe / dark matter / galaxies: kinematics and dynamics
© S. Saga et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
In the concordance model of largescale structure formation, the matter content of the Universe is dominated by collisionless cold dark matter (CDM) following VlasovPoisson equations (Peebles 1982, 1984; Blumenthal et al. 1984). The cold nature of initial conditions implies that the dark matter distribution is concentrated on a 3D sheet evolving in 6D phasespace. At the shellcrossing, that is, in places where the phasespace sheet first selfintersects, the seeds of the first dark matter haloes are created. In these regions, the fluid enters a multistream regime during which violent relaxation takes place (LyndenBell 1967) to form primordial CDM haloes. Numerical investigations suggest that dark matter haloes formed during this process initially have a powerlaw density profile ρ ∝ r^{α}with α ≈ –1.5 (Moutarde et al. 1991; Diemand et al. 2005; Ishiyama 2014; Angulo et al. 2017; Ogiya & Hahn 2018; Delos et al. 2018; Colombi 2021). During their subsequent evolution, which includes successive mergers, dark matter haloes relax to the wellknown universal NavarroFrenkWhite (NFW) profile (Navarro et al. 1996, 1997).
From an analytical point of view, many approaches have been proposed to describe the results of simulations and the different steps of dark matter halo history, relying, for example, on entropy maximisation (LyndenBell 1967; Hjorth & Williams 2010; Carron & Szapudi 2013; Pontzen & Governato 2013), selfsimilarity (Fillmore & Goldreich 1984; Bertschinger 1985; Henriksen & Widrow 1995; Sikivie et al. 1997; Zukin & Bertschinger 2010b,a; Alard 2013), or, more recently, a postcollapse perturbative treatment (Colombi 2015; Taruya & Colombi 2017; Rampf et al. 2021). However, because of the highly nonlinear complex processes taking place during the postcollapse phase, the formation and evolution of dark matter haloes remain a subject of debate, and a consistent framework that explains the different phases of halo history remains to be proposed.
The early growth of largescale structures up to the first shellcrossing, on the other hand, is well understood thanks to perturbation theory. Indeed, restricting ourselves to the early phase of structure formation (i.e. to the singlestream regime before the shellcrossing), we can employ perturbation theory as long as fluctuations in the density field remain small (see e.g. Bernardeau et al. 2002, and references therein for a review). Lagrangian perturbation theory (LPT; e.g. Shandarin & Zeldovich 1989; Bouchet et al. 1992, 1995; Buchert 1992; Buchert & Ehlers 1993; Bernardeau 1994) uses the displacement field as a small quantity in the expansion of the equations of motion. Firstorder LPT corresponds to the classic Zel’dovich approximation (Zel’dovich 1970), and in 1D space it is an exact solution until the shellcrossing (Novikov 1969). Because the Lagrangian description follows elements of fluid along the motion, the Zel’dovich approximation and higherorder LPT provide us with a rather accurate description of the largescale matter distribution, even in the nonlinear regime, shortly after the shellcrossing. The families of singularities that form at the shellcrossing and after have been examined in detail within the Lagrangian dynamics framework (Arnold et al. 1982; Hidding et al. 2014; Feldbrugge et al. 2018), and the structure of cosmological systems at the shellcrossing has been investigated for specific initial conditions (Novikov 1969; Rampf & Frisch 2017; Saga et al. 2018; Rampf 2019) and for random initial conditions (Rampf & Hahn 2021).
As described above, investigations into the early stages of largescale structure formation represent a key to understanding the bottomup scenario underlying the CDM model. Practically, perturbation theory under the singlestream assumption, which is strictly valid only during the early stages of the dynamical evolution, provides the basis for predicting statistics of the largescale structure distribution and has been successfully confronted with Nbody simulations and observations (see e.g. Bernardeau et al. 2002, for a review). In order to incorporate the effects of multistreaming at small scales, introducing effective fluid equations with a nonvanishing stress tensor in the dark matter fluids has recently been proposed, the socalled effective field theory of largescale structure (Baumann et al. 2012; Carrasco et al. 2012; Hertzberg 2014; Baldauf et al. 2015). Although this approach needs parameters in the stress tensor to be calibrated with Nbody simulations, it has attracted much attention and has recently been applied to real observational datasets to derive cosmological parameter constraints (Ivanov et al. 2020; dAmico et al. 2020). Accordingly, the understanding of multistreaming effects from first principles, even at an early stage, would provide significant insights into the precision theoretical modelling of largescale structure.
The aim of this paper is to extend the investigations of Saga et al. (2018), hereafter STC18. In that work, we studied the phasespace structure of protohaloes at the shellcrossing with simplified initial conditions composed of three crossed sine waves, following in the footsteps of Moutarde et al. (1991, 1995). Thanks to a detailed comparison between LPT predictions up to the tenth order and numerical simulations performed with the stateoftheart VlasovPoisson solver ColDICE (Sousbie & Colombi 2016), we explicitly showed that the convergence of the LPT series slows down when going from quasi1D to triaxialsymmetric initial conditions, where a spiky structure in phasespace appears. Using a global fitting form, we were able to formally extrapolate the LPT solution to the infinite order, obtaining remarkable agreement with the simulations, even at the shellcrossing. In the present work, we thoroughly examine the structure further, both at and shortly after the shellcrossing of these systems, by considering the 2D case as well. Although sine wave initial conditions are restrictive, they are to a large extent representative of high peaks of a smooth random Gaussian field (see e.g. Bardeen et al. 1986) and should provide fruitful guidance in the general case. Furthermore, initial conditions expressed only in terms of sine waves considerably simplify analytical calculations while still allowing for an insightful exploration of the pre and postcollapse dynamics.
The prominent features shortly after the shellcrossing are, for instance, the appearance of caustics and, in the multistream region delimited by these caustics, nontrivial vorticity patterns (Doroshkevich 1973; Chernin 1993; Pichon & Bernardeau 1999; Hahn et al. 2015). Thanks to a description of the phasespace sheet at the quadratic level in ColDICE, we can measure these quantities in highresolution Vlasov simulations, in particular the vorticity, with unprecedented accuracy. To perform theoretical predictions shortly after the shellcrossing, we use a simple dynamical approximation based on ballistic motion applied to the state of the system described by highorder LPT solutions at the shellcrossing. While convergence with the perturbation order of LPT remains a complex subject of investigation (Zheligovsky & Frisch 2014; Rampf et al. 2015; Michaux et al. 2021), it seems to take place at least up to the shellcrossing, not only for the threesinewave configurations we aim to examine (STC18) but also for more general, random initial conditions (Rampf & Hahn 2021), although the effects of cutoffs on power spectra remain to be investigated further in the latter case. Therefore, as long as the backreaction from multistream flows on postcollapse dynamics remains negligible, the approximation of the dynamics we propose here should work shortly after the shellcrossing. This is the first step towards a proper analytical description of postcollapse dynamics in 6D phasespace.
This paper is organised as follows. In Sect. 2, we begin by introducing the basics of LPT and its recurrence relations, as well as its applications to initial conditions given by linear superpositions of trigonometric functions, such as the initial conditions of our sine waves. In Sect. 3, the important features of the Vlasov solver ColDICE are briefly summarised. We also address some aspects of measurement uncertainties in simulation data. In Sect. 4, we examine the phasespace structure and radial profiles at the shellcrossing with a comparison of analytic predictions with simulations, in the framework of singularity theory. In Sect. 5, the structure shortly after the shellcrossing is explored using the ballistic approximation by examining phasespace diagrams, caustics, and density and vorticity fields. Again, analytical predictions are compared with the Vlasov runs. Finally, Sect. 6 summaries the main results of this article. To complement the main text, Appendices A–D provide details on the explicit form of highorder LPT solutions for the initial conditions of the sine waves up to the fifth order, the measurements in the Vlasov simulations, the predictions from quasi1D LPT, and the expected structure of singularities at the collapse time, respectively.
Throughout the paper, we use the following units: a box size L = 1 and an inverse of the Hubble parameter at present time, , for the dimensions of length and time, respectively. In particular, the Lagrangian coordinate, Eulerian coordinate, velocity, and vorticity are explicitly expressed as q/L, x/L, υ/(LH_{0}), and ω/Η_{0} = (L∇) × (υ/(LH_{0})), respectively.
2 Lagrangian Perturbation Theory
The Lagrangian equation of motion of a fluid element is given by (e.g. Peebles 1980) (1)
where the quantities x, a, and H(t) = a^{−1}da/dt are the Eulerian comoving position, the scale factor of the Universe, and Hubble parameter, respectively. The derivative operator ∇_{x} = ∂/∂x is a spatial derivative in Eulerian space. The Newton gravitational potential ϕ(x) is related to the matter density contrast with being the background mass density, through the Poisson equation: (2)
In this framework, the velocity of each mass element is given by υ = a dx/dt.
Taking the divergence and curl of Eq. (1) with respect to Eulerian coordinates, Eqs. (1) and (2) can be expressed by the equivalent set of equations: (3) (4)
where the dot represents the Lagrangian derivative of time, d/dt.
In the Lagrangian approach, for each mass element, the Eulerian position x at the time of interest t is related to the Lagrangian coordinate (initial position) q through the displacement field Ψ(q, t) by (5)
The velocity field is expressed as υ(q, t) = a dΨ/dt. In the single flow regime (i.e. before the first shellcrossing time t_{sc}), mass conservation implies , which leads to (6)
where the quantity J = det J_{ij} is the Jacobian of the matrix J_{ij} defined by (7)
and its inverse is given as (8)
While Eq. (6) is well defined only until the first shellcrossing time t_{sc}, that is, the first occurrence of J = 0, Eqs. (7) and (8) can be formally used beyond t_{sc}, as long as they are expressed in terms of Lagrangian coordinates (except that Eq. (8) might become singular).
2.1 Recursion Relation
In LPT, the displacement field, Ψ, is the fundamental building block and is considered to be a small quantity. It can be systematically formally expanded as follows, (9)
In what follows, we assume that the fastest growing modes dominate. They are known to be given to a very good approximation by (see e.g. Bernardeau et al. 2002, and references therein) (10)
where the purely timedependent function D_{+}(t) corresponds to the linear growth factor. The velocity field is then given by (11)
where function f(t) ≡ d1n D_{+}/d In a corresponds to the logarithmic derivative of the growth factor. We note that the analyses performed in subsequent sections will consider the Einsteinde Sitter cosmology, that is, a total matter density parameter Ω_{m} = 1 and a cosmological constant density parameter Ω_{Λ} = 0. In this case, one simply has D_{+} ∝ a and f = 1.
During early phases of structure formation, the Universe approaches Einsteinde Sitter cosmology and f is approximately given by . With the further approximation , as implicitly assumed in all the subsequent calculations (see e.g. Peebles 1980; Bernardeau et al. 2002; Matsubara 2015), by substituting Eqs. (5)–(9) into Eqs. (3) and (4), one obtains simple recurrence formulae for the longitudinal and transverse parts of the motion (Rampf 2012; Zheligovsky & Frisch 2014; Rampf et al. 2015; Matsubara 2015): (12) (13)
in the 3D case. Here and in what follows, we adopt the Einstein summation convention when the equation includes repeated lowercase Latin letters, and the subscripts i, j,… take 1, 2 or 3 (1 or 2) in the 3D (2D) case. In the above, we defined Ψ_{i,j}·≡ ∂Ψ_{i}/∂q_{j}, and the tensors ɛ_{ij} and are, respectively, the 2D and 3D LeviCivita symbols. The symbol stands for a differential operator: (16)
In obtaining the recursion relations, we used the following formulae of linear algebra: (17)
for the 3D case.
In the fastestgrowing mode regime (10), which will be assumed in the LPT calculations of this work, the time dependence in the recursion relations simplifies and one obtains (19) (20) (21) (22)
for the 3D case.
2.2 Twosinewave and Threesinewave Initial Conditions
Throughout this paper, we focus on specific initial conditions: two or three sine waves in a periodic box covering the interval [–L/2, L/2[: (23)
The parameters ϵ_{i} < 0 with ϵ_{x} ≥ ϵ_{y}≥ ϵ_{z} quantify the linear amplitudes of the sine waves in each direction. The initial time, t_{ini}, is chosen such that D_{+}(t_{ini})ϵ_{i} ≤ 0.012 ≪ 1, which makes the fastest growing mode approximation very accurate, as shown by STC18. In this framework, the dependence on ϵ_{i} of the dynamics is reduced to a function of the ratios ϵ_{2D} = ϵ_{y}/ϵ_{x} and ϵ_{3D} = (ϵ_{y}/ϵ_{x},ϵ_{Z}/ϵ_{x}), respectively, for two and threesinewave initial conditions. These ratios will therefore be the quantities of relevance to define our initial conditions. In this setting, the initial density field, given by δ^{ini} ≃ –∇_{q}. Ψ^{ini} = D_{+}(t_{ini})Σ_{i} ϵ_{i} cos (2π/Lq_{i}), presents a small peak at the origin, and mass elements subsequently infall towards the central overdense region. With the proper choice of ϵ_{i}, this initial overdensity can asymptotically match any peak of a smooth random Gaussian field (see e.g. Bardeen et al. 1986), which actually makes the threesinewave initial conditions quite generic, hence providing many insights into the dynamics during the early stages of protohalo formation.
Naturally, this very symmetrical setup remains unrealistic, with a tidal environment restricted to periodic replica, but has the advantage of belonging to the family of linear superpositions of trigonometric functions, here simple sine functions, which considerably simplifies LPT calculations, as described below. Initial conditions that only involve linear superpositions of trigonometric functions include exact Fourier transforms, so are in principle very general. Furthermore, with only a few Fourier modes, one can theoretically account for more realistic initial conditions with proper tidal environments and mergers, while keeping the analytical description still very tractable. However, as far as we are concerned, applications beyond two or three sine waves are left for future work.
For initial conditions given by linear superpositions of trigonometric functions, it is trivial to see that all the terms of the perturbation series are also expressed in the same way: (24) (25)
where the nthorder scalar coefficients and the nthorder vector coefficients are obtained recursively by calculating the righthandside of Eqs. (21) and (22), which depends on lowerorder terms, starting from the n = 1 coefficients determined by Eq. (23).
By imposing the conditions and , one can build up the perturbative solutions and using simple algebraic manipulations involving coefficients and : (26) (27)
These solutions lead to the nthorder displacement field given by . Using this solution for Ψ^{(n)} as well as Eqs. (21) and (22), we subsequently construct the source of the (n + 1)thorder derivatives and . By repeating the above operation together with the recursive relations, one can derive, in principle, arbitrary highorder LPT solutions. In Appendix A, we present the explicit forms of the LPT solutions up to the fifth order derived in this way^{1}.
The above prescription is valid in 3D space and can be applied to the 2D case by cancelling all fluctuations along the z axis, that is, by performing 3D calculations with ϵ_{z} = 0 (i.e. ϵ_{3D} = (ϵ_{y}/ϵ_{x}, 0)) for the three sine waves case. However, one can also realise that in 2D, vector coefficients become scalars , and that the solutions take the following form: (28) (29)
Our analytical investigations can easily cover a large range of values of ϵ_{2D} and ϵ_{3D}, while the simulations, much more costly, will only focus on three configurations, as detailed in Table 1, reflecting various regimes in the dynamics: quasi 1D with ϵ_{x} ≫ ϵ_{y,z}, anisotropic with ϵ_{x} > ϵ_{y} > ϵ_{z}, and what we design by axialsymmetric, with ϵ_{x} = ϵ_{y}( = e_{z}), denoted by Q1D, ANI, and SYM, respectively^{2}.
3 Vlasovpoisson Simulations
To perform the numerical experiments, we use the public parallel Vlasov solver ColDICE (Sousbie & Colombi 2016). ColDICE follows the phasespace sheet with an adaptive tessellation of simplices, composed, in 2D and 3D, of connected triangles and connected tetrahedra, respectively. The phasespace coordinates of the vertices of the tessellation, [X(t), V(t)], follow the standard Lagrangian equations of motion, similarly as in an Nbody code, but matter is distributed linearly inside each simplex instead of being transported by the vertices. In what follows, after providing a few additional technical details on the algorithm (Sect. 3.1), we address measurement uncertainties issues (Sect. 3.2).
Parameters of the runs performed with ColDICE (Sousbie & Colombi 2016).
3.1 The Vlasov Code ColDICE
The Lagrangian coordinates defined in Sect. 2.1 correspond to the following unperturbed initial setup, (30) (31)
where t_{0} is a virtual time corresponding to a = 0 and Q is the Lagrangian coordinate of each vertex. We note that we use capital letters to distinguish between vertex coordinates and the actual coordinates of fluid elements of the phasespace sheet that they are supposed to trace. These notations are used in Appendix B.
Vertices phasespace coordinates are then perturbed using Zel’dovich motion to set actual initial conditions defined in Sect. 2.2: (32) (33)
with Ψ^{ini} given by Eq. (23).
To update the position and the velocity of each vertex, a standard secondorder predictorcorrector scheme with slowly varying time step is employed. Constraints on the value of the time step combine bounds on the relative variations of the expansion factor, classical CourantFriedrichsLewy criterion and a harmonic condition related to the local projected density. More details about the parameters used for the time step constraints in the 3D simulations are given in Colombi (2021, hereafter C21), so we do not repeat these details here. Additionally, we chose the same constraints on the time step for the 2D runs as for the 3D simulations with n_{g} = n_{s} = 512 (see below for the definitions of n_{g} and n_{s}).
The Poisson equation is solved using the fastFourier technique in a mesh of fixed resolution n_{g}. To estimate the projected density ρ(x) on this mesh, the intersection of each simplex of the phasespace sheet with each voxel or pixel of the mesh is computed exactly up to the linear order with a special raytracing technique. Once the gravitational potential is obtained on the mesh, the gravitational force is computed from the gradient of the potential using a standard fourpoint stencil, and then it is interpolated to the vertices using secondorder triangular shape cloud interpolation (see e.g. Hockney & Eastwood 1988) in order to update their velocities.
Initially, the tessellation is constructed from a regular network of vertices corresponding, respectively, to and simplices in 2D and 3D. ColDICE allows for local refinement of the simplices following criteria based on local Poincaré invariant conservation as explained more in detail in Sousbie & Colombi (2016). Values of the refinement criterion parameter I used for our 3D runs are listed in detail in C21. For completeness, following the notations of C21, we used I = 10^{−8} for the 2D runs.
To perform local refinement, the phasespace sheet is locally described at the quadratic order inside each simplex with the help of three and six additional tracers per simplex in 2D and 3D, respectively. At the dynamical times considered in this work, which correspond at most to short periods after collapse, refinement is not triggered, except for a small number of simplices in the 3D axialsymmetric simulation, SYM3SIN, so we do not deem it necessary to discuss more about refinement. However, the ability to describe the phasespace sheet at the quadratic level is important to have correct estimates of derivatives of the velocity field, in particular of the local vorticity of the mean flow.
The parameters used for all the simulations are listed in Table 1, in particular the resolution n_{g} of the mesh used to solve Poisson equation and the initial number of vertices of the tessellation, n_{s}^{3}. In Appendix B, we explain how measurements of various quantities are performed, such as the set of curves corresponding to the intersection of the phasespace sheet with the hyperplane y( = z) = 0 used in Sects. 4.2 and 5.2, the collapse time t_{sc} shown in Table 1, the radial profiles used in Sect. 4.3, as well as the caustic network, the projected density and the vorticity analysed in Sect. 5.3.
3.2 Accuracy and Related Considerations
Accuracy and possible defects of ColDICE are discussed in Sousbie & Colombi (2016) and in C21 in the 3D case. In our work, precise determination of the collapse time is critical, since we consider, in the analyses performed below, either the collapse time or a moment very shortly after it. The uncertainty on the collapse time depends on the nature of the initial conditions and the parameters that control the accuracy of the simulations, in particular the number of simplices used to represent the phasespace sheet and the spatial resolution n_{g} of the mesh used to solve Poisson equation. This latter turns out to be a very important parameter. Indeed, decreasing n_{g} augments the collapse time, as discussed in C21. In Sect. B.2, relying on measurements in ANI3SIN configurations with various spatial resolutions, we estimate that the quantity a_{sc} shown in Table 1 should be accurate at the fourth significant digit level or better for all configurations, for example a_{sc} = 0.02919 ± 10^{−4} for ANI3SIN, except for SYM3SIN, where the uncertainty on the measured a_{sc} could be a few 10^{−4}. Indeed, as just mentioned above, accuracy on the determination of the collapse time also depends on the nature of initial conditions, the quasi1D case and the triaxialsymmetric configuration being the easiest and the most difficult to deal with, respectively, as expected.
We note that accuracy regarding the collapse time also depends on the ability to determine the exact position of the shellcrossing point. Using the tessellation technique, the shellcrossing coincides with a temporal change of sign of the orientation of the simplices. In our symmetrical setup, the position of the shellcrossing point is simply the centre of the system and the determination of the collapse time can be simply performed using the intersection of the phasespace sheet with the y = z= 0 hyperplane (see Appendix B.2). However, even in a more complex case where the centre of the forming halo would be moving, one expects to be able to accurately compute the collapse time thanks to the tessellation representation as long as sampling of the phasespace sheet is sufficiently dense and, of course, if the quadratic representation inside each simplex is fully exploited. Sampling of the phasespace sheet is controlled by n_{s}, as well as the refinement parameter I limiting deviations from symplectic motion. As discussed in Sect. 3.1, for our simulations, parameter I does not influence the results except (probably marginally) for SYM3SIN. Effects related to the choice of n_{s} are discussed in Sousbie & Colombi (2016) and in C21. They should be negligible for our sine wave simulations compared with other sources of uncertainty on the collapse time for the values of n_{s} we adopted. However, we did not test the accuracy of ColDICE in more complex cases, because we do not need to.
An additional advantage of our setting is that the axes of the sine waves are aligned with the mesh used to solve Poisson equation, which helps make the simulations more accurate: the same configurations in an oblique fashion would not achieve the same level of accuracy. This is illustrated by Appendix H of Sousbie & Colombi (2016) in the single sine wave case. As shown in this Appendix, anisotropies due to misalignments between preferred directions in the representation and preferred directions in the dynamics are the source of significant accuracy loss. However, these cumulative effects should remain small for ColDICE at the early times considered in this work.
Another critical issue is to determine accurately the shape of the caustics and then the density and vorticity fields inside the multistream region. The technique employed to extract the caustics curves and surfaces from the 2D and the 3D simulations, respectively, is rather simple and robust, as briefly sketched in Appendix B.3. Based on a detection of a spatial sign change on the simplices orientation, it consists in extracting a subset of segments (in 2D) or triangles (in 3D) from the tessellation. These segments and triangles correspond to intersections between adjacent simplices. Obviously this is a rather crude way to draw the caustics, since, in Lagrangian space, it ends up into a subset of segments/triangles of a regular pattern. Therefore, the accuracy in the determination of the caustics location in Lagrangian space is at best of the order of the distance between two neighbouring vertices: L/n_{S}, where L is the size of the simulation cube. This uncertainty, which adds up to the effects discussed above, is difficult to transpose to Eulerian space due to the variations in the strain tensor. It is however probably the main source of inaccuracy on the determination of the caustics in ColDICE. A way to improve greatly on this would consist in exploiting the quadratic representation inside each simplex, which we do not do here.
Turning to the fields, the vorticity, because it is composed of derivatives, is particularly challenging to measure even with the considerable gain brought by the tessellation technique compared with more traditional representations based on particles. Appendix B.5 explains in details how we measure density and vorticity, making use, this time, of the locally quadratic representation of the phasespace sheet inside each simplex. Accuracy on the fields measurements is not discussed quantitatively, as it was not deemed necessary for the analyses performed in this work, but we notice that vorticity becomes very noisy when very close to the caustic curves or surfaces, which we keep in mind for the analyses carried out in Sect. 5.3.
4 Shellcrossing Structure
We are now in a position to study the structure of our protohaloes at the collapse time, t_{sc}, and concentrate our investigations on phasespace diagrams and radial profiles, with comparisons of LPT pushed up to the tenth order to the Vlasov runs. This section is organised as follows. First, Sect. 4.1 presents the calculation of the collapse time itself. Indeed, this quantity depends on initial conditions and perturbation order, and the ability of LPT to provide an accurate determination of t_{sc} is of prime importance. We discuss the extrapolation to the infinite order of the LPT series proposed by STC18 and generalise it to the 2D case. Second, Sect. 4.2 examines the convergence of LPT at collapse with phasespace diagrams, extending the earlier investigation of STC18 to the 2D case. For comparison, we also tested the formal extension of LPT to the infinite order and the predictions of the quasi1D approach proposed by Rampf & Frisch (2017, hereafter RF17). Finally, in Sect. 4.3, LPT predictions and their convergence are studied in terms of radial profiles of the density, the velocity as well as the pseudo phasespace density, and put into perspective in relation to singularity theory.
4.1 Shellcrossing Time
This subsection presents estimates of the expansion factor at the first shellcrossing, which we refer to as the collapse time. In the following, the expansion factor will be formally identified to a time variable, still denoted by a to contrast with actual physical time t. We compare the value of the collapse time obtained from rathorder LPT and its extrapolation to the infinite order , as described in STC18, with the value a_{sc} measured in the Vlasov runs as explained in Appendix B.2. Our approach follows closely that of STC18. It is basically intended to repeat its main steps and to supplement it with additional discussions and comparisons with Vlasov runs in the 2D case.
Using rathorder LPT predictions, we explore the sequence of the shellcrossing times, , as a function of order n up to n = 10, by solving J^{(n)} = 0 at the origin, where J^{(n)} is the Jacobian of the nthorder LPT solution. Generally, except in the pure 1D case where firstorder LPT is exact before collapse, becomes smaller with increasing perturbation order n (see e.g. RF17; Rampf & Hahn 2021, hereafter RH21, for recent works). As illustrated by Fig. 1, the shellcrossing time calculated at the nthorder with LPT is very accurately described by the following fitting form (STC18):
with e > 0. This fitting form, also used for each coordinate of the positions and velocities in Fig. 2 below, does not necessarily represent the sole choice for approximating the n dependence of the collapse time, but using the exponential of a power law might be the only way to match the convergence speed of LPT at large n, when considering quantities computed at the collapse time of each respective order. Equation (34) also implies when n ≫ 1, which might, at first glance, seem incompatible with the findings of RH21, who examined the LPT series in the slightly different context of a Gaussian random field in a periodic box. RH21 concluded from their analyses that convergence speed was asymptotically compatible with a powerlaw times an exponential of n. This would naively correspond to , while our measurements suggest values of e different from unity, ranging from, for example, e ~ 0.003 (SYM3SIN) to e ~ 0.8 (Q1D3SIN) for fitting and from e ~ 0.6–0.7 (SYM3SIN) to e ~ 1.3–1.5 (Q1D3SIN) for fitting the Eulerian position. However, in RH21, the LPT series is studied in terms of the Taylor coefficients of the displacement field expanded as a function of the linear growing mode D_{+}, while we consider the sequence of the shellcrossing times or equivalently . This additional n dependence introduced in the time variable fundamentally changes the nature of the LPT series as a function of n, which makes a direct comparison of Eq. (34) with the findings of RH21 inappropriate. Yet, it would be interesting to investigate the convergence properties of our systems seeded with sine waves using the approach of RH21.
One important thing to note is that convergence of the collapse time with order is rather slow, except in the quasi1D case, which still requires at least the third order for reaching an approximately percent level of accuracy for approximate convergence, while a much higher order is required for other configurations, especially the 3D axialsymmetric case, ϵ_{3D}= (1,1), for which convergence does not even seem to be achieved at the tenth order. Figure 1 thus demonstrates that loworder perturbation theory cannot be used to accurately estimate the collapse times (see also STC 18, in particular, Fig. 2 of this article, and RH21).
Despite these convergence issues, the extrapolated values of obtained with Eq. (34) are very accurate, as illustrated by Table 1. Indeed, the relative difference between and the value a_{sc} measured in the Vlasov simulations remains of the order of the percent or below, which suggests that the error on the extrapolated estimate of the collapse time remains small and should not affect significantly the analytical predictions of the phasespace structure and the radial profiles at the shellcrossing. Therefore, we use for the analytical predictions of the ‘exact’ collapse time in the rest of Sect. 4. As for the simulations data, we use the snapshot with the closest possible available value â_{sc} of the expansion factor to a_{sc}, as indicated in Table 1.
Fig. 1 Convergence behaviour of the shellcrossing time. Top: shellcrossing time calculated at nthorder LPT (dots), together with the fitting function given by Eq. (34) (solid lines), as a function of the perturbation order for various initial conditions, as indicated in the panel. The dashed horizontal lines present the values estimated from the simulations (see Appendix B.1), corresponding to the sixth column of Table 1. Bottom: relative error between the shellcrossing time calculated with nthorder LPT and the one obtained with the fitting formula. 
4.2 Phasespace Structure
Figure 2 shows the (x, υ_{x}) subspace of the phasespace structure at the shellcrossing by considering the intersection of the phasespace sheet with the hyperplane y = 0 for two sine waves and y = z = 0 for three sine waves. For comparison, in addition to standard LPT results, we also present the formal extension of LPT to the infinite order as developed by STC18 and described in Sect. 4.1 for the shellcrossing time, as well as the analytical prediction obtained in the context of the quasi1D approach developed by RF17, which we extended to the second order in the transverse direction (see Appendix C for details).
For quasi1D (Q1D2SIN and Q1D3SIN) and anisotropic (ANI2SIN and ANI3SIN) initial conditions, the phasespace sheet first selfintersects along the x axis, and, until then, the dynamics is similar to the pure 1D case, in which 1LPT (the Zel’dovich solution) is exact until the shellcrossing. As illustrated by Fig. 2, the LPT prediction quickly converges with perturbation order n and describes the simulation results well even at the shellcrossing, especially in the Q1D case. From visual inspection of Fig. 2, it appears that the agreement between the LPT predictions and the simulations is the worse in the vicinity of the extrema of the xcoordinate of the velocity field as a function of position x. We note that, due to the symmetry of the initial conditions, these extrema are global and not specific to the phasespace slice represented in Fig. 2.
In the axialsymmetric cases, SYM2SIN and SYM3SIN, the phasespace sheet selfintersects simultaneously along all the axes of the dynamics. Interestingly, the phasespace structure for ϵ_{2D} = 1 (SYM2SIN) is still qualitatively similar to 1D collapse and LPT again quickly converges with perturbation order to the exact solution. This is clearly not the case for the 3D axialsymmetric configuration, ϵ_{3D} = (1,1), where LPT convergence is very slow. This setup is indeed qualitatively different from other initial conditions, with the appearance of a spiky structure in phasespace (STC18) similarly as in spherical collapse. Interestingly, in the spherical case and in the Einsteinde Sitter cosmology that we consider in this paper, convergence of the LPT is likely to be lost at the collapse time (e.g. Rampf 2019), and the analyses of RH21 suggest that the velocity blows up when convergence is lost. The exact properties of the spike we observe in our numerical data, in particular, whether the velocity diverges or remains finite, and, whether if finite, the velocity is actually smooth at the fine level, remain unknown. While this spike is not present in the LPT predictions at a finite order, it is well reproduced by the formal extrapolation to the infinite order. This is a hint that convergence of the LPT series might be lost at collapse for SYM3SIN.
A question might arise whether, because of such a spike and because of their highly contrasted nature, close to axialsymmetric 3D configurations correspond to a potentially different population of protohaloes. A partial answer can be found in C21, who followed numerically the evolution of the three sinewave configurations further in the nonlinear regime. As described in C21, collapse of our protohaloes is followed by a violent relaxation phase leading to a powerlaw profile ρ(r) ∝ r^{−α}, with α ∈ [1.5,1.7] and then by relaxation to an NFWlike universal profile. After violent relaxation, C21 did not find specific signatures in the density profile nor the pseudo phasespace density for the axialsymmetric case compared with the nonaxialsymmetric ones, except that α tends to augment when going from Q1D to axialsymmetric, and that the axialsymmetric configuration is subject to significant (and expected) radial orbit instabilities.
To complete the analyses, we note that the formal extrapolation of LPT to the infinite order matches well (but not perfectly) the simulation results for all the configurations, as already found by STC18 in the 3D case. We also notice that the quasi1D approach of RF17 can describe very well the quasi1D configurations. Interestingly, at the second order in the transverse fluctuations considered here, the predictions are rather similar to those of standard fourthorder LPT, irrespective of initial conditions. Thus, as expected, the quasi1D approach becomes less accurate when the ratio ϵ_{2D} or ϵ_{3D} approaches unity, particularly in 3D as already shown by STC18.
Fig. 2 Phasespace structure for two and threesinewave initial conditions at the collapse time: Q1D2SIN (top left), ANI2SIN (centre left), SYM2SIN (bottom left), Q1D3SIN (top right), ANI3SIN (centre right), and SYM3SIN (bottom right). The intersection of the phasespace sheet with the y = 0 plane for two sine waves and y = z = 0 hyperplane for three sine waves is displayed in (x, υ_{x}) subspace. Simulation results are compared with standard LPT predictions, which are supplemented with the blue line, denoted by’ΈΧΤ’, which corresponds to the formal extrapolation to the infinite order proposed by STC18 and sketched in Sect. 4.1 for the collapse time. For completeness, the quasi1D approach (Rampf & Frisch 2017), denoted by Q1D, is also presented (see Appendix C for details). 
4.3 Radial Profiles
We now focus on radial profiles, and define Eulerian polar and spherical coordinates for two and threesinewave initial conditions by (x, y) = (rcosθ, rsinθ) and (x, y, z) = (rsinθcosϕ, rsinθ sinϕ, rcosθ), respectively. Then the angular averaged radial profiles of the density , velocity dispersion υ^{2}, radial velocity dispersion , and infall velocity υ_{r} are given by (35) (36) (37) (38)
where with r = x being the radial coordinate. In these equations, we used the angle average defined by (39)
for twosinewave initial conditions, and (40)
for threesinewave initial conditions. It is important to note that the angular coordinates θ and ϕ in the integrands are the Eulerian coordinates, and the integrands should be evaluated in terms of the Eulerian coordinate by solving the equation x = x(q, t).
To complete the analyses, we also study the pseudo phasespace density Q(r, t), defined by (41)
We additionally note that we present not only the radial profiles but also their logarithmic slopes, defined by (42)
where the quantity X stands for the radial profile under consideration.
Radial profiles at the shellcrossing are presented in Figs. 3–6. In these figures, measurements in the Vlasov code are performed by replacing each simplex of the phasespace sheet with a large ensemble of particles as explained in detail in C21 (see also Appendix B.4). Two approaches of collapse are considered. In Figs. 3 and 4, which respectively correspond to initial conditions with two and three sine waves, calculations are performed at the exact shellcrossing time, that is, the extrapolated value for LPT and the approximate value â_{sc} for the simulations. Figures 5 and 6 are analogous, but examine LPT predictions of the rath order synchronised to their own shellcrossing time , which allows us to examine in detail the structure of the singularity at collapse produced at each perturbation order. It is important to note that from now on, we do not further examine the quasi1D approach proposed by RF17. We also do not perform the formal extension to the infinite order (Eq. (34)), because calculating radial profiles in this framework was found to be too costly with the computational methods we are currently using.
The examination of these figures confirms the conclusions of the phasespace diagram analysis in Sect. 4.2. In particular, quasi1D profiles (Q1D2SIN and Q1D3SIN) are well described by LPT predictions even at low order. At the ‘exact’ collapse, considered in Figs. 3 and 4, LPT predictions generally accurately describe the outer part of the halo, where density contrasts are lower, and then deviate from the exact solution when the radius becomes smaller than some value r_{min}(n) decreasing with increasing n, where n is the perturbation order. LPT can describe arbitrarily large densities, as long as n is large enough. In the results presented here, density contrasts as large as 100 or larger can be probed accurately by LPT, depending on the order considered and on the nature of initial conditions. Again, convergence worsens when approaching 3D axialsymmetry, the worst case being, as expected, SYM3SIN with ϵ_{3D} = (1,1). It is worth noticing that velocity profiles require significantly higher LPT orders than the density profiles to achieve a comparable visual match between theory and measurements. In particular, deviations between LPT predictions and simulations arise at quite larger radii r_{min}(n) for the velocities than for the density, with differences as large as an order of magnitude. With synchronisation, performances of LPT predictions greatly improve, as expected and as illustrated by Figs. 5 and 6. All perturbation orders predict strikingly similar density profiles, in close to perfect match with the simulations, except for SYM3SIN discussed further below, while the velocity profiles still diverge slightly from each other except for quasi1D initial conditions, Q1D2SIN and Q1D3SIN.
We now turn to the logarithmic slope of various profiles shown in the bottom insert of each panel of Figs. 3–6. The nature of the singularities of CDM structures at the collapse time and beyond has been widely studied in the literature, particularly in the framework of Zel’dovich motion (e.g. Zel’dovich 1970; Arnold et al. 1982; Hidding et al. 2014; Feldbrugge et al. 2018). In Appendix D, we rederive the asymptotic properties of our protohaloes at small radii by Taylor expanding the equations of motion up to the third order in the Lagrangian coordinate q. For our sinewave configurations, whatever the order n of the perturbation order^{3}, three kinds of singularities are expected at the shellcrossing, as summarised in Table 2: the classic 1D pancake with a powerlaw profile at small radii of the form (S1) ρ(r) ∝ r^{−2/3}, valid for all the configurations expect SYM2SIN and SYM3SIN; then (S2) ρ(r) ∝ r^{−4/3} and (S3) ρ(r) ∝ r^{−2}, for SYM2SIN and SYM3SIN, respectively. On the other hand, velocities are expected to follow the same powerlaw pattern whatever initial conditions or dimensionality, with logarithmic slopes equal to 2/3, 2/3, and 1/3 respectively, for υ^{2}, , and υ_{r}, which in turn implies Q(r) ∝ r^{−7/3} for SYM2SIN, r^{−3} for SYM3SIN, and r^{−5/3} for other configurations.
Figures 3 and 4 consider LPT predictions for perturbation order n calculated at the exact theoretical collapse time and not at the individual collapse time at this perturbation order, so the shellcrossing is not reached exactly, but gets nearer as n increases. Consequently, the asymptotic slope is only approached approximately, and better so with larger n. We note that convergence is slower for velocities than for the density, especially for axialsymmetric configurations, in particular SYM3SIN. With synchronisation, as illustrated by Figs. 5 and 6, the convergence at a small radius to the prediction of singularity theory becomes clear.
Except for SYM3SIN, for which the best available snapshot is still too far from collapse, with â_{sc} > a_{sc}, one notices that all the simulated data show a close to perfect agreement between the measured slope at small radii and the predicted ones for all the radial profiles. This result is nontrivial in the sense that the Taylor expansion at small q underlying singularity theory is not necessarily valid in the actual fully nonlinear framework. Indeed, while singularity theory seems to apply to each order n for an arbitrarily large value of n in LPT framework, this does not mean that the limit to n → ∞ should converge to the same singularity. The fact that the simulation agrees with singularity theory predictions beyond the wellknown, but nonetheless somewhat trivial, 1D case, is remarkable even though naturally expected.
Several additional remarks are in order. First, we point out that in this framework if crossed sine waves are considered generic approximations of local peaks of a smooth random Gaussian field, the probability of having exactly or is null, and hence one expects the classic 1D pancake to be exclusively dominant, from a pure mathematical point of view. However, this can be unrealistic at the coarse level, where highdensity peaks can be close to spherical or filaments locally close to cylindrically symmetric. This is illustrated by Fig. 7, which examines the transition between various regimes in the vicinity of (top panels) and (bottom panels) for tenthorder LPT. Clearly, at sufficiently large radii, values of ϵ_{2D} or ϵ_{3D} close to and give similar results for the profiles, even for the logarithmic slope, which can approach that of the axialsymmetric case. The exact mathematical asymptotic behaviour is indeed reached only at very small radii. In practice, the axialsymmetric configurations cannot be ignored, which shows the limits of singularity theory if applied blindly.
Second, we note that the slope of the density profile predicted by LPT at collapse in the axialsymmetric 3D case, SYM3SIN, is different from the prediction of pure spherical collapse, ρ(r) ∝ r^{−12/7}, for an initial profile with the same asymptotic behaviour at a small radius of the form ρ^{ini} ∝ 1 − αr^{2} with α > 0 (see e.g. Nakamura 1985; Moutarde et al. 1995). Hence, the anisotropic nature of the axialsymmetric case (contained in terms beyond leading order r^{2} in the Taylor expansion of trigonometric functions) provides significantly different results from pure spherical collapse. We note, though, that this state of fact is not fully proved by our simulations measurements, which as mentioned above, are slightly beyond the actual collapse time.
Third, we can compare, in the 3D case, the logarithmic slope of the pseudo phasespace density with that seen after violent relaxation and at late stages of the evolution of dark matter haloes, Q(r) ~ r^{−β} (see e.g. Taylor & Navarro 2001; Navarro et al. 2010; Ludlow et al. 2010), which has been found to be close to the prediction of secondary infall model, β_{sph} = 1.875 (Bertschinger 1985). This is also the case of our three sine waves haloes, as analysed in detail by C21, including SYM3SIN, which shows again that this particularly symmetric setup, while being significantly more singular at collapse, with β = 3, relaxes to a state similar to more generic 3D configurations. The slope prior to collapse in the ‘generic’ cases,β = –5/3, although of the same order of β_{sph}, remains still different. Concerning the 2D case, which is different from the dynamical point of view, further analyses of the evolution of the haloes beyond collapse need to be done.
Fig. 3 Radial profiles and their logarithmic slopes for the twosinewave initial conditions at the shellcrossing time. As indicated in topleft panel, LPT predictions are given by solid lines of various colours and Vlasov simulations results are represented by red dots. From left to right, the initial conditions are given by Q1D2SIN (ϵ_{2D} = 1/6), ANI2SIN (ϵ_{2D} = 2/3), and SYM2SIN (ϵ_{2D} = 1), respectively. From top to bottom, we present the radial profiles of the normalised density , the velocity dispersion υ^{2}, the radial velocity dispersion , the infall velocity υ_{r}, and the pseudo phasespace density Q(r), respectively. We note that when plotting the logarithmic slopes in Vlasov simulations, we used the SavitzkyGolay filter implemented in savgol_filter of SciPy (Virtanen et al. 2020) to smooth the data for presentation purposes. In the logarithmic slope panels, the horizontal dashed lines correspond to the theoretical predictions at the small radii of Appendix D, as listed in Table 2. 
Fig. 4 Same as Fig. 3 but for the threesinewave initial conditions, from left to right, Q1D3SIN [ϵ_{3D} = (1/6,1/8)], ANI3SIN [ϵ_{3D} = (3/4,1/2)], and SYM3SIN [ϵ_{3D} = (1,1)], respectively. We note that in SYM3SIN, the closest snapshot from collapse we had at our disposal from our Vlasov runs, â_{sc} = 0.03190, is significantly beyond the actual shellcrossing time estimated by the method described in Sect. B.2, a_{sc} = 0.03155, which explains some discrepancies between the theory and the simulation at small radii. 
Fig. 5 Same as Fig. 3 but all analytical predictions evaluated at individual shellcrossing times computed at each perturbation order. 
Fig. 6 Same as Fig. 4 but all analytical predictions evaluated at individual shellcrossing times computed at each perturbation order. 
Summary of the logarithmic slopes at the shellcrossing obtained by singularity theory applied to LPT predictions at various orders (see Appendix D for details).
Fig. 7 Radial profiles at the shellcrossing using tenthorder LPT in the vicinity of (top panels) and (bottom panels). In the top panels, ϵ_{2D} is progressively varying in the range [0.9,1] with linear intervals of Δϵ = 0.1/64. In the bottom panels, ϵ_{3D} progressively changes from (0.97,0.97) to (1,1) with linear intervals of Δϵ = 0.03/64. From left to right, we present radial profiles of the normalised density , the velocity dispersion υ^{2}, the radial velocity dispersion the infall velocity –υ_{r}, and the pseudo phasespace density Q, respectively. In the lower inserts corresponding to logarithmic slopes, the horizontal dashed lines correspond to the values expected from singularity theory as computed in Appendix D and listed in Table 2. 
5 Structure of the System Shortly After Shellcrossing
Once the shellcrossing time is passed, the system enters into a highly nonlinear phase of violent relaxation. In this regime, Standard LPT is not applicable anymore, but it is still possible, shortly after collapse, to take into account the effects of the multistreaming flow on the force field using an LPT approach (see e.g. Colombi 2015; Taruya & Colombi 2017, for the introduction of ‘postcollapse’ perturbation theory). While the motion can still be described in a perturbative way for a short time after the shellcrossing, it is no longer, overall, a smooth function of time due to the formation of the singularity at the collapse time (Rampf et al. 2021). Still, it is reasonable to assume that very shortly after the shellcrossing, the backreaction corrections due to the multistreaming flow can be neglected, as a first approximation. However, it is important to notice that the convergence radius in time of the LPT series is finite in the generic case (e.g. Zheligovsky & Frisch 2014; Rampf et al. 2015; RH21). As illustrated by the measurement of the previous section, the convergence of the perturbation series is expected (although not proven) at least up to the collapse time, but not necessarily far beyond collapse (see also RH21). Therefore, in what follows, we shall use the LPT solution of the nth order at collapse as a starting point, and from there, the standard ballistic approximation, where velocity field is frozen, to study the structure of the system shortly after collapse.
In practice, one would like to use sufficiently large perturbation order so that convergence of LPT is achieved, and ideally the formal extrapolation to the infinite order proposed by STC18 and reintroduced in Sect. 4.1. However, this extrapolation is insufficiently accurate for the purpose of the calculations performed next, where very small time intervals after collapse are considered. Also, it is very costly from the computational point of view to exploit this extrapolation when solving the multivalued problem intrinsic to the multistreaming solutions. So from now on, we shall consider only higherorder LPT calculations and not their extrapolation to the infinite order (nor the quasi1D approximation proposed by RF17).
This section is organised as follows. In Sect. 5.1, we introduce the ballistic approximation. We also relate, in the multistream regime, the Eulerian density and velocity fields to their Lagrangian counterparts, and introduce, in this framework, the vorticity field. Indeed, after the shellcrossing, caustics form and nonzero vorticity is generated in the multistream region delimited by the caustics. In this section we also discuss the caustic network created up to the second order by our systems seeded by sine waves. Sect. 5.2 turns to actual comparisons of analytical predictions with Vlasov simulations using, similarly as in Sect. 4.2, phasespace diagrams, but shortly after collapseIn particular, we explore the limits of the ballistic approximation. Finally, Sect. 5.3 focuses on the overall structure of the multistream region in configuration space, by successively testing LPT predictions against simulations for the caustic pattern, the projected density and the vorticity fields.
5.1 Multistream Regime and Ballistic Approximation
Because the collapse time depends on the perturbation order, it only makes sense to test various perturbation orders shortly after the shellcrossing and compare them with simulations only if the collapse times are synchronised. In other words, in what follows, we consider the time for LPT of nth order and for the simulation, as listed in last column of Table 1. The value of Δa used for each sine waves configuration is given by the difference between the values of the last and sixth columns of the table. We recall that the quantity a_{sc} corresponds to the ‘true’ collapse time measured in the Vlasov runs as explained in Appendix B.2.
From the Eulerian coordinate and velocity fields given by nthorder LPT solutions at the shellcrossing time of each perturbation order, we model the Eulerian coordinate after collapse as follows: (43) (44)
while the nthorder velocity field is frozen to its value at the shellcrossing time of each perturbation order,
In the multistream region, given the Eulerian coordinate x, the solution of the equation x = x(q) has an odd number of solutions for the Lagrangian coordinate, which we denote as q_{F} labelled by the subscript F = [1,…, n_{F}]. Shortly after collapse, if ϵ_{i} ≠ ϵ_{j} for i ≠ j, there are at most three streams, n_{F} ≤ 3, in a given point of space. For symmetric cases, for example ϵ_{2D} = 1 in 2D or ϵ_{x} = ϵ_{y} in 3D, due to simultaneous shellcrossing along several axes, the number of streams can reach 9 and 27 in 2D and 3D, respectively, as discussed further in Sect. 5.2.
In what follows, we omit the time dependence for brevity. After defining the Lagrangian density and velocity fields, (45) (46)
the Eulerian density and velocity fields are expressed by the superpositions of each flow, (47) (48)
where the summation Σ_{F} is performed over all the solutions q_{F} of the equation x = x(q_{F}).
While vorticity is not expected in singlestream regions, unless already present in the initial conditions (and in that case it corresponds to a decaying mode), its generation is one of the prominent features of the shellcrossing (e.g. Pichon & Bernardeau 1999; Pueblas & Scoccimarro 2009). By taking the curl of the velocity field with respect to the Eulerian coordinate, the vorticity is given, respectively in 2D and 3D, by (49) (50)
The vorticity fields in 2D and 3D space are scalar and vector quantities, respectively. Because the acceleration derives from a gravitational potential, local vorticity on each individual fold of the phasespace sheet cancels. Only the nonlinear superposition of folds in the first term of Eqs. (49) and (50) induces nonzero vorticity, while the second term should not contribute. Using LPT however generates spurious vorticity in Eulerian space due to the truncation at a finite perturbation order (see e.g. Buchert & Ehlers 1993; Buchert 1994 and also Uhlemann et al. 2019 for a recent numerical investigation). This spurious vorticity is mainly coming from the second term in Eqs. (49) and (50), and can be especially noticeable in the singlestream region. We neglected it in our predictions by simply forcing this term to zero.
Before studying the structure of the system beyond collapse, we take the time to analytically investigate the properties of the ballistic solution for the sine waves case in the firstand secondorder LPT, in order to understand better the measurements performed in the next paragraphs. In particular, it is important to know how the solution behaves in the multistream region, of which the boundaries are given by the caustics, where J(q) = 0.
First, we start with 1LPT (Zel’dovich approximation). In the ballistic approximation, the Eulerian coordinate can be easily calculated shortly after the shellcrossing: (51)
with A = x, y, or z and stands for the growth factor at the shellcrossing time evaluated by using 1LPT. We note that in Eq. (51), the growth rate f and the expansion factor a_{sc} are evaluated at the same time as (we recall that f = 1 in the Einsteinde Sitter cosmology considered in this work). The value of is obtained in practice by solving the equation at the origin, keeping in mind that ϵ_{x}≥ϵ_{y,z}. We recall that Eq. (51) is exact in the pure ID case, that is, when ϵ_{y}= ϵ_{z} = 0. The condition for the caustics, J^{(1LPT)} (q,a) = 0, is reduced to the relation (52)
Equation (52) implies that the equations of the caustics are given by q = q_{0}, with q_{0} being a constant vector depending on time (i.e. on Δα). Therefore, the 1LPT caustics consist of ID lines in 2D and 2D planes in 3D, even at the collapse time. This configuration is actually degenerate and is not expected in realistic cases, for instance in the framework of a smooth Gaussian random field (see e.g. Pichon & Bernardeau 1999). Indeed, the regions where the first shellcrossing occurs should be composed, in the nondegenerate case (that is nonvanishing e_{;}), of a set of points. In the vicinity of such points and shortly after the shellcrossing, even in 1LPT, the equation of the caustics should correspond, at leading order in q and in the nonaxialsymmetric case, ϵ_{i}≠ϵ_{j} ,i ≠ j, to the equation of an ellipse in 2D and an ellipsoid in 3D (see e.g. Hidding et al. 2014; Feldbrugge et al. 2018). The reason for this not being the case here is due to the extremely restrictive class of initial conditions we have chosen.
Due to the degenerate nature of the sine waves initial conditions, it is therefore needed to go beyond the first order of the perturbative development of the dynamical equations to obtain a more realistic shape for the caustics. Indeed, for instance, 2LPT brings nonlinear couplings between various axes of the dynamics: (53)
where stands for the growth factor at the shellcrossing time evaluated by using 2LPT, which can be obtained as a function of ϵ_{2D} or ϵ_{3D} by solving the following secondorder polynomial: (54)
In Eq. (53), the growth rate f and the scale factor a_{sc} are evaluated at the same time as . The additional terms compared with Eq. (51) imply, shortly after collapse, nonvanishing contributions from all axes at the quadratic level (and above) in the Lagrangian coordinate q in the equation J^{(2LPT)}(q, α) = 0 (as long as none of the ϵ_{i} cancels). This leads to, in the nonaxialsymmetric case (i.e. Q1D2SIN, Q1D3SIN, ANI2SIN, and ANI3SIN), elliptic and ellipsoid shapes for the caustic curve/surface in Lagrangian space, respectively in 2D and 3D, as expected and as illustrated in the 3D case by the topleft panel of Fig. 8, while the bottomleft one shows the expected typical pancake shape in Eulerian space. Accordingly, in the analyses performed below, except for the phasespace diagrams that we examine first, LPT is examined from the second order and above. We note that, as will be studied below, the axialsymmetric cases ϵ_{2D} = 1 and ϵ_{3D} = (1,1) are a bit more complex, as they require Taylor expansion of the motion beyond the quadratic order in q to obtain a full description of the correct and more intricate topology of the caustic surfaces (e.g. the fourth order in q instead of second for SYM2SIN), as illustrated in 3D by the right panels of Fig. 8, due to simultaneous collapses along several axes. It would be beyond the scope of this article to go further into the analytic details of catastrophe theory for these very particular configurations, so we shall be content with a descriptive analysis of the LPT results in the axialsymmetric cases.
Fig. 8 3D view of the expected caustic pattern shortly after the shellcrossing for threesinewave initial conditions. The left panels correspond to the most typical pancake, here embodied by ANI3SIN, but Q1D3SIN would look similar, while the right ones show the axialsymmetric case, i.e. SYM3SIN. The calculations are performed using 2LPT along with ballistic approximation (Eq. (53)), but higherorder LPT would provide the same topology from the qualitative point of view, except that the position of the caustic surfaces would change, especially in the axialsymmetric case. 
Fig. 9 Tests of the ballistic approximation in Vlasov simulations: measured phasespace structure for twoand threesinewave initial conditions. This figure is analogous to Fig. 2, except that it shows a zoomedin view of the central part of the system shortly after collapse. We test the validity of the ballistic approximation by using two simulation snapshots slightly before (red curves, with as shown in Table 1 except for ϵ_{3D} = (1,1), for which a = 0.03103) and slightly after the collapse time (black, with a = a_{sc} + Δa). The ballistic approximation, as described in the main text, is applied to the red curves to obtain the magenta ones, to be compared directly with the exact solution, in black, given by the simulation. 
5.2 Phasespace Structure
Figures 9 and 10 display phasespace diagrams for our various sinewave systems shortly after collapse, namely the intersection of the phasespace sheet with the hyperplanes y = 0 and y = z = 0, for 2D and 3D configurations, respectively. The first figure examines the validity of the ballistic approximation directly from simulations data, while the second one compares predictions from LPT at various orders with the simulations.
Interestingly, Fig. 9 shows that the ballistic approximation can deviate significantly from the exact solution even quite shortly after collapse. While it remains precise for all configurations in 2D, with a slight worsening of the accuracy when going from the quasi1D to the axialsymmetric setup, as expected, significant or even dramatic deviations from the exact solution can be seen in 3D for the anisotropic configuration (ANI3SIN) and the axialsymmetric configuration (SYM3SIN). In the last case (bottomright panel), the intersection of the phasespace sheet with the hyperplane y = z = 0 is composed, in the multistream region, of three curves with an ‘S’ shape instead of a single one, as a result of simultaneous shellcrossing along all coordinates axes. The prediction from the ballistic approximation has then to be compared with the ‘S’ with the largest amplitude along the x axis. Similar arguments apply to the 2D case (bottomleft panel), but in this case, the intersection with the y = 0 hyperplane is composed of two curves instead of three^{4}.
The results shown in Fig. 9 must however be interpreted with caution. Strictly speaking, we aim to apply the ballistic approximation from the collapse time, which is not exactly the case in this figure. The red curves, which correspond to the snapshots used to originate the ballistic motion, all correspond to a time a_{bc} slightly before the actual collapse (namely a_{bc} = 0.03103 for SYM3SIN and , as given in Table 1 for other cases) because we do not have access to the exact moment of the actual collapse. The level of proximity to collapse is, at least from the visual point of view, the lowest for the axialsymmetric cases, particularly SYM3SIN. In this case, the ratio Δa_{eff} /a_{bc} ≃ 0.03 is the largest, where Δa_{eff} is the difference between a_{pc} ≡ a_{sc} + Δa and a_{bc} and is used to test the ballistic approximation. We also note that Δa_{eff} /a_{bc} ≃ 0.03 is also of the same order for ANI3SIN. The greater the value of Δa_{eff} /a_{bc}, the greater the deviation due to nonlinear corrections of the motions to be expected. Furthermore, we also expect the ballistic approximation to worsen from quasiID cases to fully axialsymmetric.
When examining Fig. 9 more in detail, we also note that for ANI3SIN, the central part of the ‘S ‘ shape is well reproduced by the ballistic approximation, while the tails underestimate velocities. This last defect is in fact present to some extent in all configurations except for Q1D cases, given the uncertainties on the measurements. In the 3D axialsymmetric case, not only the velocities in the tails are strongly underestimated, but also the magnitude of the position of the caustics along the x axis (local extremum of x coordinate). This suggests that the effects of acceleration in the vicinity of collapse are strong, which implies that the ballistic approximation can be applied only for a very short time (or it could be that it is not applicable, mathematically, due to the strength of the singularity).
Now we turn to the comparison of LPT predictions of various orders with the simulations, as shown in Fig. 10. We recall that the ballistic approximation is applied from the respective collapse times of each perturbation order, which obviously makes predictions of LPT artificially more realistic than it should be.
The conclusions of Sect. 4.2, where we compared LPT with simulations at collapse, still stand at the coarse level, obviously, since the time considered in Fig. 10 is nearly equal to the collapse time. When zooming on the central part of the system, we notice that the ballistic approximation applied to LPT works increasingly well with the order, as expected, especially when approaching quasiID dynamics. It can fail in the tails of the ‘S’ shape of the phasespace sheet, as a combination of limits of LPT to describe the system in the vicinity of the velocity extrema and the limits of the ballistic approximation itself just discussed above.
It is important to note at this point that the ballistic approximation is not necessarily the best choice outside the multistream region, where LPT predictions can still perform very well. This is illustrated in Fig. 10 by the cyan curves, which show tenthorder LPT results without assuming ballistic approximation; they should be compared with the dark purple curves, which correspond to tenthorder LPT with ballistic approximation. While the differences are very small given the time considered and the accuracy of the measurements, it can be seen, when just examining middle right panel of Fig. 10, that pure tenthorder LPT seems to perform better than its ballistic counterpart in the tails of the ‘S’ and nearby the local extrema of the velocity.
On the other hand, in bottom panels of Fig. 10, which correspond to axialsymmetric configurations, it is not clear at all that the cyan curves bring any improvement over the dark purple ones outside the multistream region. Here, LPT performs significantly more poorly than for other configurations, with or without assuming ballistic approximation. In this case, the topology of the phasespace structure in the multistream region (two or three curves according to the number of dimensions) is correct and the very central part of the shell corresponding to xaxis motion remains synchronised with the simulation, but except for this, the structure of the system is reproduced only qualitatively. This will be confirmed by the analyses that follow next.
Fig. 10 Phasespace structure for two and threesinewave initial conditions shortly after collapse. This figure is analogous to Fig. 2, except that we compare predictions of LPT at various orders in the ballistic approximation framework with measurements in Vlasov simulations shortly after collapse, i.e. for a = a_{sc} + Δa, as listed in Table 1 and discussed in the main text. We note that the additional cyan curve corresponds to the LPT prediction at the tenth order without using the ballistic approximation. 
Fig. 11 Caustic pattern shortly after the shellcrossing in the 2D case: comparison of LPT predictions using ballistic approximation with Vlasov runs. The left, middle, and right panels correspond, respectively, to Q1D2SIN, ANI2SIN, and SYM2SIN configurations, while the top and bottom panels show the caustics in Lagrangian and Eulerian spaces, respectively. 
5.3 Configuration Space: Caustics, Density, and Vrticity
We now enter into the heart of this work, which consists of examining in detail the structure of the system shortly after collapse in configuration space, both in Lagrangian and Eulerian coordinates. The multistream region is delimited by caustics, where density and vorticity are singular, so an accurate description of the caustic pattern represents the first test of LPT predictions. In most cases, the singularity corresponds to a simple fold of the phasespace sheet, with ρ(r) ~ r^{−1/2} (see e.g. Feldbrugge et al. 2018) and likewise for the vorticity along the direction orthogonal to the caustic (see e.g. Pichon & Bernardeau 1999). This means that a large part of the vorticity signal is generated in the vicinity of caustics. Having the correct shape of caustics is therefore fundamental because this information mostly determines the subsequent internal structure of the multistream region both for the density and the vorticity fields. Hence, in what follows, we first comment on Figs. 11 and 12, which compare, respectively in the 2D and 3D cases, the caustic pattern at various orders of LPT with the Vlasov code for our various sine waves initial conditions. Then, we turn to the normalised density and the vorticity in Figs. 13 to 18, and compare LPT predictions at the second and tenth order with the simulations. We remind the reader again that LPT predictions with the ballistic approximation are computed by synchronising the respective collapse times obtained at each perturbation order, which obviously artificially improves the performances of LPT. We cautiously note that in plotting the results below in 3D cases, the caustics outputs (Figs. 11 and 12) are shown at the intersections with the x = 0, y = 0, and z = 0 planes, while the density and vorticity outputs (Figs. 14 and 16–18) are shown in slices slightly away from the origin because the vorticity vanishes in the x = 0, y = 0, and z = 0 planes due to the symmetries in the chosen setting.
The examination of left four panels of Fig. 11 and top twelve panels of Fig. 12 provides us insights on the caustic pattern in the most typical configurations, ϵ_{i} ≠ ϵ_{j} i ≠ j, that is the nonaxialsymmetric case, both in 2D and in 3D. We note that in the 3D case, to have a clearer view of the caustic pattern, we consider the intersection of the caustic surfaces with the planes x = 0, y = 0 and z = 0 (we recall, however, that a 3D view of the caustic surfaces was given for 2LPT in Fig. 8). In the case of the simulations, as explained in Sect. 3.2 and Appendix B.3, the caustic pattern corresponds to a tessellation, that is, a set of ensembles of connected segments in 2D and of connected triangles in 3D. The intersection of a tessellation of triangles with a plane also gives a set of ensembles of connected segments. In the figures, only the vertices of these segments are represented, and the initial regular pattern of the tessellation used to represent the phasespace sheet is clearly visible in Lagrangian space.
Figures 11 and 12 confirm the results of the previous section, namely that highorder LPT provides a rather accurate description of the caustic pattern in the Q1D case, with a clear convergence to the simulation when the perturbation order increases. In the anisotropic cases, ANI2SIN and ANI3SIN, the match between LPT and simulations is not perfect, even at the tenth order, but remains still reasonably good, especially in 2D. For the axialsymmetric configurations, the convergence of LPT with order seems much slower, particularly in Eulerian space and in 3D.
While the shape of the caustic pattern is very simple in the generic cases, it becomes significantly more intricate in the axialsymmetric cases, particularly in Eulerian space. Indeed, the caustics are composed of two connected curves instead of one for ϵ_{2D} = 1, and three connected surfaces instead of one for ϵ_{3D} = (1,1). For ϵ_{2D} = 1, the shapes of simulation caustics are approximately reproduced by LPT, but we already see that convergence to the exact solution is slow. In particular, the extension of the outer caustic in Eulerian space is significantly underestimated by LPT, and it seems that increasing the perturbation order to arbitrary values is not going to improve the agreement with the simulations. The situation is much worse in 3D, particularly in Eulerian space, although the simulation measurements are very noisy, which makes comparison with theoretical predictions difficult. We know from the previous paragraph that these mismatches are at least partly attributable to limits of the ballistic approximation, along with limits of LPT to be able to describe the velocity field at collapse in the vicinity of its extrema, particularly the spike observed on the phasespace diagram in the 3D case.
Not surprisingly, the projected density measurements of Figs. 13 and 14 fully confirm the analyses of the caustic pattern. Putting aside the axialsymmetric cases, we can see that the agreement between simulations and LPT of high order is quite good, inside and outside the caustics. We note that the very high density contrasts observed in the axialsymmetric cases are due to the multiple foldings of the phasespace sheet, even in 2D, which explains the limits of the ballistic approximation, since feedback effects from the gravitational force field after collapse are expected to be orders of magnitude larger than in the generic case.
Turning to the vorticity, as shown in Figs. 15 to 18, we obtain again a good agreement between theory and measurements in the nonaxialsymmetric configurations, except that, at the exact location of the caustics, the simulation measurements are expected to be spurious. This is discussed in Appendix B.5 (see also end of Sect. 3.2), which provides details on the measurements of various fields, in particular how we exploit a quadratic description of the phasespace sheet to achieve vorticity measurements with unprecedented accuracy, but yet still limited by strong variations of the fields in the vicinity of the caustics, in particular the strong discontinuous transition between the inner part and the outer part of the multistream regions at the precise locations of the caustics.
From the qualitative point of view, the topology of the vorticity field for the generic configurations, ϵ_{i} ≠ϵ_{j}, i ≠ j, agrees perfectly with the predictions of Pichon & Bernardeau (1999) based on Zel’dovich dynamics^{5}. In 2D, the vorticity field is a scalar that can be decomposed into the four sectors inside the caustic, two of positive sign and two of negative one. In 3D, the vorticity field is a vector. Because the shellcrossing takes place along the x axis, each component of this vector field has specific properties, which are related to variations of the velocity and density field along with the phasespace sheet. In particular, it is easy to convince oneself that if collapse happens along the x axis, the strongest variations of all the fields are expected along this direction. That means, since each coordinate of the vorticity vector depends on variations in the velocity field in other coordinates, that the magnitude of ω_{x} is expected to be small. Due to symmetries, we also expect that ω_{y} and ω_{z} are, respectively, an odd function of z and y, which means that ω_{y} approaches zero when approaching the z = 0 plane, and similarly for ω_{z} when approaching the y = 0 plane, which explains the pattern of the vorticity field in Figs. 16 and 17. These symmetries impose us to perform measurements on slices slightly shifted from the centre of the system, as detailed in the caption of Fig. 14. We also note that in these figures the pattern is analogous to 2D for the y and z coordinate of the vorticity fields in the x−z and x−y subspace, respectively.
Turning to the axialsymmetric case, the agreement between theory and measurements is again only partial. Yet the high magnitude of vorticity is of the same order for theory and measurements. The structure of the vorticity field, significantly more complex than in the generic case due to multiple foldings of the phasespace sheet, is qualitatively in agreement between LPT and Vlasov code in 2D, while simulations measurements are too noisy in the 3D case to make definitive conclusions. What is clear, however, is that the outer part of the multistream region seems to be qualitatively reproduced by the theory, but its size is totally wrong.
Fig. 12 Structure of caustics shortly after the shellcrossing for threesinewave initial conditions: comparison of LPT predictions using ballistic approximation with Vlasov runs. The top six, middle six, and bottom two panels respectively correspond to Q1D3SIN, ANI3SIN, and SYM3SIN. In each group of six panels, the top and bottom lines correspond to Lagrangian and Eulerian spaces, and the intersection of the caustic surfaces with the plane are q_{x} = 0, q_{y} = 0, and q_{z} = 0, respectively, for the first, second, and third column of each group. In the bottom group of two panels, the left and right panels correspond respectively to Lagrangian and Eulerian space and only show the intersection of the caustic surface network with the plane q_{z} = 0 due to the symmetry of the system. 
Fig. 13 2D density shortly after the shellcrossing: comparison of LPT predictions using ballistic approximation with Vlasov runs. From top to bottom: Q1D2SIN, ANI2SIN, and SYM 2SIN. From left to right: 2LPT, 10LPT, and measurements in Vlasov simulations. 
Fig. 14 Slices of the projected density shortly after the shellcrossing: comparison of LPT using ballistic approximation with Vlasov runs. The left and right groups of nine panels correspond respectively to Q1D3SIN (top, middle, and bottom line of panels: x=−1.55×l0^{−4}, y = −1.17×10^{−3}, and z= −1.56 × 10^{−3} slice) and ANI3SIN (x= −5.16 × 10^{−4}, y = −2.15 × 10^{−4}, and z= −5.47 × 10^{−4} slice), while the bottom group of three panels corresponds to SYM3SIN (z= −1.17 × 10^{−5} slice). In each group of panels, the left, middle, and right columns give respectively the secondorder LPT prediction, the tenthorder LPT prediction, and the Vlasov code measurements. Due to the symmetry of the system for SYM3SIN, only one slice is shown for the bottom panels. 
6 Summary
In this paper, following in the footsteps of Saga et al. (2018), we have investigated the structure of primordial CDM haloes seeded by two or three crossed sine waves of various amplitudes at and shortly after the shellcrossing, by thoroughly comparing LPT, up to the tenth order, with highresolution VlasovPoisson simulations performed with the public Vlasov solver ColDICE. We devoted our attention first to the phasespace structure and radial profiles of the density and velocities at the shellcrossing and, second, to the phasespace structure, caustics, and density and vorticity fields shortly after the shellcrossing. In particular, measurements of unprecedented accuracy of the vorticity in the simulations were made possible by exploiting the fact that ColDICE is able to follow the phasespace sheet structure locally at the quadratic level.
We studied three qualitatively different initial conditions characterised by the amplitude of three crossed sine waves, as summarised in Table 1: quasi 1D (Q1D2SIN, Q1D3SIN), where one amplitude of the sine waves dominates over the other(s), anisotropic (ANI2SIN, ANI3SIN), where the amplitude of each wave is different but remains of the same order, and axialsymmetric (SYM2SIN, SYM3SIN), where all amplitudes are the same. In order to predict the protohalo structure shortly after the shellcrossing in an analytical way, we used the ballistic approximation, where the acceleration is neglected after the shellcrossing time.
Our main findings can be summarised as follows: – Phasespace diagrams at collapse: Except for SYM3SIN, one expects the system to build up a classic pancake singularity at the shellcrossing, with a phasespace structure along the x axis analogous to what is obtained in 1D. This pancake is reproduced by LPT well, which converges increasingly well to the exact solution when approaching quasi1D dynamics, as expected. The local extrema of the velocity field around the singularity are the locations where LPT differs most from the exact solution, underestimating the magnitude of the velocity. Convergence with the perturbation order becomes very slow when approaching 3D axialsymmetry, where spikes appear on the velocity field on each side of the singularity;
Phasespace diagrams at collapse: Except for SYM3SIN, one expects the system to build up a classic pancake singularity at the shellcrossing, with a phasespace structure along the x axis analogous to what is obtained in 1D. This pancake is reproduced by LPT well, which converges increasingly well to the exact solution when approaching quasi1D dynamics, as expected. The local extrema of the velocity field around the singularity are the locations where LPT differs most from the exact solution, underestimating the magnitude of the velocity. Convergence with the perturbation order becomes very slow when approaching 3Daxialsymmetry, where spikes appear on the velocity field on each side of the singularity;
Radial profiles at collapse: With a sufficiently high order of perturbation, LPT can reproduce arbitrarily high density contrasts at collapse, but we note a slower convergence when turning to velocity profiles. Still, the convergence of LPT is sufficiently good to probe the asymptotic logarithmic slope at the small radii expected at collapse for various profiles from singularity theory, as summarised in Table 2 and confirmed by simulation measurements. These profiles are ρ(r) ∝ r^{−2/3} for generic initial conditions (i.e. Q1D and ANI), ρ(r) ∝ r^{−4/3} for SYM2SIN, and ρ(r) ∝ r^{−2} for SYM3SIN, while velocity profiles always present the same powerlaw behaviour, for example υ^{2}(r) ∝ r^{2/3}. Confirming predictions of singularity theory at collapse and explicit convergence to asymptotic profiles at small radii was expected but is not trivial;
Synchronisation: Agreement with singularity theory predictions is made even more explicit by computing LPT profiles at the respective collapse times computed at each perturbation order. In this case, convergence to the powerlaw profile predicted at small radii by singularity theory is explicit. In fact, the radial profiles obtained from LPT with such synchronisation are strikingly alike, particularly for the density, which motivates the use of a ballistic approximation to study the motion slightly beyond collapse;
Ballistic approximation: Measurements in simulations show that the ballistic approximation provides an accurate description of the phasespace structure of the system slightly beyond collapse, except in the axialsymmetric case SYM3SIN, where the very highly contrasted nature of the system due to multiple superpositions of phasespace sheet folds introduces strong force feedback effects. We notice, however, that even in the generic pancake case, the ballistic approximation can still slightly underestimate velocities in the vicinity of the singularity. Obviously, these conclusions depend on how long this approximation is used. In this work, we considered maximum relative variations of the expansion factor of the order of three percent. For these values, the deviations from pure ballistic motion were the most significant, as expected;
Structure beyond collapse: Given the limits discussed above, we find, when considering generic, nonaxialsymmetric configurations, that LPT of sufficient order combined with the ballistic approximation provides a very good description of the structure of the system beyond the collapse time, with excellent agreement between the predicted density and the vorticity fields inside the multistream region and those measured in the Vlasov simulations. Turning to axialsymmetric cases, this agreement is only obtained at the qualitative level, even for tenthorder LPT, particularly for SYM3SIN, where the size of the outer caustics is strongly underestimated.
Obviously, the ballistic approximation is only the first step in a more complete calculation that would take into account the feedback due to the force field, as first proposed in the 1D case by Colombi (2015) and Taruya & Colombi (2017) (see also the recent work of Rampf et al. 2021), and formulated in terms of postcollapse perturbation theory. The next step is to implement an LPT approach where the small parameter is the interval of time Δa = a − a_{sc} from the collapse time and where a Taylor expansion of the phasespace sheet is performed around the singularity in terms of Lagrangian coordinates. Shortly after collapse, in the generic, nonaxialsymmetric case, the system presents an ‘S’ shape in x – υ_{x} subspace, similar to the 1D case where positions and velocities can be approximated as thirdorder polynomials of the Lagrangian coordinate q. The calculation of the force field requires a threevalue problem to be solved in the multistream region, similar to the 1D case. While technically challenging, generalisation of postcollapse perturbation theory to 2D and 3D seems possible. It might provide significant insights on nonlinear corrections due to multistream dynamics on largescale structure statistics, for example predictions of the power spectrum of the largescale matter distribution from higherorder perturbation theory.
Fig. 15 2D vorticity fields: comparison of LPT predictions with Vlasov runs. 
Fig. 16 Vorticity components shortly after the shellcrossing: comparison of LPT with the ballistic approximation with Vlasov runs for Q1D3SIN. The 2D slices are the same as those shown in the topleft group of nine panels of Fig. 14, except that the slice changes from left to right and the vorticity component changes from top to bottom. Again, on each line of three panels, 2LPT (left panel) and 1OLPT (middle panel) are compared with simulation measurements (rightpanel). 
Fig. 17 Same as Fig. 16, but for ANI3SIN. The 2D slices considered are the same as in the topright group of nine panels of Fig. 14. 
Fig. 18 Same as Fig. 16, but for SYM3SIN. Due to the symmetry of the initial conditions, only an xy slice is shown, which is the same as in bottom three panels of Fig. 14. 
Acknowledgements
This work was supported in part by JSPS Research Fellow GrantinAid Number 17J10553 (SS), MEXT/JSPS KAKENHI Grants Numbers JP17H06359, JP20H05861, and JP21H01081, JST AIP Acceleration Research grant JP20317829 (AT), ANR grant ANR13MONU0003 (SC), as well as Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, cofunded by CEA and CNES (SC). Numerical computation with ColDICE was carried out using the HPC resources of CINES (Occigen supercomputer) under the GENCI allocations c2016047568, 2017A0040407568 and 2018A0040407568. Posttreatment of ColDICE data were performed on HORIZON cluster of Institut d’Astrophysique de Paris. Calculations of theoretical predictions have made use of the Yukawa Institute Computer Facility.
Appendix A Expressions of the LPT Solutions
In this appendix, we present the LPT solutions up to the fifth order, which are obtained by solving the recursion relations given in Eqs. (21) and (22). Since higherorder solutions are straightforwardly derived in the same way, we do not explicitly show them here.
The results are partly presented in Moutarde et al. (1991) up to the third order, taking into account the contributions other than the fastest growth mode (see also Buchert et al. 1997, for the solutions including decaying modes). Here, we first explicitly show the analytical solutions of the sine waves initial conditions up to the fifth order.
As shown in Sect. 2.1, the fastest growing mode can be expanded as follows: (A.1)
For presentation purposes, we only show the xcomponents of the displacement field, that is, with the definition (ϵ_{1}, ϵ_{2}) ≡ (ϵ_{y}/ϵ_{x}, ϵ_{z}/ϵ_{x}). Given the xcomponents, the y and zcomponents can be derived by permutating all x, y, and z: (A.2) (A.3) (A.4) (A.5) (A.6)
We note that because of the underlying symmetry of the initial conditions, the xcomponents of the LPT solutions are symmetric under the exchange of y ↔ z and ϵ_{1} ↔ ϵ_{2}.
Appendix B Measurements on the Tessellation of Vlasov Simulations
This appendix explains in detail how measurements of various quantities studied in this article are performed on the output of ColDICE, which consists of a tessellation of the phasespace sheet with simplices, respectively triangles and tetrahedra in 4D and 6D phasespace.
B.1 Phasespace Diagrams
The intersection of a hypersurface of dimension D_{1} = D with a hyperplane P of dimension D_{2} = D + 1 inside a phasespace of dimension D_{3} = 2D is of dimension D_{2} + D_{2} − D_{3} = 1. In other words, this intersection corresponds in the nontrivial case to a set of curves. Therefore, in 2D, the intersection of the phasespace sheet with the hyperplane x = 0 is a set of curves, and likewise in 3D, the intersection of the phasespace sheet with the hyperplane x = y = 0. Furthermore, as discussed more in detail in C21, because the phasespace sheet remains at all times a fully connected hypersurface, this set of curves should also be fully connected, with no hanging point.
To produce a phasespace diagram, we employed an approach valid at the linear order, consisting in simply computing, for each simplex, its geometric intersection with the hyperplane P, which is empty, a point, or a segment. As a result, phasespace diagrams extracted from the Vlasov runs consist of sets of connected segments, the ends of which are plotted in Figs. 2 and 10.
B.2 Expansion Factor at Collapse: a_{sc}
An accurate estimate of the value of the expansion factor corresponding to the collapse time is essential when studying the properties of the system at the shellcrossing and shortly after. The curves generated in Appendix B.1 can be used for this purpose, at two expansion times a_{i}, i = 1, 2, supposed to be just before and just after collapse. At both these times, we consider the unique phasespace diagram segment portion S, one of the ends of which is as close as possible to the origin and the other, with abscissa x_{i}, has υ_{x} > 0; the magnitude of the other coordinate(s) of the velocity, in theory null, are as small as possible. The latter condition excludes, in the axialsymmetric case, the component(s) of the flow corresponding to simultaneous collapse along the y or/and z direction, which has/have a nonzero value of υ_{y} or/and υ_{z}. Once this segment is identified at both times, corresponding to a_{1} and a_{2}, a simple linear interpolation is used to estimate the expansion factor at collapse: (B.1)
This formula does not require ai and a2 to be just before and after the collapse time to estimate a_{sc}, but, to provide sufficiently accurate results, the times need to be sufficiently close to the actual collapse time. A first guess of the collapse time is estimated by using perturbation theory predictions extrapolated to the infinite order (STC18), as listed in 5th column of Table 1. Examination of Table 1 shows that the value predicted this way agrees extremely well with the measurements in the Vlasov runs provided by Eq. (B.1), and shown in sixth column of the table.
We note, however, that, beyond the approximate nature of the linear approximation underlying Eq. (B.1), the collapse time estimate  in addition to other estimates, as discussed in Sect. 3.2 can be significantly influenced by the force resolution, that is, the value of n_{g}. As discussed in detail in C21, decreasing force resolution delays the collapse time; hence, to have an accurate estimate of the collapse time, a sufficiently large value of n_{g} is required. We performed extensive force resolution tests for the 3D simulation with ϵ_{3D} = (3/4,1/2). Our measurements of a_{sc} with the above method provide a_{sc} = 0.02912 (in excellent agreement with the predicted value and 0.0293 respectively for n_{g} = 1024, 512 and 256. Our n_{g} = 512 simulation thus provides, for this value of ϵ_{3D}, an estimate of the collapse expansion factor accurate at approximately the 10^{−4} level, and we expect this to apply as well to the quasi1D case ϵ_{3D}= (1/6,1/8). However, for the axialsymmetric case, ϵ_{3D} =(1, 1), the measured value might still overestimate the actual one by an amount larger than ~ 10^{−4} due to the high strength of the singularity building up at the centre of the system. On the contrary, in the 2D case, given the high value of n_{g} = 2048 used to perform the simulations, we expect our estimates of collapse time to be very accurate, at an order better than the 10^{−4} level, but we did not test this explicitly.
B.3
Caustics are regions where the determinant J of the Jacobian matrix changes sign. At the linear order in the local description of the phasespace sheet, geometrically this means that the orientation of the simplex changes in configuration space, which allows one to define unambiguously regions corresponding to the intersections of simplices with J ≥ 0 and simplices with J ≤ 0, where the sign of J is directly estimated from the current orientation of the simplex with respect to the original one. This is actually performed during runtime by ColDICE, which can output caustics directly when needed. In 2D, the phasespace sheet is composed of a tessellation of triangles, and hence the caustics estimated this way are given by sets of segments, the ends of which are shown in Fig. 11. In 3D, the caustics are given by sets of triangles, the intersection of which we compute with the y = z = 0 plane to again get a set of segments; their extremities are shown in Fig. 12. Because we are using a leading order approach, the caustic lines or surfaces are not necessarily smooth but should trace accurately enough the actual caustics for the purpose of this work (see also Sect. 3.2).
B.4 Radial Profiles
To measure radial profiles in logarithmic bins, each simplex is replaced with a large number of particles as explained for the 3D case in Appendix A2 of C21. We refer to this work for the reader interested in the details of this procedure, which we straightforwardly generalised to the 2D case.
B.5 Density Field and Vorticity: From Linear to Quadratic Order
In this section we aim to compute the following five quantities: (i) the Jacobian of the transformation between initial and final positions, J(q), as defined in Eq. (7); (ii) the Lagrangian projected density, ρ_{L}(q), defined in Eq. (45); (iii) the total Eulerian projected density, which stems from the superposition of one or more folds of the phasespace sheet, as described by Eq. (47); (iv) the Eulerian velocity field, given by Eq. (48); and (v) the vorticity, ω^{2D} and ω^{3D}, defined in Eqs. (49) and (50). Again, because the acceleration derives from a potential, local vorticity on each phasespace sheet fold cancels. Only the nonlinear superposition of phasespace sheet folds contained in the first term of these equations induces nonzero vorticity, while the second term should not contribute, although we shall take it into account in our numerical calculations.
In what follows, we explain how to compute the various quantities defined above by exploiting the decomposition in simplices of the phasespace sheet by ColDICE. In Sect. B.5.1 we introduce barycentric coordinates, which are useful to define the position of any point inside each simplex. In Sect. B.5.2 we show how the barycentric coordinates can be used to perform calculations at the linear level inside each simplex, in particular partial derivatives of a function, as already discussed for instance by Hahn et al. (2015). In Sect. B.5.3 we generalise the procedure to the case when a quadratic description of the phasespace sheet is available. Finally, in Sect. B.5.4 we describe the way we sample various quantities described above on an highresolution Cartesian grid, by using proper sets of sampling particles associated with each simplex.
B.5.1 Barycentric Coordinates
Since the phasespace sheet is sampled with simplices, it is useful to define a wellknown local system of coordinates on each simplex. Given the positions X_{k}, k = 1,…, N_{s} of the simplex vertices, with N_{s} = D + 1 where D is the dimension of configuration space, and a function g(X), one can define the following linear interpolation: (B.2)
where the barycentric coordinates ξ_{k} are positive quantities verifying ξ_{k} = 1, ξ_{k’≠k} = 0, for X = X_{k} and Σ_{k} ξ_{k} = 1. When working in Lagrangian space, the space of initial positions, the phasespace sheet is flat and the linear interpolation, (B.4)
is exact. In what follows, we shall therefore use Eq. (B.4) to define barycentric coordinates. With this definition of this means that at given time t, Eq. (B.3) remains valid, but only at the linear level, since dynamical evolution of the phasespace sheet produces curvature.
B.5.2 Calculations of Various Quantities at the Linear Level
With the barycentric coordinate representation, any function defined at vertices positions can be estimated locally at the linear level using Eq. (B.2). The derivative of the function inside each simplex can thus be written (B.5)
where M is the matrix of the partial derivatives of the vector P ≡ (Σ_{k} ξ_{k}, X1, X_{2}, X_{3}) with respect to ξ_{k}, with X given by Eq. (B.3). In other words, (B.6) (B.7)
with X_{k,α} the α^{th} coordinate of k^{th} vertex with position X_{k}. At the linear level, the derivative defined by Eq. (B.5) is therefore simply constant inside each simplex.
This calculation can be performed in Lagrangian space to directly obtain the Jacobian of the transformation between the initial and present position, Eq. (7), from the Jacobian matrix, (B.8)
where W is defined similarly as M but in Lagrangian space, (B.9) (B.10)
with Q_{k,α} the α^{th} Lagrangian coordinate of k^{th} vertex with Lagrangian position Q_{k}.
This means that the linear description allows one to obtain the Jacobian J(q) and the corresponding Lagrangian density ρ_{L}(q) = 1/J(q) as constant quantities inside each simplex. Unfortunately, this is not sufficient, since Eqs. (49) and (50) involve partial derivatives of the projected density. Rigorously speaking, this means that a representation of the simplices at the quadratic level is required.
B.5.3 Calculations of Various Quantities at the Quadratic Level
Fortunately, ColDICE employs a quadratic description of the phasespace sheet for performing local refinement of the tessellation, by using tracers defined as midpoints of edges of the simplices in Lagrangian space. For instance, the tracer associated with vertices (k, k’) corresponds to barycentric coordinates ξ_{k} = ξ_{k’} = 1/2, ξ_{j}, = 0 for j ≠ k and j ≠ k’.
This allows one to define in a unique way a quadratic hypersurface inside each simplex coinciding with vertices and tracers belonging to it. We note that the secondorder representation preserves the continuity of the phasespace sheet but does not warrant its smoothness at the differential level. Exploiting this representation will nevertheless provide, in practice, sufficiently accurate measurements of the vorticity.
In the quadratic representation, one defines the following conventional shape functions N_{k}, (B.11) (B.12)
where functions K(k) and L(k) are appropriately chosen to cover all the combinations 1 ≤ K(k) ≤ L(k) ≤ Ns, and where Ns = D + is the number of simplices and Nt = D(D + 1)/2 is the number of tracers, given D the dimension of space. Then Eq. (B.2) becomes (B.13)
where g_{k} is now also defined over tracers. This is true in particular for quadratically interpolated positions, (B.14)
while corresponding Lagrangian coordinates Q are still given by the linear representation (B.4) since the phasespace sheet is initially flat, which means that matrix W does not change.
From Eq. (B.14), we obtain the matrix M^{quad} of partial derivatives of the vector P ≡ (Σ_{k} ξ_{k}, X_{1}, X_{2},X_{3}) with respect to ξ_{k}, which replaces Eqs. (B.6) and (B.7): (B.15) (B.16)
with X_{k’,α} the α^{th} coordinate of k’^{th} vertex/tracer with position X_{k’}.While, at the linear level, this matrix was constant inside each simplex, it now depends linearly on barycentric coordinates, given the quadratic nature of the shape functions N_{k}.
The calculation of the Jacobian matrix exploiting the quadratic representation then just consists in replacing matrix M with M^{quad} in Eq. (B.8), since matrix W remains unchanged, which allows us to compute J(q) in the quadratic representation. In particular, we can compute J_{k} = J(Q_{k}) at each vertex and tracer position to prepare it for the quadratic interpolation procedure given by Eq. (B.13) (but staying aware of the fact that this calculation is accurate only up to the linear order), which, in addition to the phasespace coordinates themselves, is enough to compute all the quantities we are interested in.
Once we have a scalar function g(x) defined at vertices and tracers positions, its derivative is given from Eq. (B.13) by (B.17) (B.18)
This quantity depends in a nonlinear way on barycentric position and requires matrix M^{quad} to be invertible, so to be able to have a finite estimate of the partial derivative, it is necessary to avoid the subspace occupied by the caustics (or, more exactly, the element of surface or curve approximating the actual caustics).
In particular, we estimate the partial derivatives of the projected density as follows (B.19)
using Eq. (B.18) to estimate J/ xß. Indeed, the Jacobian is a smooth function of the Lagrangian coordinate. This is not the case for function ρ_{L}(q) = 1/J(q), which presents strong variations in the vicinity of caustics. Using the Jacobian instead of the density to estimate derivatives allows us to capture better asymptotic behaviours in the vicinity of the caustics.
B.5.4 Practical Estimate: Particle Representation and Local Coarse Graining
In order to compute any quantity at position x_{0}, one has to compute the intersection of the phase space sheet with the position x = x_{0}, or to resolve inside each simplex the equation Χ(ξ) = x_{0}. While this is a straightforward linear problem in the linear representation of the phasespace sheet (Eq. B.3), it becomes nontrivial in the quadratic case (Eq. B.14), for which some iterative procedure seems necessary. Here, to simplify the approach (and also to reduce effects of aliasing and divergences near caustics), we adopt a forward point of view consisting in projecting the tessellation on a Cartesian mesh of resolution n_{ana} in configuration space, similarly as it is performed by ColDICE to compute the projected density. This means that we apply additional coarse graining over each pixel or voxel of the Cartesian mesh: Integrals of the form (B.20)
with υ_{pixel} and υ_{voxel} being respectively the area of the pixel and the volume of voxel under consideration.
While Sousbie & Colombi (2016) use actual raytracing exact to the linear order to compute such an integral, we replaced each simplex with a large number of particles, which allows us to give an account of the quadratic shape of the simplices in a very simple way. However, employing particles introduces some discrete noise. To reduce the noise, we used the cloudincell (CIC) interpolation procedure on the target Cartesian mesh (Hockney & Eastwood 1988), that is to say, each particle is replaced by a pixel (voxel) of the same small size as the pixels (voxels) of the Cartesian mesh, and we associate a weight W_{inter} with the particle proportional to the volume of intersection between the pixel (voxel) associated with the particle and the pixel (voxel) of the mesh. In addition, we do not place particles at random within the simplex, but instead refine the simplex hierarchically in an homogeneous fashion ℓ_{max} times, and then replace each subsimplex obtained this way with a particle corresponding to the centre of the simplex in the barycentric representation. The calculation of ℓ_{max} is such that discreteness effects due to the projection inside each voxel or pixel of the Cartesian mesh are kept under control: (B.22)
where D is the dimension of space, Δ is the step of the Cartesian mesh and (B.23)
with X_{k,α} the α^{th} coordinate of k^{th} vertex with position X_{k} and likewise for X_{k’,α}.
In practice, prior to CIC interpolation, each particle p samples a small element of volume δV_{p} in the spatial part of the integral of right hand side of Eq. (B.21), (B.24)
where V_{L} is the Lagrangian volume of the simplex and Jp the Jacobian (quadratically) interpolated at position of the particle. So, at the end, the contribution of each particle in right hand integral of Eq. (B.21) is equal to δV_{p} W_{inter} g(X_{p}, V_{p}), where we recall that W_{inter} is the CIC weight defined previously and where function g is estimated at particle phasespace position (X_{p}, V_{p}) using the methods described in previous paragraphs.
We note that coarse graining due to the CIC interpolation is expected to introduce defects or biases in regions where large variations are expected, that is, in the vicinity of caustics. Obviously these effects become more dramatic when differentiating quantities, so one expects the vorticity to be affected numerically nearby the caustics, as can be noted in Figs. 15, 16, 17 and 18.
Appendix C Perturbative Treatment of Quasi1D Flow
In Sect. 4.2 we test theoretical predictions relying on the perturbative treatment of a quasi1D flow as proposed by RF17. In this appendix, we briefly present this approach, focusing mainly on threesinewave initialcondition configurations. Generalisation to the two sine waves case is straightforward.
When the system is exactly 1D, the firstorder LPT solution is exact before the shellcrossing. This fact leads to an approximate treatment in 3D space. We exploit the exact solution of the 1D problem along the x axis given by Zel’dovich approximation as the unperturbed zerothorder state: (C.1)
with the subscript i = x, y for two sine waves and i = x, y, z for three sine waves. The transverse fluctuations are considered as small firstorder perturbations to this setup: (C.2)
In the case of the three sine waves, we assume the explicit form of the functions and Ψ^{(1,1)} to be (C.3) (C.4)
Then, substituting these initial conditions into the recursion relations in Eqs. (21) and (22), we obtain (C.5) (C.6)
By solving Eqs. (C.5) and (C.6), one can construct the Q1D firstorder solutions.
We further extend the Q1D treatment proposed by RF17 up to the second order in the transverse fluctuations as described below. On top of the zeroth and firstorder perturbations, we define the secondorder perturbation as (C.7)
Again, using Eqs. (21) and (22), we obtain the following recursion relation: (C.8) (C.9)
The recursion is initialised by (C.10) (C.11)
Using the above recursion relations, the perturbative expansion is then performed by keeping terms proportional to and up to the second order, n + m = 2, while keeping terms proportional to ϵ_{x}^{k} up to the tenth order, k = 10, as shown in Fig. 2. This approach provides a very accurate description of the dynamics when the system is initially quasi 1D (i.e. ϵ_{x} » ϵ_{y,z}).
The results of the Q1D solution at the first order in the transverse fluctuations are presented in RF17. In what follows we show (for the first time) the Q1D solution at the second order in the transverse fluctuations and up to the fifth order in the longitudinal direction in terms of the growth factor (i.e. Ψ^{(2,5)}). We use the same notations as in Appendix A.
For the xcomponents of the Q1D solutions, we have (C12) (C13) (C14) (C15) (C16) (C17) (C18) (C19)
We note that the ycomponents and zcomponents are symmetric under the exchange of q_{y} ↔ q_{z} and ϵ_{y} ↔ ϵ_{z} in the above expressions.
Appendix D Asymptotic Structure of the Singularity at Collapse
In this appendix, we analytically investigate the dependence on initial conditions of the logarithmic slopes of various profiles expected at the shellcrossing and small radii, as discussed in Sect. 4.3. The main results of this investigation are summarised in Table 2. Since the extension to the 2D case is obvious by simply setting ϵ_{z} = 0, we present here the structure of the singularity obtained from LPT to some order for the three sine waves case. We note that the calculations presented in this appendix are conceptually not new since singularity theory applied to cosmology is already well known (see e.g. Arnold et al. 1982; Shandarin & Zeldovich 1989; Hidding et al. 2014; Feldbrugge et al. 2018).
D.1 Taylor Expansion Around the Singularity
To facilitate the analysis, we focus on the relation between Eulerian and Lagrangian coordinates around the origin. Expanding the relation (5) between the Eulerian coordinate x and the Lagrangian coordinate q around the origin, we obtain, after neglecting O(q^{4}) and higherorder terms, (D.1)
with Ã_{ij}(t) ≡ δ_{ij} + A_{ij}(t), where A_{ij}(t) is some function of time. In these equations, we have exploited the symmetric nature of the threesinewave initial conditions, which imply that Eulerian coordinates around the shellcrossing point can be expanded in terms of odd thirdorder polynomials of the Lagrangian coordinates, that is, polynomial forms P that verify P(−q) = −P(q) (see e.g. Colombi 2015; Taruya & Colombi 2017, for the 1D case).
The coefficients A_{ij} and C_{ijkl} are expressed in terms of partial derivatives of the displacement field with respect to the Lagrangian coordinates at the origin as follows, (D.2)
Using Eq. (D.1), the Jacobian matrix is given by (D.3)
The shellcrossing time t_{sc} is obtained by solving the equation J(0, t_{sc}) = 0, that is, (D.4)
Hereafter, the dependence on the collapse time t_{sc} will be omitted in the notations.
Using Eq. (D.3), the Jacobian is given by (D.5)
Further simplifying the calculations, we noticed that, up to the tenth order of the LPT development, for an initial displacement field given by three orthogonal sine waves aligned with each coordinate axis, the matrix is diagonal and the matrix C_{ijkl} verifies (D.6)
We did not find any simple way to demonstrate that these properties stand at any order, but the fact that they are verified up to tenthorder LPT is a really convincing clue. In this case, the shellcrossing condition, Eq. (D.4), is reduced to Ã_{xx} Ã_{yy} Ã_{zz} = 0. Prior to the shellcrossing, the Eulerian coordinates read (D.7)
for A, B = x, y, z, and the Jacobian (D.8)
In this last expression, we did not exploit yet explicitly the symmetries on C_{i jkl}; matrix. On the basis of these equations, we can now investigate the slope of the density profile for the three types of singularities we consider, depending on the values of the eigenvalues of matrix Ã_{ij}.
D.2 Asymptotic Behaviour of the Profiles at Collapse
We now analyse the asymptotic behaviour of the profiles around the origin. To this end, when we consider the following scaling, x → s x, implying r → sr, we examine how the Lagrangian coordinate changes. By doing so, we can understand the behaviour of the scaling of the Jacobian in terms of the Lagrangian coordinates, and thus reveal the behaviour of the density profile at the origin.
It is important to note that this proof is somewhat simplified and ignores details on the angular dependence when performing integrals over spherical shells to obtain the radial profiles. Therefore, the proof given in this appendix is not mathematically rigorous, but it leads to the same conclusions as the exact calculations in which proper form factors are estimated. The purpose of this section is indeed to provide a simplified rephrasing of singularity theory already presented in a more rigorous fashion in other works (see e.g. Arnold et al. 1982; Shandarin & Zeldovich 1989; Hidding et al. 2014; Feldbrugge et al. 2018).
First, we examine the case 0 = Ã_{xx} ≠ Ã_{yy} ≠ Ã_{zz}, which corresponds to the following sine waves initial conditions: Q1D2SIN, Q1D3SIN, ANI2SIN, and ANI3SIN (i.e. 0 ≤ ϵ_{2D} < 1 or 0 ≤ ϵ_{3D} < 1 for i = 1, 2), as summarised in first line of Table 2. This configuration corresponds to the shellcrossing along the x axis.
When computing the radial profiles close to the centre of the system, we consider spherical shells of radius r with r^{2} = x^{2} + y^{2} + z^{2} ≪ 1. At leading order in q we have, at the shellcrossing, (D.9) (D.10) (D.11)
Applying the scaling x → sx implies (D.12)
Then, applying the same scaling in Eq. (D.9) and taking the limit s ≪ 1, we simply obtain (D.13)
After exploiting the symmetries of matrix C_{i jkl} (Eq. D.6), the Jacobian up at leading order in q is given by (D.14)
Using the scalings in Eqs. (D.12) and (D.13), we have (D.15)
Then, we have J ∝ s^{2/3} when s ≪ 1, resulting in ρ ∝ r^{−2/3}.
Second, we consider the case 0 = Ã_{xx} = Ã_{yy} ≠ Ã_{zz}, which corresponds to the initial condition SYM2SIN (i.e. ϵ_{2D} = 1 or ϵ_{3D,1} = 1 and ϵ_{3D,2} < 1), as summarised in second line of Table 2. This configuration corresponds to simultaneous shellcrossings along x and y axes. In this case, the main contribution to Eulerian coordinates reads, at the shellcrossing, (D.16) (D.17) (D.18)
where we have exploited the symmetry between x and y, which implies C_{xxxx} = C_{yyyy} = C_{1}, C_{xxyy} = C_{yyxx} = C_{2} and C_{xxzz} = C_{yyzz} = C_{3}. Applying the same reasoning as above, the scaling x → sx implies (D.19)
which leads us, when applying the scaling to Eqs. (D.16) and (D.17) and assuming s ≪ 1, to neglect the last term in these equations, hence (D.20) (D.21)
which leads to the scalings (D.22)
With Ã_{xx} = Ã_{yy} = 0, the first term in Eq. (D.8) vanishes, and the second termbecomes the leading contribution: (D.23)
Exploiting Eq. (D.6) then allows us to write this term as (D.24)
where b_{i jk} are constant coefficients. Using the scalings given by Eqs. (D.19) and (D.22), we have (D.25)
The lowest power of s in Eq. (D.25) is s^{4/3} with i = j = 1 and k = 0. Thus, we have ρ ∝ r^{−4/3}.
Third, we consider the case Ã_{xx} = Ã_{yy} = Ã_{zz} = 0, which corresponds in the three sine waves case to the SYM3SIN (i.e. ϵ_{3D} = (1,1)), as summarised in third line of Table 2. This configuration corresponds to simultaneous shellcrossings along all the axes. The approach is analogous to the previous case. The Eulerian coordinates can be written as follows: (D.26) (D.27) (D.28)
with, again C_{1} = C_{xxxx} = C_{yyyy} = C_{zzzz}, C_{2} = C_{iij j} for i ≠ j. Equations (D.26)–(D.28) directly lead to the scaling of q: (D.29)
Now, the first and second terms in Eq. (D.8) vanish, and we have (D.30)
which reads, after exploitation of the symmetries (D.6), (D.31)
Using the scaling given by Eq. (D.29), we have (D.32)
The lowest power of s in Eq. (D.32) is s^{2}. Thus, we have ρ ∝ r^{−2}.
Finally, we focus on the velocity and pseudo phasespace density profiles. Using Eq. (D.1), the velocity field up to leading order in the Lagrangian coordinate is given by (D.33)
where the dot denotes derivative with respect to the cosmic time t. As can be seen from this relation, the velocity field is proportional to the Lagrangian coordinate, irrespective of initial conditions, because the shellcrossing condition Ã_{ij}(t_{sc}) = 0 does not imply Since, according to the calculations performed above, the leading contribution from the Lagrangian vector always come from the scaling q_{i} → s^{1/3}q_{i} when s ≪ 1, it is fairly easy to convince oneself that the velocity profiles are given, at small radii, by υ^{2} ∝ r^{2/3}, , and −υ_{r} ∝ r^{1/3}, from which we can infer as well the logarithmic slope of the pseudo phasespace density Q(r) through Eq. (41).
References
 Alard, C. 2013, MNRAS, 428, 340 [Google Scholar]
 Angulo, R.E., Hahn, O., Ludlow, A.D. & Bonoli, S. 2017, MNRAS, 471, 4687 [NASA ADS] [CrossRef] [Google Scholar]
 Arnold, V.I., Shandarin, S.F., & Zeldovich, I.B. 1982, Geophys. Astrophys. Fluid Dyn., 20, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Baldauf, T., Mercolli, L., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 123007 [NASA ADS] [CrossRef] [Google Scholar]
 Bardeen, J.M., Bond, J.R., Kaiser, N., & Szalay, A.S. 1986, ApJ, 304, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, J. Cosmology Astropart. Phys., 2012, 051 [CrossRef] [Google Scholar]
 Bernardeau, F. 1994, ApJ, 427, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, Phys. Rep., 3671 [Google Scholar]
 Bertschinger, E. 1985, ApJS, 58, 39 [Google Scholar]
 Blumenthal, G.R., Faber, S.M., Primack, J.R., & Rees, M.J. 1984, Nature, 311, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Bouchet, F.R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Bouchet, F.R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575 [NASA ADS] [Google Scholar]
 Buchert, T. 1992, MNRAS, 254, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Buchert, T. 1994, MNRAS, 267, 811 [NASA ADS] [Google Scholar]
 Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375 [Google Scholar]
 Buchert, T., Karakatsanis, G., Klaffl, R., & Schiller, P. 1997, A&A, 318, 1 [NASA ADS] [Google Scholar]
 Carrasco, J.J.M., Hertzberg, M.P., & Senatore, L. 2012, J. High Energy Phys., 2012, 82 [CrossRef] [Google Scholar]
 Carron, J., & Szapudi, I. 2013, MNRAS, 432, 3161 [Google Scholar]
 Chernin, A.D. 1993, A&A, 267, 315 [NASA ADS] [Google Scholar]
 Colombi, S. 2015, MNRAS, 446, 2902 [Google Scholar]
 Colombi, S. 2021, A&A, 647, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 dAmico, G., Gleyzes, J., Kokron, N., et al. 2020, J. Cosmology Astropart. Phys., 2020, 005 [CrossRef] [Google Scholar]
 Delos, M.S., Erickcek, A.L., Bailey, A.P., & Alvarez, M.A. 2018, Phys. Rev. D, 97, 041303 [NASA ADS] [CrossRef] [Google Scholar]
 Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389 [Google Scholar]
 Doroshkevich, A.G. 1973, Astrophys. Lett., 14, 11 [Google Scholar]
 Feldbrugge, J., van de Weygaert, R., Hidding, J., & Feldbrugge, J. 2018, J. Cosmology Astropart. Phys., 2018, 027 [CrossRef] [Google Scholar]
 Fillmore, J.A., & Goldreich, P. 1984, ApJ, 281, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hahn, O., Angulo, R.E., & Abel, T. 2015, MNRAS, 454, 3920 [CrossRef] [Google Scholar]
 Henriksen, R.N. & Widrow, L.M. 1995, MNRAS, 276, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Hertzberg, M.P. 2014, Phys. Rev. D, 89, 043521 [NASA ADS] [CrossRef] [Google Scholar]
 Hidding, J., Shandarin, S.E., & van de Weygaert, R. 2014, MNRAS, 437, 3442 [CrossRef] [Google Scholar]
 Hjorth, J., & Williams, L.L.R. 2010, ApJ, 722, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R.W., & Eastwood, J.W. 1988, Computer simulation using particles [CrossRef] [Google Scholar]
 Ishiyama, T. 2014, ApJ, 788, 27 [Google Scholar]
 Ivanov, M.M., Simonovic, M., & Zaldarriaga, M. 2020, J. Cosmology Astropart. Phys., 2020, 042 [CrossRef] [Google Scholar]
 Ludlow, A.D., Navarro, J.F., Springel, V., et al. 2010, MNRAS, 406, 137 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 1967, MNRAS, 136, 101 [Google Scholar]
 Matsubara, T. 2015, Phys. Rev. D, 92, 023534 [Google Scholar]
 Michaux, M., Hahn, O., Rampf, C., & Angulo, R.E. 2021, MNRAS, 500, 663 [Google Scholar]
 Moutarde, F., Alimi, J.M., Bouchet, F.R., Pellat, R., & Ramani, A. 1991, ApJ, 382, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Moutarde, F., Alimi, J.M., Bouchet, F.R., & Pellat, R. 1995, ApJ, 441, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Nakamura, T. 1985, Prog. Theor. Phys., 73, 657 [CrossRef] [Google Scholar]
 Navarro, J.F., Frenk, C.S., & White, S.D.M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J.F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Novikov, E.A. 1969, Sov. J. Exp. Theor. Phys., 30, 512 [NASA ADS] [Google Scholar]
 Ogiya, G., & Hahn, O. 2018, MNRAS, 473, 4339 [Google Scholar]
 Peebles, P.J.E. 1980, The largescale structure of the universe (Princeton University Press) [Google Scholar]
 Peebles, P.J.E. 1982, ApJ, 263, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P.J.E. 1984, ApJ, 277, 470 [NASA ADS] [CrossRef] [Google Scholar]
 Pichon, C., & Bernardeau, F. 1999, A&A, 343, 663 [NASA ADS] [Google Scholar]
 Pontzen, A., & Governato, F. 2013, MNRAS, 430, 121 [Google Scholar]
 Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C. 2012, J. Cosmology Astropart. Phys., 2012, 004 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C. 2019, MNRAS, 484, 5223 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C., & Frisch, U. 2017, MNRAS, 471, 671 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C., & Hahn, O. 2021, MNRAS, 501, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Rampf, C., Villone, B. & Frisch, U. 2015, MNRAS, 452, 1421 [CrossRef] [Google Scholar]
 Rampf, C., Frisch, U., & Hahn, O. 2021, MNRAS, 505, L90 [NASA ADS] [CrossRef] [Google Scholar]
 Saga, S., Taruya, A., & Colombi, S. 2018, Phys. Rev. Lett., 121, 241302 [Google Scholar]
 Shandarin, S.F., & Zeldovich, Y.B. 1989, Rev. Mod. Phys., 61, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Sikivie, P., Tkachev, I.L., & Wang, Y. 1997, Phys. Rev. D, 56, 1863 [NASA ADS] [CrossRef] [Google Scholar]
 Sousbie, T., & Colombi, S. 2016, J. Comput. Phys., 321, 644 [Google Scholar]
 Taruya, A., & Colombi, S. 2017, MNRAS, 470, 4858 [Google Scholar]
 Taylor, J.E., & Navarro, J.F. 2001, ApJ, 563, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlemann, C., Rampf, C., Gosenca, M., & Hahn, O. 2019, Phys. Rev. D, 99, 083524 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T.E., et al. 2020, Nature Methods, 17, 261 [CrossRef] [Google Scholar]
 Zel’dovich, Y.B. 1970, A&A, 5, 84 [Google Scholar]
 Zheligovsky, V., & Frisch, U. 2014, J. Fluid Mech., 749, 404 [Google Scholar]
 Zukin, P., & Bertschinger, E. 2010a, Phys. Rev. D, 82, 104044 [Google Scholar]
 Zukin, P., & Bertschinger, E. 2010b, Phys. Rev. D, 82, 104045 [Google Scholar]
The designations named here, Q1D3SIN, ANI3SIN, and SYM3SIN, are the same as Q1DS, ASYSb, and SYMS, used in Saga et al. (2018), respectively.
When examining the bottom panels of Figs. 9 (black curves) and 10, we see that solving the equation x(q) = x_{0} in the multistream region does not seem to have the expected number of solutions, that is, 3, 9, or 27 (only the bottomright panel has an expected number of solutions, i.e., 27). But this is because we are sitting exactly at the centre of the system and symmetries imply superposition of solutions in the subspace y = 0 in 2D and y = z = 0 in 3D. In particular, as discussed in Sect. 5.1, firstorder LPT has each coordinate axis totally decoupled from the others, which means that only one ‘S’ shape, which corresponds respectively, in the (x, υ_{x}) subspace, to 3 and 9 superposed ‘S’ shapes, is visible as the orange curve on bottom panels of Fig. 10. On the other hand, higherorder LPT induces nonlinear couplings between various axes of the dynamics, which implies that we have, respectively, only 2 and 3 visible ‘S’ shapes on the bottomleft and right panels of Figs. 9 (black curves) and 10, both for the simulations and LPT predictions at order n ≥ 2 among those the remaining 3−2=1 and 9−3 = 6 curves overlap with their visible ‘S’ counterpart in (x, υ_{x}) subspace.
We recall, however, that the crossed sine waves configurations we consider are degenerate with respect to 1LPT, as discussed at the end of Sect. 5.1, so the reasoning of Pichon & Bernardeau (1999) requires 2LPT to be applicable in this particular case.
All Tables
Summary of the logarithmic slopes at the shellcrossing obtained by singularity theory applied to LPT predictions at various orders (see Appendix D for details).
All Figures
Fig. 1 Convergence behaviour of the shellcrossing time. Top: shellcrossing time calculated at nthorder LPT (dots), together with the fitting function given by Eq. (34) (solid lines), as a function of the perturbation order for various initial conditions, as indicated in the panel. The dashed horizontal lines present the values estimated from the simulations (see Appendix B.1), corresponding to the sixth column of Table 1. Bottom: relative error between the shellcrossing time calculated with nthorder LPT and the one obtained with the fitting formula. 

In the text 
Fig. 2 Phasespace structure for two and threesinewave initial conditions at the collapse time: Q1D2SIN (top left), ANI2SIN (centre left), SYM2SIN (bottom left), Q1D3SIN (top right), ANI3SIN (centre right), and SYM3SIN (bottom right). The intersection of the phasespace sheet with the y = 0 plane for two sine waves and y = z = 0 hyperplane for three sine waves is displayed in (x, υ_{x}) subspace. Simulation results are compared with standard LPT predictions, which are supplemented with the blue line, denoted by’ΈΧΤ’, which corresponds to the formal extrapolation to the infinite order proposed by STC18 and sketched in Sect. 4.1 for the collapse time. For completeness, the quasi1D approach (Rampf & Frisch 2017), denoted by Q1D, is also presented (see Appendix C for details). 

In the text 
Fig. 3 Radial profiles and their logarithmic slopes for the twosinewave initial conditions at the shellcrossing time. As indicated in topleft panel, LPT predictions are given by solid lines of various colours and Vlasov simulations results are represented by red dots. From left to right, the initial conditions are given by Q1D2SIN (ϵ_{2D} = 1/6), ANI2SIN (ϵ_{2D} = 2/3), and SYM2SIN (ϵ_{2D} = 1), respectively. From top to bottom, we present the radial profiles of the normalised density , the velocity dispersion υ^{2}, the radial velocity dispersion , the infall velocity υ_{r}, and the pseudo phasespace density Q(r), respectively. We note that when plotting the logarithmic slopes in Vlasov simulations, we used the SavitzkyGolay filter implemented in savgol_filter of SciPy (Virtanen et al. 2020) to smooth the data for presentation purposes. In the logarithmic slope panels, the horizontal dashed lines correspond to the theoretical predictions at the small radii of Appendix D, as listed in Table 2. 

In the text 
Fig. 4 Same as Fig. 3 but for the threesinewave initial conditions, from left to right, Q1D3SIN [ϵ_{3D} = (1/6,1/8)], ANI3SIN [ϵ_{3D} = (3/4,1/2)], and SYM3SIN [ϵ_{3D} = (1,1)], respectively. We note that in SYM3SIN, the closest snapshot from collapse we had at our disposal from our Vlasov runs, â_{sc} = 0.03190, is significantly beyond the actual shellcrossing time estimated by the method described in Sect. B.2, a_{sc} = 0.03155, which explains some discrepancies between the theory and the simulation at small radii. 

In the text 
Fig. 5 Same as Fig. 3 but all analytical predictions evaluated at individual shellcrossing times computed at each perturbation order. 

In the text 
Fig. 6 Same as Fig. 4 but all analytical predictions evaluated at individual shellcrossing times computed at each perturbation order. 

In the text 
Fig. 7 Radial profiles at the shellcrossing using tenthorder LPT in the vicinity of (top panels) and (bottom panels). In the top panels, ϵ_{2D} is progressively varying in the range [0.9,1] with linear intervals of Δϵ = 0.1/64. In the bottom panels, ϵ_{3D} progressively changes from (0.97,0.97) to (1,1) with linear intervals of Δϵ = 0.03/64. From left to right, we present radial profiles of the normalised density , the velocity dispersion υ^{2}, the radial velocity dispersion the infall velocity –υ_{r}, and the pseudo phasespace density Q, respectively. In the lower inserts corresponding to logarithmic slopes, the horizontal dashed lines correspond to the values expected from singularity theory as computed in Appendix D and listed in Table 2. 

In the text 
Fig. 8 3D view of the expected caustic pattern shortly after the shellcrossing for threesinewave initial conditions. The left panels correspond to the most typical pancake, here embodied by ANI3SIN, but Q1D3SIN would look similar, while the right ones show the axialsymmetric case, i.e. SYM3SIN. The calculations are performed using 2LPT along with ballistic approximation (Eq. (53)), but higherorder LPT would provide the same topology from the qualitative point of view, except that the position of the caustic surfaces would change, especially in the axialsymmetric case. 

In the text 
Fig. 9 Tests of the ballistic approximation in Vlasov simulations: measured phasespace structure for twoand threesinewave initial conditions. This figure is analogous to Fig. 2, except that it shows a zoomedin view of the central part of the system shortly after collapse. We test the validity of the ballistic approximation by using two simulation snapshots slightly before (red curves, with as shown in Table 1 except for ϵ_{3D} = (1,1), for which a = 0.03103) and slightly after the collapse time (black, with a = a_{sc} + Δa). The ballistic approximation, as described in the main text, is applied to the red curves to obtain the magenta ones, to be compared directly with the exact solution, in black, given by the simulation. 

In the text 
Fig. 10 Phasespace structure for two and threesinewave initial conditions shortly after collapse. This figure is analogous to Fig. 2, except that we compare predictions of LPT at various orders in the ballistic approximation framework with measurements in Vlasov simulations shortly after collapse, i.e. for a = a_{sc} + Δa, as listed in Table 1 and discussed in the main text. We note that the additional cyan curve corresponds to the LPT prediction at the tenth order without using the ballistic approximation. 

In the text 
Fig. 11 Caustic pattern shortly after the shellcrossing in the 2D case: comparison of LPT predictions using ballistic approximation with Vlasov runs. The left, middle, and right panels correspond, respectively, to Q1D2SIN, ANI2SIN, and SYM2SIN configurations, while the top and bottom panels show the caustics in Lagrangian and Eulerian spaces, respectively. 

In the text 
Fig. 12 Structure of caustics shortly after the shellcrossing for threesinewave initial conditions: comparison of LPT predictions using ballistic approximation with Vlasov runs. The top six, middle six, and bottom two panels respectively correspond to Q1D3SIN, ANI3SIN, and SYM3SIN. In each group of six panels, the top and bottom lines correspond to Lagrangian and Eulerian spaces, and the intersection of the caustic surfaces with the plane are q_{x} = 0, q_{y} = 0, and q_{z} = 0, respectively, for the first, second, and third column of each group. In the bottom group of two panels, the left and right panels correspond respectively to Lagrangian and Eulerian space and only show the intersection of the caustic surface network with the plane q_{z} = 0 due to the symmetry of the system. 

In the text 
Fig. 13 2D density shortly after the shellcrossing: comparison of LPT predictions using ballistic approximation with Vlasov runs. From top to bottom: Q1D2SIN, ANI2SIN, and SYM 2SIN. From left to right: 2LPT, 10LPT, and measurements in Vlasov simulations. 

In the text 
Fig. 14 Slices of the projected density shortly after the shellcrossing: comparison of LPT using ballistic approximation with Vlasov runs. The left and right groups of nine panels correspond respectively to Q1D3SIN (top, middle, and bottom line of panels: x=−1.55×l0^{−4}, y = −1.17×10^{−3}, and z= −1.56 × 10^{−3} slice) and ANI3SIN (x= −5.16 × 10^{−4}, y = −2.15 × 10^{−4}, and z= −5.47 × 10^{−4} slice), while the bottom group of three panels corresponds to SYM3SIN (z= −1.17 × 10^{−5} slice). In each group of panels, the left, middle, and right columns give respectively the secondorder LPT prediction, the tenthorder LPT prediction, and the Vlasov code measurements. Due to the symmetry of the system for SYM3SIN, only one slice is shown for the bottom panels. 

In the text 
Fig. 15 2D vorticity fields: comparison of LPT predictions with Vlasov runs. 

In the text 
Fig. 16 Vorticity components shortly after the shellcrossing: comparison of LPT with the ballistic approximation with Vlasov runs for Q1D3SIN. The 2D slices are the same as those shown in the topleft group of nine panels of Fig. 14, except that the slice changes from left to right and the vorticity component changes from top to bottom. Again, on each line of three panels, 2LPT (left panel) and 1OLPT (middle panel) are compared with simulation measurements (rightpanel). 

In the text 
Fig. 17 Same as Fig. 16, but for ANI3SIN. The 2D slices considered are the same as in the topright group of nine panels of Fig. 14. 

In the text 
Fig. 18 Same as Fig. 16, but for SYM3SIN. Due to the symmetry of the initial conditions, only an xy slice is shown, which is the same as in bottom three panels of Fig. 14. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.