# The gravitational force field of proto-pancakes

Shohei Saga, Stéphane Colombi, Atsushi Taruya

► To cite this version:

Shohei Saga, Stéphane Colombi, Atsushi Taruya. The gravitational force field of proto-pancakes. Astronomy and Astrophysics - A&A, 2023, 678, pp.A168. 10.1051/0004-6361/202346968. hal-04122084

HAL Id: hal-04122084

https://hal.science/hal-04122084

Submitted on 10 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The gravitational force field of proto-pancakes

Shohei Saga<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />✉, Stéphane Colombi<fnref n="3" />, and Atsushi Taruya<fnref n="5" />,<fnref n="6" />

Received 22 May 2023 / Accepted 9 September 2023

## ABSTRACT

It is well known that the first structures that form from small fluctuations in a self-gravitating, collisionless, and initially smooth cold dark matter (CDM) fluid are pancakes. We studied the gravitational force generated by such pancakes just after shell crossing and have found a simple analytical formula for the force along the collapse direction, which can be applied to both the single- and multi-stream regimes. We tested the formula on the early growth of CDM proto-haloes seeded by two or three crossed sine waves. Adopting the high-order Lagrangian perturbation theory (LPT) solution as a proxy for the dynamics, we confirm that our analytical prediction agrees well with the exact solution computed via a direct resolution of the Poisson equation, as long as the local caustic structure remains sufficiently one-dimensional. These results are further confirmed by comparisons of the LPT predictions performed this way to measurements in Vlasov simulations performed with the public code COLDICE. We also show that the component of the force orthogonal to the collapse direction preserves its single-stream nature – it does not change qualitatively before or after the collapse – allowing sufficiently high-order LPT acceleration to be used to approximate it accurately as long as the LPT series converges. As expected, solving the Poisson equation on the density field generated with LPT displacement provides a more accurate force than the LPT acceleration itself, as a direct consequence of the faster convergence of the LPT series for the positions than for the accelerations. This may provide a clue as to how we can improve standard LPT predictions. Our investigations represent a very needed first step in the study of gravitational dynamics in the multi-stream regime analytically: we estimate, at the leading order in time and space, the proper backreaction on the gravitational field inside the pancakes.

**Key words.** gravitation – galaxies: kinematics and dynamics – dark matter

## 1. Introduction

Cold dark matter (CDM) is widely believed to dominate the matter content of the Universe and is microscopically modelled as a self-gravitating collisionless fluid that obeys the Vlasov-Poisson equations (Peebles 1982, 1984; Blumenthal et al. 1984). Due to its initially virtually null local velocity dispersion, the CDM phase-space distribution function can be described as a 3D sheet evolving in 6D phase-space. This sheet originally represents a single-stream flow, but as a consequence of the evolution under self-gravity, it can at some point self-intersect in configuration space. Such shell crossings mark the formation of singularities of various kinds, in particular pancake-like structures accompanied by apparent divergences of the density field (see e.g. Zel'dovich 1970; Arnold et al. 1982; Shandarin & Zeldovich 1989; Gouda & Nakamura 1988; Melott & Shandarin 1989; Hidding et al. 2014; Feldbrugge et al. 2018). After the first shell crossings, the sheet repeatedly self-interacts and folds to form intricate multi-stream structures, which include filaments and dark matter haloes. Although numerical simulations have revealed a number of details regarding the dynamical history of dark matter, it is still difficult to develop an analytical theory capable of predicting the entire

growth history of these structures in a fully self-consistent way due to the highly non-linear processes involved in multi-stream dynamics. An accurate description of the early stages of the evolution of multi-stream regions is fundamental to understanding these processes, and this must go through the calculation of the gravitational force field sourced by pancakes, which is the object of this article.

While we primarily focus on perturbation theory (PT) and multi-stream dynamics around the first shell-crossing time, it is also interesting to note that pancakes represent, at the coarse level (that is, after a certain level of smoothing to remove small-scale fluctuations in the matter distribution), fundamental building blocks of large-scale structures. They lie at the frontiers of voids in the galaxy distribution, and their presence can influence local Hubble flow estimates. Combined with voids, pancakes could therefore be used as powerful dynamical tools for studying, for example, the origin of the Hubble constant tension between the late and early Universe (see e.g. Verde et al. 2019; Riess 2020; Perivolaropoulos & Skara 2022), as discussed in, for example, Gurzadyan et al. (2023).

To approximate dark matter dynamics on large scales, one usually relies on PT (see e.g. Bernardeau et al. 2002, for a review). This theory has been widely used to predict, in the

<footnote n="1">¹ Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="2">² Institute for Advanced Research, Nagoya University, Furo-cho Chikusa-ku, Nagoya 464-8601, Japan
e-mail: shohei.saga@yukawa.kyoto-u.ac.jp</footnote>

<footnote n="3">³ Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Chikusa-ku, Nagoya 464-8602, Japan</footnote>

<footnote n="4">⁴ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris, CNRS, 5 place Jules Janssen, 92190 Meudon, France</footnote>

<footnote n="5">⁵ Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="6">⁶ Kavli Institute for the Physics and Mathematics of the Universe (WPI), Todai institute for Advanced Study, University of Tokyo, Kashiwa, Chiba 277-8568, Japan</footnote>
weakly non-linear regime, large-scale structure statistics such as power-spectrum and higher-order correlations. Many techniques have been developed in this framework (e.g. Crocce & Scoccimarro 2006, 2008; Valageas 2007; Taruya & Hiramatsu 2008; Matsubara 2008, 2011; Bernardeau et al. 2008, 2012, 2014; Pietroni 2008; Taruya et al. 2009, 2012; Valageas et al. 2013) and applied to observational data in order to constrain cosmological models (e.g. Blake et al. 2011; Beutler et al. 2017; Zhao et al. 2019; Ivanov et al. 2020; Tröster et al. 2020). In standard PT, a single-stream flow is imposed, and the small parameter is the Eulerian density contrast. However, the single-stream approximation is valid only during the early phases of structure formation, and its relevance can be questioned in the PT formalism (see e.g. Bernardeau et al. 2014; Blas et al. 2014; Nishimichi et al. 2016; Halle et al. 2020). Beyond the single-stream approximation, it is challenging to incorporate backreactions from the multi-stream regions into the analytical predictions in generic situations (see e.g. Rampf 2021, for a recent review). A way to account for multi-streaming in large-scale structure statistics consists in using effective field theory (Baumann et al. 2012; Carrasco et al. 2012; Hertzberg 2014; Baldauf et al. 2015). While effective field theory can provide meaningful constraints on cosmological models from observational data for practical applications (e.g. Ivanov et al. 2020; d'Amico et al. 2020), it involves free parameters that need to be calibrated with *N*-body simulations. For a more rigorous treatment beyond the phenomenological approach, one solid approach would be to consider high-order velocity moments of Vlasov equations that account for multi-streaming. Recently, Garny et al. (2023a,b) succeeded in making the problem analytically tractable, greatly facilitating the convergence of PT predictions. However, the approach they developed still depends on simplifying assumptions, and there is room for improvement. Accurately incorporating multi-stream effects into statistical predictions of large-scale structure is thus now considered an important factor for properly extracting cosmological information from observations.

Pancakes generally correspond to the first stages of multi-stream evolution; filaments and then dark matter haloes represent the next stages when shell crossing subsequently occurs along transverse directions of motion. The importance of accounting for multi-stream flows has therefore also been recognised in the formation process of proto-haloes that supposedly develop monolithically during an early violent relaxation phase (Lynden-Bell 1967) and subsequently, in the CDM scenario, merge hierarchically to form larger haloes with a universal density profile (Navarro et al. 1996, 1997). Although numerical investigations have revealed important properties of proto-haloes, for example their power-law density profile, $\rho(r) \sim r^{-\alpha}$, with a logarithmic slope $\alpha \approx 1.5$ (Moutarde et al. 1991; Diemand et al. 2005; Ishiyama 2014; Angulo et al. 2017; Ogiya & Hahn 2018; Delos et al. 2018a,b; Colombi 2021; Delos & White 2022, 2023; White 2022), there is no exact analytical theory that accounts for multi-stream dynamics inside dark matter haloes, despite the multiple approaches to the problem, for instance, self-similarity (Fillmore & Goldreich 1984; Bertschinger 1985; Henriksen & Widrow 1995; Sikivie et al. 1997; Yano & Gouda 1998; Yano et al. 2004; Zukin & Bertschinger 2010a,b; Alard 2013) or entropy maximisation (Lynden-Bell 1967; Hjorth & Williams 2010; Carron & Szapudi 2013; Pontzen & Governato 2013).

One way to push the dynamics beyond shell crossing while preserving total mass conservation consists in using a Lagrangian approach that follows the motion of matter

elements as functions of their initial position. Lagrangian perturbation theory (LPT), according to which the displacement field is the small parameter, has been widely employed to accurately describe the large-scale matter distribution in the quasi-linear regime, even beyond shell crossing (e.g. Zel'dovich 1970; Shandarin & Zeldovich 1989; Bouchet et al. 1992, 1995; Buchert 1992; Buchert & Ehlers 1993; Bernardeau 1994). Unfortunately, LPT quickly becomes inaccurate after shell crossing because it does not correctly account for the force feedback inside multi-stream regions. Recently, some progress has been achieved in this regard in the 1D case, which corresponds to the simplest dynamical setup, that is, a pure pancake that reduces to the interaction between infinite parallel planes. In 1D, linear LPT is exact prior to shell crossing (Novikov 1969), but subsequent multi-stream evolution still does not have a fully general analytical solution. However, it is possible to derive some approximate, asymptotically exact solutions just after shell crossing based on LPT but extended beyond collapse and correctly taking the backreaction of the gravitational force inside the multi-stream region into account (Colombi 2015; Taruya & Colombi 2017; Rampf et al. 2021). This article aims to generalise this post-collapse PT approach to 3D by calculating the gravitational force in proto-pancakes with the idea of developing a general analytical treatment of the early stages of multi-stream motion.

Contrary to the 1D case, one of the difficulties in developing a 3D post-collapse theory is the absence of an exact solution, even before shell crossing. This solution can be approached using sufficiently high-order LPT, but as the system approaches the first shell crossing, the perturbative treatment worsens and its convergence speed decreases (e.g. Zheligovsky & Frisch 2014; Rampf et al. 2015, 2022; Rampf & Hahn 2021, for work on the convergence radius). Convergence speed strongly depends on the nature of the initial conditions. It is facilitated when approaching quasi-1D initial conditions (e.g. Rampf & Frisch 2017; Saga et al. 2018) and is made more difficult when approaching axisymmetric configurations or spherical symmetry (e.g. Saga et al. 2018; Rampf 2019). In the simplified approach considered in the present work, we focus on pancakes seeded by a locally symmetric displacement field, a restrictive but still relatively generic setup that applies, for instance, to high peaks of a Gaussian random field Bardeen et al. (1986). To test our analytical predictions, we considered systems evolving from the initial crossed sine wave conditions, which we already studied at collapse in Saga et al. (2018) and slightly beyond shell crossing in Saga et al. (2022, hereafter, STC). For these systems, we find that LPT is able, at a sufficiently high order, to provide an accurate description of the density distribution around the shell crossing; therefore, they can be used to test our analytical predictions for the gravitational force field.

This paper is organised as follows. In Sect. 2 we introduce the basic equations of motion in the Lagrangian description and the LPT framework. We also discuss important aspects of the gravitational force calculation. In Sect. 3 we provide a simple analytical formula for the force inside a pancake in a case where the causal structure is seeded by a locally symmetric displacement. In Sect. 4 our analytical predictions are tested in systems with initial sine wave conditions, using high-order LPT as a proxy for the dynamics. The approximation of the force field is compared to the exact solution of the Poisson equation, in which the density field is generated by the LPT motion. This is followed in Sect. 5 by comparisons of the analytical predictions to the force field directly measured in Vlasov-Poisson simulations performed with the public code ColdICE (Sousbie & Colombi 2016).
Finally, Sect. 6 is devoted to the summary of our main findings. To supplement the main text, Appendices A–C, respectively test the validity of the series expansion at the third order of the displacement field used in Sect. 3, the accuracy of our Green function approach used to compute the force in the theoretical calculations, and finally the method we use to measure the force field in the ColdICE simulations.

## 2. Basic equations and gravitational force

We now introduce the basic equations in the Lagrangian description (Sect. 2.1) and discuss important aspects of the gravitational force calculation in the framework of LPT (Sect. 2.2).

### 2.1. Lagrangian dynamical framework

We considered the Lagrangian equation of motion of a fluid element at the Eulerian comoving position, $\mathbf{x}$, in the presence of gravity (e.g. Peebles 1980):

$$ \frac{d^2\mathbf{x}}{dt^2} + 2H\frac{d\mathbf{x}}{dt} = -\frac{1}{a^2}\nabla_x\phi(\mathbf{x}), \qquad (1) $$

where the quantities $a$, $H(t) = a^{-1}da/dt$, and $\phi(\mathbf{x})$ are the scale factor of the Universe, the Hubble parameter, and the Newton gravitational potential, respectively, and the operator $\nabla_x = \partial/\partial\mathbf{x}$ is the spatial gradient in Eulerian space. The gravitational potential is related to the matter density contrast $\delta(\mathbf{x}) = \rho(\mathbf{x})/\bar{\rho}-1$ with $\bar{\rho}$ the background mass density, through the Poisson equation

$$ \nabla_x^2\phi(\mathbf{x}) = 4\pi G a^2 \bar{\rho} \delta(\mathbf{x}). \qquad (2) $$

The Lagrangian description relates the initial, Lagrangian position, $\mathbf{q}$, of each mass element to the Eulerian position at time $t$, $\mathbf{x}(\mathbf{q}, t)$, through the introduction of the displacement field, $\mathbf{\Psi}(\mathbf{q}, t)$. In this framework, the Eulerian position $\mathbf{x}$ and velocity $\mathbf{v}$ are given by

$$ \mathbf{x}(\mathbf{q}, t) = \mathbf{q} + \mathbf{\Psi}(\mathbf{q}, t), \qquad (3) $$

$$ \mathbf{v} = a \frac{d\mathbf{\Psi}}{dt}. \qquad (4) $$

Assuming homogeneous initial density, $\rho(\mathbf{q})/\bar{\rho} = 1$, mass conservation reads $d^3\mathbf{q} = (1 + \delta(\mathbf{x})) d^3\mathbf{x}$ before the first shell-crossing time, $t_{sc}$. Hence,

$$ 1 + \delta(\mathbf{x}) = \frac{1}{J}, \qquad (5) $$

where the quantity $J = \det J_{ij}$ is the Jacobian of the matrix $J_{ij}$ defined by

$$ J_{ij}(\mathbf{q}, t) = \frac{\partial x_i(\mathbf{q}, t)}{\partial q_j} = \delta_{ij} + \Psi_{i,j}(\mathbf{q}, t). \qquad (6) $$

The first occurrence of $J = 0$ determines the first shell-crossing time, $t_{sc}$.

Until the first shell crossing, we can employ a perturbative treatment to predict the fluid motion, namely LPT (e.g. Zel'dovich 1970; Shandarin & Zeldovich 1989; Bouchet et al. 1992, 1995; Buchert 1992; Buchert & Ehlers 1993; Bernardeau 1994). In LPT, the displacement field, $\mathbf{\Psi}$, is considered a small quantity, which is systematically expanded as

$$ \mathbf{\Psi}(\mathbf{q}, t) = \sum_{n=1}^{\infty} \mathbf{\Psi}^{(n)}(\mathbf{q}, t). \qquad (7) $$

Assuming that the fastest growing modes dominate, the perturbative solution is approximated quite well by the form (see e.g. Bernardeau et al. 2002, and references therein)

$$ \mathbf{\Psi}^{(n)}(\mathbf{q}, t) = D_{+}^{n}(t) \mathbf{\Psi}^{(n)}(\mathbf{q}), \qquad (8) $$

with the time-dependent function $D_+(t)$ being the linear growth factor. Substituting Eq. (8) into Eq. (4), the velocity field is given by

$$ \mathbf{v}(\mathbf{q}, t) = a H f \sum_{n=1}^{\infty} n D_{+}^{n}(t) \mathbf{\Psi}^{(n)}(\mathbf{q}), \qquad (9) $$

where we define the linear growth rate by $f(t) \equiv d \ln D_+/d \ln a$.

In Eqs. (8) and (9) the nth order displacement $\mathbf{\Psi}^{(n)}$ is computed recursively by exploiting Eqs. (1) and (2) using well-known algebraic techniques (see e.g. Rampf 2012; Zheligovsky & Frisch 2014; Rampf et al. 2015; Matsubara 2015). Specific expressions for the three-sine-wave case examined below are given in STC and we do not consider it necessary to repeat them here. The Lagrangian framework we use in this article relies in practice on PT as a proxy of the dynamics. However, the analytic expressions computed in Sect. 3 are more general in the sense that they apply (asymptotically, that is, just beyond shell crossing) to any non-degenerate displacement field locally symmetric and Taylor expandable up to the third order in the Lagrangian position around the singularity of interest.

### 2.2. Gravitational force

To facilitate resolution of the Poisson Eq. (2), we decomposed it into two independent equations:

$$ \nabla_x^2 \phi_p = 4\pi G \bar{\rho} a^2 (1 + \delta(\mathbf{x})), \qquad (10) $$

$$ \nabla_x^2 \bar{\phi} = -4\pi G \bar{\rho} a^2. \qquad (11) $$

Thus,

$$ \phi(\mathbf{x}) = \phi_p(\mathbf{x}) + \bar{\phi}(\mathbf{x}). \qquad (12) $$

The first and second terms of the right-hand side of this equation are respectively the gravitational potential coming from the total density $\bar{\rho}(1 + \delta(\mathbf{x}))$ and the negative background density $-\bar{\rho}$ as counter-term. After solving Eqs. (10) and (11), the gravitational acceleration, also abusively referred to here as force $\mathbf{F}(\mathbf{x}) = -\nabla_x\phi(\mathbf{x})$, is given by

$$ \mathbf{F}(\mathbf{x}) = 4\pi G a^2 \bar{\rho} \left( \int \frac{d^d \mathbf{x}'}{2^{d-1}\pi} \frac{(1 + \delta(\mathbf{x}'))(\mathbf{x} - \mathbf{x}')}{|\mathbf{x} - \mathbf{x}'|^d} + \frac{1}{d} \mathbf{x} \right), \qquad (13) $$

$$ = 4\pi G a^2 \bar{\rho} \left( \int \frac{d^d \mathbf{q}}{2^{d-1}\pi} \frac{\mathbf{x} - \mathbf{x}'(\mathbf{q})}{|\mathbf{x} - \mathbf{x}'(\mathbf{q})|^d} + \frac{1}{d} \mathbf{x} \right), \qquad (14) $$

where $d$ represents the dimension of space ($d = 2$ or $3$ considered in this paper). In the second line, we used mass conservation in $d$-dimensional space: $(1 + \delta(\mathbf{x}')) d^d\mathbf{x}' = d^d\mathbf{q}$. The second term in parentheses represents the gravitational force arising from the background density in Eq. (11).

Thanks to the change of variable $\mathbf{x} \to \mathbf{q}$, the expression of the force in Eq. (14) does not depend explicitly on the density contrast nor does it require a detailed knowledge of the multi-stream structure, which is, regardless of its complexity, implicitly contained in the function $\mathbf{x}(\mathbf{q})$. This considerably facilitates the numerical calculation of the integrals without the need to solve a multi-value problem. We note that this property
was exploited before to compute the gravity field from caus-
tic rings (see e.g. Sikivie 1998, 1999; Charmousis et al. 2003;
Onemli 2006; Natarajan & Sikivie 2006, 2007; Duffy & Sikivie
2008; Onemli & Sikivie 2009; Chakrabarty & Sikivie 2018;
Tam 2012), and the calculations we perform here are analogous.

While integral (14) seems simpler to estimate than integral
(13), because it is performed in Lagrangian space, it remains
challenging to compute it analytically. The main goal of the
present work is to find explicit expressions that approximate it.

Another way to estimate the gravitational acceleration con-
sists in simply computing the second time derivative of the
Lagrangian displacement (provided that it is sourced by pure
gravity). Consider the presence of $n_S$ streams at Eulerian posi-
tion $\mathbf{x}$. In this case, the acceleration can be formally written as
the local average of the time derivatives of velocities over all the
streams weighted by the density of each stream:

$$
F(\mathbf{x}) = \frac{\sum_{i=1}^{n_S} \rho_i(\mathbf{x}) \Gamma_i(\mathbf{x})}{\sum_{i=1}^{n_S} \rho_i(\mathbf{x})}, \quad (15)
$$

$$
\Gamma_i(\mathbf{x}) \equiv \frac{d(a v_i(\mathbf{x}))}{d t}, \qquad (16)
$$

where the quantities $\rho_i$ and $a v_i$ stand for the density and peculiar velocity of $i$th stream, respectively. If one has access to the exact solution of the dynamics, this expression is somewhat trivial since, in this case, $\Gamma_i(\mathbf{x}) = \Gamma_j(\mathbf{x})$, $i \neq j$. On the other hand, if the accelerations in Eq. (15) are given by second time derivatives of the LPT displacement computed at some order, the force in Eq. (15) does not generally agree with Eq. (14) applied to the same displacement field, even in the single-stream regime. Indeed, the LPT solutions are derived not by directly solving the Poisson equation as in Eq. (14) but by perturbatively solving the Lagrangian equations of motion.

Equation (14), which is strongly non-linear in essence as it
can account accurately of multi-streams, acts as a re-summation
of the LPT acceleration: it is expected to provide a more accu-
rate prediction of the gravitational force field than the second
time derivative of the LPT displacement. However, as long as
the LPT series converges, we expect the higher-order LPT accel-
eration to converge to the force given by Eq. (14) in the single-
stream regime, and this property will turn out to be useful even in
the multi-stream regime when estimating the gravitational force
orthogonal to the shell-crossing direction (coplanar with the pan-
cake) by using Eq. (15).

**3. Analytical predictions for the gravitational force**

In this section we aim to compute the gravitational force shortly
after the first shell crossing in 3D space. As detailed in Sect. 3.1,
we restricted our computations to the formation of a symmet-
ric pancake seeded by a locally axisymmetric motion. The cal-
culation of the component of the force along the shell-crossing
direction is the most challenging. However, after Taylor expand-
ing the Lagrangian displacement field around the singularity just
after shell crossing, it turns out to be very similar to the pure
1D case already treated in Gurevich & Zybin (1995), Colombi
(2015), Taruya & Colombi (2017), and Rampf et al. (2021). In
particular it involves the resolution of a three-value problem
related to the three flows inside the proto-pancake, as detailed in
Sect. 3.2. The expression for the force along the shell-crossing
direction is given in Sect. 3.3. In this subsection we also argue
that the force field in the transverse direction should not be sig-
nificantly affected by the multi-stream nature of the flow, which

will allow us to estimate it directly as the second time derivative
of the displacement estimated with high-order LPT.

3.1. Main assumptions

In what follows, the calculations were all performed in 3D space,
but the extension to 2D is made straightforward by ignoring or
setting to zero all the contributions that depend on z. We also
supposed that the first shell crossing takes place at the origin,
**q** = **x** = **0**, along the x-axis direction, and also that the system
exhibits locally axisymmetric dynamics. This setup, illustrated
by Fig. 1, seemingly appears to be very particular, but locally
represents the expected motion around high peaks of a Gaussian
random field (see e.g. Bardeen et al. 1986).

Axisymmetric dynamics translates as follows on the dis-
placement, $\Psi(\mathbf{q})$:

$$
\begin{align}
\Psi_x(q_x, q_y, q_z) &= \Psi_x(q_x, -q_y, q_z) = \Psi_x(q_x, q_y, -q_z) \nonumber \\
&= -\Psi_x(-q_x, q_y, q_z), \tag{17}
\end{align}
$$

$$
\begin{align}
\Psi_y(q_x, q_y, q_z) &= \Psi_y(q_x, q_y, -q_z) = \Psi_y(-q_x, q_y, q_z) \nonumber \\
&= -\Psi_y(q_x, -q_y, q_z), \tag{18}
\end{align}
$$

$$
\begin{align}
\Psi_z(q_x, q_y, q_z) &= \Psi_z(-q_x, q_y, q_z) = \Psi_z(q_x, -q_y, q_z) \nonumber \\
&= -\Psi_z(q_x, q_y, -q_z). \tag{19}
\end{align}
$$

Here and hereafter, we omit the time dependence in the nota-
tions. Expanding these functions around the origin **q** = 0, we
have

$$
\Psi_x(\mathbf{q}) = \sum_{i,j,k=0} \frac{\psi_{2i+1\,2j\,2k}}{(2i+1)!(2j)!(2k)!} q_x^{2i+1} q_y^{2j} q_z^{2k}, \quad (20)
$$

$$
\Psi_y(\mathbf{q}) = \sum_{i,j,k=0} \frac{\psi_{2i\,2j+1\,2k}}{(2i)!(2j+1)!(2k)!} q_x^{2i} q_y^{2j+1} q_z^{2k}, \quad (21)
$$

$$
\Psi_z(\mathbf{q}) = \sum_{i,j,k=0} \frac{\psi_{2i\,2j\,2k+1}}{(2i)!(2j)!(2k+1)!} q_x^{2i} q_y^{2j} q_z^{2k+1}, \quad (22)
$$

with $\psi_{ijk}$ being some functions of time. Substituting
Eqs. (20)–(22) into Eq. (3), and neglecting $O(q^4)$ and higher-
order terms, we obtain

$$
x(\mathbf{q}) \approx (1 + \psi_{100}) q_x + \frac{1}{2} (\psi_{120} q_y^2 + \psi_{102} q_z^2) q_x + \frac{1}{6} \psi_{300} q_x^3, \quad (23)
$$

$$
y(\mathbf{q}) \approx (1 + \psi_{010}) q_y + \frac{1}{2} (\psi_{012} q_z^2 + \psi_{210} q_x^2) q_y + \frac{1}{6} \psi_{030} q_y^3, \quad (24)
$$

$$
z(\mathbf{q}) \approx (1 + \psi_{001}) q_z + \frac{1}{2} (\psi_{201} q_x^2 + \psi_{021} q_y^2) q_z + \frac{1}{6} \psi_{003} q_z^3. \quad (25)
$$

While this local representation of the motion is minimal,
it remains accurate shortly after collapse as illustrated by
Appendix A.

We now write the necessary conditions that the coefficients in
Eqs. (23)–(25) must satisfy for a pancake to exist near the origin.
It is important to note that these conditions do not necessarily
imply that a halo subsequently forms, this would require more
restrictive constraints.

Since we considered a system in which first shell crossing
just took place along the *x* direction at **q** = **0**, we imposed

$$
\frac{\partial x(\mathbf{0})}{\partial q_x} \equiv -h = 1 + \psi_{100} < 0, \quad 0 < h \ll 1, \qquad (26)
$$

$$
\frac{\partial y(\mathbf{0})}{\partial q_y} = 1 + \psi_{010} > 0, \qquad (27)
$$

$$
\frac{\partial z(\mathbf{0})}{\partial q_z} = 1 + \psi_{001} > 0. \qquad (28)
$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic representation of caustic structure shortly after shell crossing. (A) Left: schematic representation of a 3D Eulerian caustic shortly after shell crossing along the x-axis direction, together with a 2D slice with y = y₀ (light red plane). Right: Intersection of the caustic surface with the slice (red curve), also shown in the left panel. (B) Schematic representation of the x component of the Lagrangian coordinate, qₓ, as a function of the Eulerian coordinate, x, for fixed (y₀, z₀) (solid blue line). Given x₀, the solution of the three-value problem x(qₓ) = x₀ is given by qₓ,n(x₀, y₀, z₀) for n = 1, 2, and 3, as in Eq. (41).</figcaption></figure>

Additional constraints can be obtained from the expression of the Jacobian determinant, $J = \det \partial\mathbf{x}/\partial\mathbf{q}$ at the leading order in $\mathbf{q}$, which reads

$$J \approx \frac{1}{2}(1 + \psi_{010})(1 + \psi_{001})(-2h + \psi_{120}q_y^2 + \psi_{102}q_z^2 + \psi_{300}q_x^2). \quad (29)$$

From catastrophe theory (see e.g. Hidding et al. 2014; Feldbrugge et al. 2018, and references therein), the caustic surface should be an ellipsoid outside of which the Jacobian determinant must have a positive value $J > 0$, so we imposed

$$\psi_{120} > 0, \quad \psi_{102} > 0, \quad \psi_{300} > 0. \quad (30)$$

Finally, the smallness of the h parameter induces an additional simplification of Eqs. (23)–(25) if one supposes that Lagrangian coordinates are restricted to lie in the neighbourhood of the pancake, namely $q_x \sim q_y \sim q_z \sim O(h^{1/2})$ from Eq. (29). In this case, one realises that, when examining Eqs. (23)–(25),

$$x(\mathbf{q}) = O(h^{3/2}), \quad (31)$$

$$y(\mathbf{q}) \approx (1 + \psi_{010})q_y + O(h^{3/2}), \quad (32)$$

$$z(\mathbf{q}) \approx (1 + \psi_{001})q_z + O(h^{3/2}), \quad (33)$$

which implies $|x| \sim O(h^{3/2}) \ll |y|, |z| \sim O(h^{1/2})$, a signature of the pancake nature of the system: the extension of the caustic region along the x-axis is asymptotically infinitely smaller than that along the other axes in the limit $h \to 0$, as illustrated by panel A of Fig. 1. Accordingly, inside and in the vicinity of the caustic region, we can ignore the higher-order terms in the expressions of y and z, and reduce Eqs. (23)–(25) to

$$x(\mathbf{q}) \approx (1 + \psi_{100})q_x + \frac{1}{2}(\psi_{120}q_y^2 + \psi_{102}q_z^2)q_x + \frac{1}{6}\psi_{300}q_x^3, \quad (34)$$

$$y(\mathbf{q}) \approx (1 + \psi_{010})q_y, \quad (35)$$

$$z(\mathbf{q}) \approx (1 + \psi_{001})q_z. \quad (36)$$

Equations (34)–(36) represent our starting point to derive the x component of the gravitational force inside a pancake.

## 3.2. The three-value problem

Despite its apparent simplicity, Eq. (14) is not easily exploitable when it comes to estimate the gravitational force analytically, even with as simple expressions as Eqs. (34)–(36). Indeed, although the multi-stream nature of the flow does not appear explicitly in the integral (14), we see in Sect. 3.3 that the calculation of the x component of the gravitational force still requires solving the three-value problem implicit in Eqs. (34)–(36), that is, finding $\mathbf{q}$ given $\mathbf{x}$.

From Eqs. (35) and (36), we trivially obtain

$$q_y = \frac{y}{1 + \psi_{010}}, \quad q_z = \frac{z}{1 + \psi_{001}}. \quad (37)$$

The calculation of $q_x(\mathbf{x})$ is more complex because it requires solving the following cubic equation, as illustrated by panel B of Fig. 1,

$$q_x^3 + 3A(y, z)q_x + 2B(x) = 0, \quad (38)$$

where we define

$$A(y, z) = \frac{1}{\psi_{300}} \left( -2h + \psi_{120} \left( \frac{y}{1 + \psi_{010}} \right)^2 + \psi_{102} \left( \frac{z}{1 + \psi_{001}} \right)^2 \right), \quad (39)$$

$$B(x) = -\frac{3x}{\psi_{300}}. \quad (40)$$

The roots of cubic Eq. (38) are given by (see e.g. Abramowitz & Stegun 1972)

$$q_{x,n}(\mathbf{x}) = \omega^{n-1} \frac{-A(y, z)}{(\sqrt{D(\mathbf{x})} - B(\mathbf{x}))^{1/3}} + \omega^{4-n} (\sqrt{D(\mathbf{x})} - B(\mathbf{x}))^{1/3}, \quad (41)$$

for $n = 1, 2$, and 3. Here, the factor $\omega$ is one of the complex cubic roots of unity (i.e. $\omega = (-1 \pm i\sqrt{3})/2$), which are the solutions of $\omega^2 + \omega + 1 = 0$. We define the discriminant $D(\mathbf{x})$, which determines the properties of the roots (41), by (e.g. Abramowitz & Stegun 1972)

$$D(\mathbf{x}) = (A(y, z))^3 + (B(x))^2. \quad (42)$$
We first note that the equation $D(\mathbf{x}) = 0$ defines the caustics surfaces in Eulerian space:

$$ \left[ -2h + \psi_{120} \left( \frac{y}{1 + \psi_{010}} \right)^2 + \psi_{102} \left( \frac{z}{1 + \psi_{001}} \right)^2 \right]^3 + 9\psi_{300}x^2 = 0. \quad (43) $$

According to the sign of $D(\mathbf{x})$, the solutions of the cubic equation, Eq. (38), can then be classified as follows. If $D(\mathbf{x}) < 0$, we are inside the multi-stream region delimited by the caustic surfaces: the cubic equation has three real solutions, with $q_{x,2} < q_{x,3} < q_{x,1}$ and $q_{x,2} < 0 < q_{x,1}$. If $D(\mathbf{x}) > 0$, we are outside the multi-stream region: only $q_{x,1}$ or $q_{x,2}$ is real, and the other solutions are a complex conjugate to each other.

As a final note, from Eq. (43), the maximum values of the Eulerian coordinates of caustics along the x-, y-, and z-axes are given by

$$ x_{\max} = \sqrt{\frac{8h^3}{9\psi_{300}}}, \quad (44) $$

$$ y_{\max} = (1 + \psi_{010}) \sqrt{\frac{2h}{\psi_{120}}}, \quad (45) $$

$$ z_{\max} = (1 + \psi_{001}) \sqrt{\frac{2h}{\psi_{102}}}, \quad (46) $$

Once again, $x_{\max} \ll y_{\max}$, $z_{\max}$ in the limit $h \to 0$.

### 3.3. The force field inside a pancake

We now present the main results of this paper. Focusing on the dynamics just after the first shell crossing along the x direction, we first derive an analytical formula for the x component of the gravitational force after collapse. Next, we discuss orthogonal components and argue that they can be directly approximated by the high-order LPT acceleration.

By zooming in very closely on the multi-stream region, we notice that the structure of the density field becomes almost 1D when $h \to 0$. This follows from the fact the size of the pancake along the x-axis is much smaller than along orthogonal directions (see Eqs. (31)–(33) and also panel A of Fig. 1). Consequently, given the Eulerian position $(x, y_0, z_0)$ for $y_0$ and $z_0$ fixed, solving the Poisson Eq. (10), is asymptotically reduced to a 1D problem along the x-axis, as schematically shown in panel B of Fig. 1. In this case, we can neglect local variations of the density along the y and z directions, and Eq. (10) reduces to the following 1D Poisson equation for the x coordinate of the gravitational acceleration,

$$ \begin{aligned} \frac{d^2 \phi_p(x, y_0, z_0)}{dx^2} &= 4\pi G \bar{\rho} a^2 (1 + \delta(x, y_0, z_0)), \\ &= 4\pi G \bar{\rho} a^2 \frac{1}{(1 + \psi_{010})(1 + \psi_{001})} \left| \frac{\partial x}{\partial q_x} \right|^{-1}, \end{aligned} \quad (47) $$

where, in the second equality, we have used Eq. (5) with Eqs. (34)–(36).

To solve this equation, we followed in the footsteps of Colombi (2015), Taruya & Colombi (2017), and Rampf et al. (2021) and employed a Green function approach to derive the x component of the force in the multi-stream region, $F_x(\mathbf{x}_0) \approx$

$-d\phi_p(\mathbf{x}_0)/dx$ with $\mathbf{x}_0 = (x_0, y_0, z_0)$ (inside the pancake, $x \sim h^{3/2}$, and the background contribution $\phi$ is negligible):

$$ \begin{aligned} F_x(\mathbf{x}_0) &\approx \int dx \frac{4\pi G \bar{\rho} a^2}{(1 + \psi_{010})(1 + \psi_{001})} \sum_{n=0,1,2} \left| \frac{\partial x}{\partial q_{x,n}} \right|^{-1} \frac{1}{2} [\Theta_H(x - x_0) - \Theta_H(x_0 - x)], \\ &= \int dq_x \frac{4\pi G \bar{\rho} a^2}{(1 + \psi_{010})(1 + \psi_{001})} \frac{1}{2} \left[ \Theta_H \left( x \left( q_x, \frac{y_0}{1 + \psi_{010}}, \frac{z_0}{1 + \psi_{001}} \right) - x_0 \right) \right. \\ &\qquad \left. - \Theta_H \left( x_0 - x \left( q_x, \frac{y_0}{1 + \psi_{010}}, \frac{z_0}{1 + \psi_{001}} \right) \right) \right], \\ &= - \frac{4\pi G \bar{\rho} a^2}{(1 + \psi_{010})(1 + \psi_{001})} [q_{x,1}(\mathbf{x}_0) + q_{x,2}(\mathbf{x}_0) - q_{x,3}(\mathbf{x}_0)] \end{aligned} \quad (48) $$

where $\Theta_H$ represents the Heaviside step function, and the quantities $q_{x,n}$ are the three-value problem solutions given in Eq. (41). In the first line of Eq. (48), we have dropped the negligible contributions in the 1D Green function that are not proportional to the Heaviside step function.

Interestingly, one can also write, asymptotically,

$$ F_x(\mathbf{x}_0) \approx \Gamma_x[q_{x,1}(\mathbf{x}_0)] + \Gamma_x[q_{x,2}(\mathbf{x}_0)] - \Gamma_x[q_{x,3}(\mathbf{x}_0)], \quad (49) $$

where $\Gamma_x$, defined in Eq. (16), is the x coordinate of the LPT acceleration (computed by ignoring shell crossing). This approximation will be tested in Sect. 5.

Equation (48) applies only to the multi-stream region. However, as mentioned above, in the single-stream region, either $q_{x,1}$ or $q_{x,2}$ is real while the two other values are complex conjugate to each other. Taking the real part of right hand side of Eq. (48), the unphysical contribution arising from the complex conjugate solutions vanishes, leaving only the physical contribution from one real solution<fnref n="1" />:

$$ F_x(\mathbf{x}_0) \approx - \frac{4\pi G \bar{\rho} a^2}{(1 + \psi_{010})(1 + \psi_{001})} \mathrm{Re} \left[ q_{x,1}(\mathbf{x}_0) + q_{x,2}(\mathbf{x}_0) - q_{x,3}(\mathbf{x}_0) \right]. \quad (51) $$

Here, 'Re' denotes the real part.

Equation (51) is the most important formula in our work. Together with the solution of the cubic Eq. (41), it consists of a simple algebraic form for the x component of the force in the vicinity of the pancake, inside and outside it. This formula is the natural extension to 3D of the 1D calculations of Gurevich & Zybin (1995), Colombi (2015), Taruya & Colombi (2017), and Rampf et al. (2021).

The components $F_y$ and $F_z$ of the force orthogonal to the shell-crossing direction, hence coplanar with the pancake, should, on the other hand, remain quite insensitive to the effects of shell crossing, as we show with numerical tests in following sections. This means that if the LPT displacement field is entirely sourced by gravity, its second time derivative should provide a very good approximation of the local gravitational acceleration, even slightly beyond shell crossing. As a consequence, the acceleration directly derived from the high-order LPT solution (with no Taylor expansion at the third order in $\mathbf{q}$ space) should provide a good approximation for $F_y$ and $F_z$. Of course this is true only if one remains within the convergence radius of the LPT series, which is finite in the general case (e.g. Zheligovsky & Frisch 2014; Rampf et al. 2015; Rampf & Hahn 2021). To overcome the three-value problem, the acceleration

$$ F_x(\mathbf{x}_0) = - \frac{4\pi G \bar{\rho} a^2}{(1 + \psi_{010})(1 + \psi_{001})} \left[ q_{x,1}(\mathbf{x}_0) + q_{x,2}(\mathbf{x}_0) - (q_{x,3}(\mathbf{x}_0))^* \right] \quad (50) $$

where the asterisk denotes the complex conjugate.

<footnote n="1">¹ Equivalently, we can rewrite Eq. (48) as follows:</footnote>
can be obtained with the weighted average (15) over the different flows inside the pancake, even though this is not absolutely necessary, since it is expected, in the vicinity of the pancake, that $F_y[x(q_1)] \approx F_y[x(q_2)] \approx F_y[x(q_3)]$ and similarly for the z component of the force: this is due to the fact that the vectors $q_1, q_2$ and $q_3$ are nearly equal.

## 4. Examination of the accuracy of the formulas

In this section we test the validity and the accuracy of our prescription for calculating the post-collapse force. To this end, we assume Einstein-de Sitter cosmology and consider simplified initial conditions composed of two or three crossed sine waves, following the early investigations of Moutarde et al. (1991, 1995) and our previous works (Saga et al. 2018). We present these initial conditions in Sect. 4.1, where we also give practical details on the way the analyses are performed. Section 4.2 focuses on the force along the shell-crossing direction. We test, at different orders of the perturbative development, the accuracy of Eq. (51) against the exact solution obtained by direct resolution of Poisson equation. This is followed in Sect. 4.3 by a comparison, in the transverse direction, of the LPT acceleration to the exact force. All these analyses are performed very shortly after collapse to make sure that assumptions of Sect. 3.1 are verified. One indeed expects increasing discrepancies with time between the exact solution and the approximations of the dynamics, as examined in Sect. 4.4.

Throughout this section, we use the LPT solution as a proxy for the 'exact' Eulerian position. To calculate the exact force field, we simply injected the LPT solution into Eq. (14) and numerically performed the 2D or 3D integrals<fnref n="2" />. The comparisons to simulations performed in Sect. 5 will show that this proxy turns to provide a quite accurate description of the true displacements at the times we consider.

### 4.1. Two- and three-sine-wave initial conditions: Practical implementation

We considered initial conditions seeded by two or three crossed sine waves in a periodic box [-L/2, L/2] in which the initial displacement field at initial time $t_{ini}$ is expressed by (see STC)

$$ \Psi_i^{ini}(q, t_{ini}) = \frac{L}{2\pi} D_+(t_{ini}) \epsilon_i \sin\left(\frac{2\pi}{L} q_i\right), \quad (52) $$

with $\epsilon_i < 0$ and $|\epsilon_x| \ge |\epsilon_y| \ge |\epsilon_z|$. The initial density field, $\delta^{ini} \approx -\nabla_q \cdot \Psi^{ini} = D_+(t_{ini}) \sum_i |\epsilon_i| \cos(2\pi/L q_i)$, presents a small peak at the origin. Subsequently, mass elements fall towards the overdense central region, and shell crossing takes place at the origin. The initial time, $t_{ini}$, is set to satisfy $D_+(t_{ini})|\epsilon_i| \le 0.012 \ll 1$, so that the fastest growing mode approximation is accurate (Saga et al. 2018). We note that in the Einstein-de Sitter universe, the growth factor is simply proportional to the scale factor and we have $f = 1$. Hence, we hereafter use the scale factor to describe the time rather than $D_+$.

In this setup, the dynamics is determined by the ratios $\epsilon_{2D} = \epsilon_y/\epsilon_x$ and $\epsilon_{3D} = (\epsilon_y/\epsilon_x, \epsilon_z/\epsilon_x)$, respectively for two- and three-sine-wave initial conditions (STC). We considered three qualita-

tively different initial conditions, as detailed in Table 1: quasi-1D with $|\epsilon_x| \gg |\epsilon_{y,z}|$ (Q1D), anisotropic with $|\epsilon_x| > |\epsilon_y| > |\epsilon_z|$ (ANI), and axial-symmetric with $|\epsilon_x| = |\epsilon_y| = |\epsilon_z|$ (SYM). The Q1D and ANI initial conditions are the primary targets of our analyses because first shell crossing occurs only along the x direction, satisfying the assumptions of Sect. 3.1. On the other hand, for the SYM case, since shell crossing takes place simultaneously along two or three axial directions, the prescriptions proposed in Sect. 3.3 to approximate the force field becomes improper. However, the SYM case might provide insight into rare primordial haloes and remains mathematically interesting. In this case, we explore only numerically the force field by comparing predictions of higher-order LPT combined with the Green function approach (Eq. (14)) to measurements in the simulations (see Sect. 5).

Sine wave initial conditions (52) correspond to a low-order trigonometric polynomial, which makes LPT calculations relatively cheap and allows us to reach the 40th and 15th order for the 2D and 3D case, respectively (see Sect. 2 in STC for the explicit procedure)<fnref n="3" />. Having obtained the higher-order LPT solutions, we analytically predicted the force along the shell-crossing direction using Eq. (51), by expanding the LPT solutions around the origin up to the third order in terms of Lagrangian coordinates to have the expressions for the coefficients $\psi_{ijk}$ in Eqs. (20)–(22) hence in formula (51). Also, we analytically predicted the force orthogonal to the shell-crossing direction by directly using the LPT acceleration (without Taylor expanding it in $q$) and averaging it (numerically) over different flows at same Eulerian position $x$ using Eq. (15).

Because the collapse time, traced here by the value of the expansion factor $a_{sc}^{(n)}$, decreases when perturbative order $n$ augments (see e.g. STC), LPT predictions are in practice synchronised to their own respective shell-crossing times. Then the system is evolved beyond collapse by a very small amount of time:

$$ a^{(n)} = a_{sc}^{(n)} (1 + \Delta), \quad (53) $$

where $\Delta$ is a small parameter, such that the assumptions of Sect. 3.1 are valid, in particular $h \approx \Delta |\partial v_x / \partial q_x / (aH)|_{q=0, t=t_{sc}} \ll 1$. In practice, we take $\Delta = 0.001$ in Sects. 4.2 and 4.3, while it is allowed to be larger in Sect. 4.4. As a final note, we truncated the integral in Eq. (14) to a finite interval $[-q_{max}, q_{max}]$ in each dimension, setting $q_{max} = L/2$ for $F_x$ and $q_{max} = 10L$ for $F_{y,z}$ to insure convergence of the integral, as studied in Appendix B.

In the rest of the paper, we use the following units: a box size $L = 1$ and an inverse of the Hubble parameter at the present time $H_0 = 1$ for the dimensions of length and time, respectively. We also present the normalised force $\tilde{F}(x)$ given by

$$ \tilde{F}(x) = \frac{F(x)}{4\pi G a^2 \bar{\rho}}. \quad (54) $$

### 4.2. Force along the shell-crossing direction, $F_x$

Figure 2 displays, just after shell crossing, the x component of the force $F_x$ as a function of $x$ in the vicinity of the pancake for various values of $y$. The perturbative order increases from left to right, while initial conditions are given, from top to bottom, by Q1D-2SIN, ANI-2SIN, Q1D-3SIN, and ANI-3SIN.

<footnote n="2">² In numerically integrating Eq. (14), we used the NIntegrate function in Mathematica with the options {MaxRecursion -> 100 (10000), PrecisionGoal -> 5 (8), Method -> {Automatic, "SymbolicProcessing" -> 0}} for $F_x$ ($F_{y,z}$). In this instance, the computation of the results depicted in Figs. 2 or 3, on an 8-core CPU laptop, takes several hours.</footnote>

<footnote n="3">³ The high-order LPT solutions can be provided upon request as a Mathematica notebook.</footnote>
<table><caption>**Table 1.** Parameters of the runs performed with COIDICE (Sousbie & Colombi 2016).</caption><thead><tr><th>Designation</th><th>ε<sub>2D</sub> or ε<sub>3D</sub></th><th>ε<sub>x</sub></th><th>n<sub>g</sub></th><th>n<sub>s</sub></th><th>a<sup>LPT</sup><sub>sc</sub></th><th>a<sub>sc</sub></th><th>a<sub>sim</sub></th><th>Δ<sub>sim</sub> ≡ (a<sub>sim</sub> - a<sub>sc</sub>)/a<sub>sc</sub></th></tr></thead><tbody><tr><td>Quasi 1D</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Q1D-2SIN</td><td>1/6</td><td>-18</td><td>2048</td><td>2048</td><td>0.05279</td><td>0.05285</td><td>0.05402</td><td>0.02219</td></tr><tr><td>Q1D-3SIN</td><td>(1/6, 1/8)</td><td>-24</td><td>512</td><td>256</td><td>0.03814</td><td>0.03832</td><td>0.03907</td><td>0.01974</td></tr><tr><td>Anisotropic</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ANI-2SIN</td><td>2/3</td><td>-18</td><td>2048</td><td>2048</td><td>0.04534</td><td>0.04545</td><td>0.04601</td><td>0.01230</td></tr><tr><td>ANI-3SIN</td><td>(3/4, 1/2)</td><td>-24</td><td>512</td><td>512</td><td>0.02912</td><td>0.02919</td><td>0.03003</td><td>0.02890</td></tr><tr><td>Axial-symmetric</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SYM-2SIN</td><td>1</td><td>-18</td><td>2048</td><td>2048</td><td>0.04087</td><td>0.04090</td><td>0.04101</td><td>0.002717</td></tr><tr><td>SYM-3SIN</td><td>(1, 1)</td><td>-18</td><td>512</td><td>512</td><td>0.03236</td><td>0.03155</td><td>0.03201</td><td>0.01446</td></tr></tbody></table>

**Notes.** The first column indicates the designation of the run. The second column corresponds to the relative amplitudes of the initial sine waves, namely, $ε_{2D} = ε_y/ε_x$ and $ε_{3D} = (ε_y/ε_x, ε_z/ε_x)$ for two and three sine waves, respectively. The third column gives the value of $ε_x$. The fourth and fifth columns indicate the resolutions of the numerical simulations: the spatial resolution of the grid used to solve the Poisson equation, and the spatial resolution of the mesh of vertices used to construct the initial tessellation, respectively. In the sixth column, $a_{sc}^{LPT}$ is the shell-crossing time estimated by 40LPT for the two-sine-wave initial conditions and by 15LPT for the three-sine-wave initial conditions. The seventh, eighth, and ninth columns indicate the scale factor $a_{sc}$ at shell crossing measured in the simulations (see Appendix A1 in STC), the scale factor of the output time used for the analyses performed shortly after collapse in the simulations, and the fractional difference between $a_{sc}$ and $a_{sim}$, respectively.

We first notice the perfect agreement between the analytical formula (51) (solid curves) and the exact solution (14) (dashed curves), which fully justifies the mathematical relevance of the formalism developed in Sect. 3. In particular, the sharp transition between the single-stream and the multi-stream region is perfectly described by Eq. (51), with the variations with y (and z, not shown here) accounted for correctly. The discontinuity observed on the derivative of the x component of the force field is a typical signature of the presence of caustics, as found previously in the 1D case (see e.g. Gurevich & Zybin 1995; Colombi 2015), as well as in caustic rings (see e.g. Fig. 8 in Sikivie 1999).

Interestingly, the results do not strongly depend on the perturbation order in Fig. 2. This is because both the analytical formula (51) and the numerical calculation (14) rely on LPT displacement, which shows fast convergence behaviour with LPT order after synchronisation of shell-crossing times (see Fig. 10 in STC).

Obviously, our approach works because parameter $Δ$ in Eq. (53) is very small, which makes the pancake quasi-1D, independently of the values of $ε_{2D}$ and $ε_{3D}$ as long as $ε_y, ε_z < ε_x$. We see in Sect. 4.4 how the accuracy of the description deteriorates with increasing $Δ$, that is, with increasing time interval after collapse.

### 4.3. Force orthogonal to the shell-crossing direction

We now turn to the component $F_y$ of the force orthogonal to the shell-crossing direction (hence coplanar with the pancake) inside and in the vicinity of the pancake, as shown in Fig. 3 very shortly after collapse (note that the plot of the z component, $F_z$, in the 3D case, would give results very similar to $F_y$, so is not shown here). As already mentioned in Sect. 3.3, contrary to $F_x$, the exact gravitational force field given by Eq. (14), displayed as dashed curves, is not significantly affected by the presence of caustics, with a weak dependence on $x$, $F_y ∝ y$ and $F_z ∝ z$ close to the origin, while $F_y$ and $F_z$ are locally even with respect to $z$ and $y$, respectively (not shown on the figures). In particular we find, by comparing $F_{y,z}$ just before and after collapse ($Δ = ±0.001$), that it hardly changes during this period of time. These results suggest that shell crossing along the x direction does not strongly affect the dynamics along other axes. This property allows us to still use the high-order LPT solution, that

is, the acceleration computed as the second time derivative of the LPT displacement, to describe the y and z components of the force as long as the LPT series converges, as shown by the solid curves in Fig. 3 after averaging other the streams according to Eq. (15). The solid curves in Fig. 3 converge to the exact solution when the perturbative order, $n$, increases, except the 15th order is still insufficient for ANI-3SIN. Convergence is eased when approaching quasi-1D initial conditions (Q1D-2SIN and Q1D-3SIN), thanks to the much faster convergence of the LPT series in this case. We also checked in Eq. (15) that the LPT accelerations computed inside each stream are nearly the same.

Another important finding is that the dashed curves are only weakly dependent on LPT order, which is a consequence of synchronisation of LPT predictions with their own collapse time. We see in Sect. 5 that they actually agree very well with the measurements in COIDICE simulations. As shown by STC, the density field, which stems from synchronised LPT displacements and sources the force field, is very well described by LPT predictions even of relatively low order. On the other hand, the solid curves in Fig. 3 show that the LPT acceleration converges much more slowly than the LPT displacement, even after synchronisation. This is a rather obvious consequence of the fact that time derivatives give more weight than cumulative quantities to deviations from the exact solution, since these last ones increase with time. Similarly, STC found that convergence of LPT for velocities fields was slower than for density fields. We thus note that Eq. (14) provides, at fixed LPT order, a much more accurate way to compute the gravitational acceleration than the second time derivative of the LPT displacement, as already mentioned in last paragraph of Sect. 2.2. In other words, it provides a re-summation of the gravitational force that could be used to improve on LPT predictions even prior to shell crossing (i.e. even in the absence of synchronisation; see Fig. 6 in Sect. 5).

### 4.4. Time dependence

In the analyses performed in the two previous sections, we made sure that the assumptions of Sect. 3.1 were valid by taking a very small value of $Δ = 0.001$ in Eq. (53). In Fig. 4 we investigate how the results change when $Δ$ increases, up to the values $Δ_{sim}^{LPT} = Δ_{sim}(a_{sc}/a_{sc}^{LPT})$, with $Δ_{sim}$ listed in Table 1, that are used for measurements of Sect. 5 in COIDICE simulations.
<figure><img src="image_3.png" /></figure>

**Fig. 2.** *x* component of force as a function of *x* inside and in the vicinity of a pancake seeded by three sine waves in the context of LPT dynamics, with $-1.8 \le x/x_{max} \le 1.8$, where $x_{max}$ is the maximum extension of the caustic region along the *x*-axis, as given by Eq. (44) in our approximate formalism. The output time is set to $\Delta = 0.001$ with synchronisation (see Eq. (53)). In each panel, the curves of various colours correspond to different values of *y* (and *z* = 0 in the 3D case), namely $y/y_{max} = 0$ (orange), 0.75 (green), and 1.5 (red, outside the multi-stream region), where $y_{max}$ is the maximum extension of the caustic region along the *y*-axis (Eq. (45)). From top to bottom, we consider Q1D-2SIN, ANI-2SIN, Q1D-3SIN, and ANI-3SIN. From left to right, the predictions are made with 5, 20 (10), and 40LPT (15LPT) for two-sine-wave (three-sine-wave) initial conditions. The dashed and solid lines respectively represent the ‘exact’ force given by Eq. (14) and the analytic prediction (51).

We first focus on $F_x$ by examining left panels of Fig. 4. As expected, the analytical prediction (51) deviates more and more from the exact solution as $\Delta$ increases. For a given value of $\Delta$, for example $\Delta = 0.01$, deviations are slightly larger for ANI cases than Q1D cases, in agreement with intuition. Yet, even for non Q1D initial conditions, the local quasi-1D nature of the dynamics dominates if $\Delta$ is small enough. Obviously our approach has its limits, as illustrated by bottom right panel of the left part of Fig. 4. Importantly, deviations between dashes and solid lines in the left part of Fig. 4 are not related to performances of high-order LPT but instead, to our truncation to the third order of the Taylor expansion in Eqs. (20)–(22), as illustrated by Fig. A.1. This also explains why deviations also increase with the value of $y$, when passing from the green to the red curve in the left panels of Fig. 4. We note that it would be possible to correct the solutions (41) of the three-value problem by including perturbatively

<figure><img src="image_4.png" /></figure>

**Fig. 3.** Same as Fig. 2 but for the *y* component of the force, $F_y$, in the range $-1.8 \le y/y_{max} \le 1.8$ for $x/x_{max} = 0$ (orange), 0.75 (green), and 1.5 (red). The dashed and solid lines represent, respectively, the force given in Eq. (14) and the LPT acceleration given in Eq. (15). Due to the weak dependence of $F_y$ on the presence of the caustic and the very small range of values of $x$ considered, the curves for $x/x_{max} = 0, 0.75$, and 1.5 nearly perfectly overlap.

higher-order terms in the Taylor expansion in *q*, which would most certainly improve the agreement between the solid and the dashed curve in the ANI-3SIN case, even for $\Delta = 0.0289$. Of course, such a procedure cannot apply to arbitrarily large $\Delta$. It would moreover become pointless since LPT becomes a worse prescription of the dynamics as $\Delta$ increases. Indeed LPT has a finite convergence radius in time and does not take into account the feedback of the gravitational force inside the caustics, which is the objective we have in mind when computing the force field.

Next, we focused on $F_y$ and examined the right panels of Fig. 4^4, which confirm what we find in the previous section for $\Delta = 0.001$: except for the ANI-3SIN case, the LPT acceleration given by Eq. (15) agrees well with the exact solution (14) at all times considered. It stays pretty insensitive to the existence of the multi-stream region, as long as the pancake remains very thin and preserves the 1D nature of the local dynamics, which facilitates convergence of the LPT series. This is obviously not the case of ANI-3SIN, where $F_y$ presents significant variations with $x$ that increase with time, especially for the exact solution. We note that the good agreement between theory and exact

<footnote n="4">⁴ The *z* coordinate of the force, not shown here, would give very similar results.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Time dependence of the x and y components of the force for the various initial sine wave conditions considered in this article, as indicated in the panels. In each row, the output time, traced by the parameter $\Delta$ in Eq. (53), increases from left to right, as indicated in each panel, with $\Delta = \Delta_{sim} (a_{sc}/a_{sc}^{LPT})$ in the right panels, where $\Delta_{sim}$ corresponds to the output time in the simulations discussed in Sect. 5 (see also Table 1). Note that in this last case, the value indicated in the figure is not $\Delta$ but $\Delta_{sim}$, but the differences are small. The exact solution given by Eq. (14) applied to the LPT displacement (dashes) is compared to the analytical formula (51) for $F_x$ and with the LPT acceleration (15) for $F_y$ (solid lines). The two top and two bottom rows use 40LPT and 15LPT, respectively, to compute the displacement field and its acceleration. In the two bottom rows, we assume $z=0$.</figcaption></figure>

solution is obviously greatly facilitated by 40LPT in 2D, while only 15LPT is used in the 3D cases. Still, although the Q1D-3SIN case uses only the 15LPT solution for the acceleration, theory and exact solutions still agree very well with each other. On the other hand, in the ANI-3SIN case, the LPT acceleration deviates significantly from the exact force given by Eq. (14) and this even before shell crossing, while the 15LPT displacement itself remains a very good approximation of the true displacement measured in simulations, even for $\Delta = \Delta_{sim}^{LPT}$, as shown below in Sect. 5.

**5. Comparison to simulations**

In the previous sections, we tested the accuracy of force field calculations by using LPT as a proxy of the dynamics. We now rely on actual measurements in Vlasov-Poisson simulations to test LPT itself, since it is used for computing the coefficients in Eqs. (34)–(36) that lead to the asymptotic expression (51) as

well as the gravitational acceleration as a second time derivative of the LPT displacement in Eqs. (49) and (15).

The simulations we used to conduct the analyses were performed by Colombi (2021) and STC with the public Vlasov code CO1DICE (Sousbie & Colombi 2016). This solver directly follows the evolution of a self-gravitating 3D (2D) phase-space sheet in 6D (4D) phase-space with an adaptive tessellation of tetrahedra (triangles). The initial configuration uses a regular pattern of $n_s$ vertices with null velocities to construct the tessellation. At the beginning of the simulations, this pattern is perturbed with Zel'dovich motion following Eq. (52). During runtime, Poisson equation is solved in CO1DICE on a mesh of fixed resolution $n_g$. Table 1 indicates the values of $n_s$ and $n_g$ adopted for our runs, as well as the expansion factor $a_{sim}$ of the snapshot used for the analyses slightly beyond first shell crossing. Other technical details on the simulations we used are already provided in Colombi (2021) and STC, so we do not repeat them here. In Appendix C, we explain how we measured the force field from the tessellation.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison between simulations and analytical predictions for $F_x$ (left panels) and $F_y$ (right panels) slightly beyond first shell-crossing time. From top to bottom, we consider Q1D-2SIN, ANI-2SIN, Q1D-3SIN, and ANI-3SIN initial conditions, respectively. For each initial condition, we present the force for three different values of y and/or z, as indicated in each panel. The largest values of y and z are outside the multi-stream region. The dots stand for the measurements in the CO1DICE simulations. The solid green curves in the left panels give the theoretical prediction for $F_x$ from the analytical formula (51). The solid blue lines correspond to the theoretical predictions for $F_x$ and $F_y$ obtained from the LPT acceleration using Eq. (49) (only in the multi-stream region) and Eq. (15), respectively. The dashed orange lines represent the force field resulting from solving the Poisson Eq. (14), when using the LPT displacement instead of the supposedly exact positions of particles in the CO1DICE runs.</figcaption></figure>

To compare LPT predictions to measurements in simulations, we proceeded as in Sect. 4 and synchronised high-order LPT solutions (here, 40LPT for the 2D case and 15LPT for the 3D case) as follows. First, LPT solutions were evolved to their own shell-crossing time, designated in Table 1 by $a_{sc}^{\text{LPT}}$ (namely $a_{sc}^{(40)}$ and $a_{sc}^{(15)}$ for 2D and 3D, respectively, in the notations of Sect. 4). Next, time was advanced up to expansion factor $a = a_{sc}^{\text{LPT}} + (a_{\text{sim}} - a_{sc})$, where $a_{sc}$ and $a_{\text{sim}}$ are respectively the shell-crossing time measured by STC in the simulations and the expansion factor of the snapshot we considered for the analyses, as listed in Table 1. We note that the output times of the simulations, which are different for each initial condition, are, in terms of relative expansion factor values, closer to the shell-crossing time in the following order: ANI-2SIN, Q1D-3SIN, Q1D-2SIN, and ANI-3SIN (see the last column in Table 1).

Figure 5 compares analytical predictions to measurements in the simulations of the force field for the various sine wave initial conditions introduced in Sect. 4.1, except the axisymmetric ones,

which are discussed further below. The first thing to notice when examining this figure is the excellent agreement between the measurements (black dots) and the numerical solution of Poisson equation based on the Green function approach, Eq. (14), when applied to the density field sourced by the LPT displacement after synchronisation (orange dashes). This confirms the earlier investigations of STC on the density field, which compared post-collapse predictions based on the ballistic approximation, where the velocity field is frozen at collapse time, to the same simulations data. At the times considered here, the ballistic approximation should not differ much from LPT predictions pushed beyond shell crossing as performed in the present work. Such a good agreement between LPT and simulations fully validates the conclusions of the analyses of Sect. 4. Thus, in Fig. 5 we observe the same matches and discrepancies between the solid green (left panels) or blue (right panels) curves and the black dots or orange dashes, respectively, that we can discern in Fig. 4 between the solid curves and the dashes. This confirms
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Same as the right panels of Fig. 5, but without synchronisation to collapse time at each LPT order, i.e. for a given initial condition, all the curves are calculated at the same expansion factor, $a_{\text{sim}}$. In addition, several LPT orders are considered, as indicated in the right panels, to illustrate the convergence speed.</figcaption></figure>

again the validity of our theoretical predictions when two conditions are met: (i) sufficiently short period of time after collapse, so that our local Taylor expansion of the displacement field remains accurate enough to estimate the component of the force along the shell-crossing direction (orthogonal to the pancake) with Eq. (51), and (ii) sufficiently high LPT order for the LPT displacement to be accurate around shell-crossing time, as well as its second time derivative used in Eq. (15) as an approximation for the component of the force orthogonal to the direction of collapse (coplanar with the pancake). Obviously, these two conditions are not met for ANI-3SIN as already discussed in detail in Sect. 4, but are facilitated in the Q1D cases.

Interestingly, Eq. (49) (solid blue curves in the left panels of Fig. 5) approximates very well the $x$ component of the force, except again for ANI-3SIN. Summing up correctly the LPT accelerations in the multi-stream regions using either Eq. (49) or Eq. (15) also provides, shortly after shell crossing, a very good self-consistent approximation of all the components of the force field provided that the LPT series is converged.

The synchronisation process plays an important role in the analytical prediction of the force along the direction of shell crossing, $F_x$, because collapse time changes with LPT order and the effect of shell crossing is dramatic on $F_x$. However, this is not really the case for the orthogonal component – coplanar

with the pancake – that pretty much preserves the single-stream behaviour of the motion, as already discussed in Sect. 4.3. We double-checked this property by reproducing the right panels of Fig. 5 on Fig. 6 but without synchronisation of LPT predictions. For the 40th LPT order and the 15th LPT order respectively in 2D and 3D, we indeed distinguish no difference between the two figures. Additionally, various LPT orders are considered on Fig. 6. Interestingly, the LPT accelerations (solid lines) show slow convergence, while the forces computed with the Green function method from the LPT displacements (Eq. (14), dashed lines) show very fast convergence, even for ANI-3SIN. In other words, as mentioned in Sect. 4.3, Eq. (14) represents a re-summation of LPT that provides a much more accurate description of the gravitational acceleration than the second time derivative of the LPT displacement. This property might turn very useful for future applications.

Before closing this section, we examine in Fig. 7 axial-symmetric initial conditions, SYM-2SIN and SYM-3SIN. In this case, shell crossing occurs simultaneously along all the axes of the dynamics, so we expect qualitatively different behaviour from the Q1D and ANI cases. Strictly speaking, axial-symmetric configurations have zero weight from the statistical point of view but can in practice still be present at the coarse level, for example a very high peak in random Gaussian fields that are expected to be rounder (see e.g. Bardeen et al. 1986).

The axisymmetric case does not correspond to the superposition of 2 and 3 pancakes in the 2D and 3D cases, respectively (see Gurevich & Zybin 1995, for analytical predictions of the gravitational potential under these assumptions). Indeed, the caustic structure stemming from simultaneous shell crossings along several directions is convoluted, as shown by STC, which translates into multiple discontinuities of the derivative of the force field along each axis. For instance, in the top-left panel of Fig. 7 there are four sharp transition points on the force field instead of the two in the top-left panel of Fig. 5. Even if the analytical predictions discussed in Sect. 3 are irrelevant for SYM cases, it is still possible to use the Green function approach in Eq. (14) with the high-order LPT displacement field used as a proxy for the dynamics, as indicated by the dashed curves on Fig. 7. While 40LPT successfully reproduces the post-collapse gravitational field in the 2D case, 15LPT is not accurate enough. Turning to the 3D case, 15LPT is not even good from the qualitative point of view, as already noticed for STC for the density field itself. However, SYM-3SIN is further away from shell crossing than SYM-2SIN in terms of the parameter $\Delta$ as shown in last column of Table 1, which explains partly the worse performance of LPT predictions for SYM-3SIN compared to SYM-2SIN. Another and more obvious reason for this lies in the much higher contrasts expected in SYM-3SIN compared to SYM-2SIN, which can introduce non negligible feedback effects from the multi-stream region, even very shortly after shell crossing, as already extensively discussed in STC.

## 6. Summary

In this article we have studied the gravitational force field generated by pancakes shortly after collapse. Restricting ourselves to the case where the displacement field sourcing the pancake is locally symmetric, we derived approximations for the force field by combining fundamentals of catastrophe theory and high-order LPT that we carefully validated on systems seeded by two- or three-sine-wave initial conditions. Our analyses include comparisons of theoretical predictions to measurements in Vlasov simulations performed with the public code COLDICE. The important
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Fig. 5 but for the SYM-2SIN (top panels) and SYM-3SIN (bottom panels) initial conditions and the fact that there are no blue or green curves, because our analytical recipes apply only to shell crossing occurring along one direction. Only the *x* component of the force is shown in Fig. 7 due to the symmetric nature of the system. In addition, the 15LPT prediction from Eq. (14) is shown as a dashed light orange curve in the 2D case.</figcaption></figure>

assumption we made is that the time lapse after collapse, quantified here by the parameter $\Delta$ given in Eq. (53), is very short, allowing one to assume that the multi-stream region is locally infinitely thin. The main results of our work can be summarised as follows:

(i) The calculation of the component $F_x$, the gravitational force aligned with the direction of the shell crossing (that is, orthogonal to the pancake), comes down to solving a three-value problem that reduces to the resolution of a third-order polynomial in the limit $\Delta \ll 1$. This process is very analogous, not surprisingly, to the 1D case and leads to explicit expressions in terms of Eulerian coordinates, Eqs. (37), (39)–(41), and (51). Calculation of the various time-dependent coefficients intervening in the expression of $F_x$ relies on a Taylor expansion of the LPT displacement at the third order in the Lagrangian position, while staying within the convergence radius of the displacement expanded as a high-order series in time.

(ii) The component $F_y$ (and $F_z$), the gravitational force orthogonal to the direction of shell crossing (that is, in the same plane as the pancake), is rather insensitive to the presence of caustics. It can thus be predicted with the usual LPT acceleration, that is, the second time derivative of the LPT displacement, as long as the LPT acceleration converges as a series expansion in time. While Eq. (15) can be used to perform averaging over several streams, it is not absolutely needed. However, the three-value problem still needs to be solved, either numerically as we did in this work or approximately using Eqs. (37), (39)–(41).

(iii) Equation (51) is asymptotically equivalent to Eq. (49), with the *x* component of the acceleration given by the second time derivative of the LPT displacement, which in turn makes the approach fully consistent with point (ii).

(iv) A much higher LPT order is needed for the acceleration to converge than for the displacement. This is a natural consequence of the properties of the LPT series, of which the convergence speed decreases with successive time derivatives (see e.g. Rampf et al. 2023). However, one can numerically solve the Poisson equation using the Green function approach embodied by integral Eq. (14) based on

the knowledge of the LPT displacement to reach – as is obvious – a convergence level that is as good as for the displacement. Hence, Eq. (14) acts as a re-summation procedure of LPT for accurately computing the gravitational field.

(v) While quasi-1D initial conditions facilitate convergence of PT predictions, our calculations also apply to pancakes seeded by peaks with a general ellipsoid shape and remain accurate as long as the pancake remains very flat, but become inaccurate in the extreme case where collapse takes place simultaneously along all the axes of the dynamics.

In this work we have used LPT as an approximation of the dynamics shortly beyond collapse, but there are limitations to such an approach due to the finite radius of convergence with respect to time of the perturbative series. A safer point of view, adopted and tested successfully by STC on our sine wave cases, consists in using a ballistic approximation, where velocities (and in the present case, LPT accelerations) are frozen from collapse. While for the sake simplicity we did not adopt it here, it would be easy to implement this ballistic approximation to compute the gravitational force field of pancakes instead of using the procedure described in points (i) and (ii). The advantage of the ballistic approach is it allows us to overcome the problem of LPT convergence, which is known to be guaranteed in practice up to collapse (see e.g. Rampf & Hahn 2021, STC).

The calculations presented in this work can in principle be generalised to an arbitrary (smooth and non-degenerate) displacement field. Clearly, the case of a displacement field deriving from a potential in Lagrangian space (i.e. vorticity-free in Lagrangian space) can be treated easily. While not fully representative of the exact dynamics, since zero vorticity prior to shell crossing is only expected in Eulerian space, a vorticity-free displacement in Lagrangian space remains a very good approximation up to the second order in the LPT framework, which includes the Zel'dovich approximation. In the case of a displacement field deriving from a potential, it is easy to see that, with the proper combination of affine transformations in both Lagrangian and Eulerian space, one can obtain equations similar to Eqs. (34)–(36), but with additional $q_x^2q_y$, $q_x^2q_z$, $q_y^2$, and $q_z^2$ terms in the right-hand side of Eq. (34), which makes solving the three-value problem and the force field calculation slightly more involved. We postpone more general analyses that do not impose the local symmetries (18)–(19) to future work.

The calculation of the force field generated in the vicinity of a proto-pancake represents the first step towards an accurate treatment of the dynamics in the multi-stream regime. It involves computing corrections to the motion inside the pancakes due to the force backreaction from the multi-stream regions, thereby extending the 1D calculations of Colombi (2015), Taruya & Colombi (2017), and Rampf et al. (2021) to the 3D case. Of course, this approach remains limited, as it is expected to work only shortly after shell crossing, though it might be possible to combine it with an adaptive smoothing procedure to accurately predict large-scale structure statistics, such as the power spectrum, even in the non-linear regime (see e.g. Taruya & Colombi 2017; Halle et al. 2020, for the 1D case). Shell crossing can also locally take place along other axes of the dynamics, leading to the formation of proto-filaments and protohaloes. This is also followed by violent relaxation, which is a quick folding of the phase-space sheet in multiple directions. Understanding in detail the early evolution of proto-pancakes remains crucial to understanding how multi-flow dynamics is initiated and how this affects the early evolution of the statistics of the large-scale matter distribution.
**Acknowledgements.** We thank C. Rampf for reading and commenting on the draft. S.S. is supported by JSPS Overseas Research Fellowships. This work was supported in part by MEXT/JSPS KAKENHI Grant Numbers JP20H05861 and JP21H01081 (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 COIDICE was carried out using the HPC resources of CINES (Occigen supercomputer) under the GENCI allocations c2016047568, 2017-A0040407568 and 2018-A0040407568. Post-treatment of COIDICE data were performed on HORI-ZON cluster of Institut d'Astrophysique de Paris.

## References

Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions, National Bureau of Standards Applied Mathematics Series, 55

Alard, C. 2013, MNRAS, 428, 340

Angulo, R. E., Hahn, O., Ludlow, A. D., & Bonoli, S. 2017, MNRAS, 471, 4687

Arnold, V. I., Shandarin, S. F., & Zeldovich, I. B. 1982, Geophys. Astrophys. Fluid Dyn., 20, 111

Baldauf, T., Mercolli, L., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 123007

Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15

Baumann, D., Nicolis, A., Senatore, L., & Zaldarriaga, M. 2012, J. Cosmol. Astropart. Phys., 2012, 051

Bernardeau, F. 1994, ApJ, 427, 51

Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1

Bernardeau, F., Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 78, 103521

Bernardeau, F., Crocce, M., & Scoccimarro, R. 2012, Phys. Rev. D, 85, 123519

Bernardeau, F., Taruya, A., & Nishimichi, T. 2014, Phys. Rev. D, 89, 023502

Bertschinger, E. 1985, ApJS, 58, 39

Beutler, F., Seo, H.-J., Saito, S., et al. 2017, MNRAS, 466, 2242

Blake, C., Brough, S., Colless, M., et al. 2011, MNRAS, 415, 2876

Blas, D., Garny, M., & Konstandin, T. 2014, J. Cosmol. Astropart. Phys., 2014, 010

Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Nature, 311, 517

Bouchet, F. R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5

Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575

Buchert, T. 1992, MNRAS, 254, 729

Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375

Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, J. High Energy Phys., 2012, 82

Carron, J., & Szapudi, I. 2013, MNRAS, 432, 3161

Chakrabarty, S. S., & Sikivie, P. 2018, Phys. Rev. D, 98, 103009

Charmousis, C., Onemli, V., Qiu, Z., & Sikivie, P. 2003, Phys. Rev. D, 67, 103502

Colombi, S. 2015, MNRAS, 446, 2902

Colombi, S. 2021, A&A, 647, A66

Crocce, M., & Scoccimarro, R. 2006, Phys. Rev. D, 73, 063520

Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533

d'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 005

Delos, M. S., Erickcek, A. L., Bailey, A. P., & Alvarez, M. A. 2018a, Phys. Rev. D, 97, 041303

Delos, M. S., Erickcek, A. L., Bailey, A. P., & Alvarez, M. A. 2018b, Phys. Rev. D, 98, 063527

Delos, M. S., & White, S. D. M. 2022, ArXiv e-prints [arXiv:2209.11237]

Delos, M. S., & White, S. D. M. 2023, MNRAS, 518, 3509

Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389

Duffy, L. D., & Sikivie, P. 2008, Phys. Rev. D, 78, 063508

Feldbrugge, J., van de Weygaert, R., Hidding, J., & Feldbrugge, J. 2018, J. Cosmol. Astropart. Phys., 2018, 027

Fillmore, J. A., & Goldreich, P. 1984, ApJ, 281, 1

Garny, M., Laxhuber, D., & Scoccimarro, R. 2023a, Phys. Rev. D, 107, 063539

Garny, M., Laxhuber, D., & Scoccimarro, R. 2023b, Phys. Rev. D, 107, 063540

Gouda, N., & Nakamura, T. 1988, Prog. Theor. Phys., 79, 765

Gurevich, A. V., & Zybin, K. P. 1995, Phys.-Usp., 38, 687

Gurzadyan, V. G., Fimin, N. N., & Chechetkin, V. M. 2023, A&A, 672, A95

Halle, A., Nishimichi, T., Taruya, A., Colombi, S., & Bernardeau, F. 2020, MNRAS, 499, 1769

Henriksen, R. N., & Widrow, L. M. 1995, MNRAS, 276, 679

Hernquist, L., Bouchet, F. R., & Suto, Y. 1991, ApJS, 75, 231

Hertzberg, M. P. 2014, Phys. Rev. D, 89, 043521

Hidding, J., Shandarin, S. F., & van de Weygaert, R. 2014, MNRAS, 437, 3442

Hjorth, J., & Williams, L. L. R. 2010, ApJ, 722, 851

Ishiyama, T. 2014, ApJ, 788, 27

Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, J. Cosmol. Astropart. Phys., 2020, 042

Lynden-Bell, D. 1967, MNRAS, 136, 101

Matsubara, T. 2008, Phys. Rev. D, 77, 063530

Matsubara, T. 2011, Phys. Rev. D, 83, 083518

Matsubara, T. 2015, Phys. Rev. D, 92, 023534

Melott, A. L., & Shandarin, S. F. 1989, ApJ, 343, 26

Moutarde, F., Alimi, J. M., Bouchet, F. R., Pellat, R., & Ramani, A. 1991, ApJ, 382, 377

Moutarde, F., Alimi, J. M., Bouchet, F. R., & Pellat, R. 1995, ApJ, 441, 10

Natarajan, A., & Sikivie, P. 2006, Phys. Rev. D, 73, 023510

Natarajan, A., & Sikivie, P. 2007, Phys. Rev. D, 76, 023505

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nishimichi, T., Bernardeau, F., & Taruya, A. 2016, Phys. Lett. B, 762, 247

Novikov, E. A. 1969, Sov. J. Exp. Theor. Phys., 30, 512

Ogiya, G., & Hahn, O. 2018, MNRAS, 473, 4339

Onemli, V., & Sikivie, P. 2009, Phys. Lett. B, 675, 279

Onemli, V. K. 2006, Phys. Rev. D, 74, 123010

Peebles, P. J. E. 1980, The large-scale Structure of the Universe (Princeton University Press)

Peebles, P. J. E. 1982, ApJ, 263, L1

Peebles, P. J. E. 1984, ApJ, 277, 470

Perivolaropoulos, L., & Skara, F. 2022, New A Rev., 95, 101659

Pietroni, M. 2008, J. Cosmol. Astropart. Phys., 2008, 036

Pontzen, A., & Governato, F. 2013, MNRAS, 430, 121

Rampf, C. 2012, J. Cosmol. Astropart. Phys., 2012, 004

Rampf, C. 2019, MNRAS, 484, 5223

Rampf, C. 2021, Rev. Mod. Plasma Phys., 5, 10

Rampf, C., & Frisch, U. 2017, MNRAS, 471, 671

Rampf, C., & Hahn, O. 2021, MNRAS, 501, L71

Rampf, C., Villone, B., & Frisch, U. 2015, MNRAS, 452, 1421

Rampf, C., Frisch, U., & Hahn, O. 2021, MNRAS, 505, L90

Rampf, C., Frisch, U., & Hahn, O. 2022, Phys. Rev. Fluids, 7, 104610

Rampf, C., Saga, S., Taruya, A., & Colombi, S. 2023, ArXiv e-prints [arXiv:2303.12832]

Riess, A. G. 2020, Nat. Rev. Phys., 2, 10

Saga, S., Taruya, A., & Colombi, S. 2018, Phys. Rev. Lett., 121, 241302

Saga, S., Taruya, A., & Colombi, S. 2022, A&A, 664, A3

Shandarin, S. F., & Zeldovich, Y. B. 1989, Rev. Mod. Phys., 61, 185

Sikivie, P. 1998, Phys. Lett. B, 432, 139

Sikivie, P. 1999, Phys. Rev. D, 60, 063501

Sikivie, P., Tkachev, I. I., & Wang, Y. 1997, Phys. Rev. D, 56, 1863

Sousbie, T., & Colombi, S. 2016, J. Comput. Phys., 321, 644

Tam, H. 2012, ArXiv e-prints [arXiv:1205.1260]

Taruya, A., & Colombi, S. 2017, MNRAS, 470, 4858

Taruya, A., & Hiramatsu, T. 2008, ApJ, 674, 617

Taruya, A., Nishimichi, T., Saito, S., & Hiramatsu, T. 2009, Phys. Rev. D, 80, 123503

Taruya, A., Bernardeau, F., Nishimichi, T., & Codis, S. 2012, Phys. Rev. D, 86, 103528

Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10

Valageas, P. 2007, A&A, 465, 725

Valageas, P., Nishimichi, T., & Taruya, A. 2013, Phys. Rev. D, 87, 083522

Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891

White, S. D. M. 2022, MNRAS, 517, L46

Yano, T., & Gouda, N. 1998, ApJS, 118, 267

Yano, T., Koyama, H., Buchert, T., & Gouda, N. 2004, ApJS, 151, 185

Zel'dovich, Y. B. 1970, A&A, 5, 84

Zhao, G.-B., Wang, Y., Saito, S., et al. 2019, MNRAS, 482, 3497

Zheligovsky, V., & Frisch, U. 2014, J. Fluid Mech., 749, 404

Zukin, P., & Bertschinger, E. 2010a, Phys. Rev. D, 82, 104045

Zukin, P., & Bertschinger, E. 2010b, Phys. Rev. D, 82, 104044
# Appendix A: Series expansion of the displacement field

<figure><img src="image_9.png" /></figure>

**Fig. A.1.** x coordinate as a function of $q_x$ for Q1D-3SIN (top panels) and ANI-3SIN (bottom panels) initial conditions at various output times, with $\Delta = 0.005$ (left), 0.01 (centre), and $\Delta_{sim} = a_{sim}/a_{sc} - 1$, as given in the right column of Table 1 (right). The solid and dashed curves correspond to $q_y = 0$ and $q_y \approx a_{y}^{\text{max}}$, where $q_{y}^{\text{max}}$ is the maximum extension of the caustic in Lagrangian space along the $q_y$-axis, while $q_z = 0$. We present the results predicted by the pure 15LPT solution (black) and its Taylor expansion expression in the Lagrangian coordinate up to up to $O(q^3)$ (orange) and up to $O(q^5)$ (blue).

In Sect. 3.1 we derive an analytical expression for computing the x component of gravitational force, starting from the Taylor expansion of the displacement field with respect to the Lagrangian coordinate up to the third order in **q**, given in Eqs. (23)–(25). Figure A.1 tests the accuracy of Eqs. (23)–(25) in $q_x$-x space for our non-axisymmetric three-sine-wave setups. Clearly, at the dynamical times considered in this paper, a third-order approximation for the displacement field is sufficient in the Q1D case, but significant deviations can be seen in the bottom-right panel of the figure. As expected, these deviations increase when moving away from the origin of coordinates, hence with increasing, $|x|$, $|y|$ (and $|z|$ in 3D). They explain, at least partly, the discrepancies between theory and exact solution in Fig. 4. Another significant source of error is the approximation of y and z coordinates at the linear order in **q** (Eqs. 35 and 36), which

we do not examine here. We note that, as long as higher-order corrections remain small, it would be possible to correct perturbatively the three-value solution derived from the third-order Taylor expansion to include higher-order contributions. Here, we notice that the fifth order already brings very significant improvements, though we should also check what happens in other planes.

# Appendix B: Calculation of the force: Truncation of integral (14)

Throughout this article, we compare our analytical predictions for the force to direct calculations of the integral (14) over a finite interval, [-$q_{max}$, $q_{max}$]. Our systems seeded by initial sine wave conditions are periodic, and hence the problem comes down to computing the gravitational contribution of each matter element inside the simulation volume of size $L$, plus all the periodic replica. Many techniques exist in the literature to perform such a calculation quickly, for instance Ewald summation (see e.g. Hernquist et al. 1991). Our brute force technique allows us to have rather accurate estimates of pairwise interactions between nearby elements of mass, which is necessary given the very thin nature of the caustic structures, but fails to account for the contribution of all remote elements, including replica, which can make the force field calculation inaccurate.

To validate the choice of our somewhat simplistic technique, Fig. B.1 presents the **x** (top panels) and **y** components (bottom panels) of the force for our (non-axisymmetric) sine wave initial conditions, restricting to the local neighbourhood of the pancake. All the output times are the same as in Fig. 5, which compares analytical predictions to simulations. We see that the **x** component of the force at the scales of interest is dominated by the influence of the caustics and is thus quite insensitive to the replica given these small values of **x**, so a truncation of the integral (14) at $q_{max} = L/2$ is enough for all the initial conditions (including SYM-2SIN and SYM-3SIN, not shown in Fig. 5). On the other hand, the **y** component of the force can be significantly affected by the replicas for Q1D initial conditions, but from lower panels of Fig. 5, $q_{max} = 5L/2$ seems sufficient for adequate convergence of the integral (14). In practice, we adopted the values listed in Table C.1 for the measurements in the COLDICE simulations, and $q_{max} = L/2$ for $F_x$, $q_{max} = 20L/2$ for $F_{y,z}$, for other calculations using the LPT displacement in Eq. (14).
<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Tests of the choice of the bounds of integral (14). The x and y coordinates of the force are plotted, respectively, in the top and bottom panels, for various choices of the integration range [-q<sub>max</sub>, q<sub>max</sub>] with q<sub>max</sub> ranging from L/2 to 20L/2, as indicated in the upper-left panel. From left to right, we consider Q1D-2SIN, ANI-2SIN, Q1D-3SIN, and ANI-3SIN initial conditions. In the top panels, various values of y are considered, as indicated in the second panel of the top row, while the bottom panels only assume x = 0 since other values of x would not lead to significantly different results. We also set z = 0 in the 3D cases considered in the two top-right and two bottom-right panels. Note that the output time chosen in this figure is the same as in Fig. 5. We see that all the curves corresponding to different values of q<sub>max</sub> are superimposed on each other, except for q<sub>max</sub> = L/2 for the Q1D cases.</figcaption></figure>

Appendix C: Force measurement in simulations

In COLDICE, the dark matter distribution is represented with
an adaptive tessellation of simplices, that is, an ensemble of
connected triangles and tetrahedra, respectively in 2D and 3D.
To compute the force field F̃(x) at a given point of space, we
employed a direct approach consisting in replacing each simplex
with a set of particles. To do so, if needed, each simplex was
refined isotropically ℓmax times. At the end, each sub-simplex
was replaced with a single particle lying at the barycentre of
the its vertices. To obtain the highest accuracy, the refinement
process exploits the quadratic nature of each simplex thanks to
additional tracers used in COLDICE. To avoid divergences due to
the singular nature of the force field induced by a point particle
distribution, we introduced additional softening as follows:

$$
\tilde{F}(\mathbf{x}) = \frac{m}{2^{d-2} \times 2\pi} \frac{\mathbf{x} - \mathbf{x}_0}{(|\mathbf{x} - \mathbf{x}_0|^2 + \varepsilon^2)^{d/2}}, \quad (C.1)
$$

where, *d* = 2 or 3 is the dimension of space, *m* and *x*<sub>0</sub> are respec-
tively the (normalised) mass and the position of the particle. This
equation does not account for the background correction propor-
tional to **x** in Eq. (14), which has to be added at the end of the
calculation. Depending on the coordinate of the force field con-
sidered, a set of periodic replica of each particle with positions

$$
x_0 + (i, j, k) L \tag{C.2}
$$

can contribute, with integers *i*, *j*, *k* ∈ [-*n*<sub>rep</sub>, *n*<sub>rep</sub>], which is equiv-
alent to restricting the integral (14) to the interval [-*q*<sub>max</sub>, *q*<sub>max</sub>],
with *q*<sub>max</sub> = (*n*<sub>rep</sub> + 1/2)*L*. The values of *q*<sub>max</sub> we adopted are
listed in Table C.1. They are large enough according to the con-
vergence tests discussed in Appendix B.

Because we examined the system just after shell crossing, it
is important to have an accurate description of the phase-space
sheet inside a region containing the caustics. This region was

chosen to be a thin rectangular parallelepiped covering the inter-
vals x ∈ [-x<sub>max</sub>, x<sub>max</sub>], y ∈ [-y<sub>max</sub>, y<sub>max</sub>], and z ∈ [-z<sub>max</sub>, z<sub>max</sub>]
in each dimension (except for the last one, which is in 3D only).
Outside this region, each simplex was replaced with a single par-
ticle and inside it, each simplex was adaptively refined l<sub>max</sub> ≥ 1
times, with a value of l<sub>max</sub> depending on initial conditions.

The various parameters introduced in this appendix are listed
in Table C.1, which completes Table 1 of the main text. W note, for
instance, that ε = x<sub>max</sub>/64. Because of the high mass resolution of
the 2D simulations (n<sub>s</sub> = 2048), is was enough to take ℓ<sub>max</sub> = 2 for
Q1D-2SIN, ANI-2SIN and SYM-2SIN. Furthermore, due to the
lower cost of the force calculation in 2D, many periodic replica
could be used, n<sub>rep</sub> = 30. On the other hand, the 3D case is
much more involved computationally, which imposes us to adopt
a much smaller value of n<sub>rep</sub>, yet still large enough according to
the tests performed in Appendix B. The most delicate case was
Q1D-3SIN, where ℓ<sub>max</sub> = 4 was necessary to have (barely) suf-
ficiently accurate measurements of the force field. There are two
reasons for this. First, in this simulation, the initial tessellation has
a lower number of simplices, 6 × n<sub>s</sub><sup>3</sup>, with n<sub>s</sub> = 256, compared to
n<sub>s</sub> = 512 for ANI-3SIN and SYM-3SIN, which already imposes
at least one additional refinement level in Q1D-3SIN compared to
the two other cases. Second, the discrete nature of the sampling of
the phase-space sheet we adopted can introduce systematic biases
to the force field that are stronger in Q1D settings due to the fact
that particles, which initially form a regular mesh, tend to cluster
together much more along the x direction inside and in the vicinity
of the caustics than in orthogonal direction(s). This induces arti-
ficial fluctuations to the force field, especially on axes orthogonal
to x (see e.g. middle panel of the right group in third row of Fig. 5),
that would not appear with a proper, smooth representation of the
phase-space sheet. While it would be possible (but not trivial) to
implement a more optimal calculation of the force field, we did
not find it necessary for our analyses.
<table><caption>**Table C.1.** Parameters involved in the measurement of the force field in the Vlasov simulations, as described in Appendix C.</caption><thead><tr><th>Designation</th><th>q<sub>max</sub></th><th>x<sub>max</sub></th><th>y<sub>max</sub></th><th>z<sub>max</sub></th><th>l<sub>max</sub></th><th>ε</th></tr></thead><tbody><tr><td><em>Quasi 1D</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Q1D-2SIN</td><td>L/2 for F<sub>x</sub>, 30.5L otherwise</td><td>0.00075</td><td>0.15</td><td></td><td>2</td><td>1.1718750 × 10<sup>-5</sup></td></tr><tr><td>Q1D-3SIN</td><td>L/2 for F<sub>x</sub>, 5L/2 otherwise</td><td>0.0006</td><td>0.15</td><td>0.20</td><td>4</td><td>9.375 × 10<sup>-6</sup></td></tr><tr><td><em>Anisotropic</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ANI-2SIN</td><td>L/2 for F<sub>x</sub>, 30.5L otherwise</td><td>0.0004</td><td>0.025</td><td></td><td>2</td><td>6.25 × 10<sup>-6</sup></td></tr><tr><td>ANI-3SIN</td><td>L/2 for F<sub>x</sub>, 7L/2 otherwise</td><td>0.002</td><td>0.0275</td><td>0.07</td><td>2 for q<sub>max</sub> = L/2<br>1 for q<sub>max</sub> = 7L/2</td><td>3.125 × 10<sup>-5</sup></td></tr><tr><td><em>Axial-symmetric</em></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SYM-2SIN</td><td>L/2</td><td>0.00008</td><td>0.00008</td><td></td><td>2</td><td>1.25 × 10<sup>-6</sup></td></tr><tr><td>SYM-3SIN</td><td>L/2</td><td>0.0015</td><td>0.0015</td><td>0.0015</td><td>0</td><td>2.34375 × 10<sup>-5</sup></td></tr></tbody></table>

**Notes.** From left to right: designation of the run, truncation parameter $q_{max}$, coordinates $x_{max}$, $y_{max}$ and $z_{max}$ of the upper right corner of the rectangular region where higher sampling of simplices is performed using $l_{max}$ successive refinements as indicated in the next column; finally, $\epsilon$ is the softening parameter of the force field (Eq. C.1).