Open Access
Issue
A&A
Volume 664, August 2022
Article Number A3
Number of page(s) 31
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202142756
Published online 03 August 2022

© S. Saga et al. 2022

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe-to-Open model. Subscribe to A&A to support open access publication.

1 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, 1984; Blumenthal et al. 1984). The cold nature of initial conditions implies that the dark matter distribution is concentrated on a 3D sheet evolving in 6D phase-space. At the shell-crossing, that is, in places where the phase-space sheet first self-intersects, the seeds of the 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 dark matter haloes formed during this process initially 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 (NFW) profile (Navarro et al. 1996, 1997).

From an analytical point of view, many approaches have been proposed to describe the results of simulations and the different steps of dark matter halo history, relying, for example, on entropy maximisation (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; Rampf et al. 2021). However, because of the highly non-linear complex processes taking place during the post-collapse phase, the formation and evolution of dark matter haloes remain a subject of debate, and a consistent framework that explains the different phases of halo history remains to be proposed.

The early growth of large-scale structures up to the first shell-crossing, on the other hand, is well understood thanks to perturbation theory. Indeed, restricting ourselves to the early phase of structure formation (i.e. to the single-stream regime before the shell-crossing), we can employ perturbation theory as long as fluctuations in the density field remain small (see e.g. Bernardeau et al. 2002, and references therein for a review). Lagrangian perturbation theory (LPT; e.g. Shandarin & Zeldovich 1989; Bouchet et al. 1992, 1995; Buchert 1992; Buchert & Ehlers 1993; Bernardeau 1994) uses the displacement field as a small quantity in the expansion of the equations of motion. First-order LPT corresponds to the classic Zel’dovich approximation (Zel’dovich 1970), and in 1D space it is an exact solution until the shell-crossing (Novikov 1969). Because the Lagrangian description follows elements of fluid along the motion, the Zel’dovich approximation and higher-order LPT provide us with a rather accurate description of the large-scale matter distribution, even in the non-linear regime, shortly after the shell-crossing. The families of singularities that form at the shell-crossing and after have been examined in detail within the Lagrangian dynamics framework (Arnold et al. 1982; Hidding et al. 2014; Feldbrugge et al. 2018), and the structure of cosmological systems at the 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 (Rampf & Hahn 2021).

As described above, investigations into the early stages of large-scale structure formation represent a key to understanding 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 large-scale structure distribution and has been successfully confronted with N-body simulations and observations (see e.g. Bernardeau et al. 2002, for a review). In order to incorporate the effects of multi-streaming at small scales, introducing effective fluid equations with a non-vanishing stress tensor in the dark matter fluids has recently been proposed, 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 recently been applied to real observational datasets to derive cosmological parameter constraints (Ivanov et al. 2020; dAmico 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 that work, we studied the phase-space structure of protohaloes at the shell-crossing with simplified initial conditions composed of three crossed sine waves, following in the footsteps of Moutarde et al. (1991, 1995). Thanks to a detailed comparison between LPT predictions up to the tenth order and numerical simulations performed with the state-of-the-art Vlasov-Poisson solver ColDICE (Sousbie & Colombi 2016), we explicitly showed that the convergence of the LPT series slows down when going from quasi-1D to triaxial-symmetric initial conditions, where a spiky structure in phase-space appears. Using a global fitting form, we were able to formally extrapolate the LPT solution to the infinite order, obtaining remarkable agreement with the simulations, even at the shell-crossing. In the present work, we thoroughly examine the structure further, both at and shortly after the shell-crossing of these systems, by considering the 2D case as well. Although sine wave initial conditions are restrictive, they are to a large extent representative of high peaks of a smooth random Gaussian field (see e.g. Bardeen et al. 1986) and should provide fruitful guidance in the general case. Furthermore, initial conditions expressed only in terms of sine waves considerably simplify analytical calculations while still allowing for an insightful exploration of the pre- and post-collapse dynamics.

The prominent features shortly after the shell-crossing are, for instance, the appearance of caustics and, in the multi-stream region delimited by these caustics, non-trivial vorticity patterns (Doroshkevich 1973; Chernin 1993; Pichon & Bernardeau 1999; Hahn et al. 2015). Thanks to a description of the phase-space 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 the 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 the shell-crossing. While convergence with the perturbation order of LPT remains a complex subject of investigation (Zheligovsky & Frisch 2014; Rampf et al. 2015; Michaux et al. 2021), it seems to take place at least up to the shell-crossing, not only for the three-sine-wave configurations we aim to examine (STC18) but also for more general, random initial conditions (Rampf & Hahn 2021), although the effects of cutoffs on power spectra remain to be investigated further in the latter case. Therefore, as long as the back-reaction from multi-stream flows on post-collapse dynamics remains negligible, the approximation of the dynamics we propose here should work shortly after the shell-crossing. This is the first step towards a proper analytical description of post-collapse dynamics in 6D phase-space.

This paper is organised as follows. In Sect. 2, we begin by introducing the basics of LPT and its recurrence relations, as well as its applications to initial conditions given by linear superpositions of trigonometric functions, such as the initial conditions of our sine waves. In Sect. 3, the important features of the Vlasov solver ColDICE are briefly summarised. We also address some aspects of measurement uncertainties in simulation data. In Sect. 4, we examine the phase-space structure and radial profiles at the shell-crossing with a comparison of analytic predictions with simulations, in the framework of singularity theory. In Sect. 5, the structure shortly after the shell-crossing is explored using the ballistic approximation by examining phase-space diagrams, caustics, and density and vorticity fields. Again, analytical predictions are compared with the Vlasov runs. Finally, Sect. 6 summaries the main results of this article. To complement the main text, Appendices AD provide details on the explicit form of high-order LPT solutions for the initial conditions of the sine waves up to the fifth order, the measurements in the Vlasov simulations, the predictions from quasi-1D LPT, and the expected structure of singularities at the collapse time, respectively.

Throughout the paper, we use the following units: a box size L = 1 and an inverse of the Hubble parameter at present time, , for the dimensions of length and time, respectively. In particular, the Lagrangian coordinate, Eulerian coordinate, velocity, and vorticity are explicitly expressed as q/L, x/L, υ/(LH0), and ω/Η0 = (L) × (υ/(LH0)), respectively.

2 Lagrangian Perturbation Theory

The Lagrangian equation of motion of a fluid element is given by (e.g. Peebles 1980) (1)

where the quantities x, a, and H(t) = a−1da/dt are the Eulerian comoving position, the scale factor of the Universe, and Hubble parameter, respectively. The derivative operator x = ∂/∂x is a spatial derivative in Eulerian space. The Newton gravitational potential ϕ(x) is related to the matter density contrast with being the background mass density, through the Poisson equation: (2)

In this framework, the velocity of each mass element is given by υ = a dx/dt.

Taking the divergence and curl of Eq. (1) with respect to Eulerian coordinates, Eqs. (1) and (2) can be expressed by the equivalent set of equations: (3) (4)

where the dot represents the Lagrangian derivative of time, d/dt.

In the Lagrangian approach, for each mass element, the Eulerian position x at the time of interest t is related to the Lagrangian coordinate (initial position) q through the displacement field Ψ(q, t) by (5)

The velocity field is expressed as υ(q, t) = a dΨ/dt. In the single flow regime (i.e. before the first shell-crossing time tsc), mass conservation implies , which leads to (6)

where the quantity J = det Jij is the Jacobian of the matrix Jij defined by (7)

and its inverse is given as (8)

While Eq. (6) is well defined only until the first shell-crossing time tsc, that is, the first occurrence of J = 0, Eqs. (7) and (8) can be formally used beyond tsc, as long as they are expressed in terms of Lagrangian coordinates (except that Eq. (8) might become singular).

2.1 Recursion Relation

In LPT, the displacement field, Ψ, is the fundamental building block and is considered to be a small quantity. It can be systematically formally expanded as follows, (9)

In what follows, we assume that the fastest growing modes dominate. They are known to be given to a very good approximation by (see e.g. Bernardeau et al. 2002, and references therein) (10)

where the purely time-dependent function D+(t) corresponds to the linear growth factor. The velocity field is then given by (11)

where function f(t) ≡ d1n D+/d In a corresponds to the logarithmic derivative of the growth factor. We note that the analyses performed in subsequent sections will consider the 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.

During early phases of structure formation, the Universe approaches Einstein-de Sitter cosmology and f is approximately given by . With the further approximation , as implicitly assumed in all the subsequent calculations (see e.g. Peebles 1980; Bernardeau et al. 2002; Matsubara 2015), by substituting Eqs. (5)(9) into Eqs. (3) and (4), one obtains simple recurrence formulae for the longitudinal and transverse parts of the motion (Rampf 2012; Zheligovsky & Frisch 2014; Rampf et al. 2015; Matsubara 2015): (12) (13)

in the 2D case, and (14) (15)

in the 3D case. Here and in what follows, we adopt the Einstein summation convention when the equation includes repeated lowercase Latin letters, and the subscripts i, j,… take 1, 2 or 3 (1 or 2) in the 3D (2D) case. In the above, we defined Ψi,j·≡ ∂Ψi/∂qj, and the tensors ɛij and are, respectively, the 2D and 3D Levi-Civita symbols. The symbol stands for a differential operator: (16)

In obtaining the recursion relations, we used the following formulae of linear algebra: (17)

for the 2D case, and (18)

for the 3D 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 (19) (20) (21) (22)

for the 3D case.

2.2 Two-sine-wave and Three-sine-wave Initial Conditions

Throughout this paper, we focus on specific initial conditions: two or three sine waves in a periodic box covering the interval [–L/2, L/2[: (23)

The parameters ϵi < 0 with |ϵx| ≥ |ϵy|≥ |ϵz| quantify the linear amplitudes of the sine waves in each direction. The initial time, tini, is chosen such that D+(tini)|ϵ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 = ϵyx and ϵ3D = (ϵyxZx), respectively, for two- and three-sine-wave initial conditions. These ratios will therefore be the quantities of relevance to define our initial conditions. In this setting, the initial density field, given by δini ≃ –q. Ψini = D+(tiniii| cos (2π/Lqi), presents a small peak at the origin, and mass elements subsequently infall towards the central overdense region. With the proper choice of ϵi, this initial overdensity can asymptotically match any peak of a smooth random Gaussian field (see e.g. Bardeen et al. 1986), which actually makes the three-sine-wave 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: (24) (25)

where the nth-order scalar coefficients and the nth-order vector coefficients 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 Eq. (23).

By imposing the conditions and , one can build up the perturbative solutions and using simple algebraic manipulations involving coefficients and : (26) (27)

These solutions lead to the nth-order displacement field given by . Using this solution for Ψ(n) as well as Eqs. (21) and (22), we subsequently construct the source of the (n + 1)th-order derivatives and . 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 the fifth order derived in this way1.

The above prescription is valid in 3D space and can be applied to the 2D case by cancelling all fluctuations along the z axis, that is, by performing 3D calculations with ϵz = 0 (i.e. ϵ3D = (ϵyx, 0)) for the three sine waves case. However, one can also realise that in 2D, vector coefficients become scalars , and that the solutions take the following form: (28) (29)

Our analytical investigations can easily cover a large range of values of ϵ2D and ϵ3D, while the simulations, much more costly, will only focus on three configurations, as detailed in Table 1, reflecting various regimes in the dynamics: quasi 1D with |ϵx| ≫ |ϵy,z|, anisotropic with |ϵx| > |ϵy| > |ϵz|, and what we design by axial-symmetric, with |ϵx| = |ϵy|( = |ez|), denoted by Q1D, ANI, and SYM, respectively2.

3 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 2D and 3D, 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 (Sect. 3.1), we address measurement uncertainties issues (Sect. 3.2).

Table 1

Parameters of the runs performed with ColDICE (Sousbie & Colombi 2016).

3.1 The Vlasov Code ColDICE

The Lagrangian coordinates defined in Sect. 2.1 correspond to the following unperturbed initial set-up, (30) (31)

where t0 is a virtual time corresponding to a = 0 and Q is the Lagrangian coordinate of each vertex. We note that we use capital letters to distinguish between vertex coordinates and the actual coordinates of fluid elements of the 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 Sect. 2.2: (32) (33)

with Ψini given by Eq. (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 ng = ns = 512 (see below for the definitions of ng and ns).

The Poisson equation is solved using the fast-Fourier technique in a mesh of fixed resolution ng. To estimate the projected density ρ(x) on this mesh, the intersection of each simplex of the phase-space sheet with each voxel or pixel of the mesh is computed exactly up to the 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 four-point 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 vertices corresponding, respectively, to and simplices in 2D and 3D. ColDICE allows for local refinement of the simplices following criteria based on local Poincaré invariant conservation as explained more in detail in Sousbie & Colombi (2016). Values of the refinement criterion parameter I used for our 3D runs are listed in detail in C21. For completeness, following the notations of C21, we used I = 10−8 for the 2D runs.

To perform local refinement, the phase-space sheet is locally described at the quadratic order inside each simplex with the help of three and six additional tracers per simplex in 2D and 3D, respectively. At the dynamical times considered in this work, which correspond at most to short periods after collapse, refinement is not triggered, except for a small number of simplices in the 3D 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 ng of the mesh used to solve Poisson equation and the initial number of vertices of the tessellation, ns3. 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 Sects. 4.2 and 5.2, the collapse time tsc shown in Table 1, the radial profiles used in Sect. 4.3, as well as the caustic network, the projected density and the vorticity analysed in Sect. 5.3.

3.2 Accuracy and Related Considerations

Accuracy and possible defects of ColDICE are discussed in Sousbie & Colombi (2016) and in C21 in the 3D case. In our work, precise determination of the collapse time is critical, since we consider, in the analyses performed below, either the collapse time or a moment very shortly after it. The uncertainty on the collapse time depends on the nature of the initial conditions and the parameters that control the accuracy of the simulations, in particular the number of simplices used to represent the phasespace sheet and the spatial resolution ng of the mesh used to solve Poisson equation. This latter turns out to be a very important parameter. Indeed, decreasing ng augments the collapse time, as discussed in C21. In Sect. B.2, relying on measurements in ANI-3SIN configurations with various spatial resolutions, we estimate that the quantity asc shown in Table 1 should be accurate at the fourth significant digit level or better for all configurations, for example asc = 0.02919 ± 10−4 for ANI-3SIN, except for SYM-3SIN, where the uncertainty on the measured asc could be a few 10−4. Indeed, as just mentioned above, accuracy on the determination of the collapse time also depends on the nature of initial conditions, the quasi-1D case and the triaxial-symmetric configuration being the easiest and the most difficult to deal with, respectively, as expected.

We note that accuracy regarding the collapse time also depends on the ability to determine the exact position of the shell-crossing point. Using the tessellation technique, the 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 centre of the system and the determination of the 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 centre of the forming halo would be moving, one expects to be able to accurately compute the collapse time thanks to the tessellation representation as long as sampling of the 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 ns, as well as the refinement parameter I limiting deviations from symplectic motion. As discussed in Sect. 3.1, for our simulations, parameter I does not influence the results except (probably marginally) for SYM-3SIN. Effects related to the choice of ns are discussed in Sousbie & Colombi (2016) and in C21. They should be negligible for our sine wave simulations compared with other sources of uncertainty on the collapse time for the values of ns we adopted. However, we did not test the accuracy of ColDICE in more complex cases, because we do not need to.

An additional advantage of our setting is that the axes of the sine waves are aligned with the mesh used to solve Poisson equation, which helps make the simulations more accurate: the same configurations in an oblique fashion would not achieve the same level of accuracy. This is illustrated by Appendix H of Sousbie & Colombi (2016) in the single sine wave case. As shown in this Appendix, anisotropies due to misalignments between preferred directions in the representation and preferred directions in the dynamics are the source of significant accuracy loss. However, these cumulative effects should remain small for ColDICE at the early times considered in this work.

Another critical issue is to determine accurately the shape of the caustics and then the density and vorticity fields inside the 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: L/nS, where L is the size of the simulation cube. This uncertainty, which adds up to the effects discussed above, is difficult to transpose to Eulerian space due to the variations in the strain tensor. It is however probably the main source of inaccuracy on the determination of the caustics in ColDICE. A way to improve greatly on this would consist in exploiting the quadratic representation inside each simplex, which we do not do here.

Turning to the fields, the vorticity, because it is composed of derivatives, is particularly challenging to measure even with the considerable gain brought by the tessellation technique compared with more traditional representations based on particles. Appendix B.5 explains in details how we measure density and vorticity, making use, this time, of the locally quadratic representation of the 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 or surfaces, which we keep in mind for the analyses carried out in Sect. 5.3.

4 Shell-crossing Structure

We are now in a position to study the structure of our protohaloes at the collapse time, tsc, and concentrate our investigations on phase-space diagrams and radial profiles, with comparisons of LPT pushed up to the tenth order to the Vlasov runs. This section is organised as follows. First, Sect. 4.1 presents the calculation of the collapse time itself. Indeed, this quantity depends on initial conditions and perturbation order, and the ability of LPT to provide an accurate determination of tsc is of prime importance. We discuss the extrapolation to the infinite order of the LPT series proposed by STC18 and generalise it to the 2D case. Second, Sect. 4.2 examines the convergence of LPT at collapse with phase-space diagrams, extending the earlier investigation of STC18 to the 2D case. For comparison, we also tested the formal extension of LPT to the infinite order and the predictions of the quasi-1D approach proposed by Rampf & Frisch (2017, hereafter RF17). Finally, in Sect. 4.3, LPT predictions and their convergence are studied in terms of radial profiles of the density, the velocity as well as the pseudo phase-space density, and put into perspective in relation to singularity theory.

4.1 Shell-crossing Time

This subsection presents estimates of the expansion factor at the first shell-crossing, which we refer to as the collapse time. In the following, the expansion factor will be formally identified to a time variable, still denoted by a to contrast with actual physical time t. We compare the value of the collapse time obtained from rath-order LPT and its extrapolation to the infinite order , as described in STC18, with the value asc measured in the Vlasov runs as explained in Appendix B.2. Our approach follows closely that of STC18. It is basically intended to repeat its main steps and to supplement it with additional discussions and comparisons with Vlasov runs in the 2D case.

Using rath-order LPT predictions, we explore the sequence of the shell-crossing times, , as a function of order n up to n = 10, by solving J(n) = 0 at the origin, where J(n) is the Jaco-bian of the nth-order LPT solution. Generally, except in the pure 1D case where first-order LPT is exact before collapse, becomes smaller with increasing perturbation order n (see e.g. RF17; Rampf & Hahn 2021, hereafter RH21, for recent works). As illustrated by Fig. 1, the shell-crossing time calculated at the nth-order with LPT is very accurately described by the following fitting form (STC18):

(34)

with e > 0. This fitting form, also used for each coordinate of the positions and velocities in Fig. 2 below, does not necessarily represent the sole choice for approximating the n dependence of the collapse time, but using the exponential of a power law might be the only way to match the convergence speed of LPT at large n, when considering quantities computed at the collapse time of each respective order. Equation (34) also implies when n ≫ 1, which might, at first glance, seem incompatible with the findings of RH21, who examined the LPT series in the slightly different context of a Gaussian random field in a periodic box. RH21 concluded from their analyses that convergence speed was asymptotically compatible with a power-law times an exponential of n. This would naively correspond to , while our measurements suggest values of e different from unity, ranging from, for example, e ~ 0.003 (SYM-3SIN) to e ~ 0.8 (Q1D-3SIN) for fitting and from 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 the shell-crossing times or equivalently . This additional n dependence introduced in the time variable fundamentally changes the nature of the LPT series as a function of n, which makes a direct comparison of Eq. (34) with the findings of RH21 inappropriate. Yet, it would be interesting to investigate the convergence properties of our systems seeded with sine waves using the approach of RH21.

One important thing to note is that convergence of the collapse time with order is rather slow, except in the quasi-1D case, which still requires at least the third order for reaching an approximately percent level of accuracy for approximate convergence, while a much higher order is required for other configurations, especially the 3D axial-symmetric case, ϵ3D= (1,1), for which convergence does not even seem to be achieved at the tenth order. Figure 1 thus demonstrates that low-order perturbation theory cannot be used to accurately estimate the collapse times (see also STC 18, in particular, Fig. 2 of this article, and RH21).

Despite these convergence issues, the extrapolated values of obtained with Eq. (34) are very accurate, as illustrated by Table 1. Indeed, the relative difference between and the value asc measured in the Vlasov simulations remains of the order of the percent or below, which suggests that the error on the extrapolated estimate of the collapse time remains small and should not affect significantly the analytical predictions of the phase-space structure and the radial profiles at the shell-crossing. Therefore, we use for the analytical predictions of the ‘exact’ collapse time in the rest of Sect. 4. As for the simulations data, we use the snapshot with the closest possible available value âsc of the expansion factor to asc, as indicated in Table 1.

thumbnail Fig. 1

Convergence behaviour of the shell-crossing time. Top: shell-crossing time calculated at nth-order LPT (dots), together with the fitting function given by Eq. (34) (solid lines), as a function of the perturbation order for various initial conditions, as indicated in the panel. The dashed horizontal lines present the values estimated from the simulations (see Appendix B.1), corresponding to the sixth column of Table 1. Bottom: relative error between the shell-crossing time calculated with nth-order LPT and the one obtained with the fitting formula.

4.2 Phase-space Structure

Figure 2 shows the (x, υx) subspace of the phase-space structure at the shell-crossing by considering the intersection of the phase-space sheet with the hyperplane y = 0 for two sine waves and y = z = 0 for three sine waves. For comparison, in addition to standard LPT results, we also present the formal extension of LPT to the infinite order as developed by STC18 and described in Sect. 4.1 for the shell-crossing time, as well as the analytical prediction obtained in the context of the quasi-1D approach developed by RF17, which we extended to the second order in the transverse direction (see Appendix C for details).

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 the x axis, and, until then, the dynamics is similar to the pure 1D case, in which 1LPT (the Zel’dovich solution) is exact until the shell-crossing. As illustrated by Fig. 2, the LPT prediction quickly converges with perturbation order n and describes the simulation results well even at the 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. We 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 1D collapse and LPT again quickly converges with perturbation order to the exact solution. This is clearly not the case for the 3D 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 (STC18) 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 the collapse time (e.g. Rampf 2019), and the analyses of RH21 suggest that the velocity blows up when convergence is lost. The exact properties of the spike we observe in our numerical data, in particular, whether the velocity diverges or remains finite, and, whether if finite, the velocity is actually smooth at the fine level, remain unknown. While this spike is not present in the LPT predictions at a finite order, it is well reproduced by the formal extrapolation to the infinite order. This is a hint that convergence of the LPT series might be lost at collapse for SYM-3SIN.

A question might arise whether, because of such a spike and because of their highly contrasted nature, close to axial-symmetric 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 with the non-axial-symmetric 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, we note that the formal extrapolation of LPT to the infinite order matches well (but not perfectly) the simulation results for all the configurations, as already found by STC18 in the 3D case. We also notice that the quasi-1D 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 fourth-order LPT, irrespective of initial conditions. Thus, as expected, the quasi-1D approach becomes less accurate when the ratio ϵ2D or ϵ3D approaches unity, particularly in 3D as already shown by STC18.

thumbnail Fig. 2

Phase-space structure for two- and three-sine-wave initial conditions at the collapse time: Q1D-2SIN (top left), ANI-2SIN (centre left), SYM-2SIN (bottom left), Q1D-3SIN (top right), ANI-3SIN (centre right), and SYM-3SIN (bottom right). 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, υx) subspace. Simulation results are compared with standard LPT predictions, which are supplemented with the blue line, denoted by’ΈΧΤ’, which corresponds to the formal extrapolation to the infinite order proposed by STC18 and sketched in Sect. 4.1 for the collapse time. For completeness, the quasi-1D approach (Rampf & Frisch 2017), denoted by Q1D, is also presented (see Appendix C for details).

4.3 Radial Profiles

We now focus on radial profiles, and define Eulerian polar and spherical coordinates for two- and three-sine-wave initial conditions by (x, y) = (rcosθ, rsinθ) and (x, y, z) = (rsinθcosϕ, rsinθ sinϕ, rcosθ), respectively. Then the angular averaged radial profiles of the density , velocity dispersion υ2, radial velocity dispersion , and infall velocity r are given by (35) (36) (37) (38)

where with r = |x| being the radial coordinate. In these equations, we used the angle average defined by (39)

for two-sine-wave initial conditions, and (40)

for three-sine-wave initial conditions. It is important to note that the angular coordinates θ and ϕ in the integrands are the Eulerian coordinates, and the integrands should be evaluated in terms of the Eulerian coordinate by solving the equation x = x(q, t).

To complete the analyses, we also study the pseudo phase-space density Q(r, t), defined by (41)

We additionally note that we present not only the radial profiles but also their logarithmic slopes, defined by (42)

where the quantity X stands for the radial profile under consideration.

Radial profiles at the shell-crossing are presented in Figs. 36. 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 respectively correspond to initial conditions with two and three sine waves, calculations are performed at the exact shell-crossing time, that is, the extrapolated value for LPT and the approximate value âsc for the simulations. Figures 5 and 6 are analogous, but examine LPT predictions of the rath order synchronised to their own shell-crossing time , which allows us to examine in detail the structure of the singularity at collapse produced at each perturbation order. It is important to note that from now on, we do not further examine the quasi-1D approach proposed by RF17. We also do not perform the formal extension to the infinite order (Eq. (34)), because calculating radial profiles in this framework was found to be too costly with the computational methods we are currently using.

The examination of these figures confirms the conclusions of the phase-space diagram analysis in Sect. 4.2. In particular, quasi-1D profiles (Q1D-2SIN and Q1D-3SIN) are well described by LPT predictions even at low order. At the ‘exact’ collapse, considered in Figs. 3 and 4, LPT predictions generally accurately describe the outer part of the halo, where density contrasts are lower, and then deviate from the exact solution when the radius becomes smaller than some value rmin(n) decreasing with increasing n, where n is the perturbation order. LPT can describe arbitrarily large densities, as long as n is large enough. In the results presented here, density contrasts as large as 100 or larger can be probed accurately by LPT, depending on the order considered and on the nature of initial conditions. Again, convergence worsens when approaching 3D axial-symmetry, the worst case being, as expected, SYM-3SIN with ϵ3D = (1,1). It is worth noticing that velocity profiles require significantly higher LPT orders than the density profiles to achieve a comparable visual match between theory and measurements. In particular, deviations between LPT predictions and simulations arise at quite larger radii rmin(n) for the velocities than for the density, with differences as large as an order of magnitude. With synchronisation, performances of LPT predictions greatly improve, as expected and as illustrated by Figs. 5 and 6. All perturbation orders predict strikingly similar density profiles, in close to perfect match with the simulations, except for SYM-3SIN discussed further below, while the velocity profiles still diverge slightly from each other except for quasi-1D 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. 36. The nature of the singularities of CDM structures at the collapse time and beyond has been widely studied in the literature, particularly in the framework of Zel’dovich motion (e.g. Zel’dovich 1970; Arnold et al. 1982; Hidding et al. 2014; Feldbrugge et al. 2018). In Appendix D, we re-derive the asymptotic properties of our protohaloes at small radii by Taylor expanding the equations of motion up to the third order in the Lagrangian coordinate q. For our sine-wave configurations, whatever the order n of the perturbation order3, three kinds of singularities are expected at the shell-crossing, as summarised in Table 2: the classic 1D 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 υ2, , and 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.

Figures 3 and 4 consider LPT predictions for perturbation order n calculated at the exact theoretical collapse time and not at the individual collapse time at this perturbation order, so the shell-crossing is not reached exactly, but gets nearer as n increases. Consequently, the asymptotic slope is only approached approximately, and better so with larger n. We note that convergence is slower for velocities than for the density, especially for axial-symmetric 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 from collapse, with âsc > asc, 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 non-linear 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, we point out that in this framework if crossed sine waves are considered generic approximations of local peaks of a smooth random Gaussian field, the probability of having exactly or is null, and hence one expects the classic 1D pancake to be exclusively dominant, from a pure mathematical point of view. However, this can be unrealistic at the coarse level, where 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 (top panels) and (bottom panels) for tenth-order LPT. Clearly, at sufficiently large radii, values of ϵ2D or ϵ3D close to and give similar results for the profiles, even for the logarithmic slope, which can approach that of the 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 − αr2 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 r2 in the Taylor expansion of trigonometric functions) provides significantly different results from pure spherical collapse. We note, though, that this state of fact is not fully proved by our simulations measurements, which as mentioned above, are slightly beyond the actual collapse time.

Third, we can compare, in the 3D case, the logarithmic slope of the pseudo phase-space density with that seen after violent relaxation and at late stages of the evolution of dark matter haloes, Q(r) ~ rβ (see e.g. Taylor & Navarro 2001; Navarro et al. 2010; Ludlow et al. 2010), which has been found to be close to the prediction of secondary infall model, βsph = 1.875 (Bertschinger 1985). This is also the case of our three sine waves haloes, as analysed in detail by C21, including 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 3D configurations. The slope prior to collapse in the ‘generic’ cases,β = –5/3, although of the same order of βsph, remains still different. Concerning the 2D case, which is different from the dynamical point of view, further analyses of the evolution of the haloes beyond collapse need to be done.

thumbnail Fig. 3

Radial profiles and their logarithmic slopes for the two-sine-wave initial conditions at the shell-crossing time. As indicated in top-left panel, LPT predictions are given by solid lines of various colours and Vlasov simulations results are represented by red dots. From left to right, the initial conditions are given by 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 υ2, the radial velocity dispersion , the infall velocity -υr, and the pseudo phase-space density Q(r), respectively. We 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 slope panels, the horizontal dashed lines correspond to the theoretical predictions at the small radii of Appendix D, as listed in Table 2.

thumbnail Fig. 4

Same as Fig. 3 but for the three-sine-wave initial conditions, from left to right, Q1D-3SIN [ϵ3D = (1/6,1/8)], ANI-3SIN [ϵ3D = (3/4,1/2)], and SYM-3SIN [ϵ3D = (1,1)], respectively. We note that in SYM-3SIN, the closest snapshot from collapse we had at our disposal from our Vlasov runs, âsc = 0.03190, is significantly beyond the actual shell-crossing time estimated by the method described in Sect. B.2, asc = 0.03155, which explains some discrepancies between the theory and the simulation at small radii.

thumbnail Fig. 5

Same as Fig. 3 but all analytical predictions evaluated at individual shell-crossing times computed at each perturbation order.

thumbnail Fig. 6

Same as Fig. 4 but all analytical predictions evaluated at individual shell-crossing times computed at each perturbation order.

Table 2

Summary of the logarithmic slopes at the shell-crossing obtained by singularity theory applied to LPT predictions at various orders (see Appendix D for details).

thumbnail Fig. 7

Radial profiles at the shell-crossing using tenth-order LPT in the vicinity of (top panels) and (bottom panels). In the top panels, ϵ2D is progressively varying in the range [0.9,1] with linear intervals of Δϵ = 0.1/64. In the bottom panels, ϵ3D progressively changes from (0.97,0.97) to (1,1) with linear intervals of Δϵ = 0.03/64. From left to right, we present radial profiles of the normalised density , the velocity dispersion υ2, the radial velocity dispersion the infall velocity –υr, and the pseudo 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.

5 Structure of the System Shortly After Shell-crossing

Once the 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 in a perturbative way for a short time after the shell-crossing, it is no longer, overall, a smooth function of time due to the formation of the singularity at the collapse time (Rampf et al. 2021). Still, it is reasonable to assume that very shortly after the 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 the collapse time, but not necessarily far beyond collapse (see also RH21). Therefore, in what follows, we shall use the LPT solution of the nth order at collapse as a starting point, and from there, the standard ballistic approximation, where velocity field is frozen, to study the structure of the system shortly after collapse.

In practice, one would like to use sufficiently large perturbation order so that convergence of LPT is achieved, and ideally the formal extrapolation to the infinite order proposed by STC18 and reintroduced in Sect. 4.1. However, this extrapolation is insufficiently accurate for the purpose of the calculations performed next, where very small time intervals after collapse are considered. Also, it is very costly from the computational point of view to exploit this extrapolation when solving the multi-valued problem intrinsic to the multi-streaming solutions. So from now on, we shall consider only higher-order LPT calculations and not their extrapolation to the infinite order (nor the quasi-1D approximation proposed by RF17).

This section is organised as follows. In Sect. 5.1, we introduce the ballistic approximation. We also relate, in the multi-stream regime, the Eulerian density and velocity fields to their Lagrangian counterparts, and introduce, in this framework, the vorticity field. Indeed, after the 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 the second order by our systems seeded by sine waves. Sect. 5.2 turns to actual comparisons of analytical predictions with Vlasov simulations using, similarly as in Sect. 4.2, phase-space diagrams, but shortly after collapseIn particular, we explore the limits of the ballistic approximation. Finally, Sect. 5.3 focuses on the overall structure of the multi-stream region in configuration space, by successively testing LPT predictions against simulations for the caustic pattern, the projected density and the vorticity fields.

5.1 Multi-stream Regime and Ballistic Approximation

Because the collapse time depends on the perturbation order, it only makes sense to test various perturbation orders shortly after the shell-crossing and compare them with simulations only if the collapse times are synchronised. In other words, in what follows, we consider the time for LPT of nth order and for the simulation, as listed in last column of Table 1. The value of Δa used for each sine waves configuration is given by the difference between the values of the last and sixth columns of the table. We recall that the quantity asc 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 the shell-crossing time of each perturbation order, we model the Eulerian coordinate after collapse as follows: (43) (44)

while the nth-order velocity field is frozen to its value at the shell-crossing time of each perturbation order,

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, which we denote as qF labelled by the subscript F = [1,…, nF]. Shortly after collapse, if ϵi ≠ ϵj for ij, there are at most three streams, nF ≤ 3, in a given point of space. For symmetric cases, for example ϵ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 in 2D and 3D, respectively, as discussed further in Sect. 5.2.

In what follows, we omit the time dependence for brevity. After defining the Lagrangian density and velocity fields, (45) (46)

the Eulerian density and velocity fields are expressed by the superpositions of each flow, (47) (48)

where the summation ΣF is performed over all the solutions qF of the equation x = x(qF).

While vorticity is not expected in single-stream regions, unless already present in the initial conditions (and in that case it corresponds to a decaying mode), its generation is one of the prominent features of the 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 2D and 3D, by (49) (50)

The vorticity fields in 2D and 3D space are scalar and vector quantities, respectively. Because the acceleration derives from a gravitational potential, local vorticity on each individual fold of the phase-space sheet cancels. Only the non-linear superposition of folds in the first term of Eqs. (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 a finite perturbation order (see e.g. Buchert & Ehlers 1993; Buchert 1994 and also Uhlemann et al. 2019 for a recent numerical investigation). This spurious vorticity is mainly coming from the second term in Eqs. (49) and (50), and can be especially noticeable in the 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 second-order 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 calculated shortly after the shell-crossing: (51)

with A = x, y, or z and stands for the growth factor at the shell-crossing time evaluated by using 1LPT. We note that in Eq. (51), the growth rate f and the expansion factor asc are evaluated at the same time as (we recall that f = 1 in the Einstein-de Sitter cosmology considered in this work). The value of is obtained in practice by solving the equation at the origin, keeping in mind that |ϵx|≥|ϵy,z|. We recall that Eq. (51) is exact in the pure ID case, that is, when ϵy= ϵz = 0. The condition for the caustics, J(1LPT) (q,a) = 0, is reduced to the relation (52)

Equation (52) implies that the equations of the caustics are given by q = q0, with q0 being a constant vector depending on time (i.e. on Δα). Therefore, the 1LPT caustics consist of ID lines in 2D and 2D planes in 3D, even at the collapse time. This configuration is actually degenerate and is not expected in realistic cases, for instance in the framework of a smooth Gaussian random field (see e.g. Pichon & Bernardeau 1999). Indeed, the regions where the first shell-crossing occurs should be composed, in the non-degenerate case (that is non-vanishing e;), of a set of points. In the vicinity of such points and shortly after the shell-crossing, even in 1LPT, the equation of the caustics should correspond, at leading order in q and in the non-axial-symmetric case, ϵiϵj ,ij, 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: (53)

where 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: (54)

In Eq. (53), the growth rate f and the scale factor asc are evaluated at the same time as . The additional terms compared with 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, α) = 0 (as long as none of the ϵi cancels). This leads to, in the non-axial-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 3D 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 the second order and above. We note that, as will be studied below, the axial-symmetric cases ϵ2D = 1 and ϵ3D = (1,1) are a bit more complex, as they require Taylor expansion of the motion beyond the quadratic order in q to obtain a full description of the correct and more intricate topology of the caustic surfaces (e.g. the fourth order in q instead of second for SYM-2SIN), as illustrated in 3D by the right panels of Fig. 8, due to simultaneous collapses along several axes. It would be beyond the scope of this article to go further into the analytic details of catastrophe theory for these very particular configurations, so we shall be content with a descriptive analysis of the LPT results in the axial-symmetric cases.

thumbnail Fig. 8

3D view of the expected caustic pattern shortly after the shell-crossing for three-sine-wave initial conditions. The left panels correspond to the most typical pancake, here embodied by ANI-3SIN, but Q1D-3SIN would look similar, while the right ones show the axial-symmetric case, i.e. 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.

thumbnail Fig. 9

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

5.2 Phase-space Structure

Figures 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 with the simulations.

Interestingly, Fig. 9 shows that the ballistic approximation can deviate significantly from the exact solution even quite shortly after collapse. While it remains precise for all configurations in 2D, with a slight worsening of the accuracy when going from the 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 phase-space sheet with the hyperplane y = z = 0 is composed, in the multi-stream region, of three curves with an ‘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 with the ‘S’ with the largest amplitude along the 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 three4.

The results shown in Fig. 9 must however be interpreted with caution. Strictly speaking, we aim to apply the ballistic approximation from the collapse time, which is not exactly the case in this figure. The red curves, which correspond to the snapshots used to originate the ballistic motion, all correspond to a time abc slightly before the actual collapse (namely abc = 0.03103 for SYM-3SIN and , as given in Table 1 for other cases) because we do not have access to the exact moment of the actual collapse. The level of proximity to collapse is, at least from the visual point of view, the lowest for the axial-symmetric cases, particularly SYM-3SIN. In this case, the ratio Δaeff /abc ≃ 0.03 is the largest, where Δaeff is the difference between apcasc + Δa and abc and is used to test the ballistic approximation. We also note that Δaeff /abc ≃ 0.03 is also of the same order for ANI-3SIN. The greater the value of Δaeff /abc, the greater the deviation due to non-linear corrections of the motions to be expected. Furthermore, we also expect the ballistic approximation to worsen from quasi-ID 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 the x axis (local extremum of x coordinate). This suggests that the effects of acceleration in the vicinity of collapse are strong, which implies that the ballistic approximation can be applied only for a very short time (or it could be that it is not applicable, mathematically, due to the strength of the singularity).

Now we turn to the comparison of LPT predictions of various orders with the simulations, as shown in Fig. 10. We recall that the ballistic approximation is applied from the respective collapse times of each perturbation order, which obviously makes predictions of LPT artificially more realistic than it should be.

The conclusions of Sect. 4.2, where we compared LPT with simulations at collapse, still stand at the coarse level, obviously, since the time considered in Fig. 10 is nearly equal to the collapse time. When zooming on the central part of the system, we notice that the ballistic approximation applied to LPT works increasingly well with the order, as expected, especially when approaching quasi-ID dynamics. It can fail in the tails of the ‘S’ shape of the phase-space sheet, as a combination of limits of LPT to describe the system in the vicinity of the velocity extrema and the limits of the ballistic approximation itself just discussed above.

It is important to note at this point that the ballistic approximation is not necessarily the best choice outside the multi-stream region, where LPT predictions can still perform very well. This is illustrated in Fig. 10 by the cyan curves, which show tenth-order LPT results without assuming ballistic approximation; they should be compared with the dark purple curves, which correspond to tenth-order 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 tenth-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.

thumbnail Fig. 10

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

thumbnail Fig. 11

Caustic pattern shortly after the shell-crossing in the 2D case: comparison of LPT predictions using ballistic approximation with Vlasov runs. The left, middle, and right panels correspond, respectively, to Q1D-2SIN, ANI-2SIN, and SYM-2SIN configurations, while the top and bottom panels show the caustics in Lagrangian and Eulerian spaces, respectively.

5.3 Configuration Space: Caustics, Density, and Vrticity

We now enter into the heart of this work, which consists of examining in detail the structure of the system shortly after collapse in configuration space, both in Lagrangian and Eulerian coordinates. The 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 with the Vlasov code for our various sine waves initial conditions. Then, we turn to the normalised density and the vorticity in Figs. 13 to 18, and compare LPT predictions at the second and tenth order with the simulations. We remind the reader again that LPT predictions with the ballistic approximation are computed by synchronising the respective collapse times obtained at each perturbation order, which obviously artificially improves the performances of LPT. We cautiously note that in plotting the results below in 3D cases, the caustics outputs (Figs. 11 and 12) are shown at the intersections with the x = 0, y = 0, and z = 0 planes, while the density and vorticity outputs (Figs. 14 and 1618) 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 ij, that is the non-axial-symmetric case, both in 2D and in 3D. We note that in the 3D case, to have a clearer view of the caustic pattern, we consider the intersection of the caustic surfaces with the planes x = 0, y = 0 and z = 0 (we recall, however, that a 3D view of the caustic surfaces was given for 2LPT in Fig. 8). In the case of the simulations, as explained in Sect. 3.2 and Appendix B.3, the caustic pattern corresponds to a tessellation, that is, a set of ensembles of connected segments in 2D and of connected triangles in 3D. The intersection of a tessellation of triangles with a plane also gives a set of ensembles of connected segments. In the figures, only the vertices of these segments are represented, and the initial regular pattern of the tessellation used to represent the phase-space sheet is clearly visible in Lagrangian space.

Figures 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 the tenth 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 significantly underestimated by LPT, and it seems that increasing the perturbation order to arbitrary values is not going to improve the agreement with the simulations. The situation is much worse in 3D, particularly in Eulerian space, although the simulation measurements are very noisy, which makes comparison with theoretical predictions difficult. We know from the previous paragraph that these mismatches are at least partly attributable to limits of the ballistic approximation, along with limits of LPT to be able to describe the velocity field at collapse in the vicinity of its extrema, particularly the spike observed on the 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. We 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 Sect. 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, ij, agrees perfectly with the predictions of Pichon & Bernardeau (1999) based on Zel’dovich dynamics5. In 2D, the vorticity field is a scalar that can be decomposed into the four sectors inside the caustic, two of positive sign and two of negative one. In 3D, the vorticity field is a vector. Because the shell-crossing takes place along the x axis, each component of this vector field has specific properties, which are related to variations of the velocity and density field along with the phase-space sheet. In particular, it is easy to convince oneself that if collapse happens along the x axis, the strongest variations of all the fields are expected along this direction. That means, since each coordinate of the vorticity vector depends on variations in the velocity field in other coordinates, that the magnitude of ωx is expected to be small. Due to symmetries, we also expect that ωy and ωz are, respectively, an odd function of z and y, which means that ωy approaches zero when approaching the z = 0 plane, and similarly for ωz when approaching the y = 0 plane, which explains the pattern of the vorticity field in Figs. 16 and 17. These symmetries impose us to perform measurements on slices slightly shifted from the centre of the system, as detailed in the caption of Fig. 14. We also note that in these figures the pattern is analogous to 2D for the y and z coordinate of the vorticity fields in the xz and xy 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.

thumbnail Fig. 12

Structure of caustics shortly after the shell-crossing for three-sine-wave initial conditions: comparison of LPT predictions using ballistic approximation with Vlasov runs. The top six, middle six, and bottom two panels respectively correspond to Q1D-3SIN, ANI-3SIN, and SYM-3SIN. In each group of six panels, the top and bottom lines correspond to Lagrangian and Eulerian spaces, and the intersection of the caustic surfaces with the plane are qx = 0, qy = 0, and qz = 0, respectively, for the first, second, and third column of each group. In the bottom group of two panels, the left and right panels correspond respectively to Lagrangian and Eulerian space and only show the intersection of the caustic surface network with the plane qz = 0 due to the symmetry of the system.

thumbnail Fig. 13

2D density shortly after the shell-crossing: comparison of LPT predictions using ballistic approximation with Vlasov runs. From top to bottom: Q1D-2SIN, ANI-2SIN, and SYM 2SIN. From left to right: 2LPT, 10LPT, and measurements in Vlasov simulations.

thumbnail Fig. 14

Slices of the projected density shortly after the shell-crossing: comparison of LPT using ballistic approximation with Vlasov runs. The left and right groups of nine panels correspond respectively to Q1D-3SIN (top, middle, and bottom line of panels: x=−1.55×l0−4, y = −1.17×10−3, and z= −1.56 × 10−3 slice) and ANI-3SIN (x= −5.16 × 10−4, y = −2.15 × 10−4, and z= −5.47 × 10−4 slice), while the bottom group of three panels corresponds to SYM-3SIN (z= −1.17 × 10−5 slice). In each group of panels, the left, middle, and right columns give respectively the second-order LPT prediction, the tenth-order LPT prediction, and the Vlasov code measurements. Due to the symmetry of the system for SYM-3SIN, only one slice is shown for the bottom panels.

6 Summary

In this paper, following in the footsteps of Saga et al. (2018), we have investigated the structure of primordial CDM haloes seeded by two or three crossed sine waves of various amplitudes at and shortly after the shell-crossing, by thoroughly comparing LPT, up to the tenth order, with 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 the shell-crossing and, second, to the phase-space structure, caustics, and density and vorticity fields shortly after the 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 the phase-space sheet structure locally at the quadratic level.

We studied three qualitatively different initial conditions characterised by the amplitude of three crossed sine waves, as summarised in Table 1: quasi 1D (Q1D-2SIN, Q1D-3SIN), where one amplitude of the sine waves dominates over the other(s), anisotropic (ANI-2SIN, ANI-3SIN), where the amplitude of each wave is different but remains of the same order, and axial-symmetric (SYM-2SIN, SYM-3SIN), where all amplitudes are the same. In order to predict the protohalo structure shortly after the shell-crossing in an analytical way, we used the 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 the shell-crossing, with a phase-space structure along the x axis analogous to what is obtained in 1D. This pancake is reproduced by LPT well, which converges increasingly well to the exact solution when approaching quasi-1D dynamics, as expected. The local extrema of the velocity field around the singularity are the locations where LPT differs most from the exact solution, underestimating the magnitude of the velocity. Convergence with the perturbation order becomes very slow when approaching 3D axial-symmetry, where spikes appear on the velocity field on each side of the singularity;

  • Phase-space diagrams at collapse: Except for SYM-3SIN, one expects the system to build up a classic pancake singularity at the shell-crossing, with a phase-space structure along the x axis analogous to what is obtained in 1D. This pancake is reproduced by LPT well, which converges increasingly well to the exact solution when approaching quasi-1D dynamics, as expected. The local extrema of the velocity field around the singularity are the locations where LPT differs most from the exact solution, underestimating the magnitude of the velocity. Convergence with the perturbation order becomes very slow when approaching 3Daxial-symmetry, where spikes appear on the velocity field on each side of the singularity;

  • Radial profiles at collapse: With a sufficiently high order of perturbation, LPT can reproduce arbitrarily high density contrasts at collapse, but we note a slower convergence when turning to velocity profiles. Still, the convergence of LPT is sufficiently good to probe the asymptotic logarithmic slope at the small radii expected at collapse for various profiles from singularity theory, as summarised in Table 2 and confirmed by simulation measurements. These profiles are ρ(r) ∝ r−2/3 for generic initial conditions (i.e. Q1D and ANI), ρ(r) ∝ r−4/3 for SYM-2SIN, and ρ(r) ∝ r−2 for SYM-3SIN, while velocity profiles always present the same power-law behaviour, for example υ2(r) ∝ r2/3. Confirming predictions of singularity theory at collapse and explicit convergence to asymptotic profiles at small radii was expected but is not trivial;

  • Synchronisation: Agreement with singularity theory predictions is made even more explicit by computing LPT profiles at the respective collapse times computed at each perturbation order. In this case, convergence to the 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 the 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 tenth-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 in a more complete calculation that would take into account the feedback due to the force field, as first proposed in the 1D case by Colombi (2015) and Taruya & Colombi (2017) (see also the recent work of Rampf et al. 2021), and formulated in terms of post-collapse perturbation theory. The next step is to implement an LPT approach where the small parameter is the interval of time Δa = aasc from the 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-axial-symmetric case, the system presents an ‘S’ shape in x – υx subspace, similar to the 1D case where positions and velocities can be approximated as third-order polynomials of the Lagrangian coordinate q. The calculation of the force field requires a three-value problem to be solved in the multi-stream region, similar to the 1D case. While technically challenging, generalisation of post-collapse perturbation theory to 2D and 3D seems possible. It might provide significant insights on non-linear corrections due to multi-stream dynamics on large-scale structure statistics, for example predictions of the power spectrum of the large-scale matter distribution from higher-order perturbation theory.

thumbnail Fig. 15

2D vorticity fields: comparison of LPT predictions with Vlasov runs.

thumbnail Fig. 16

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

thumbnail Fig. 17

Same as Fig. 16, but for ANI-3SIN. The 2D slices considered are the same as in the top-right group of nine panels of Fig. 14.

thumbnail Fig. 18

Same as Fig. 16, but for SYM-3SIN. Due to the symmetry of the initial conditions, only an x-y slice is shown, which is the same as in bottom three panels of Fig. 14.

Acknowledgements

This work was supported in part by JSPS Research Fellow Grant-in-Aid Number 17J10553 (SS), MEXT/JSPS KAKENHI Grants Numbers JP17H06359, JP20H05861, and JP21H01081, JST AIP Acceleration Research grant JP20317829 (AT), ANR grant ANR-13-MONU-0003 (SC), as well as Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES (SC). Numerical computation with ColDICE was carried out using the HPC resources of CINES (Occigen supercomputer) under the GENCI allocations c2016047568, 2017-A0040407568 and 2018-A0040407568. Post-treatment of ColDICE data were performed on HORIZON cluster of Institut d’Astrophysique de Paris. Calculations of theoretical predictions have made use of the Yukawa Institute Computer Facility.

Appendix A Expressions of the LPT Solutions

In this appendix, we present the LPT solutions up to the fifth order, which are obtained by solving the recursion relations given in Eqs. (21) and (22). Since higher-order solutions are straightforwardly derived in the same way, we do not explicitly show them here.

The results are partly presented in Moutarde et al. (1991) up to the third order, taking into account the contributions other than the fastest growth mode (see also Buchert et al. 1997, for the solutions including decaying modes). Here, we first explicitly show the analytical solutions of the sine waves initial conditions up to the fifth order.

As shown in Sect. 2.1, the fastest growing mode can be expanded as follows: (A.1)

For presentation purposes, we only show the x-components of the displacement field, that is, with the definition (ϵ1, ϵ2) ≡ (ϵyx, ϵz/ϵx). Given the x-components, the y- and z-components can be derived by permutating all x, y, and z: (A.2) (A.3) (A.4) (A.5) (A.6)

We note that because of the underlying symmetry of the initial conditions, the x-components of the LPT solutions are symmetric under the exchange of yz and ϵ1ϵ2.

Appendix B Measurements on the Tessellation of Vlasov Simulations

This appendix explains in detail how measurements of various quantities studied in this article are performed on the output of ColDICE, which consists of a tessellation of the phase-space sheet with simplices, respectively triangles and tetrahedra in 4D and 6D phase-space.

B.1 Phase-space Diagrams

The intersection of a hypersurface of dimension D1 = D with a hyperplane P of dimension D2 = D + 1 inside a phase-space of dimension D3 = 2D is of dimension D2 + D2D3 = 1. In other words, this intersection corresponds in the non-trivial case to a set of curves. Therefore, in 2D, the intersection of the phase-space 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 employed an approach valid at the linear order, consisting in simply computing, for each simplex, its geometric intersection with the hyperplane P, which is empty, a point, or a segment. As a result, phase-space diagrams extracted from the Vlasov runs consist of sets of connected segments, the ends of which are plotted in Figs. 2 and 10.

B.2 Expansion Factor at Collapse: asc

An accurate estimate of the value of the expansion factor corresponding to the collapse time is essential when studying the properties of the system at the shell-crossing and shortly after. The curves generated in Appendix B.1 can be used for this purpose, at two expansion times ai, 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, one of the ends of which is as close as possible to the origin and the other, with abscissa xi, has υx > 0; the magnitude of the other coordinate(s) of the velocity, in theory null, are as small as possible. The latter condition excludes, in the axial-symmetric case, the component(s) of the flow corresponding to simultaneous collapse along the y or/and z direction, which has/have a non-zero value of υy or/and υz. Once this segment is identified at both times, corresponding to a1 and a2, a simple linear interpolation is used to estimate the expansion factor at collapse: (B.1)

This formula does not require ai and a2 to be just before and after the collapse time to estimate asc, but, to provide sufficiently accurate results, the times need to be sufficiently close to the actual collapse time. A first guess of the collapse time is estimated by using perturbation theory predictions extrapolated to the infinite order (STC18), as listed in 5th column of Table 1. Examination of Table 1 shows that the value predicted this way agrees extremely well with the measurements in the Vlasov runs provided by Eq. (B.1), and shown in sixth column of the table.

We note, however, that, beyond the approximate nature of the linear approximation underlying Eq. (B.1), the collapse time estimate - in addition to other estimates, as discussed in Sect. 3.2 -can be significantly influenced by the force resolution, that is, the value of ng. As discussed in detail in C21, decreasing force resolution delays the collapse time; hence, to have an accurate estimate of the collapse time, a sufficiently large value of ng is required. We performed extensive force resolution tests for the 3D simulation with ϵ3D = (3/4,1/2). Our measurements of asc with the above method provide asc = 0.02912 (in excellent agreement with the predicted value and 0.0293 respectively for ng = 1024, 512 and 256. Our ng = 512 simulation thus provides, for this value of ϵ3D, an estimate of the collapse expansion factor accurate at approximately the 10−4 level, and we expect this to apply as well to the quasi-1D 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 centre of the system. On the contrary, in the 2D case, given the high value of ng = 2048 used to perform the simulations, we expect our estimates of collapse time to be very accurate, at an order better than the 10−4 level, but we did not test this explicitly.

B.3

Caustics are regions where the determinant J of the Jacobian matrix changes sign. At the linear order in the local description of the 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, and hence the caustics estimated this way are given by sets of segments, the ends of which are shown in Fig. 11. In 3D, the caustics are given by sets of triangles, the intersection of which we compute with the y = z = 0 plane to again get a set of segments; their extremities are shown in Fig. 12. Because we are using a leading order approach, the caustic lines or surfaces are not necessarily smooth but should trace accurately enough the actual caustics for the purpose of this work (see also Sect. 3.2).

B.4 Radial Profiles

To measure radial profiles in logarithmic bins, each simplex is replaced with a large number of particles as explained for the 3D case in Appendix A2 of C21. We refer to this work for the reader interested in the details of this procedure, which we straightforwardly generalised to the 2D case.

B.5 Density Field and Vorticity: From Linear to Quadratic Order

In this section we aim to compute the following five quantities: (i) the Jacobian of the transformation between initial and final positions, J(q), as defined in Eq. (7); (ii) the Lagrangian projected density, ρL(q), defined in Eq. (45); (iii) the total Eulerian projected density, which stems from the superposition of one or more folds of the phase-space sheet, as described by Eq. (47); (iv) the Eulerian velocity field, given by Eq. (48); and (v) the vorticity, ω2D and ω3D, defined in Eqs. (49) and (50). Again, because the acceleration derives from a potential, local vorticity on each phase-space sheet fold cancels. Only the non-linear superposition of phase-space sheet folds contained in the first term of these equations induces non-zero 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 Sect. B.5.1 we introduce barycentric coordinates, which are useful to define the position of any point inside each simplex. In Sect. B.5.2 we show how the barycentric coordinates can be used to perform calculations at the linear level inside each simplex, in particular partial derivatives of a function, as already discussed for instance by Hahn et al. (2015). In Sect. B.5.3 we generalise the procedure to the case when a quadratic description of the phase-space sheet is available. Finally, in Sect. 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 with 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 Xk, k = 1,…, Ns of the simplex vertices, with Ns = D + 1 where D is the dimension of configuration space, and a function g(X), one can define the following linear interpolation: (B.2)

for (B.3)

where the barycentric coordinates ξk are positive quantities verifying ξk = 1, ξk’k = 0, for X = Xk and Σk ξk = 1. When working in Lagrangian space, the space of initial positions, the phase-space sheet is flat and the linear interpolation, (B.4)

is exact. In what follows, we shall therefore use Eq. (B.4) to define barycentric coordinates. With this definition of this means that at given time t, Eq. (B.3) remains valid, but only at the linear level, since dynamical evolution of the 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 Eq. (B.2). The derivative of the function inside each simplex can thus be written (B.5)

where M is the matrix of the partial derivatives of the vector P ≡ (Σk ξk, X1, X2, X3) with respect to ξk, with X given by Eq. (B.3). In other words, (B.6) (B.7)

with Xk,α the αth coordinate of kth vertex with position Xk. At the linear level, the derivative defined by Eq. (B.5) is therefore simply constant inside each simplex.

This calculation can be performed in Lagrangian space to directly obtain the Jacobian of the transformation between the initial and present position, Eq. (7), from the Jacobian matrix, (B.8)

where W is defined similarly as M but in Lagrangian space, (B.9) (B.10)

with Qk,α the αth Lagrangian coordinate of kth vertex with Lagrangian position Qk.

This means that the linear description allows one to obtain the Jacobian J(q) and the corresponding Lagrangian density ρL(q) = 1/|J(q)| as constant quantities inside each simplex. Unfortunately, this is not sufficient, since Eqs. (49) and (50) involve partial derivatives of the projected density. Rigorously speaking, this means that a representation of the simplices at the quadratic level is required.

B.5.3 Calculations of Various Quantities at the Quadratic Level

Fortunately, ColDICE employs a quadratic description of the 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 with vertices (k, k’) corresponds to barycentric coordinates ξk = ξk’ = 1/2, ξj, = 0 for jk and jk’.

This allows one to define in a unique way a quadratic hyper-surface inside each simplex coinciding with vertices and tracers belonging to it. We 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 Nk, (B.11) (B.12)

where functions K(k) and L(k) are appropriately chosen to cover all the combinations 1 ≤ K(k) ≤ L(k) ≤ Ns, and where Ns = D + is the number of simplices and Nt = D(D + 1)/2 is the number of tracers, given D the dimension of space. Then Eq. (B.2) becomes (B.13)

where gk is now also defined over tracers. This is true in particular for quadratically interpolated positions, (B.14)

while corresponding Lagrangian coordinates Q are still given by the linear representation (B.4) since the phase-space sheet is initially flat, which means that matrix W does not change.

From Eq. (B.14), we obtain the matrix Mquad of partial derivatives of the vector P ≡ (Σk ξk, X1, X2,X3) with respect to ξk, which replaces Eqs. (B.6) and (B.7): (B.15) (B.16)

with Xk’,α the αth coordinate of k’th vertex/tracer with position Xk’.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 Nk.

The calculation of the Jacobian matrix exploiting the quadratic representation then just consists in replacing matrix M with Mquad in Eq. (B.8), since matrix W remains unchanged, which allows us to compute J(q) in the quadratic representation. In particular, we can compute Jk = J(Qk) at each vertex and tracer position to prepare it for the quadratic interpolation procedure given by Eq. (B.13) (but staying aware of the fact that this calculation is accurate only up to the linear order), which, in addition to the 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 Eq. (B.13) by (B.17) (B.18)

This quantity depends in a non-linear way on barycentric position and requires matrix Mquad to be invertible, so to be able to have a finite estimate of the partial derivative, it is necessary to avoid the subspace occupied by the caustics (or, more exactly, the element of surface or curve approximating the actual caustics).

In particular, we estimate the partial derivatives of the projected density as follows (B.19)

using Eq. (B.18) to estimate J/ xß. Indeed, the Jacobian is a smooth function of the Lagrangian coordinate. This is not the case for function ρL(q) = 1/|J(q)|, which presents strong variations in the vicinity of caustics. Using the Jacobian instead of the density to estimate derivatives allows us to capture better asymptotic behaviours in the vicinity of the caustics.

B.5.4 Practical Estimate: Particle Representation and Local Coarse Graining

In order to compute any quantity at position x0, one has to compute the intersection of the phase space sheet with the position x = x0, or to resolve inside each simplex the equation Χ(ξ) = x0. While this is a straightforward linear problem in the linear representation of the phase-space sheet (Eq. B.3), it becomes non-trivial in the quadratic case (Eq. B.14), for which some iterative procedure seems necessary. Here, to simplify the approach (and also to reduce effects of aliasing and divergences near caustics), we adopt a forward point of view consisting in projecting the tessellation on a Cartesian mesh of resolution nana in configuration space, similarly as it is performed by ColDICE to compute the projected density. This means that we apply additional coarse graining over each pixel or voxel of the Cartesian mesh: Integrals of the form (B.20)

become (B.21)

with υpixel and υvoxel being respectively the area of the pixel and the volume of voxel under consideration.

While Sousbie & Colombi (2016) use actual ray-tracing exact to the linear order to compute such an integral, we replaced each simplex with a large number of particles, which allows us to give an account of the quadratic shape of the simplices in a very simple way. However, employing particles introduces some discrete noise. To reduce the noise, we used the cloud-in-cell (CIC) interpolation procedure on the target Cartesian mesh (Hockney & Eastwood 1988), that is to say, each particle is replaced by a pixel (voxel) of the same small size as the pixels (voxels) of the Cartesian mesh, and we associate a weight Winter with the particle proportional to the volume of intersection between the pixel (voxel) associated with the particle and the pixel (voxel) of the mesh. In addition, we do not place particles at random within the simplex, but instead refine the simplex hierarchically in an homogeneous fashion max times, and then replace each 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 or pixel of the Cartesian mesh are kept under control: (B.22)

where D is the dimension of space, Δ is the step of the Cartesian mesh and (B.23)

with Xk,α the αth coordinate of kth vertex with position Xk and likewise for Xk’,α.

In practice, prior to CIC interpolation, each particle p samples a small element of volume δVp in the spatial part of the integral of right hand side of Eq. (B.21), (B.24)

where VL is the Lagrangian volume of the simplex and Jp the Jacobian (quadratically) interpolated at position of the particle. So, at the end, the contribution of each particle in right hand integral of Eq. (B.21) is equal to δVp Winter g(Xp, Vp), where we recall that Winter is the CIC weight defined previously and where function g is estimated at particle phase-space position (Xp, Vp) using the methods described in previous paragraphs.

We note that coarse graining due to the CIC interpolation is expected to introduce defects or biases in regions where large variations are expected, that is, in the vicinity of caustics. Obviously these effects become more dramatic when differentiating quantities, so one expects the vorticity to be affected numerically nearby the caustics, as can be noted in Figs. 15, 16, 17 and 18.

Appendix C Perturbative Treatment of Quasi-1D Flow

In Sect. 4.2 we test theoretical predictions relying on the pertur-bative treatment of a quasi-1D flow as proposed by RF17. In this appendix, we briefly present this approach, focusing mainly on three-sine-wave initial-condition configurations. Generalisation to the two sine waves case is straightforward.

When the system is exactly 1D, the first-order LPT solution is exact before the shell-crossing. This fact leads to an approximate treatment in 3D space. We exploit the exact solution of the 1D problem along the x axis given by Zel’dovich approximation as the unperturbed zeroth-order state: (C.1)

with the subscript i = x, y for two sine waves and i = x, y, z for three sine waves. The transverse fluctuations are considered as small first-order perturbations to this set-up: (C.2)

In the case of the three sine waves, we assume the explicit form of the functions and Ψ(1,1) to be (C.3) (C.4)

Then, substituting these initial conditions into the recursion relations in Eqs. (21) and (22), we obtain (C.5) (C.6)

By solving Eqs. (C.5) and (C.6), one can construct the Q1D firstorder solutions.

We further extend the Q1D treatment proposed by RF17 up to the second order in the transverse fluctuations as described below. On top of the zeroth- and first-order perturbations, we define the second-order perturbation as (C.7)

Again, using Eqs. (21) and (22), we obtain the following recursion relation: (C.8) (C.9)

The recursion is initialised by (C.10) (C.11)

Using the above recursion relations, the perturbative expansion is then performed by keeping terms proportional to and up to the second order, n + m = 2, while keeping terms proportional to ϵxk up to the tenth order, k = 10, as shown in Fig. 2. This approach provides a very accurate description of the dynamics when the system is initially quasi 1D (i.e. |ϵx| » |ϵy,z).

The results of the Q1D solution at the first order in the transverse fluctuations are presented in RF17. In what follows we show (for the first time) the Q1D solution at the second order in the transverse fluctuations and up to the fifth order in the longitudinal direction in terms of the growth factor (i.e. Ψ(2,5)). We use the same notations as in Appendix A.

For the x-components of the Q1D solutions, we have (C12) (C13) (C14) (C15) (C16) (C17) (C18) (C19)

We note that the y-components and z-components are symmetric under the exchange of qyqz and ϵyϵz in the above expressions.

Appendix D Asymptotic Structure of the Singularity at Collapse

In this appendix, we analytically investigate the dependence on initial conditions of the logarithmic slopes of various profiles expected at the shell-crossing and small radii, as discussed in Sect. 4.3. The main results of this investigation are summarised in Table 2. Since the extension to the 2D case is obvious by simply setting ϵz = 0, we present here the structure of the singularity obtained from LPT to some order for the three sine waves case. We note that the calculations presented in this appendix are conceptually not new since singularity theory applied to cosmology is already well known (see e.g. Arnold et al. 1982; Shandarin & Zeldovich 1989; Hidding et al. 2014; Feldbrugge et al. 2018).

D.1 Taylor Expansion Around the Singularity

To facilitate the analysis, we focus on the relation between Eulerian and Lagrangian coordinates around the origin. Expanding the relation (5) between the Eulerian coordinate x and the Lagrangian coordinate q around the origin, we obtain, after neglecting O(q4) and higher-order terms, (D.1)

with Ãij(t) ≡ δij + Aij(t), where Aij(t) is some function of time. In these equations, we have exploited the symmetric nature of the three-sine-wave 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 that verify P(−q) =P(q) (see e.g. Colombi 2015; Taruya & Colombi 2017, for the 1D case).

The coefficients Aij and Cijkl are expressed in terms of partial derivatives of the displacement field with respect to the Lagrangian coordinates at the origin as follows, (D.2)

Using Eq. (D.1), the Jacobian matrix is given by (D.3)

The shell-crossing time tsc is obtained by solving the equation J(0, tsc) = 0, that is, (D.4)

Hereafter, the dependence on the collapse time tsc will be omitted in the notations.

Using Eq. (D.3), the Jacobian is given by (D.5)

Further simplifying the calculations, we noticed that, up to the tenth order of the LPT development, for an initial displacement field given by three orthogonal sine waves aligned with each coordinate axis, the matrix is diagonal and the matrix Cijkl verifies (D.6)

We did not find any simple way to demonstrate that these properties stand at any order, but the fact that they are verified up to tenth-order LPT is a really convincing clue. In this case, the shell-crossing condition, Eq. (D.4), is reduced to Ãxx Ãyy Ãzz = 0. Prior to the shell-crossing, the Eulerian coordinates read (D.7)

for A, B = x, y, z, and the Jacobian (D.8)

In this last expression, we did not exploit yet explicitly the symmetries on Ci jkl; matrix. On the basis of these equations, we can now investigate the slope of the density profile for the three types of singularities we consider, depending on the values of the eigenvalues of matrix Ãij.

D.2 Asymptotic Behaviour of the Profiles at Collapse

We now analyse the asymptotic behaviour of the profiles around the origin. To this end, when we consider the following scaling, xs x, implying rsr, we examine how the Lagrangian coordinate changes. By doing so, we can understand the behaviour of the scaling of the Jacobian in terms of the Lagrangian coordinates, and thus reveal the behaviour of the density profile at the origin.

It is important to note that this proof is somewhat simplified and ignores details on the angular dependence when performing integrals over spherical shells to obtain the radial profiles. Therefore, the proof given in this appendix is not mathematically rigorous, but it leads to the same conclusions as the exact calculations in which proper form factors are estimated. The purpose of this section is indeed to provide a simplified rephrasing of singularity theory already presented in a more rigorous fashion in other works (see e.g. Arnold et al. 1982; Shandarin & Zeldovich 1989; Hidding et al. 2014; Feldbrugge et al. 2018).

First, we examine the case 0 = ÃxxÃyyÃzz, which corresponds to the following sine waves initial conditions: Q1D-2SIN, Q1D-3SIN, ANI-2SIN, and ANI-3SIN (i.e. 0 ≤ ϵ2D < 1 or 0 ≤ ϵ3D < 1 for i = 1, 2), as summarised in first line of Table 2. This configuration corresponds to the shell-crossing along the x axis.

When computing the radial profiles close to the centre of the system, we consider spherical shells of radius r with r2 = x2 + y2 + z2 ≪ 1. At leading order in q we have, at the shell-crossing, (D.9) (D.10) (D.11)

Applying the scaling xsx implies (D.12)

Then, applying the same scaling in Eq. (D.9) and taking the limit s ≪ 1, we simply obtain (D.13)

After exploiting the symmetries of matrix Ci jkl (Eq. D.6), the Jacobian up at leading order in q is given by (D.14)

Using the scalings in Eqs. (D.12) and (D.13), we have (D.15)

Then, we have Js2/3 when s ≪ 1, resulting in ρ ∝ r−2/3.

Second, we consider the case 0 = Ãxx = ÃyyÃzz, which corresponds to the initial condition SYM-2SIN (i.e. ϵ2D = 1 or ϵ3D,1 = 1 and ϵ3D,2 < 1), as summarised in second line of Table 2. This configuration corresponds to simultaneous shell-crossings along x and y axes. In this case, the main contribution to Eulerian coordinates reads, at the shell-crossing, (D.16) (D.17) (D.18)

where we have exploited the symmetry between x and y, which implies Cxxxx = Cyyyy = C1, Cxxyy = Cyyxx = C2 and Cxxzz = Cyyzz = C3. Applying the same reasoning as above, the scaling xsx implies (D.19)

which leads us, when applying the scaling to Eqs. (D.16) and (D.17) and assuming s ≪ 1, to neglect the last term in these equations, hence (D.20) (D.21)

which leads to the scalings (D.22)

With Ãxx = Ãyy = 0, the first term in Eq. (D.8) vanishes, and the second termbecomes the leading contribution: (D.23)

Exploiting Eq. (D.6) then allows us to write this term as (D.24)

where bi jk are constant coefficients. Using the scalings given by Eqs. (D.19) and (D.22), we have (D.25)

The lowest power of s in Eq. (D.25) is s4/3 with i = j = 1 and k = 0. Thus, we have ρr−4/3.

Third, we consider the case Ãxx = Ãyy = Ãzz = 0, which corresponds in the three sine waves case to the SYM-3SIN (i.e. ϵ3D = (1,1)), as summarised in third line of Table 2. This configuration corresponds to simultaneous shell-crossings along all the axes. The approach is analogous to the previous case. The Eulerian coordinates can be written as follows: (D.26) (D.27) (D.28)

with, again C1 = Cxxxx = Cyyyy = Czzzz, C2 = Ciij j for ij. Equations (D.26)(D.28) directly lead to the scaling of q: (D.29)

Now, the first and second terms in Eq. (D.8) vanish, and we have (D.30)

which reads, after exploitation of the symmetries (D.6), (D.31)

Using the scaling given by Eq. (D.29), we have (D.32)

The lowest power of s in Eq. (D.32) is s2. Thus, we have ρ ∝ r−2.

Finally, we focus on the velocity and pseudo phase-space density profiles. Using Eq. (D.1), the velocity field up to leading order in the Lagrangian coordinate is given by (D.33)

where the dot denotes derivative with respect to the cosmic time t. As can be seen from this relation, the velocity field is proportional to the Lagrangian coordinate, irrespective of initial conditions, because the shell-crossing condition Ãij(tsc) = 0 does not imply Since, according to the calculations performed above, the leading contribution from the Lagrangian vector always come from the scaling qis1/3qi when s ≪ 1, it is fairly easy to convince oneself that the velocity profiles are given, at small radii, by υ2r2/3, , and −υrr1/3, from which we can infer as well the logarithmic slope of the pseudo phase-space density Q(r) through Eq. (41).

References

  1. Alard, C. 2013, MNRAS, 428, 340 [Google Scholar]
  2. Angulo, R.E., Hahn, O., Ludlow, A.D. & Bonoli, S. 2017, MNRAS, 471, 4687 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arnold, V.I., Shandarin, S.F., & Zeldovich, I.B. 1982, Geophys. Astrophys. Fluid Dyn., 20, 111 [NASA ADS] [CrossRef] [Google Scholar]
  4. Baldauf, T., Mercolli, L., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 123007 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bardeen, J.M., Bond, J.R., Kaiser, N., & Szalay, A.S. 1986, ApJ, 304, 15 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, J. Cosmology Astropart. Phys., 2012, 051 [CrossRef] [Google Scholar]
  7. Bernardeau, F. 1994, ApJ, 427, 51 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bernardeau, F., Colombi, S., Gaztanaga, E., & Scoccimarro, R. 2002, Phys. Rep., 3671 [Google Scholar]
  9. Bertschinger, E. 1985, ApJS, 58, 39 [Google Scholar]
  10. Blumenthal, G.R., Faber, S.M., Primack, J.R., & Rees, M.J. 1984, Nature, 311, 517 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bouchet, F.R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bouchet, F.R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575 [NASA ADS] [Google Scholar]
  13. Buchert, T. 1992, MNRAS, 254, 729 [NASA ADS] [CrossRef] [Google Scholar]
  14. Buchert, T. 1994, MNRAS, 267, 811 [NASA ADS] [Google Scholar]
  15. Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375 [Google Scholar]
  16. Buchert, T., Karakatsanis, G., Klaffl, R., & Schiller, P. 1997, A&A, 318, 1 [NASA ADS] [Google Scholar]
  17. Carrasco, J.J.M., Hertzberg, M.P., & Senatore, L. 2012, J. High Energy Phys., 2012, 82 [CrossRef] [Google Scholar]
  18. Carron, J., & Szapudi, I. 2013, MNRAS, 432, 3161 [Google Scholar]
  19. Chernin, A.D. 1993, A&A, 267, 315 [NASA ADS] [Google Scholar]
  20. Colombi, S. 2015, MNRAS, 446, 2902 [Google Scholar]
  21. Colombi, S. 2021, A&A, 647, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. dAmico, G., Gleyzes, J., Kokron, N., et al. 2020, J. Cosmology Astropart. Phys., 2020, 005 [CrossRef] [Google Scholar]
  23. Delos, M.S., Erickcek, A.L., Bailey, A.P., & Alvarez, M.A. 2018, Phys. Rev. D, 97, 041303 [NASA ADS] [CrossRef] [Google Scholar]
  24. Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389 [Google Scholar]
  25. Doroshkevich, A.G. 1973, Astrophys. Lett., 14, 11 [Google Scholar]
  26. Feldbrugge, J., van de Weygaert, R., Hidding, J., & Feldbrugge, J. 2018, J. Cosmology Astropart. Phys., 2018, 027 [CrossRef] [Google Scholar]
  27. Fillmore, J.A., & Goldreich, P. 1984, ApJ, 281, 1 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hahn, O., Angulo, R.E., & Abel, T. 2015, MNRAS, 454, 3920 [CrossRef] [Google Scholar]
  29. Henriksen, R.N. & Widrow, L.M. 1995, MNRAS, 276, 679 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hertzberg, M.P. 2014, Phys. Rev. D, 89, 043521 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hidding, J., Shandarin, S.E., & van de Weygaert, R. 2014, MNRAS, 437, 3442 [CrossRef] [Google Scholar]
  32. Hjorth, J., & Williams, L.L.R. 2010, ApJ, 722, 851 [NASA ADS] [CrossRef] [Google Scholar]
  33. Hockney, R.W., & Eastwood, J.W. 1988, Computer simulation using particles [CrossRef] [Google Scholar]
  34. Ishiyama, T. 2014, ApJ, 788, 27 [Google Scholar]
  35. Ivanov, M.M., Simonovic, M., & Zaldarriaga, M. 2020, J. Cosmology Astropart. Phys., 2020, 042 [CrossRef] [Google Scholar]
  36. Ludlow, A.D., Navarro, J.F., Springel, V., et al. 2010, MNRAS, 406, 137 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lynden-Bell, D. 1967, MNRAS, 136, 101 [Google Scholar]
  38. Matsubara, T. 2015, Phys. Rev. D, 92, 023534 [Google Scholar]
  39. Michaux, M., Hahn, O., Rampf, C., & Angulo, R.E. 2021, MNRAS, 500, 663 [Google Scholar]
  40. Moutarde, F., Alimi, J.M., Bouchet, F.R., Pellat, R., & Ramani, A. 1991, ApJ, 382, 377 [NASA ADS] [CrossRef] [Google Scholar]
  41. Moutarde, F., Alimi, J.M., Bouchet, F.R., & Pellat, R. 1995, ApJ, 441, 10 [NASA ADS] [CrossRef] [Google Scholar]
  42. Nakamura, T. 1985, Prog. Theor. Phys., 73, 657 [CrossRef] [Google Scholar]
  43. Navarro, J.F., Frenk, C.S., & White, S.D.M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
  44. Navarro, J.F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  45. Navarro, J.F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21 [NASA ADS] [CrossRef] [Google Scholar]
  46. Novikov, E.A. 1969, Sov. J. Exp. Theor. Phys., 30, 512 [NASA ADS] [Google Scholar]
  47. Ogiya, G., & Hahn, O. 2018, MNRAS, 473, 4339 [Google Scholar]
  48. Peebles, P.J.E. 1980, The large-scale structure of the universe (Princeton University Press) [Google Scholar]
  49. Peebles, P.J.E. 1982, ApJ, 263, L1 [NASA ADS] [CrossRef] [Google Scholar]
  50. Peebles, P.J.E. 1984, ApJ, 277, 470 [NASA ADS] [CrossRef] [Google Scholar]
  51. Pichon, C., & Bernardeau, F. 1999, A&A, 343, 663 [NASA ADS] [Google Scholar]
  52. Pontzen, A., & Governato, F. 2013, MNRAS, 430, 121 [Google Scholar]
  53. Pueblas, S., & Scoccimarro, R. 2009, Phys. Rev. D, 80, 043504 [NASA ADS] [CrossRef] [Google Scholar]
  54. Rampf, C. 2012, J. Cosmology Astropart. Phys., 2012, 004 [NASA ADS] [CrossRef] [Google Scholar]
  55. Rampf, C. 2019, MNRAS, 484, 5223 [NASA ADS] [CrossRef] [Google Scholar]
  56. Rampf, C., & Frisch, U. 2017, MNRAS, 471, 671 [NASA ADS] [CrossRef] [Google Scholar]
  57. Rampf, C., & Hahn, O. 2021, MNRAS, 501, L71 [NASA ADS] [CrossRef] [Google Scholar]
  58. Rampf, C., Villone, B. & Frisch, U. 2015, MNRAS, 452, 1421 [CrossRef] [Google Scholar]
  59. Rampf, C., Frisch, U., & Hahn, O. 2021, MNRAS, 505, L90 [NASA ADS] [CrossRef] [Google Scholar]
  60. Saga, S., Taruya, A., & Colombi, S. 2018, Phys. Rev. Lett., 121, 241302 [Google Scholar]
  61. Shandarin, S.F., & Zeldovich, Y.B. 1989, Rev. Mod. Phys., 61, 185 [NASA ADS] [CrossRef] [Google Scholar]
  62. Sikivie, P., Tkachev, I.L., & Wang, Y. 1997, Phys. Rev. D, 56, 1863 [NASA ADS] [CrossRef] [Google Scholar]
  63. Sousbie, T., & Colombi, S. 2016, J. Comput. Phys., 321, 644 [Google Scholar]
  64. Taruya, A., & Colombi, S. 2017, MNRAS, 470, 4858 [Google Scholar]
  65. Taylor, J.E., & Navarro, J.F. 2001, ApJ, 563, 483 [NASA ADS] [CrossRef] [Google Scholar]
  66. Uhlemann, C., Rampf, C., Gosenca, M., & Hahn, O. 2019, Phys. Rev. D, 99, 083524 [NASA ADS] [CrossRef] [Google Scholar]
  67. Virtanen, P., Gommers, R., Oliphant, T.E., et al. 2020, Nature Methods, 17, 261 [CrossRef] [Google Scholar]
  68. Zel’dovich, Y.B. 1970, A&A, 5, 84 [Google Scholar]
  69. Zheligovsky, V., & Frisch, U. 2014, J. Fluid Mech., 749, 404 [Google Scholar]
  70. Zukin, P., & Bertschinger, E. 2010a, Phys. Rev. D, 82, 104044 [Google Scholar]
  71. Zukin, P., & Bertschinger, E. 2010b, Phys. Rev. D, 82, 104045 [Google Scholar]

1

The LPT solutions up to the tenth order can be provided upon request as a Mathematica notebook.

2

The designations named here, Q1D-3SIN, ANI-3SIN, and SYM-3SIN, are the same as Q1D-S, ASY-Sb, and SYM-S, used in Saga et al. (2018), respectively.

3

Strictly speaking, our calculations are performed for n < 10, but it is reasonable to expect that the result applies to any arbitrary order.

4

When examining the bottom panels of Figs. 9 (black curves) and 10, we see that solving the equation x(q) = x0 in the multi-stream region does not seem to have the expected number of solutions, that is, 3, 9, or 27 (only the bottom-right panel has an expected number of solutions, i.e., 27). But this is because we are sitting exactly at the centre of the system and symmetries imply superposition of solutions in the subspace y = 0 in 2D and y = z = 0 in 3D. In particular, as discussed in Sect. 5.1, 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, υx) subspace, to 3 and 9 superposed ‘S’ shapes, is visible as the orange curve on bottom panels of Fig. 10. On the other hand, higher-order LPT induces non-linear 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, υx) subspace.

5

We recall, however, that the crossed sine waves configurations we consider are degenerate with respect to 1LPT, as discussed at the end of Sect. 5.1, so the reasoning of Pichon & Bernardeau (1999) requires 2LPT to be applicable in this particular case.

All Tables

Table 1

Parameters of the runs performed with ColDICE (Sousbie & Colombi 2016).

Table 2

Summary of the logarithmic slopes at the shell-crossing obtained by singularity theory applied to LPT predictions at various orders (see Appendix D for details).

All Figures

thumbnail Fig. 1

Convergence behaviour of the shell-crossing time. Top: shell-crossing time calculated at nth-order LPT (dots), together with the fitting function given by Eq. (34) (solid lines), as a function of the perturbation order for various initial conditions, as indicated in the panel. The dashed horizontal lines present the values estimated from the simulations (see Appendix B.1), corresponding to the sixth column of Table 1. Bottom: relative error between the shell-crossing time calculated with nth-order LPT and the one obtained with the fitting formula.

In the text
thumbnail Fig. 2

Phase-space structure for two- and three-sine-wave initial conditions at the collapse time: Q1D-2SIN (top left), ANI-2SIN (centre left), SYM-2SIN (bottom left), Q1D-3SIN (top right), ANI-3SIN (centre right), and SYM-3SIN (bottom right). 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, υx) subspace. Simulation results are compared with standard LPT predictions, which are supplemented with the blue line, denoted by’ΈΧΤ’, which corresponds to the formal extrapolation to the infinite order proposed by STC18 and sketched in Sect. 4.1 for the collapse time. For completeness, the quasi-1D approach (Rampf & Frisch 2017), denoted by Q1D, is also presented (see Appendix C for details).

In the text
thumbnail Fig. 3

Radial profiles and their logarithmic slopes for the two-sine-wave initial conditions at the shell-crossing time. As indicated in top-left panel, LPT predictions are given by solid lines of various colours and Vlasov simulations results are represented by red dots. From left to right, the initial conditions are given by 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 υ2, the radial velocity dispersion , the infall velocity -υr, and the pseudo phase-space density Q(r), respectively. We 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 slope panels, the horizontal dashed lines correspond to the theoretical predictions at the small radii of Appendix D, as listed in Table 2.

In the text
thumbnail Fig. 4

Same as Fig. 3 but for the three-sine-wave initial conditions, from left to right, Q1D-3SIN [ϵ3D = (1/6,1/8)], ANI-3SIN [ϵ3D = (3/4,1/2)], and SYM-3SIN [ϵ3D = (1,1)], respectively. We note that in SYM-3SIN, the closest snapshot from collapse we had at our disposal from our Vlasov runs, âsc = 0.03190, is significantly beyond the actual shell-crossing time estimated by the method described in Sect. B.2, asc = 0.03155, which explains some discrepancies between the theory and the simulation at small radii.

In the text
thumbnail Fig. 5

Same as Fig. 3 but all analytical predictions evaluated at individual shell-crossing times computed at each perturbation order.

In the text
thumbnail Fig. 6

Same as Fig. 4 but all analytical predictions evaluated at individual shell-crossing times computed at each perturbation order.

In the text
thumbnail Fig. 7

Radial profiles at the shell-crossing using tenth-order LPT in the vicinity of (top panels) and (bottom panels). In the top panels, ϵ2D is progressively varying in the range [0.9,1] with linear intervals of Δϵ = 0.1/64. In the bottom panels, ϵ3D progressively changes from (0.97,0.97) to (1,1) with linear intervals of Δϵ = 0.03/64. From left to right, we present radial profiles of the normalised density , the velocity dispersion υ2, the radial velocity dispersion the infall velocity –υr, and the pseudo 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.

In the text
thumbnail Fig. 8

3D view of the expected caustic pattern shortly after the shell-crossing for three-sine-wave initial conditions. The left panels correspond to the most typical pancake, here embodied by ANI-3SIN, but Q1D-3SIN would look similar, while the right ones show the axial-symmetric case, i.e. 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.

In the text
thumbnail Fig. 9

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

In the text
thumbnail Fig. 10

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

In the text
thumbnail Fig. 11

Caustic pattern shortly after the shell-crossing in the 2D case: comparison of LPT predictions using ballistic approximation with Vlasov runs. The left, middle, and right panels correspond, respectively, to Q1D-2SIN, ANI-2SIN, and SYM-2SIN configurations, while the top and bottom panels show the caustics in Lagrangian and Eulerian spaces, respectively.

In the text
thumbnail Fig. 12

Structure of caustics shortly after the shell-crossing for three-sine-wave initial conditions: comparison of LPT predictions using ballistic approximation with Vlasov runs. The top six, middle six, and bottom two panels respectively correspond to Q1D-3SIN, ANI-3SIN, and SYM-3SIN. In each group of six panels, the top and bottom lines correspond to Lagrangian and Eulerian spaces, and the intersection of the caustic surfaces with the plane are qx = 0, qy = 0, and qz = 0, respectively, for the first, second, and third column of each group. In the bottom group of two panels, the left and right panels correspond respectively to Lagrangian and Eulerian space and only show the intersection of the caustic surface network with the plane qz = 0 due to the symmetry of the system.

In the text
thumbnail Fig. 13

2D density shortly after the shell-crossing: comparison of LPT predictions using ballistic approximation with Vlasov runs. From top to bottom: Q1D-2SIN, ANI-2SIN, and SYM 2SIN. From left to right: 2LPT, 10LPT, and measurements in Vlasov simulations.

In the text
thumbnail Fig. 14

Slices of the projected density shortly after the shell-crossing: comparison of LPT using ballistic approximation with Vlasov runs. The left and right groups of nine panels correspond respectively to Q1D-3SIN (top, middle, and bottom line of panels: x=−1.55×l0−4, y = −1.17×10−3, and z= −1.56 × 10−3 slice) and ANI-3SIN (x= −5.16 × 10−4, y = −2.15 × 10−4, and z= −5.47 × 10−4 slice), while the bottom group of three panels corresponds to SYM-3SIN (z= −1.17 × 10−5 slice). In each group of panels, the left, middle, and right columns give respectively the second-order LPT prediction, the tenth-order LPT prediction, and the Vlasov code measurements. Due to the symmetry of the system for SYM-3SIN, only one slice is shown for the bottom panels.

In the text
thumbnail Fig. 15

2D vorticity fields: comparison of LPT predictions with Vlasov runs.

In the text
thumbnail Fig. 16

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

In the text
thumbnail Fig. 17

Same as Fig. 16, but for ANI-3SIN. The 2D slices considered are the same as in the top-right group of nine panels of Fig. 14.

In the text
thumbnail Fig. 18

Same as Fig. 16, but for SYM-3SIN. Due to the symmetry of the initial conditions, only an x-y slice is shown, which is the same as in bottom three panels of Fig. 14.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.