© S. Saga et al. 2022

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 A–D 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, $H_0^{-1}=1$, for the dimensions of length and time, respectively. In particular, the Lagrangian coordinate, Eulerian coordinate, velocity, and vorticity are explicitly expressed as q/L, x/L, υ/(LH₀), and ω/Η₀ = (L∇) × (υ/(LH₀)), respectively.

2 Lagrangian Perturbation Theory

The Lagrangian equation of motion of a fluid element is given by (e.g. Peebles 1980) $$\frac{{\text{d}}^{2}x}{\text{d}{t}^2}+2H\frac{\text{d}x}{\text{d}t}=-\frac{1}{a^2}{\nabla }_x\varphi \left(x\right),$$(1)

where the quantities x, a, and H(t) = a⁻¹da/dt are the Eulerian comoving position, the scale factor of the Universe, and Hubble parameter, respectively. The derivative operator ∇_x = ∂/∂x is a spatial derivative in Eulerian space. The Newton gravitational potential ϕ(x) is related to the matter density contrast $\delta \left(x\right)=\rho \left(x\right)/\overline{\rho }-1$ with $\overline{\rho }$ being the background mass density, through the Poisson equation: $${\nabla }_x^2\varphi \left(x\right)=4\pi G\overline{\rho }{a}^2\delta \left(x\right).$$(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: $${\nabla }_x\cdot \left(\ddot{x}+2H\dot{x}\right)=4\pi G\overline{\rho }\delta \left(x\right),$$(3) $${\nabla }_x\times \left(\ddot{x}+2H\dot{x}\right)=0,$$(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 $$x\left(q,t\right)=q+\Psi \left(q,t\right).$$(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 t_sc), mass conservation implies $\overline{\rho }{\text{d}}^{3}q=\rho \left(x\right){\text{d}}^{3}x$, which leads to $$1+\delta \left(x\right)=\frac{\rho \left(x\right)}{\overline{\rho }}=\frac{1}{J},$$(6)

where the quantity J = det J_ij is the Jacobian of the matrix J_ij defined by $$J_{ij}\left(q,t\right)=\frac{\partial x_i\left(q,t\right)}{\partial q_j}={\delta }_{ij}+{\Psi }_{i,j}\left(q,t\right),$$(7)

and its inverse is given as $$J_{ij}^{-1}=\frac{\partial q_i}{\partial x_j}.$$(8)

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

2.1 Recursion Relation

In LPT, the displacement field, Ψ, is the fundamental building block and is considered to be a small quantity. It can be systematically formally expanded as follows, $$\Psi \left(q,t\right)={\displaystyle \sum _{n=1}^{\infty }{\Psi }^{(n)}}\left(q,t\right).$$(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) $${\Psi }^{(n)}\left(q,t\right)={\Psi }^{(n)}\left(q\right),$$(10)

where the purely time-dependent function D₊(t) corresponds to the linear growth factor. The velocity field is then given by $$\upsilon \left(q,t\right)=aHf{\displaystyle \sum _{n=1}^{\infty }n{D}_{+}^n\left(t\right){\Psi }^{(n)}\left(q\right)}$$(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 $f\approx {\text{Ω}}_{\text{m}}^{3/5}$. With the further approximation $f\approx {\text{Ω}}_m^{1/2}$, 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): $$\left(\widehat{\mathcal{T}}-\frac{3}{2}\right){\Psi }_{k,k}^{\left(n\right)}=-{\epsilon }_{im}{\epsilon }_{jk}{\displaystyle \sum _{n_1+n_2=n}{\Psi }_{m,k}^{\left(n_1\right)}}\left(\widehat{\mathcal{T}}-\frac{3}{4}\right){\Psi }_{i,j}^{\left(n_2\right)},$$(12) $${\epsilon }_{ij}\widehat{\mathcal{T}}{\Psi }_{j,i}^{\left(n\right)}={\epsilon }_{ij}{\displaystyle \sum _{n_1+n_2=n}{\Psi }_{k,j}^{\left(n_1\right)}}\widehat{\mathcal{T}}{\Psi }_{k,i}^{\left(n_2\right)},$$(13)

in the 2D case, and $$\begin{array}{cc}\left(\widehat{\mathcal{T}}-\frac{3}{2}\right){\Psi }_{k,k}^{\left(n\right)}=& -{\epsilon }_{ijk}{\epsilon }_{ipq}{\displaystyle \sum _{n_1+n_2=n}{\Psi }_{j,p}^{\left(n_1\right)}}\left(\widehat{\mathcal{T}}-\frac{3}{4}\right){\Psi }_{k,q}^{\left(n_2\right)}\\ & -\frac{1}{2}{\epsilon }_{ijk}{\epsilon }_{pqr}\\ & \times {\displaystyle \sum _{n_1+n_2+n_3=n}{\Psi }_{i,p}^{\left(n_1\right)}{\Psi }_{j,q}^{\left(n_2\right)}\left(\widehat{\mathcal{T}}-\frac{1}{2}\right){\Psi }_{k,r}^{\left(n_3\right)},}\end{array}$$(14) $${\epsilon }_{ijk}\widehat{\mathcal{T}}{\Psi }_{j,k}^{\left(n\right)}=-{\epsilon }_{ijk}{\displaystyle \sum _{n_1+n_2=n}{\Psi }_{p,j}^{\left(n_1\right)}}\widehat{\mathcal{T}}{\Psi }_{p,k}^{\left(n_2\right)},$$(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/∂q_j, and the tensors ɛ_ij and are, respectively, the 2D and 3D Levi-Civita symbols. The symbol $\widehat{\mathcal{T}}$ stands for a differential operator: $$\widehat{\mathcal{T}}\equiv \frac{{\partial }^2}{\partial \ln D_{+}{}^2}+\frac{1}{2}\frac{\partial }{\partial \ln D_{+}}.$$(16)

In obtaining the recursion relations, we used the following formulae of linear algebra: $$\begin{array}{cc}J=\frac{1}{2}{\epsilon }_{ij}{\epsilon }_{kr}{J}_{ik}{J}_{jr},& J_{ij}^{-1}\end{array}=\frac{1}{J}{\epsilon }_{ik}{\epsilon }_{jr}{J}_{rk},$$(17)

for the 2D case, and $$\begin{array}{cc}J=\frac{1}{6}{\epsilon }_{ijk}{\epsilon }_{pqr}{J}_{ip}{J}_{jq}{J}_{kr},& J_{ij}^{-1}=\frac{1}{2J}{\epsilon }_{jkp}{\epsilon }_{iqr}{J}_{kq}{J}_{pr}\end{array},$$(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 $${\Psi }_{k,k}^{\left(n\right)}=-{\epsilon }_{im}{\epsilon }_{jk}{\displaystyle \sum _{n_1+n_2=n}\frac{a{n}_2^2+2n_2-3}{1\left(n-1\right)\left(2n+3\right)}{\Psi }_{m,k}^{\left(n_1\right)}{\Psi }_{i,j}^{\left(n_2\right)}},$$(19) $${\epsilon }_{ij}{\Psi }_{j.i}^{\left(n\right)}{\epsilon }_{ij}{\displaystyle \sum _{n_1+n_2=n}\frac{n_2\left(2n_2+1\right)}{n\left(2n+1\right)}{\Psi }_{k,j}^{\left(n_1\right)}{\Psi }_{k,i}^{\left(n_2\right)}},$$(20) $$\begin{array}{ccc}{\Psi }_{k,k}^{\left(n\right)}& =& -{\epsilon }_{ijk}{\epsilon }_{ipq}{\displaystyle \sum _{n_1+n_2=n}\frac{4n_2^2+2n_2-3}{2\left(n-1\right)\left(2n+3\right)}{\Psi }_{j,p}^{\left(n_1\right)}{\Psi }_{k,q}^{\left(n_2\right)}}\\ & & -{\epsilon }_{ijk}{\epsilon }_{pqr}\\ & & \times {\displaystyle \sum _{n_1+n_2+n_3=n}\frac{\left(n_3+1\right)}{2\left(n-1\right)\left(2n+3\right)}{\Psi }_{i,p}^{\left(n_1\right)}{\Psi }_{j,q}^{\left(n_2\right)}{\Psi }_{k,r}^{\left(n_3\right)},}\end{array}$$(21) $${\epsilon }_{ijk}{\Psi }_{jk}^{\left(n\right)}=-{\epsilon }_{ijk}{\displaystyle \sum _{n_1+n_2=n}\frac{n_2\left(2n_2+1\right)}{n\left(2n+1\right)}{\Psi }_{p,j}^{\left(n_1\right)}{\Psi }_{p,k}^{\left(n_2\right)}},$$(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[: $${\Psi }_i^{\text{ini}}\left(q,t_{\text{ini}}\right)=\frac{L}{2\pi }{D}_{+}\left(t_{\text{ini}}\right){ϵ}_i\sin \left(\frac{2\pi }{L}qi\right)$$(23)

The parameters ϵ_i < 0 with |ϵ_x| ≥ |ϵ_y|≥ |ϵ_z| quantify the linear amplitudes of the sine waves in each direction. The initial time, t_ini, is chosen such that D₊(t_ini)|ϵ_i| ≤ 0.012 ≪ 1, which makes the fastest growing mode approximation very accurate, as shown by STC18. In this framework, the dependence on ϵ_i of the dynamics is reduced to a function of the ratios ϵ_2D = ϵ_y/ϵ_x and ϵ_3D = (ϵ_y/ϵ_x,ϵ_Z/ϵ_x), respectively, for two- and 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 δⁱⁿⁱ ≃ –∇_q. Ψⁱⁿⁱ = D₊(t_ini)Σ_i |ϵ_i| cos (2π/Lq_i), presents a small peak at the origin, and mass elements subsequently infall towards the central overdense region. With the proper choice of ϵ_i, this initial overdensity can asymptotically match any peak of a smooth random Gaussian field (see e.g. Bardeen et al. 1986), which actually makes the 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: $${\nabla }_q\cdot {\Psi }_{\text{L}}^{\left(n\right)}={\displaystyle \sum _m{\alpha }_m^{\left(n\right)}}{\text{e}}^{\text{i\hspace{0.17em}}m\cdot q},$$(24) $${\nabla }_q\times {\Psi }_{\text{T}}^{\left(n\right)}={\displaystyle \sum _m{\beta }_m^{\left(n\right)}}{\text{e}}^{\text{i\hspace{0.17em}}m\cdot q},$$(25)

where the nth-order scalar coefficients ${\alpha }_m^{\left(n\right)}$ and the nth-order vector coefficients ${\beta }_m^{\left(n\right)}$ 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 ${\nabla }_q\times {\Psi }_{\text{L}}^{\left(n\right)}=0$ and ${\nabla }_q\cdot {\Psi }_{\text{T}}^{\left(n\right)}=0$, one can build up the perturbative solutions ${\Psi }_{\text{L}}^{\left(n\right)}$ and ${\Psi }_{\text{T}}^{\left(n\right)}$ using simple algebraic manipulations involving coefficients ${\alpha }_m^{\left(n\right)}$ and ${\beta }_m^{\left(n\right)}$: $${\Psi }_{\text{T}}^{\left(n\right)}={\displaystyle \sum _m\left(-\text{i}\right){\alpha }_m^{\left(n\right)}}\frac{m}{{\left|m\right|}^2}{\text{e}}^{\text{i}m\cdot q},$$(26) $${\Psi }_{\text{T}}^{\left(n\right)}={\displaystyle \sum _m\text{i}}\frac{m\times {\beta }_m^{\left(n\right)}}{{\left|m\right|}^2}{\text{e}}^{\text{i\hspace{0.17em}}m\cdot q}.$$(27)

These solutions lead to the nth-order displacement field given by ${\Psi }^{\left(n\right)}={\Psi }_{\text{L}}^{\left(n\right)}+{\Psi }_{\text{T}}^{\left(n\right)}$. Using this solution for Ψ⁽ⁿ⁾ as well as Eqs. (21) and (22), we subsequently construct the source of the (n + 1)th-order derivatives ${\nabla }_q\cdot {\Psi }_{\text{L}}^{\left(n+1\right)}$ and ${\nabla }_q\times {\Psi }_{\text{T}}^{\left(n+1\right)}$. 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 way¹.

The above prescription is valid in 3D space and can be applied to the 2D case by cancelling all fluctuations along the z axis, that is, by performing 3D calculations with ϵ_z = 0 (i.e. ϵ_3D = (ϵ_y/ϵ_x, 0)) for the three sine waves case. However, one can also realise that in 2D, vector coefficients ${\beta }_m^{\left(n\right)}$ become scalars ${\beta }_m^{\left(n\right)}$, and that the solutions take the following form: $${\Psi }_{\text{L}}^{\left(n\right)}={\displaystyle \sum _m\left(-\text{i}\right)}\frac{{\alpha }_m^{\left(n\right)}}{{\left|m\right|}^2}\left(\begin{array}{c}{m}_x\\ m_y\end{array}\right){\text{e}}^{\text{i\hspace{0.17em}}m\cdot q},$$(28) $${\Psi }_{\text{T}}^{\left(n\right)}={\displaystyle \sum _m\text{i}}\frac{{\beta }_m^{\left(n\right)}}{{\left|m\right|}^2}\left(\begin{array}{c}{m}_y\\ -m_x\end{array}\right){\text{e}}^{\text{i\hspace{0.17em}}m\cdot q}.$$(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|( = |e_z|), denoted by Q1D, ANI, and SYM, respectively².

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).

3.1 The Vlasov Code ColDICE

The Lagrangian coordinates defined in Sect. 2.1 correspond to the following unperturbed initial set-up, $$X\left(Q,t_0\right)\equiv Q,$$(30) $$V\left(Q,t_0\right)\equiv 0,$$(31)

where t₀ 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: $$X\left(Q,t_{\text{ini}}\right)=Q+{\Psi }^{\text{ini}}\left(Q,t_{\text{ini}}\right),$$(32) $$V\left(Q,t_{\text{ini}}\right)=aHf{\Psi }^{\text{ini}}\left(Q,t_{\text{ini}}\right),$$(33)

with Ψⁱⁿⁱ 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 n_g = n_s = 512 (see below for the definitions of n_g and n_s).

The Poisson equation is solved using the fast-Fourier technique in a mesh of fixed resolution n_g. To estimate the projected density ρ(x) on this mesh, the intersection of each simplex of the phase-space sheet with each voxel 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 $n_{\text{s}}^D$ vertices corresponding, respectively, to $2n_{\text{s}}^2$ and $6n_{\text{s}}^3$ 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⁻⁸ 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 n_g of the mesh used to solve Poisson equation and the initial number of vertices of the tessellation, n_s³. In Appendix B, we explain how measurements of various quantities are performed, such as the set of curves corresponding to the intersection of the phase-space sheet with the hyperplane y( = z) = 0 used in Sects. 4.2 and 5.2, the collapse time t_sc shown in Table 1, the radial profiles used in Sect. 4.3, as well as the caustic network, the projected density and the vorticity analysed in Sect. 5.3.

3.2 Accuracy and Related Considerations

Accuracy and possible defects of ColDICE are discussed in Sousbie & Colombi (2016) and in C21 in the 3D case. In our work, precise determination of the collapse time is critical, since we consider, in the analyses performed below, either the collapse time or a moment very shortly after it. The uncertainty on the collapse time depends on the nature of the initial conditions and the parameters that control the accuracy of the simulations, in particular the number of simplices used to represent the phasespace sheet and the spatial resolution n_g of the mesh used to solve Poisson equation. This latter turns out to be a very important parameter. Indeed, decreasing n_g augments the collapse time, as discussed in C21. In Sect. B.2, relying on measurements in ANI-3SIN configurations with various spatial resolutions, we estimate that the quantity a_sc shown in Table 1 should be accurate at the fourth significant digit level or better for all configurations, for example a_sc = 0.02919 ± 10⁻⁴ for ANI-3SIN, except for SYM-3SIN, where the uncertainty on the measured a_sc could be a few 10⁻⁴. 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 n_s, as well as the refinement parameter I limiting deviations from symplectic motion. As discussed in Sect. 3.1, for our simulations, parameter I does not influence the results except (probably marginally) for SYM-3SIN. Effects related to the choice of n_s are discussed in Sousbie & Colombi (2016) and in C21. They should be negligible for our sine wave simulations compared with other sources of uncertainty on the collapse time for the values of n_s we adopted. However, we did not test the accuracy of ColDICE in more complex cases, because we do not need to.

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

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

Turning to the fields, the vorticity, because it is composed of derivatives, is particularly challenging to measure even with the considerable gain brought by the tessellation technique compared with more traditional representations based on particles. Appendix B.5 explains in details how we measure density and vorticity, making use, this time, of the locally quadratic representation of the 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, t_sc, 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 t_sc is of prime importance. We discuss the extrapolation to the infinite order of the LPT series proposed by STC18 and generalise it to the 2D case. Second, Sect. 4.2 examines the convergence of LPT at collapse with 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 $a_{\text{sc}}^{\left(n\right)}$ of the collapse time obtained from rath-order LPT and its extrapolation to the infinite order $a_{\text{sc}}^{\infty }$, as described in STC18, with the value a_sc measured in the Vlasov runs as explained in Appendix B.2. Our approach follows closely that of STC18. It is basically intended to repeat its main steps and to supplement it with additional discussions and comparisons with Vlasov runs in the 2D case.

Using rath-order LPT predictions, we explore the sequence of the shell-crossing times, $a_{\text{sc}}^{\left(n\right)}$, as a function of order n up to n = 10, by solving J⁽ⁿ⁾ = 0 at the origin, where J⁽ⁿ⁾ 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, $a_{\text{sc}}^{\left(n\right)}$ 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):

$$a_{\text{sc}}^{\left(n\right)}-a_{\text{sc}}^{\infty }+{\left(b+c\exp \left[d{n}^e\right]\right)}^{-1},$$(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 $a_{\text{sc}}^{\left(n\right)}$ of each respective order. Equation (34) also implies $a_{\text{sc}}^{\left(n\right)}-a_{\text{sc}}^{\left(n-1\right)}~\text{exp}\left(-\left|d\right|n^e\right)$ 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 $a_{\text{sc}}^{\left(n\right)}-a_{\text{sc}}^{\left(n-1\right)}~n^{\kappa }\text{\hspace{0.17em}exp}\left(-\eta n\right)\left(\text{i.e.\hspace{0.17em}\hspace{0.17em}}e=1\right)$, 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 $a_{\text{sc}}^{\left(n\right)}$ 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 $a_{\text{sc}}^{\left(n\right)}$ or equivalently $D_{+}^{\left(n\right)}$. 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 $a_{\text{sc}}^{\infty }$ obtained with Eq. (34) are very accurate, as illustrated by Table 1. Indeed, the relative difference between $a_{\text{sc}}^{\infty }$ and the value a_sc measured in the Vlasov simulations remains of the order of the percent or below, which suggests that the error on the extrapolated estimate of the collapse time remains small and should not affect significantly the analytical predictions of the phase-space structure and the radial profiles at the shell-crossing. Therefore, we use $a_{\text{sc}}^{\infty }$ for the analytical predictions of the ‘exact’ collapse time in the rest of Sect. 4. As for the simulations data, we use the snapshot with the closest possible available value â_sc of the expansion factor to a_sc, as indicated in Table 1.

Fig. 1 Convergence behaviour of the 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.

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 $\rho /\overline{\rho }$, velocity dispersion υ², radial velocity dispersion $v_r^2$, and infall velocity -υ_r are given by $$\rho \left(r,t\right)/\overline{\rho }={\langle J^{-1}\left(q,t\right)\rangle }_{\Omega },$$(35) $${\upsilon }^2\left(r,t\right)=\frac{{\langle J^{-1}\left(q,t\right){\upsilon }^2\left(q,t\right)\rangle }_{\Omega }}{{\langle J^{-1}\left(q,t\right)\rangle }_{\Omega }},$$(36) $${\upsilon }_r^2\left(r,t\right)=\frac{{\langle J^{-1}\left(q,t\right){\left(\upsilon \left(q,t\right)\cdot \widehat{x}\right)}^2\rangle }_{\Omega }}{{\langle J^{-1}\left(q,t\right)\rangle }_{\Omega }},$$(37) $${\upsilon }_v\left(r,t\right)=\frac{{\langle J^{-1}\left(q,t\right)\left(\upsilon \left(q,t\right)\cdot \widehat{x}\right)\rangle }_{\Omega }}{{\langle J^{-1}\left(q,t\right)\rangle }_{\Omega }},$$(38)

where $\widehat{x}=x/r$ with r = |x| being the radial coordinate. In these equations, we used the angle average defined by $${\langle f\left(q\right)\rangle }_{\Omega }={\displaystyle \int \frac{\text{d}\theta }{4\pi }f\left(q\right){|}_{x=x\left(q,t\right)},}$$(39)

for two-sine-wave initial conditions, and $${\langle f\left(q\right)\rangle }_{\Omega }={\displaystyle \int \frac{\text{sin}\theta \text{d}\theta \text{d}\varphi }{4\pi }f\left(q\right){|}_{x=x\left(q,t\right)},}$$(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 $$Q\left(r,t\right)=\frac{\rho \left(r,t\right)/\overline{\rho }}{{\upsilon }^3\left(r,t\right)}.$$(41)

We additionally note that we present not only the radial profiles but also their logarithmic slopes, defined by $$n\equiv \frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}X}{\text{d\hspace{0.17em}ln\hspace{0.17em}}r},$$(42)

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

Radial profiles at the shell-crossing are presented in Figs. 3–6. In these figures, measurements in the Vlasov code are performed by replacing each simplex of the phase-space sheet with a large ensemble of particles as explained in detail in C21 (see also Appendix B.4). Two approaches of collapse are considered. In Figs. 3 and 4, which 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 $a_{\text{sc}}^{\infty }$ 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 $a_{\text{sc}}^{\left(n\right)}$, 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 r_min(n) decreasing with increasing n, where n is the perturbation order. LPT can describe arbitrarily large densities, as long as n is large enough. In the results presented here, density contrasts as large as 100 or larger can be probed accurately by LPT, depending on the order considered and on the nature of initial conditions. Again, convergence worsens when approaching 3D 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 r_min(n) for the velocities than for the density, with differences as large as an order of magnitude. With synchronisation, performances of LPT predictions greatly improve, as expected and as illustrated by Figs. 5 and 6. All perturbation orders predict strikingly similar density profiles, in close to perfect match with the simulations, except for 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. 3–6. The nature of the singularities of CDM structures at the collapse time and beyond has been widely studied in the literature, particularly in the framework of Zel’dovich motion (e.g. Zel’dovich 1970; Arnold et al. 1982; Hidding et al. 2014; Feldbrugge et al. 2018). In Appendix D, we 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 order³, 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⁻², for SYM-2SIN and SYM-3SIN, respectively. On the other hand, velocities are expected to follow the same power-law pattern whatever initial conditions or dimensionality, with logarithmic slopes equal to 2/3, 2/3, and 1/3 respectively, for υ², $v_r^2$, and -υ_r, which in turn implies Q(r) ∝ r^−7/3 for SYM-2SIN, r⁻³ 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 $v_r^2$ and not at the individual collapse time $a_{\text{sc}}^{\left(n\right)}$ 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 > a_sc, one notices that all the simulated data show a close to perfect agreement between the measured slope at small radii and the predicted ones for all the radial profiles. This result is non-trivial in the sense that the Taylor expansion at small |q| underlying singularity theory is not necessarily valid in the actual fully 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 ${ϵ}_{\text{2D}}={ϵ}_{\text{2D}}^{\text{SYM}}\equiv 1$ or ${ϵ}_{\text{3D}}={ϵ}_{\text{3D}}^{\text{SYM}}\equiv \left(1,1\right)$ 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 ${ϵ}_{\text{2D}}^{\text{SYM}}$ (top panels) and ${ϵ}_{\text{3D}}^{\text{SYM}}$ (bottom panels) for tenth-order LPT. Clearly, at sufficiently large radii, values of ϵ_2D or ϵ_3D close to ${ϵ}_{\text{2D}}^{\text{SYM}}$ and ${ϵ}_{\text{3D}}^{\text{SYM}}$ 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 ρⁱⁿⁱ ∝ 1 − αr² with α > 0 (see e.g. Nakamura 1985; Moutarde et al. 1995). Hence, the anisotropic nature of the axial-symmetric case (contained in terms beyond leading order r² 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.

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 $\rho /\overline{\rho }$, the velocity dispersion υ2, the radial velocity dispersion $v_r^2$, 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.

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.

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

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

Fig. 7 Radial profiles at the shell-crossing using tenth-order LPT in the vicinity of ${ϵ}_{\text{2D}}\simeq {ϵ}_{\text{2D}}^{\text{SYM}}\equiv 1$ (top panels) and ${ϵ}_{\text{3D}}\simeq {ϵ}_{\text{3D}}^{\text{SYM}}\equiv \left(1,1\right)$ (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 $\rho /\overline{\rho }$, the velocity dispersion υ2, the radial velocity dispersion $v_r^2$ 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 $a_{\text{sc}}^{\left(n\right)}$ 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 $a_{\text{pc}}^{\left(n\right)}=a_{\text{sc}}^{\left(n\right)}+\text{Δ}a$ for LPT of nth order and $a_{\text{pc}}^{\text{SIM}}=a_{\text{sc}}+\text{Δ}a$ for the simulation, as listed in last column of Table 1. The value of Δa used for each sine waves configuration is given by the difference between the values of the last and sixth columns of the table. We recall that the quantity a_sc corresponds to the ‘true’ collapse time measured in the Vlasov runs as explained in Appendix B.2.

From the Eulerian coordinate and velocity fields given by nth-order LPT solutions at the shell-crossing time $a_{\text{sc}}^{\left(n\right)}$ of each perturbation order, we model the Eulerian coordinate after collapse as follows: $${x\left(q,a\right)=x^{\left(n\right)}\left(q,a_{\text{sc}}^{\left(n\right)}\right)+\frac{\partial x^{\left(n\right)}}{\partial a}|}_{a=a_{\text{sc}}^{\left(n\right)}}\Delta a,$$(43) $${=x^{\left(n\right)}\left(q,a_{\text{sc}}^{\left(n\right)}\right)+\frac{{\upsilon }^{\left(n\right)}}{a^{2}H\left(a\right)}|}_{a=a_{\text{sc}}^{\left(n\right)}}\Delta a,$$(44)

while the nth-order velocity field is frozen to its value at the shell-crossing time of each perturbation order, $v\left(q,a\right)=v^{\left(n\right)}\left(q,a_{\text{sc}}^{\left(n\right)}\right)$

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 q_F labelled by the subscript F = [1,…, n_F]. Shortly after collapse, if ϵ_i ≠ ϵ_j for i ≠ j, there are at most three streams, n_F ≤ 3, in a given point of space. For symmetric cases, for example ϵ_2D = 1 in 2D or ϵ_x = ϵ_y in 3D, due to simultaneous 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, $${\rho }_{\text{L}}\left(q\right)={\left|\det J\left(q\right)\right|}^{-1},$$(45) $${\upsilon }_{\text{L}}\left(q\right)=a\frac{\text{d}\Psi \left(q\right)}{\text{d}t},$$(46)

the Eulerian density and velocity fields are expressed by the superpositions of each flow, $$\rho \left(x\right)={\displaystyle \sum _F{\rho }_{\text{L}}}\left(q_F\right),$$(47) $$\upsilon \left(x\right)=\frac{{\displaystyle \sum {}_F{\rho }_{\text{L}}}\left(q_F\right){\upsilon }_{\text{L}}\left(q_F\right)}{{\displaystyle \sum {}_F{\rho }_{\text{L}}}\left(q_F\right)},$$(48)

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

While vorticity is not expected in 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 $$\begin{array}{ccc}{\omega }^{2\text{D}}\left(x\right)& =& {\epsilon }_{ij}\frac{\partial {\upsilon }_j\left(x\right)}{\partial x_i}\\ & =& {\epsilon }_{ij}\frac{1}{\rho \left(x\right)}{\displaystyle \sum _F{J}_{mi}^{-1}\left(q_F\right)\frac{\partial {\rho }_{\text{L}}\left(q_F\right)}{\partial q_{F,m}}}\left[{\upsilon }_{\text{L,}j}\left(q_F\right)-{\upsilon }_i\left(x\right)\right]\\ & =& {\epsilon }_{ij}\frac{1}{\rho \left(x\right)}{\displaystyle \sum _F{J}_{mi}^{-1}\left(q_F\right)\frac{\partial {\upsilon }_{\text{L,}j}\left(q_F\right)}{\partial q_{F,m}}}{\rho }_{\text{L}}\left(q_F\right),\end{array}$$(49) $$\begin{array}{ccc}{\omega }_i^{{}^{2\text{D}}}\left(x\right)& =& {\epsilon }_{ijk}\frac{\partial {\upsilon }_k\left(x\right)}{\partial x_j},\\ & =& {\epsilon }_{ijk}\frac{1}{\rho \left(x\right)}{\displaystyle \sum _F{J}_{mj}^{-1}\left(q_F\right)\frac{\partial {\rho }_{\text{L}}\left(q_F\right)}{\partial q_{F,m}}}\left[{\upsilon }_{\text{L,}k}\left(q_F\right)-{\upsilon }_k\left(x\right)\right]\\ & =& {\epsilon }_{ijk}\frac{1}{\rho \left(x\right)}{\displaystyle \sum _F{J}_{mj}^{-1}\left(q_F\right)\frac{\partial {\upsilon }_{\text{L,}k}\left(q_F\right)}{\partial q_{F,m}}}{\rho }_{\text{L}}\left(q_F\right).\end{array}$$(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: $$A^{\left(1\text{LPT}\right)}\left(q,a\right)=q_A+\frac{L}{2\pi }{D}_{+,\text{sc}}^{\left(1\text{LPT}\right)}{ϵ}_A\left(1+f\frac{\Delta a}{a_{sc}}\right)\sin \left(\frac{2\pi }{L}{q}_A\right),$$(51)

with A = x, y, or z and $D_{\text{+,sc}}^{\left(\text{1LPT}\right)}=1/{ϵ}_x$ 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 a_sc are evaluated at the same time as $D_{\text{+,sc}}^{\left(\text{1LPT}\right)}$ (we recall that f = 1 in the Einstein-de Sitter cosmology considered in this work). The value of $D_{\text{+,sc}}^{\left(\text{1LPT}\right)}$ is obtained in practice by solving the equation $\partial x^{\left(\text{1LPT}\right)}/\partial q_x=0$ 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 $$\cos \left(\frac{2\pi }{L}{q}_i\right)={\left(1+f\frac{\Delta a}{a_{\text{sc}}}\right)}^{-1}.$$(52)

Equation (52) implies that the equations of the caustics are given by q = q₀, with q₀ 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 ,i ≠ j, to the equation of an ellipse in 2D and an ellipsoid in 3D (see e.g. Hidding et al. 2014; Feldbrugge et al. 2018). The reason for this not being the case here is due to the extremely restrictive class of initial conditions we have chosen.

Due to the degenerate nature of the sine waves initial conditions, it is therefore needed to go beyond the first order of the perturbative development of the dynamical equations to obtain a more realistic shape for the caustics. Indeed, for instance, 2LPT brings non-linear couplings between various axes of the dynamics: $$\begin{array}{ccc}{A}^{\left(2\text{LPT}\right)}\left(q,a\right)& =& q_A+\frac{L}{2\pi }{D}_{+,\text{sc}}^{\left(2\text{LPT}\right)}{ϵ}_A\left(1+f\frac{\Delta a}{a_{\text{sc}}}\right)sin\left(\frac{2\pi }{L}{q}_A\right)\\ & & -\frac{3}{14}\frac{L}{2\pi }{\left(D_{+,\text{sc}}^{\left(2\text{LPT}\right)}\right)}^2{\displaystyle \sum _B\left(1-{\delta }_{AB}\right){ϵ}_A{ϵ}_B}\\ & & \times \left(1+2f\frac{\Delta a}{a_{sc}}\right)\sin \left(\frac{2\pi }{L}{q}_A\right)\cos \left(\frac{2\pi }{L}{q}_{\text{B}}\right),\end{array}$$(53)

where $D_{\text{+,sc}}^{\left(\text{2LPT}\right)}$ 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: $$1+D_{+,\text{sc}}^{\left(2\text{LPT}\right)}{ϵ}_x-\frac{3}{14}{\left(D_{+,\text{sc}}^{\left(2\text{LPT}\right)}\right)}^2{ϵ}_x\left({ϵ}_y+{ϵ}_z\right)=0.$$(54)

In Eq. (53), the growth rate f and the scale factor a_sc are evaluated at the same time as $D_{\text{+,sc}}^{\left(\text{2LPT}\right)}$. 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.

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.

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 $a={\widehat{a}}_{\text{sc}}$ 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 three⁴.

The results shown in Fig. 9 must however be interpreted with caution. Strictly speaking, we aim to apply the ballistic approximation from the collapse time, which is not exactly the case in this figure. The red curves, which correspond to the snapshots used to originate the ballistic motion, all correspond to a time a_bc slightly before the actual collapse (namely a_bc = 0.03103 for SYM-3SIN and $a_{\text{bc}}={\widehat{a}}_{\text{sc}}$, 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 Δa_eff /a_bc ≃ 0.03 is the largest, where Δa_eff is the difference between a_pc ≡ a_sc + Δa and a_bc and is used to test the ballistic approximation. We also note that Δa_eff /a_bc ≃ 0.03 is also of the same order for ANI-3SIN. The greater the value of Δa_eff /a_bc, 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.

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.

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 16–18) are shown in slices slightly away from the origin because the vorticity vanishes in the x = 0, y = 0, and z = 0 planes due to the symmetries in the chosen setting.

The examination of left four panels of Fig. 11 and top twelve panels of Fig. 12 provides us insights on the caustic pattern in the most typical configurations, ϵ_i ≠ ϵ_j i ≠ j, that is the non-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, i ≠ j, agrees perfectly with the predictions of Pichon & Bernardeau (1999) based on Zel’dovich dynamics⁵. 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 x−z and x−y subspace, respectively.

Turning to the axial-symmetric case, the agreement between theory and measurements is again only partial. Yet the high magnitude of vorticity is of the same order for theory and measurements. The structure of the vorticity field, significantly more complex than in the generic case due to multiple foldings of the phase-space sheet, is qualitatively in agreement between LPT and Vlasov code in 2D, while simulations measurements are too noisy in the 3D case to make definitive conclusions. What is clear, however, is that the outer part of the multi-stream region seems to be qualitatively reproduced by the theory, but its size is totally wrong.

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.

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.

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⁻² for SYM-3SIN, while velocity profiles always present the same power-law behaviour, for example υ²(r) ∝ r^2/3. Confirming predictions of singularity theory at collapse and explicit convergence to asymptotic profiles at small radii was expected but is not trivial;

• Synchronisation: Agreement with singularity theory predictions is made even more explicit by computing LPT profiles at the respective collapse times computed at each perturbation order. In this case, convergence to the 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 = a − a_sc 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.

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

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).

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.

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: $$\text{Ψ}\left(q,t\right)={\displaystyle \sum _{n=1}^{\infty }{D}_{+}^n\left(t\right){\text{Ψ}}^{\left(n\right)}}\left(q\right)$$(A.1)

For presentation purposes, we only show the x-components of the displacement field, that is, ${\Psi }_x^{\left(n\right)}\left(q\right)$ with the definition (ϵ₁, ϵ₂) ≡ (ϵ_y/ϵ_x, ϵ_z/ϵ_x). Given the x-components, the y- and z-components can be derived by permutating all x, y, and z: $${\text{Ψ}}_x^{\left(1\right)}=\frac{{ϵ}_x}{2\pi }\sin \left(2\pi q_x\right),$$(A.2) $${\text{Ψ}}_x^{\left(2\right)}=\frac{3{ϵ}_x^2}{2\pi }\left[{ϵ}_1\cos \left(2\pi q_y\right)+{ϵ}_2\cos \left(2\pi q_z\right)\right]sin\left(2\pi q_x\right),$$(A.3) $$\begin{array}{l}{\text{Ψ}}_x^{\left(3\right)}=\frac{{ϵ}_x^3}{2520\pi }[78\cos \left(2\pi q_x\right)+\left({ϵ}_1\cos \left(2\pi q_y\right)+{ϵ}_2\cos \left(2\pi q_z\right)\right),\hfill \\ +160{ϵ}_1{ϵ}_2\cos \left(2\pi q_y\right)\cos \left(2\pi q_z\right)-3{ϵ}_1^2\cos \left(4\pi q_y\right)-3{ϵ}_2^2\cos \left(4\pi q_z\right)+75\left({ϵ}_1^2+{ϵ}_2^2\right)]\sin \left(2\pi q_x\right),\hfill \end{array}$$(A.4) $$\begin{array}{l}{\text{Ψ}}_x^{\left(4\right)}-\frac{{ϵ}_x^4}{776100\pi }[{ϵ}_2\cos \left(2\pi q_z\right)\left(4242\cos \left(4\pi q_x\right)+28550{ϵ}_1^1\cos \left(4\pi q_y\right)+208850{ϵ}_1^2+57015{ϵ}_2^2+89166\right)\hfill \\ +60\cos \left(4\pi q_x\right)\left(6010{ϵ}_1{ϵ}_2\cos \left(2\pi q_y\right)\cos \left(2\pi q_z\right)+1274{ϵ}_1^2\cos \left(4\pi q_y\right)+1274{ϵ}_2^2\cos \left(4\pi q_z\right)+2039\left({ϵ}_1^2+{ϵ}_2^2\right)\right)\hfill \\ +2{ϵ}_1\cos \left(4\pi q_y\right)\left(2121\cos \left(4\pi q_x\right)-9303{ϵ}_1^2\cos \left(4\pi q_y\right)+1475{ϵ}_2^2\cos \left(4\pi q_z\right)+33159{ϵ}_1^2+104425{ϵ}_2^2+44583\right)\hfill \\ -9030{ϵ}_2^3\cos \left(6\pi q_z\right)]\sin \left(2\pi q_x\right),\hfill \end{array}$$(A.5) $$\begin{array}{l}{\text{Ψ}}_x^{\left(5\right)}=\frac{{ϵ}_x^5}{36793476720000\pi }[895050\cos \left(2\pi q_x\right)\left(181296{ϵ}_1{ϵ}_2\cos \left(2\pi q_y\right)\cos \left(2\pi q_z\right)+41657\left({ϵ}_1^2+{ϵ}_2^2\right)\right)\hfill \\ \begin{array}{l}+560226137900{ϵ}_1^2{ϵ}_2\cos \left(2\pi q_x\right)\cos \left(4\pi q_y\right)\cos \left(2\pi q_z\right)+560226137900{ϵ}_1{ϵ}_2^2\cos \left(2\pi q_x\right)\cos \left(4\pi q_y\right)\cos \left(4\pi q_z\right)\\ +594423768510{ϵ}_1^3\cos \left(2\pi q_x\right)\cos \left(2\pi q_y\right)-18642271230{ϵ}_1^3\cos \left(2\pi q_x\right)\cos \left(6\pi q_y\right)\end{array}\hfill \\ +57401651835{ϵ}_1^2\cos \left(4\pi \left(q_x-q_y\right)\right)+57401651835{ϵ}_1^2\cos \left(4\pi \left(q_x+q_y\right)\right)+100617966020{ϵ}_1{ϵ}_2^2\cos \left(2\pi q_x\right)\cos \left(2\pi q_y\right)\hfill \\ +106526974554{ϵ}_1\cos \left(2\pi q_x\right)\cos \left(2\pi q_y\right)-6828912090{ϵ}_1\cos \left(6\pi q_x\right)\cos \left(2\pi q_y\right)+1006179766020{ϵ}_1^3{ϵ}_2\cos \left(2\pi q_x\right)\cos \left(2\pi q_z\right)\hfill \\ +57401651835{ϵ}_2^2\cos \left(4\pi \left(q_x-q_z\right)\right)+594423768510{ϵ}_2^3\cos \left(2\pi q_x\right)\cos \left(2\pi q_z\right)-18642271230{ϵ}_2^{3}cos\left(2\pi q_x\right)cos\left(6\pi q_z\right)\hfill \\ +57401651835{ϵ}_2^2\cos \left(4\pi \left(q_x+q_z\right)\right)+106526974554{ϵ}_2\cos \left(2\pi q_x\right)\cos \left(2\pi q_z\right)-6828912090{ϵ}_2\cos \left(6\pi q_x\right)\cos \left(2\pi q_z\right)\hfill \\ +46659033850{ϵ}_1^2{ϵ}_2^2\cos \left(4\pi \left(q_y-q_z\right)\right)+46659033850{ϵ}_1^2{ϵ}_2^2\cos \left(4\pi \left(q_y+q_z\right)\right)503168866200{ϵ}_1^3{ϵ}_2\cos \left(4\pi q_y\right)\cos \left(4\pi q_z\right)\hfill \\ -54358176600{ϵ}_1^3{ϵ}_2\cos \left(6\pi q_y\right)\cos \left(2\pi q_z\right)+503168866200{ϵ}_1{ϵ}_2^3\cos \left(2\pi q_y\right)\cos \left(2\pi q_z\right)-54358176600{ϵ}_1{ϵ}_2^3\cos \left(2\pi q_y\right)\cos \left(6\pi q_z\right)\hfill \\ +877320054000{ϵ}_1{ϵ}_2\cos \left(2\pi q_y\right)\cos \left(2\pi q_z\right)+115709716980{ϵ}_1^2{ϵ}_2^2\cos \left(2\pi q_y\right)-35064073044{ϵ}_1^4\cos \left(4\pi q_y\right)\hfill \\ -1607701095{ϵ}_1^4\cos \left(8\pi q_y\right)+241408089450{ϵ}_1^2\cos \left(4\pi q_y\right)+115709716980{ϵ}_1^2{ϵ}_2^2\cos \left(4\pi q_z\right)-35064073044{ϵ}_2^4\cos \left(4\pi q_z\right)\hfill \\ +2414082809450{ϵ}_2^2\cos \left(4\pi q_z\right)-1607701095{ϵ}_2^4\cos \left(8\pi q_z\right)\hfill \\ +546975\left({ϵ}_1^2\left(940700{ϵ}_2^2+443058\right)+171045{ϵ}_1^4+21{ϵ}_2^2\left(8145{ϵ}_2^2+21098\right)\right)]\sin \left(2\pi q_x\right).\hfill \\ \hfill \end{array}$$(A.6)