Cold dark matter protohalo structure around collapse: Lagrangian cosmological perturbation theory versus Vlasov simulations

We explore the structure around shell-crossing 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 10th order to high-resolution 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 locally the phase-space sheet at the quadratic level. To test LPT predictions beyond shell-crossing, we employ a ballistic approximation, which assumes that the velocity field is frozen just after shell-crossing. In the generic case, where the amplitudes of the sine waves are all different, high-order LPT predictions match very well the exact solution, even beyond collapse. As expected, convergence slows down when going from quasi-1D dynamics where one wave dominates over the two others, to the axial-symmetric configuration, where all the amplitudes of the waves are equal. It is also noticed that LPT convergence is slower when considering velocity related 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 axial-symmetric configurations, that are still correct from the qualitative point of view, but where multiple foldings of the phase-space sheet produce very high density contrasts, hence a strong backreaction of the gravitational force.


Introduction
In the concordance model of large-scale structure formation, the matter content of the Universe is dominated by collisionless cold dark matter (CDM) following Vlasov-Poisson equations (Peebles 1982(Peebles , 1984Blumenthal et al. 1984). The cold nature of initial conditions implies that the dark matter distribution is concentrated on a three-dimensional sheet evolving in six-dimensional phase-space. At shell crossing, that is in places where the phase-space sheet first self-intersects, the seeds of first dark matter haloes are created. In these regions, the fluid enters a multi-stream regime during which violent relaxation takes place (Lynden-Bell 1967) to form primordial CDM haloes. Numerical investigations suggest that first dark matter haloes formed during this process have a power-law 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 well known universal Navarro-Frenk-White profile (hereafter NFW, Navarro et al. 1996Navarro et al. , 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 haloes history, relying, for example, on entropy maximisation (Lynden-Bell 1967;Hjorth & Williams 2010;Carron & Szapudi 2013;Pontzen & Governato 2013), self-similarity (Fillmore & Goldreich 1984;Bertschinger 1985;Henriksen & Widrow 1995;Sikivie et al. 1997;Zukin & Bertschinger 2010b,a;Alard 2013), or more recently, a post-collapse perturbative treatment (Colombi 2015;Taruya & Colombi 2017;. 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 explaining the different phases of halo history remains yet to be proposed. The early growth of large scale structures up to first shellcrossing, on the other hand, is well understood thanks to perturbation theory. Indeed, restricting to the early phase of structure formation, i.e., to the single-stream regime before shell-crossing, we can employ perturbation theory as long as fluctuations in the density field remain small (see, e.g., Bernardeau et al. 2002, for a review and references therein). Lagrangian perturbation theory (LPT) (e.g., Shandarin & Zeldovich 1989;Bouchet et al. 1992;Buchert 1992;Buchert & Ehlers 1993;Bouchet et al. 1995;Bernardeau 1994) uses the displacement field as a small quantity in the expansion of the equations of motion. First-order LPT corresponds to the classic Zel'dovich approximation (Zel'dovich 1970), and in one-dimensional space, it is an exact solution un-Article number, page 1 of 31 arXiv:2111.08836v2 [astro-ph.CO] 15 Aug 2023 A&A proofs: manuscript no. ms_collapse til shell-crossing (Novikov 1969). Because the Lagrangian description follows elements of fluid along the motion, Zel'dovich approximation and higher-order LPT provide us with a rather accurate description of the large scale matter distribution, even in the nonlinear regime, shortly after shell crossing. The families of singularities that form at shell-crossing and after have been examined in detail in the Lagrangian dynamics framework (Arnold et al. 1982;Hidding et al. 2014;Feldbrugge et al. 2018), and the structure of cosmological systems at shell-crossing has been investigated for specific initial conditions (Novikov 1969;Rampf & Frisch 2017;Saga et al. 2018;Rampf 2019) and for random initial conditions .
As described above, investigations into the early stages of large-scale structure formation represent a key to understand the bottom-up scenario underlying the CDM model. Practically, perturbation theory under the single-stream 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 N-body simulations and observations (see, e.g., Bernardeau et al. 2002, for a review). Recently, in order to incorporate the effects of multi-streaming at small scales, it has been proposed to introduce effective fluid equations with a non-vanishing stress tensor in the dark matter fluids, the so-called effective field theory of large-scale 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 N-body simulations, it has attracted much attention, and has been recently applied to real observational datasets to derive cosmological parameter constraints (Ivanov et al. 2020;d'Amico et al. 2020). Accordingly, the understanding of multi-streaming effects from first principles, even at an early stage, would provide significant insights into the precision theoretical modelling of large-scale structure.
The aim of this paper is to extend the investigations of Saga et al. (2018), hereafter STC18. In this work, we studied the phase-space structure of protohaloes at shell-crossing with simplified initial conditions composed of three crossed sine waves, following the footsteps of Moutarde et al. (1991Moutarde et al. ( , 1995. Thanks to a detailed comparison between LPT predictions up to 10th order to numerical simulations performed with the state-of-theart Vlasov-Poisson solver ColDICE (Sousbie & Colombi 2016), we explicitly showed that the convergence of the LPT series slows down when going from quasi-one dimensional 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 infinite order, obtaining a remarkable agreement with the simulations, even at shellcrossing. At present, we further thoroughly examine the structure both at and shortly after shell-crossing of these systems by considering as well the two-dimensional case. 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 on the general case. Furthermore, initial conditions expressed only in terms of sine waves considerably simplify analytical calculations while still allowing one an insightful exploration of the pre-and post-collapse dynamics.
The prominent features shortly after shell-crossing are, for instance, the appearance of caustics, and in the multistream region delimited by these latter, non-trivial 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 high-resolution Vlasov simulations, in particular the vorticity, with unprecedented accuracy. To perform theoretical predictions shortly after shell-crossing, we use a simple dynamical approximation based on ballistic motion applied to the state of the system described by high-order LPT solutions at shell-crossing. While convergence with 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 shell crossing, not only for the three sine waves configurations we aim to examine (STC18), but also for more general, random initial conditions , although the effects of cut-offs on power-spectra remain to be investigated furthermore in the latter case. Therefore, as long as the back-reaction from multi-stream flows on post-collapse dynamics remains negligible, shortly after shell-crossing, the approximation of the dynamics we propose here should work. This is the first step toward a proper analytical description of postcollapse dynamics in 6D phase-space.
This paper is organized as follows. In Sec. 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 our sine waves initial conditions. In Sec. 3, the important features of the Vlasov solver ColDICE are briefly summarised. We also address some aspects of measurements uncertainties in simulation data. In Sec. 4, we examine the phase-space structure and radial profiles at shellcrossing with a comparison of analytic predictions to simulations, in the framework of singularity theory. In Sec. 5, the structure shortly after shell-crossing is explored using the ballistic approximation by examining phase-space diagrams, caustics, density and vorticity fields. Again, analytical predictions are compared to the Vlasov runs. Finally, Sec. 6 summaries the main results of this article. To complement the main text, Appendices A, B, C and D provide details on the explicit form of high-order LPT solutions for the sine waves initial conditions up to 5th order, the measurements in the Vlasov simulations, the predictions from quasi-1D LPT, and the expected structure of singularities at collapse time, respectively.
Throughout the paper, we will use the following units: a box size L = 1 and an inverse of the Hubble parameter at present time H −1 0 = 1 for the dimensions of length and time, respectively. In particular, the Lagrangian/Eulerian coordinates, velocity, and vorticity will be explicitly expressed as q/L, x/L, v/(L H 0 ), and ω/H 0 = (L ∇) × (v/(L H 0 )), respectively.

Lagrangian Perturbation Theory
The Lagrangian equation of motion of a fluid element is given by (e.g., Peebles 1980) 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 δ(x) = ρ(x)/ρ − 1 withρ being the background mass density, through the Poisson equation: ∇ 2 x ϕ(x) = 4πGρ a 2 δ(x).
In this framework, the velocity of each mass element is given by v = 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: 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 x(q, t) = q + Ψ(q, t).
The velocity field is expressed as v(q, t) = a dΨ/dt. In the single flow regime, i.e., before first shell-crossing time t sc , mass conservation impliesρ d 3 q = ρ(x) d 3 x, which leads to where the quantity J = det J i j is the Jacobian of the matrix J i j defined by and its inverse is given as While equation (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).

Recursion relation
In Lagrangian perturbation theory, the displacement field, Ψ, is the fundamental building block which is considered as a small quantity. It can be systematically formally expanded as follows, In what follows, we shall assume that the fastest growing modes dominate. They are well known to be given to a very good approximation by (see, e.g., Bernardeau et al. 2002, and references therein) where the purely time-dependent function D + (t) corresponds to the linear growth factor. The velocity field is then given by where function f (t) ≡ d ln D + /d ln a corresponds to the logarithmic derivative of the growth factor. Note that the analyses performed in subsequent sections will consider the Einstein-de 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.
Here and in what follows, we adopt the Einstein summation convention when the equation includes repeated lowercase Latin letters, and the subsrcipts i, j, .. take 1, 2 or 3 (1 or 2) in the three-dimensional (two-dimensional) case. In the above, we defined Ψ i, j ≡ ∂Ψ i /∂q j , and the tensors ε i j and ε i jk are, respectively, the two-dimensional and threedimensional Levi-Civita symbols. The symbolT stands for a differential operator: In obtaining the recursion relations, we used the following formulae of linear algebra: for the two-dimensional case, and for the three-dimensional case.
In the fastest-growing 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 ε i j Ψ (n) j,i = ε i j n 1 +n 2 =n n 2 (2n 2 + 1) n(2n + 1) Article number, page 3 of 31 A&A proofs: manuscript no. ms_collapse for the two-dimensional case, and for the three-dimensional case.

Two sine and three sine waves 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[: 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 three sine waves 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 , presents a small peak at the origin, and mass elements subsequently infall toward 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 three sine waves initial conditions quite generic, hence providing many insights into the dynamics during the early stages of protohalo formation. Naturally, this very symmetrical set-up 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: where the nth-order scalar coefficients α (n) m and the nth-order vector coefficients β (n) m are obtained recursively by calculating the right-hand-side of Eqs. (21) and (22), which depends on lower order terms, starting from the n = 1 coefficients determined by equation (23). By imposing the conditions ∇ q × Ψ (n) L = 0 and ∇ q · Ψ (n) T = 0, one can build up the perturbative solutions Ψ (n) L and Ψ (n) T using simple algebraic manipulations involving coefficients α (n) m and β (n) m : These solutions lead to the nth-order displacement field given by Ψ (n) = Ψ (n) L + Ψ (n) T . Using this solution for Ψ (n) as well as Eqs. (21) and (22), we subsequently construct the source of the (n + 1)th-order derivatives ∇ q · Ψ (n+1) L and ∇ q × Ψ (n+1) T . By repeating the above operation together with the recursive relations, one can derive, in principle, arbitrary high-order LPT solutions. In Appendix A, we present the explicit forms of the LPT solutions up to 5th order derived in this way 1 .
The above prescription is valid in three-dimensional space, and can be applied to the 2D case by cancelling all fluctuations along 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 realize that in two dimensions, vector coefficients β (n) m become scalars β (n) m , and that the solutions take the following form: 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 one-dimensional with |ϵ x | ≫ |ϵ y,z |, anisotropic with |ϵ x | > |ϵ y | > |ϵ z |, and what we design by axial-symmetric, with |ϵ x | = |ϵ y |(= |ϵ z |), denoted by Q1D, ANI and SYM, respectively 2 .

Vlasov-Poisson simulations
To perform the numerical experiments, we use the public parallel Vlasov solver ColDICE (Sousbie & Colombi 2016). ColDICE follows the phase-space sheet with an adaptive tessellation of simplices, composed, in 2 and 3 dimensions, of connected triangles and connected tetrahedra, respectively. The phase-space coordinates of the vertices of the tessellation, [X(t), V(t)], follow the standard Lagrangian equations of motion, similarly as in an N-body 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 (Sec. 3.1), we address measurement uncertainties issues (Sec. 3.2).  (Sousbie & Colombi 2016). The first column indicates the designation of the run. The second column corresponds to the relative amplitudes of the initial sine waves, namely, ϵ 2D = ϵ y /ϵ x and ϵ 3D = (ϵ y /ϵ x , ϵ z /ϵ x ) for two and three sine waves, respectively. The third column gives the value of ϵ x . The fourth and fith columns, respectively, indicate the spatial resolution n g of the grid used to solve the Poisson equation, and the spatial resolution n s of the mesh of vertices used to construct the initial tessellation (see Sec. 3 for details). The sixth column indicates the scale factor a ∞ sc at shell-crossing estimated by LPT extrapolated to infinite order (Saga et al. 2018), while the seventh one provides the value a sc measured in the Vlasov runs (see Appendix B.1) and which is actually used in Sec. 4.2. The eighth column indicates the valueâ sc of the expansion factor of the closest available simulation snapshot to collapse time for the comparisons performed in section 4. Finally, the last column indicates the value of the expansion factor used for the analyses performed beyond collapse time in Sec 5. Note: with the normalisation of the scale factor used here, the value of a sc has to be multiplied by |ϵ x | to obtain a more intuitive collapse time that would be equal to unity in the purely one dimensional case.

The Vlasov code ColDICE
The Lagrangian coordinates defined in Sec. 2.1 correspond to the following unperturbed initial set-up, where t 0 is a virtual time corresponding to a = 0, Q is the Lagrangian coordinate of each vertex. Note that we use capital letters to distinguish between vertices coordinates and actual coordinates of fluid elements of the phase-space sheet that they are supposed to trace. These notations are used in Appendix B. Vertices phase-space coordinates are then perturbed using Zel'dovich motion to set actual initial conditions defined in Sec. 2.2: with Ψ ini given by equation (23).
To update the position and the velocity of each vertex, a standard second order predictor-corrector 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 Courant-Friedrichs-Lewy 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 ).
Poisson equation is solved using the Fast-Fourier-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 phase-space sheet with each voxel/pixel of the mesh is computed exactly up to linear order with a special ray-tracing technique. Once the gravitational potential is obtained on the mesh, the gravitational force is computed from the gradient of the potential using a standard 4 points stencil, and then it is interpolated to the vertices using second-order 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 n D s vertices corresponding, respectively to 2n 2 s and 6n 3 s simplices in 2 and 3 dimensions. 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 phase-space sheet is locally described at quadratic order inside each simplex with the help of 3 and 6 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 axial-symmetric simulation, SYM-3SIN, so we do not deem it necessary to discuss more about refinement. However, the ability to describe the phase-space 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 3 s . In Appendix B, we explain how measurements of various quantities are performed, such as the set of curves corresponding to the intersection of the phase-space sheet with the hyperplane y(= z) = 0 used in Secs. 4.2 and 5.2, the collapse time t sc shown in Table 1, the radial profiles used in Sec. 4.3, as well as the caustic network, the projected density and the vorticity analysed in Sec. 5.3.

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 collapse time is critical, since we consider, in the analyses performed below, either collapse time or a moment very shortly after it. The uncertainty on collapse time depends on the nature of initial conditions and the parameters controlling the accuracy of the simulations, in partic-ular the number of simplices used to represent the phase-space 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 collapse time, as discussed in C21. In Sec. B.2, relying on measurements in ANI-3SIN 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, e.g. a sc = 0.02919 ± 10 −4 for ANI-3SIN, except for SYM-3SIN, where the uncertainty on the measured a sc could be of a few 10 −4 . Indeed, as just mentioned above, accuracy on the determination of collapse time also depends on the nature of initial conditions, the quasi-1D case and the triaxial-symmetric configuration being the easiest and the most difficult to deal with, respectively, as expected.
Note that accuracy on collapse time also depends on the ability to determine the exact position of the shell-crossing point. Using the tessellation technique, shell crossing coincides with a temporal change of sign of the orientation of the simplices. In our symmetrical set-up, the position of the shell-crossing point is simply the center of the system and the determination of collapse time can be simply performed using the intersection of the phase-space sheet with the y = z = 0 hyperplane (see Appendix B.2). However, even in a more complex case where the center of the forming halo would be moving, one expects to be able to compute collapse time accurately thanks to the tessellation representation as long as sampling of the phase-space sheet is sufficiently dense and, of course, if the quadratic representation inside each simplex is fully exploited. Sampling of the phase-space sheet is controlled by n s , as well as the refinement parameter I limiting deviations from symplectic motion. As discussed in Sec. 3.1, for our simulations, parameter I does not influence the results except (probably marginally) for SYM-3SIN. 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 to other sources of uncertainty on 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 making 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 multi-stream 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, i.e. 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 to 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 phase-space 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/surfaces, which we will keep in mind for the analyses carried out in Sec. 5.3.

Shell crossing structure
We are now in a position to study the structure of our protohaloes at collapse time, t sc , and concentrate our investigations on phase-space diagrams and radial profiles, with comparisons of LPT pushed up to 10th order to the Vlasov runs. This section is organized as follows. First, Sec. 4.1 presents the calculation of 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 infinite order of the LPT series proposed by STC18 and generalize it to the 2D case. Second, Sec. 4.2 examines the convergence of LPT at collapse with phase-space diagrams, extending the earlier investigation of STC18 to the 2D case. For comparison, we also test the formal extension of LPT to infinite order and the predictions of the quasi one-dimensional approach proposed by Rampf & Frisch (2017, hereafter RF17). Finally, in Sec. 4.3, LPT predictions and their convergence are studied in terms of radial profiles of the density, the velocity as well as the pseudo phase-space density, and put into perspective in relation to singularity theory.

Shell-crossing time
This subsection presents estimates of the expansion factor at first shell-crossing to which we refer to as 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 a (n) sc of collapse time obtained from nthorder LPT and its extrapolation to infinite order a ∞ sc , as described in STC18, to 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 to Vlasov runs in the 2D case.
Using nth-order LPT predictions, we explore the sequence of shell-crossing times, a (n) sc , 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 nth-order LPT solution. Generally, except in the pure onedimensional case where first-order LPT is exact before collapse, a (n) sc becomes smaller with increasing perturbation order n (see, e.g., RF17, , hereafter RH21, for recent works). As illustrated by Fig. 1, the shell-crossing time calcu-lated at nth-order 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 collapse time, but using the exponential of a powerlaw might be the only way to match the convergence speed of LPT at large n, when considering quantities computed at collapse time a (n) sc of each respective order. Equation (34) also implies a (n) sc −a (n−1) sc ∼ exp(−|d| n e ) when n ≫ 1, which might, at first sight, 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 power-law times an exponential of n. This would naively correspond to a (n) sc − a (n−1) sc ∼ n κ exp(−η n), i.e. e = 1, while our measurements suggest values of e different from unity, ranging from e.g. e ∼ 0.003 (SYM-3SIN) to e ∼ 0.8 (Q1D-3SIN) for fitting a (n) sc and from e.g. e ∼ 0.6-0.7 (SYM-3SIN) to e ∼ 1.3-1.5 (Q1D-3SIN) 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 shell-crossing times a (n) sc , or equivalently D (n) + . 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 equation (34) to 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 collapse time with order is rather slow, except in the quasi-1D case, which still requires at least third order for reaching a ∼percent level of accuracy for approximate convergence, while much higher order is required for other configurations, especially the 3D axialsymmetric case, ϵ 3D = (1, 1), for which convergence does not seem to be achieved even at 10th order. Fig. 1 thus demonstrates that low-order perturbation theory cannot be used to estimate accurately collapse times (see also STC18, in particular, Fig. 2 of this article, and RH21).
Despite these convergence issues, the extrapolated values of a ∞ sc obtained with equation (34) are very accurate, as illustrated by Table 1. Indeed, the relative difference between a ∞ sc 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 collapse time remains small and should not affect significantly the analytical predictions of the phasespace structure and the radial profiles at shell-crossing. Therefore, we will use a ∞ sc for the analytical predictions of "exact" collapse time in the remaining section 4. As for the simulations data, we will use the snapshot with the closest possible available valueâ sc of the expansion factor to a sc , as indicated in Table 1. Fig. 2 shows the (x, v x ) subspace of the phase-space structure at shell-crossing 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 infinite order as developed by STC18 and described in Sec. 4.1 for the shell-crossing time, as well as the analytical prediction obtained in the context of the quasi one-dimensional approach developed by RF17, that we extended to second order in the transverse direction (see Appendix C for details).

Phase-space structure
For quasi-1D (Q1D-2SIN and Q1D-3SIN) and anisotropic (ANI-2SIN and ANI-3SIN) initial conditions, the phase-space sheet first self-intersects along x-axis, and, until then, the dynamics is similar to the pure one-dimensional case, in which 1LPT (the Zel'dovich solution) is exact until shell-crossing. As illustrated by Fig. 2, the LPT prediction quickly converges with perturbation order n and describes the simulation results well even at shell-crossing, 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 x-coordinate of the velocity field as a function of position x. Note that, due to the symmetry of the initial conditions, these extrema are global and not specific to the phase-space slice represented in Fig. 2.
In the axial-symmetric cases, SYM-2SIN and SYM-3SIN, the phase-space sheet self-intersects simultaneously along all the axes of the dynamics. Interestingly, the phase-space structure for ϵ 2D = 1 (SYM-2SIN) is still qualitatively similar to one-dimensional collapse and LPT again quickly converges with perturbation order to the exact solution. This is clearly not the case for the three-dimensional axial-symmetric configuration, ϵ 3D = (1, 1), where LPT convergence is very slow. This set-up is indeed qualitatively different from other initial conditions, with the appearance of a spiky structure in phase-space (SCT18) similarly as in spherical collapse. Interestingly, in the spherical case and in the Einstein-de Sitter cosmology that we consider in this paper, convergence of the LPT is likely to be lost at collapse time The intersection of the phase-space sheet with the y = 0 plane for two sine waves and y = z = 0 hyper-plane for three sine waves is displayed in (x, v x ) subspace. Simulation results are compared to standard LPT predictions, that are supplemented with the blue line, denoted by "EXT", which corresponds to the formal extrapolation to infinite order proposed by STC18 and sketched in Sec. 4.1 for the collapse time. For completeness, the quasi one-dimensional approach (Rampf & Frisch 2017), denoted by Q1D, is also presented (see Appendix C for details).
(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 finite order, it is well reproduced by the formal extrapolation to infinite order. This is a hint that convergence of the LPT series might be lost at collapse for SYM-3SIN.
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 sine-wave configurations further in the non-linear regime. As described in C21, collapse of our protohaloes is followed by a violent relaxation phase leading to a power-law profile ρ(r) ∝ r −α , with α ∈ [1.5, 1.7] and then by relaxation to an NFW like universal profile. After violent relaxation, C21 did not find specific signatures in the density profile nor the pseudo phase-space density for the axial-symmetric case compared to the non-axialsymmetric ones, except that α tends to augment when going from Q1D to axial-symmetric, and that the axial-symmetric configuration is subject to significant (and expected) radial orbit instabilities.
To complete the analyses, note that the formal extrapolation of LPT to infinite order matches well (but not perfectly) the simulation results for all the configurations, as already found by SCT18 in the three-dimensional case. We also notice that the quasi one-dimensional approach of RF17 can describe very well the quasi-1D configurations. Interestingly, at the second order in the transverse fluctuations considered here, the predictions are rather similar to those of standard 4th-order LPT, irrespective of initial conditions. Thus, as expected, the quasi one-dimensional approach becomes less accurate when the ratio ϵ 2D or ϵ 3D approaches unity, particularly in 3D as already shown by STC18.

Radial profiles
We now focus on radial profiles, and define Eulerian polar and spherical coordinates for two and three sine waves initial conditions by (x, y) = (r cos θ, r sin θ) and (x, y, z) = (r sin θ cos ϕ, r sin θ sin ϕ, r cos θ), respectively. Then the angular averaged radial profiles of the density ρ/ρ, velocity dispersion v 2 , radial velocity dispersion v 2 r , and infall velocity −v r are given by wherex = x/r with r = |x| being the radial coordinate. In these equations, we used the angle average defined by for two sine waves initial conditions, and for three sine waves 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 To complete the analyses, we also study the pseudo phasespace density Q(r, t), defined by Note in addition that we present not only the radial profiles but also their logarithmic slopes defined by where the quantity X stands for the radial profile under consideration. Radial profiles at shell-crossing are presented in Figs. 3-6. In these figures, measurements in the Vlasov code are performed by replacing each simplex of the phase-space 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 correspond respectively to two and three sine waves initial conditions, calculations are performed at "exact" shell-crossing time, that is the extrapolated value a ∞ sc for LPT and the approximate valueâ sc for the simulations. Figs. 5 and 6 are analogous, but examine LPT predictions of nth-order synchronised to their own shell-crossing time a (n) sc , 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 examine further the quasi onedimensional approach proposed by RF17. We also do not perform the formal extension to 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 phase-space diagram analysis in Sec. 4.2. In particular, quasi-1D profiles (Q1D-2SIN and Q1D-3SIN) are well described by LPT predictions even at low order. At "exact" collapse, considered in Figs. 3 and 4, LPT predictions generally accurately to right, the initial conditions are given by Q1D-2SIN (ϵ 2D = 1/6), ANI-2SIN (ϵ 2D = 2/3), and SYM-2SIN (ϵ 2D = 1), respectively. From top to bottom, we present the radial profiles of the normalised density ρ/ρ, the velocity dispersion v 2 , the radial velocity dispersion v 2 r , the infall velocity −v r , and the pseudo pseudo phase-space density Q(r), respectively. Note that when plotting the logarithmic slopes in Vlasov simulations, we used the Savitzky-Golay filter implemented in savgol_filter of SciPy (Virtanen et al. 2020) to smooth the data, for presentation purposes. In the logarithmic slopes panels, the horizontal dashed lines correspond to the theoretical predictions at small radii of Appendix D as listed in Table 2.  1)], respectively. Note that in SYM-3SIN, the closest snapshot from collapse we had at disposal from our Vlasov runs,â sc = 0.03190, is significantly beyond actual shell-crossing time estimated by the method described in Sec. B.2, a sc = 0.03155, which explains some discrepancies between the theory and the simulation at small radii.
can be probed accurately by LPT, depending on the order considered and on the nature of initial conditions. Again, convergence worsens when approaching three-dimensional axial-symmetry, the worse case being, as expected, SYM-3SIN with ϵ 3D = (1, 1). It is worth noticing that velocity profiles require significantly higher LPT order than density 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 or- ders predict strikingly similar density profiles, in close to perfect match with the simulations, except for SYM-3SIN discussed further below, while the velocity profiles still diverge slightly from each other except for quasi one-dimensional initial conditions, Q1D-2SIN and Q1D-3SIN.
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 cold dark matter structures at 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 third order in the Lagrangian coordinate q. For our sine-wave configurations, whatever the order n of the pertur-S. Saga et al.: Cold dark matter structure around collapse bation order, 3 three kinds of singularities are expected at shellcrossing, as summarised in Table 2: the classic one-dimensional pancake with a power-law profile at small radii of the form (S1) ρ(r) ∝ r −2/3 , valid for all the configurations expect SYM-2SIN and SYM-3SIN; then (S2) ρ(r) ∝ r −4/3 and (S3) ρ(r) ∝ r −2 , for SYM-2SIN and SYM-3SIN, respectively. On the other hand, velocities are expected to follow the same power-law pattern whatever initial conditions or dimensionality, with logarithmic slopes equal to 2/3, 2/3, and 1/3 respectively, for v 2 , v 2 r , and −v r , which in turn implies Q(r) ∝ r −7/3 for SYM-2SIN, r −3 for SYM-3SIN, and r −5/3 for other configurations.
Figs. 3 and 4 consider LPT predictions for perturbation order n calculated at "exact" theoretical collapse time a ∞ sc and not at individual collapse time a (n) sc at this perturbation order, so 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. Note that convergence is slower for velocities than for the density, especially for axialsymmetric configurations, in particular SYM-3SIN. 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 SYM-3SIN for which the best available snapshot is still too far off 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 non-trivial 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 well known, but nonetheless somewhat trivial, 1D case, is remarkable even though naturally expected.
Several additional remarks are in order. First, if crossed sine waves are considered as generic approximations of local peaks of a smooth random Gaussian field, we also point out in this framework that the probability of having exactly ϵ 2D = ϵ SYM 2D ≡ 1 or ϵ 3D = ϵ SYM 3D ≡ (1, 1) is null, hence one expects the classic one-dimensional pancake to be exclusively dominant, from the pure mathematical view. However, this can be unrealistic at the coarse level, where high-density 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 ϵ SYM 2D (top panels) and ϵ SYM 3D (bottom panels) for 10th-order LPT. Clearly, at sufficiently large radii, values of ϵ 2D or ϵ 3D close to ϵ SYM 2D and ϵ SYM 3D give similar results for the profiles, even for the logarithmic slope which can approach that of the axial-symmetric case. The exact mathematical asymptotic behaviour is indeed reached only at very small radii. In practice, the axial-symmetric 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 axial-symmetric 3D case, SYM-3SIN, 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 axial-symmetric case (contained in terms beyond leading order r 2 in the Taylor expansion of trigonometric functions) provides significantly different results from pure spherical collapse. Note, though, that this state of fact is not fully proved by our simulations measurements, which as mentioned above, are slightly beyond actual collapse time.
Third, we can compare, in the three-dimensional case, the logarithmic slope of the pseudo phase-space density to 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 SYM-3SIN, which shows again that this particularly symmetric set-up, while being significantly more singular at collapse, with β = 3, relaxes to a state similar to more generic threedimensional configurations. The slope prior to collapse in the "generic" cases, β = −5/3, although of the same order of β sph , Sine wave amplitudes Designation Table 2. Summary of the logarithmic slopes at shell-crossing obtained by singularity theory applied to LPT predictions at various orders (see Appendix D for details). The first and second columns indicate the relative amplitudes of the initial sine waves and the designation of the runs, respectively. The third through seventh columns show the logarithmic slopes of the normalised density, velocity dispersion, radial velocity dispersion, infall velocity, and pseudo phase-space density, respectively. Note thus that the logarithmic slopes do not depend, as expected, on the perturbation order of LPT. ≡ (1, 1) (bottom panels). On top panels, ϵ 2D is progressively varying in the range [0.9, 1] with linear intervals of ∆ϵ = 0.1/64. On 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 v 2 , the radial velocity dispersion v 2 r , the infall velocity −v r , and the pseudo phase-space 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. 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.

Structure of the system shortly after shell crossing
Once shell-crossing time is passed, the system enters into a highly non-linear 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 multi-streaming flow on the force field using an LPT approach (see, e.g., Colombi 2015;Taruya & Colombi 2017, for the introduction of "post-collapse" perturbation theory). While the motion can still be described again for a short while in a perturbative way from shell-crossing, it is not, overall, a smooth function of time anymore due to the formation of the singularity at collapse time . Still, it is reasonable to assume that very shortly after shell crossing, the back-reaction corrections due to the multi-streaming 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 collapse time, but not necessarily far beyond collapse (see also RH21). Therefore, in what follows, we shall use the LPT solution of 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 infinite order proposed by STC18 and Article number, page 12 of 31 S. Saga et al.: Cold dark matter structure around collapse reintroduced in Sec. 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 multi-streaming solutions. So from now on, we shall consider only higher-order LPT calculations and not their extrapolation to infinite order (nor the quasi-1D approximation proposed by RF17).
This section is organized as follows. In Sec. 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 shell-crossing, caustics form and non-zero vorticity is generated in the multi-stream region delimited by the caustics. In this section, we also discuss the caustic network created up to second order by our systems seeded by sine waves. Sec. 5.2 turns to actual comparisons of analytical predictions to Vlasov simulations using, similarly as in Sec. 4.2, phase-space diagrams, but shortly after collapse. In particular, we explore the limits of the ballistic approximation. Finally, Sec. 5.3 focuses on the overall structure of the multi-stream region in configuration space, by successively testing LPT predictions against simulations for the caustic pattern, the projected density and the vorticity fields.

Multi-stream regime and ballistic approximation
Because collapse time a (n) sc depends on perturbation order, it only makes sense to test various perturbation orders shortly after shell-crossing and compare them to simulations only if collapse times are synchronised. In other words, in what follows, we consider the time a (n) pc = a (n) sc + ∆a for LPT of nth-order and a SIM pc = a sc + ∆a 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 6th columns of the table. Remind 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 nth-order LPT solutions at shell-crossing time a (n) sc of each perturbation order, we model the Eulerian coordinate after collapse as follows: while the nth-order velocity field is frozen to its value at shellcrossing time of each perturbation order, v(q, a) = v (n) (q, a (n) sc ). In the multi-stream region, given the Eulerian coordinate x, the solution of the equation x = x(q) has an odd number of solutions for the Lagrangian coordinate, that we denote 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, e.g. ϵ 2D = 1 in 2D or ϵ x = ϵ y in 3D, due to simultaneous shell-crossing along several axes, the number of streams can reach 9 and 27, respectively in two and three dimensions, as discussed further in Sec. 5.2.
In what follows, we omit the time dependence for brevity. After defining the Lagrangian density and velocity fields: the Eulerian density and velocity fields are expressed by the superpositions of each flow: 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 single-stream regions, unless already present in the initial conditions (and in the latter case, it corresponds to a decaying mode), its generation is one of the prominent features of shell-crossing (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 2 and 3 dimensions, by The vorticity fields in two-and three-dimensional space are scalar and vector quantities, respectively. Because the acceleration derives from a gravitational potential, local vorticity on each individual fold of the phase-space sheet cancels. Only the nonlinear superposition of folds in the first term of equations (49) and (50) induces non-zero vorticity, while the second term should not contribute. Using LPT however generates spurious vorticity in Eulerian space due to the truncation at 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 single-stream 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 first-and 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 multi-stream 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 Article number, page 13 of 31 A&A proofs: manuscript no. ms_collapse calculated shortly after shell-crossing: with A = x, y, or z and D (1LPT) +,sc = −1/ϵ x stands for the growth factor at shell-crossing time evaluated by using 1LPT. Note that in equation (51), the growth rate f and the expansion factor a sc are evaluated at the same time as D (1LPT)

+,sc
(remind that f = 1 in the Einstein-de Sitter cosmology considered in this work). The value of D (1LPT) +,sc is obtained in practice by solving the equation ∂x (1LPT) /∂q x = 0 at the origin, keeping in mind that |ϵ x | ≥ |ϵ y,z |. Remind that equation (51) is exact in the pure 1D case, that is when ϵ y = ϵ z = 0. The condition for the caustics, J (1LPT) (q, a) = 0, is reduced to the relation: Eq. (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 ∆a. Therefore the 1LPT caustics consist of onedimensional lines in 2D and two-dimensional planes in 3D, even at 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 first shell-crossing occurs should be composed, in the non-degenerate case (that is non-vanishing ϵ i ), of a set of points. In the vicinity of such points and shortly after shell-crossing, even in 1LPT, the equation of the caustics should correspond, at leading order in q and in the nonaxial-symmetric 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 non-linear couplings between various axes of the dynamics: where D (2LPT) +,sc stands for the growth factor at the shell-crossing time evaluated by using 2LPT, which can be obtained as a function of ϵ 2D or ϵ 3D by solving the following second order polynomial: In Eq. (53) , the growth rate f and the scale factor a sc are evaluated at the same time as D (2LPT) +,sc . The additional terms compared to Eq. (51) imply, shortly after collapse, non-vanishing contributions from all axes at the quadratic level (and above) in the Lagrangian coordinate q in the equation J (2LPT) (q, a) = 0   Fig. 8. Three-dimensional view of the expected caustic pattern shortly after shell-crossing for three sine waves initial conditions. The lefts panels correspond to the most typical pancake, here embodied by ANI-3SIN, but Q1D-3SIN would look alike, while the right ones show the axial-symmetric case, that is SYM-3SIN. The calculations are performed using 2LPT along with ballistic approximation (Eq. 53), but higher-order 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 axial-symmetric case.
(as long as none of the ϵ i cancels). This leads to, in the nonaxial-symmetric case (i.e., Q1D-2SIN, Q1D-3SIN, ANI-2SIN and ANI-3SIN), elliptic and ellipsoid shapes for the caustic curve/surface in Lagrangian space, respectively in 2D and 3D, as expected and as illustrated in the three-dimensional case by the top left panel of Fig 8, while the bottom left one shows the expected typical pancake shape in Eulerian space. Accordingly, in the analyses performed below, except for the phase-space diagrams that we examine first, LPT is examined from 2nd order and above. 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 quadratic order in q to obtain a full description of the correct and more intricate topology of the caustic surfaces (e.g. 4th order in q instead of second for SYM-2SIN), as illustrated in 3D by right panels of Fig. 8, due to simultaneous collapses along several axes. It would be beyond the scope of this article to go further in analytic details of catastrophe theory for these very particular configurations, so we shall be content with a descriptive analysis of the LPT results in this latter case.

Phase-space structure
Figs. 9 and 10 display phase-space diagrams for our various sine-wave systems shortly after collapse, namely the intersection of the phase-space 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 to 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 config- Tests of the ballistic approximation in Vlasov simulations: measured phase-space structure for two and three sine waves initial conditions. This figure is analogous to Fig. 2, except that it shows a zoom on 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 a =â sc as shown in Table 1 except for ϵ 3D = (1, 1), for which a = 0.03103) and slightly after 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 to the "exact" solution, in black, given by the simulation.
urations in 2D, with a slight worsening of the accuracy when going from the quasi-1D to the axial-symmetric set-up, as expected, significant or even dramatic deviations from the exact solution can be seen in 3D for the anisotropic configuration (ANI-3SIN) and the axial-symmetric configuration (SYM-3SIN). In the last case (bottom right panel), the intersection of the phasespace sheet with the hyperplane y = z = 0 is composed, in the multi-stream region, of three curves with a "S" shape instead of a single one, as a result of simultaneous shell-crossing along all coordinates axes. The prediction from the ballistic approximation has then to be compared to the "S" with the largest amplitude along x-axis. Similar arguments apply to the 2D case (bottom left panel), but in this case, the intersection with the y = 0 hyperplane is composed of two curves instead of three. 4 4 Note that when examining the bottom panels of Figs. 9 (black curves) and 10, solving the equation x(q) = x 0 in the multi-stream region does not seem to have the expected number of solutions, that is either 3, 9 or 27 (only in the bottom right panel for the latter case). 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 Sec. 5.1, first-order LPT has each coordinate axis totally decoupled from the others, which means that only one "S" shape, which corresponds respectively, in the (x, v x ) subspace, The results shown in Fig. 9 must however be interpreted with caution. Indeed, strictly speaking, we aim to apply the ballistic approximation from 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 actual collapse (namely a bc = 0.03103 for SYM-3SIN and a bc =â sc as given in Table 1 for other cases), because we do not have exactly access to this latter. The level of proximity to collapse is, at least from the visual point of view, the lowest for the axial-symmetric cases, particularly SYM-3SIN. It is also for this latter case that 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. Note also that ∆a eff /a bc ≃ 0.03 is also of the same order for ANI-3SIN. The greater the value of ∆a eff /a bc , the greater the deviation due to 3 and 9 superposed "S" shapes, is visible as the orange curve on bottom panels of Fig. 10. On the other hand, higher-order 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 bottom left 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, v x ) subspace.  Fig. 10. Phase-space structure for two and three sine waves 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 to measurements in Vlasov simulations shortly after collapse, that is for a = a sc + ∆a as listed in Table 1 and discussed in the main text. Note that the additional cyan curve corresponds to the LPT prediction at 10th-order without using the ballistic approximation. to nonlinear corrections of the motions to be expected. Furthermore, we also expect the ballistic approximation to worsen from quasi-1D cases to fully axial-symmetric. When examining Fig. 9 more in detail, we also note that for ANI-3SIN, 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 axial-symmetric case, not only the velocities in the tails are strongly underestimated, but also the magnitude of the position of the caustics along 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 to the simulations, as shown in Fig. 10. Remind 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 Sec. 4.2, where we compared LPT to simulations at collapse, still stand at the coarse level, obviously, since the time considered in Fig. 10 is nearly equal to 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 quasi-1D dynamics. It can fail in the tails of the "S" shape of the phase-space sheet, as a combination of limits of LPT to de-scribe 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 multi-stream region, where LPT predictions can still perform very well. This is illustrated on Fig. 10 by the cyan curves which show 10thorder LPT results without assuming ballistic approximation, to be compared to the dark purple curves which correspond to 10thorder 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 10th-order 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 axial-symmetric configurations, it is not clear at all that the cyan curves bring any improvement over the dark purple ones outside the multi-stream region. Here, LPT performs significantly more poorly than for other configurations, with or without assuming ballistic approximation. In this case, the topology of the phase-space structure in the multi-stream region (two or three curves according to the number of dimensions) is correct and the very central part of the shell corresponding to x-axis 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.

Configuration space: caustics, density, vorticity
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 multi-stream 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 phase-space 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 multi-stream 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 to 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 2nd and 10th order to 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. Note cautiously 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 non-axialsymmetric case, both in 2D and in 3D. Note that in the 3D case, to have a more clear view of the caustic pattern, we consider the intersection of the caustic surfaces with the planes x = 0, y = 0 and z = 0 (remind 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 Sec. 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. On the figures, only the vertices of these segments are represented, and the initial regular pattern of the tessellation used to represent the phase-space sheet is clearly visible in Lagrangian space.
Figs. 11 and 12 confirm the results of the previous section, namely that high order 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, ANI-2SIN and ANI-3SIN, the match between LPT and simulations is not perfect, even at 10th order, but remains still reasonably good, especially in 2D. For the axial-symmetric 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 axial-symmetric 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 quite underestimated by LPT and it seems that increasing 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, partic- ularly the spike observed on the phase-space 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 axial-symmetric cases, we can see that the agreement between simulations and LPT of high order is quite good, inside and outside the caustics. Note that the very high density contrasts observed in the axial-symmetric cases are due to the multiple foldings of the phase-space 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 non-axial-symmetric 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 Sec. 3.2), which provides details on the measurements of various fields, in particular how we exploit a quadratic description of the phase-space 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 multi-stream 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 which 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 shell crossing takes place along 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 phase-space sheet. In particular, it is easy to convince oneself, that if collapse happens along xaxis, the strongest variations of all the fields are expected along this direction. That means, since each coordinate of the vorticity vector depends on variations of 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 on 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. Note as well, in these figures, that 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 axial-symmetric 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 phase-space 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 multi-stream region seems to be qualitatively reproduced by the theory, but its size is totally wrong.

Summary
In this paper, following the footsteps of Saga et al. (2018), we have investigated the structure of primordial cold dark matter (CDM) haloes seeded by two or three crossed sine waves of various amplitudes at and shortly after shell-crossing, by comparing thoroughly up to 10th-order Lagrangian perturbation theory (LPT) to high-resolution Vlasov-Poisson simulations performed with the public Vlasov solver ColDICE. We devoted our attention first to the phase-space structure and radial profiles of the density and velocities at shell-crossing, and, second, to the phase-space structure, caustics, density and vorticity fields shortly after shell-crossing. 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 locally the phase-space sheet structure 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:  ballistic approximation, where the acceleration is neglected after the shell-crossing time.
Our main findings can be summarised as follows: -Phase-space diagrams at collapse: except for SYM-3SIN, one expects the system to build up a classic pancake singularity at shell-crossing, with a phase-space structure along x-axis analogous to what is obtained in one dimension. This pancake is well reproduced by LPT which converges increasingly well to the exact solution when approaching quasi onedimensional dynamics, as expected. The local extrema of the velocity field around the singularity are the locations where LPT differs most from the exact solution, by underestimat- ing the magnitude of the velocity. Convergence with perturbation order becomes very slow when approaching threedimensional axial-symmetry, where spikes appear on the velocity field on each side of the singularity. -Radial profiles at collapse: with a sufficiently high order of the 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 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, that is Q1D and ANI, and ρ(r) ∝ r −4/3 for SYM-2SIN, ρ(r) ∝ r −2 for SYM-3SIN, while velocities profiles always present the same power-law behaviour, e.g. v 2 (r) ∝ r 2/3 . Confirming predictions of singularity theory at collapse and explicit convergence to asymptotic profiles at small radii is expected but non-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 power-law 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 phase-space structure of the system slightly beyond collapse, except in the axial-symmetric case SYM-3SIN, where the very highly contrasted nature of the system due to multiple superpositions of phase-space 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, non-axial-symmetric configurations, that LPT of sufficient order combined with the ballistic approximation provides a very good description of the structure of the system beyond collapse time, with excellent agreement between the predicted density and the vorticity fields inside the multi-stream region and those measured in the Vlasov simulations. Turning to axial-symmetric cases, this agreement is only obtained at the qualitative level, even for 10th-order LPT, particularly for SYM-3SIN, where the size of the outer caustics is strongly underestimated.
Obviously, the ballistic approximation is only the first step for a more complete calculation taking into account the feedback due to the force field, as first proposed in the one dimensional case by Colombi (2015); Taruya & Colombi (2017) (see also the recent work of , and formulated in terms of "post-collapse perturbation theory". The next step is indeed to implement a Lagrangian perturbation theory approach where the small parameter is the interval of time ∆a = a − a sc from collapse time and where a Taylor expansion of the phase-space sheet is performed around the singularity in terms of Lagrangian coordinates. Shortly after collapse, in the generic, non-axialsymmetric case, the system presents a "S" shape in x − v x subspace, similar to the 1D case where position and velocities can be approximated as third-order polynomials of the Lagrangian coordinate q. The calculation of the force field requires then to solve a 3 value problem in the multi-stream region, similarly as in the 1D case. While technically challenging, generalisation of post-collapse perturbation theory to 2D and 3D seems possible. It might bring strong insights on nonlinear corrections due to multi-stream dynamics on large scale structure statistics, e.g. predictions of the power spectrum of the large-scale matter distribution from higher order perturbation theory.
words, this intersection corresponds in the non-trivial 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 phase-space sheet with the hyperplane x = y = 0. Furthermore, as discussed more in detail in C21, because the phase-space 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 phase-space diagram, we employ an approach valid at linear order, consisting in simply computing, for each simplex, its geometric intersection with the hyperplane P, which is either empty, a point or a segment. As a result, phase-space diagrams extracted from the Vlasov runs consist of sets of connected segments of which the ends 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 collapse time is essential when studying the properties of the system at shell crossing 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 phase-space diagram segment portion S of which one of the ends is as close as possible to the origin, while the other one, with abscissa x i , has v x > 0 and the magnitude of other coordinate(s) of the velocity, in theory null, as small as possible. The latter condition excludes, in the axial-symmetric case, the component(s) of the flow corresponding to simultaneous collapse along the y or/and z direction, which has/have non-zero value of v y or/and v z . Once this segment is identified at both times corresponding to a 1 and a 2 , a simple linear interpolation is used to estimate expansion factor at collapse: This formula does not require to be just before and after collapse time to estimate a sc , but, to provide sufficiently accurate results, it is obviously needed to consider times sufficiently close to actual collapse time. A first guess of collapse time is estimated by using perturbation theory predictions extrapolated to infinite order (STC18), as listed in 5th column of Table 1. Examination of Table 1 shows that the value a ∞ sc predicted this way agrees extremely well with the measurements in the Vlasov runs provided by equation (B.1), and shown in 6th column of the table.
Note however that, beyond the approximate nature of the linear approximation underlying equation (B.1), collapse time estimate can, among others, as discussed in the main text in Sec. 3.2, be significantly influenced by force resolution, that is the value of n g . As discussed in detail in C21, decreasing force resolution delays collapse time, hence, to have an accurate estimate of this latter, it is required to have a sufficiently large value of n g . 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 a ∞ sc = 0.02911), 0.0292 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 collapse expansion factor accurate at approximately the 10 −4 level and we expect this to apply as well to the quasi one-dimensional case ϵ 3D = (1/6, 1/8). However, for the axial-symmetric 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 center 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
Caustics are regions where the determinant J of the Jacobian matrix changes sign. At linear order in the local description of the phase-space 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 phase-space sheet is composed of a tessellation of triangles, hence the caustics estimated this way are given by sets of segments of which the ends are shown in Fig. 11. In 3D, the caustics are given by sets of triangles, of which we compute the intersection with the y = z = 0 plane to get again a set of segments of which the 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 Sec. 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, that 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 quantities: the jacobian of the transformation between initial and final positions, J(q) as defined in equation (7); the Lagrangian projected density ρ L (q) defined in equation (45); the total Eulerian projected density, which stems from the superposition of one or more folds of the phase-space sheet, as described by equation (47); the Eulerian velocity field, given by equation (48); the vorticity, ω 2D and ω 3D , defined in equations (49) and (50). Again, because the acceleration derives from a potential, local vorticity on each phase-space sheet fold cancels.
Only the nonlinear superposition of phase-space 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 phase-space sheet by ColDICE. In § B.5.1, we introduce barycentric coordinates, which are useful to define the position of any point inside each simplex. In § 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 § B.5.3, we generalize the procedure to the case when a quadratic description of the phase-space sheet is available. Finally, in § B.5.4, we describe the way we sample various quantities described above on an high resolution cartesian grid, by using proper sets of sampling particles associated to each simplex.

B.5.1. Barycentric coordinates
Since the phase-space sheet is sampled with simplices, it is useful to define a well known 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 dimension of configuration space, and a function g(X), one can define the following linear interpolation 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 phase-space sheet is flat and the linear interpolation In what follows, we shall therefore use equation (B.4) to define barycentric coordinates. With this definition of ξ k , this means that at given time t, equation (B.3) remains valid, but only at the linear level, since dynamical evolution of the phase-space 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 equation (B.2). The derivative of the function inside each simplex can thus be written where M is the matrix of the partial derivatives of the vector P ≡ ( k ξ k , X 1 , X 2 , X 3 ) with respect to ξ k , with X given by equation (B.3). In other words, with X k,α the α th coordinate of k th vertex with position X k . At the linear level, the derivative defined by equation (B.5) is therefore simply constant inside each simplex. This calculation can be performed in Lagrangian space to obtain directly the Jacobian of the transformation between initial and present position, equation (7), from the jacobian matrix where W is defined similarly as M but in Lagrangian space, 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 equations (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 phase-space sheet for performing local refinement of the tessellation, by using tracers defined as mid-points of edges of the simplices in Lagrangian space. For instance, the tracer associated to 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. Note that the second order representation preserves the continuity of the phase-space 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 , where functions K(k) and L(k) are appropriately chosen to cover all the combinations 1 ≤ K(k) < L(k) ≤ N s , and where N s = D + 1 is the number of simplices and N t = D(D + 1)/2 is the number of tracers, given D the dimension of space. Then equation (B.2) becomes where g k is now also defined over tracers. This is true in particular for quadratically interpolated positions while corresponding Lagrangian coordinates Q are still given by the linear representation (B.4) since the phase-space sheet is initially flat, which means that matrix W does not change.
From equation (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 equations (B.6) and (B.7): 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 by M quad in equation (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 equation (B.13) (but staying aware of the fact that this calculation is accurate only up to linear order), which, in addition to phase-space 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 equation (B.13) by This quantity depends in a non-linear 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/curve approximating the actual caustics).
In particular, we estimate the partial derivatives of the projected density as follows using equation (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(ξ) = x 0 . While this is a straightforward linear problem in the linear representation of the phase-space sheet (equation B.3), it becomes non-trivial in the quadratic case (equation 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/voxel of the cartesian mesh: Integrals of the form with v pixel/voxel the area/volume of the pixel/voxel under consideration. While Sousbie & Colombi (2016) use actual ray-tracing exact to linear order to compute such an integral, we replace 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 use cloud-in-cell (CIC) interpolation procedure on the target cartesian mesh (Hockney & Eastwood 1988), that is each particle is replaced by a voxel/pixel of the same (small) size as the voxels/pixels of the cartesian mesh, and associate a weight W inter to the particle proportional to the volume of intersection between the voxel/pixel associated with the particle and the voxel/pixel 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 sub-simplex 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/pixel of the cartesian mesh are kept under control: where D is the dimension of space, ∆ is the step of the cartesian mesh and S = max k k ′ ,α |X k,α − X k ′ ,α |, (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 equation (B.21), where V L is the Lagrangian volume of the simplex and J p the jacobian (quadratically) interpolated at position of the particle. So, at the end, the contribution of each particle in right hand integral of equation (B.21) is equal to δV p W inter g(X p , V p ), where we remind that W inter is the CIC weight defined previously and where function g is estimated at particle phase-space position (X p , V p ) using the methods described in previous paragraphs. Note that coarse graining due to the CIC interpolation is expected to introduce defects/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 quasi one-dimensional flow
In Sec. 4.2, we test theoretical predictions relying on the perturbative treatment of a quasi one-dimensional flow as proposed by RF17. In this appendix, we briefly present this approach, focusing mainly on three sine waves initial conditions configurations. Generalisation to the two sine waves case is straightforward.
When the system is exactly one dimensional, the first-order LPT solution is exact before shell-crossing. This fact leads to an approximate treatment in three-dimensional space. We exploit the exact solution of the one-dimensional problem along x-axis given by Zel'dovich approximation as the unperturbed zerothorder state: 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 first-order perturbations to this set-up: Ψ (1) (q, t) = ∞ n=1 D n + (t) Ψ (1,n) (q).
The results of the Q1D solution at first order in the transverse fluctuations are presented in RF17. Here, we show (for the first time) the Q1D solution at second second-order in the transverse fluctuations and up to 5th-order in the longitudinal direction in terms of the growth factor, i.e., Ψ (2,5) as follows. We use the same notations as in Appendix A.

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 shell-crossing and small radii, as discussed in Sec. 4.3.
The main results of this investigation are summarised in Table 2.
Since the extension to the two-dimensional 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. 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 higher order terms, x i (q, t) =Ã ia (t)q a + C iabc (t) q a q b q c , (D.1) withÃ i j (t) ≡ δ i j + A i j (t), where A i j (t) is some function of time. In these equations, we have exploited the symmetric nature of the three sine waves initial conditions, which imply that Eulerian coordinates around the shell-crossing point can be expanded in terms of odd third-order polynomials of the Lagrangian coordinates, that is polynomial forms P verifying P(−q) = −P(q) (see, e.g., Colombi 2015;Taruya & Colombi 2017, for the onedimensional case). The coefficients A i j and C i jkl are expressed in terms of partial derivatives of the displacement field with respect to the La-