# Post-collapse perturbation theory in 1D cosmology – beyond shell-crossing

Atsushi Taruya, Stéphane Colombi

► To cite this version:

Atsushi Taruya, Stéphane Colombi. Post-collapse perturbation theory in 1D cosmology – beyond shell-crossing. Monthly Notices of the Royal Astronomical Society, 2017, 470 (4), pp.4858-4884. 10.1093/mnras/stx1501. hal-01645194

HAL Id: hal-01645194

https://hal.science/hal-01645194v1

Submitted on 26 Aug 2025

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
# Post-collapse perturbation theory in 1D cosmology – beyond shell-crossing

Atsushi Taruya<fnref n="1" />,<fnref n="2" />* and Stéphane Colombi<fnref n="1" />,<fnref n="3" />

Accepted 2017 June 14. Received 2017 June 14; in original form 2017 February 6

## ABSTRACT

We develop a new perturbation theory (PT) treatment that can describe gravitational dynamics of large-scale structure after shell-crossing in the one-dimensional cosmological case. Starting with cold initial conditions, the motion of matter distribution follows at early stages the single-stream regime, which can, in one dimension, be described exactly by the first-order Lagrangian perturbation, i.e. the Zel'dovich solution. However, the single-stream flow no longer holds after shell-crossing and a proper account of the multistream flow is essential for post-collapse dynamics. In this paper, extending previous work by Colombi, we present a perturbative description for the multistream flow after shell-crossing in a cosmological setup. In addition, we introduce an adaptive smoothing scheme to deal with the bulk properties of phase-space structures. The filtering scales in this scheme are linked to the next-crossing time in the post-collapse region, estimated from our PT calculations. Our PT treatment combined with adaptive smoothing is illustrated in several cases. Predictions are compared to simulations and we find that post-collapse PT with adaptive smoothing reproduces the power spectrum and phase-space structures remarkably well even at small scales, where Zel'dovich solution substantially deviates from simulations.

**Key words:** methods: analytical – methods: numerical – large-scale structure of Universe.

## 1 INTRODUCTION

It is currently admitted that processes of structure formation in the Universe are mainly dominated at large scale by an invisible component called dark matter. Although the microscopic origin of dark matter is still unclear, it is macroscopically described as a self-gravitating collisionless fluid following the collisionless Boltzmann or Vlasov equation in a cosmological background,

$$ \left[ \frac{\partial}{\partial t} + \frac{\mathbf{v}}{a} \cdot \frac{\partial}{\partial \mathbf{x}} - \frac{1}{a} \nabla \Phi \cdot \frac{\partial}{\partial \mathbf{v}} \right] f(\mathbf{x}, \mathbf{v}, t) = 0, \quad (1) $$

supplemented with the Poisson equation for the Newton potential $\Phi$,

$$ \frac{1}{a^2} \nabla^2 \Phi = 4\pi G \left[ \frac{1}{a^3} \int d^3 \mathbf{v} f(\mathbf{x}, \mathbf{v}, t) - \bar{\rho}_m \right], \quad (2) $$

where $f(\mathbf{x}, \mathbf{v}, t)$ is the phase-space density at comoving position $\mathbf{x}$, peculiar velocity $\mathbf{v}$ and time $t$, $a$ is the expansion factor of the Universe, $\Phi$ is the gravitational potential and $\bar{\rho}_m$ is the average dark matter density.

In the standard picture of structure formation, dark matter was initially cold, i.e. with a virtually null local velocity dispersion, so the six-dimensional phase-space distribution is effectively reduced to a three-dimensional hypersurface and this remains true at all times thanks to the Hamiltonian nature of the system. At early times, dark matter thus follows the single-stream flow regime, with a velocity field $\bar{\mathbf{v}}$ uniquely determined as a function of position, and its evolution is that of a pressure-less fluid with a phase-space distribution function given by

$$ f(\mathbf{x}, \mathbf{v}, t) = \bar{\rho}_m a^3 \{1 + \delta(\mathbf{x}, t)\} \delta_D[\mathbf{v} - \bar{\mathbf{v}}(\mathbf{x}, t)], \quad (3) $$

where $\delta$ is the density contrast of the dark matter distribution, initially of very small magnitude, as well as $\bar{\mathbf{v}}$.

Substituting equation (3) into equations (1) and (2) and taking the zeroth and first velocity moments yields the Eulerian formulation of large-scale structure dynamics

$$ \frac{\partial \delta}{\partial t} + \frac{1}{a} \nabla \{1 + \delta\} \bar{\mathbf{v}} = 0, \quad (4) $$

$$ \frac{\partial \bar{\mathbf{v}}}{\partial t} + \frac{1}{a} \{\bar{\mathbf{v}} \cdot \nabla\} \bar{\mathbf{v}} = -\frac{1}{a} \nabla \Phi, \quad (5) $$

<footnote n="1">¹Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="2">²Kayli Institute for the Physics and Mathematics of the Universe (Kayli IPMU, WPI), Todai Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, CNRS UMR 7095 and UPMC, 98bis bd Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: ataruya@yukawa.kyoto-u.ac.jp</footnote>
$$ \frac{1}{a^2} \nabla^2 \Phi = 4\pi G \bar{\rho}_m \delta. \qquad (6) $$

Although the single-stream flow is, strictly speaking, valid only during the early phase of structure formation, the above equations have been shown in practice to describe accurately non-linear mode coupling in the weakly non-linear regime and provide a solid basis for perturbation theory (PT) calculations to predict statistical quantities of large-scale structure, such as the power spectrum or the two-point correlation function of the matter distribution (see, e.g. Bernardeau et al. 2002).

With the advent of cosmological observations aiming at precisely mapping galaxy clustering at large scales, techniques have been developed to improve on the slow convergence of the perturbative expansion as well as to accelerate higher order calculations (e.g. Crocce & Scoccimarro 2008; Bernardeau, Crocce & Scoccimarro 2008; Taruya et al. 2009; Taruya et al. 2012, and references therein). Improved calculations involving the next-to-next-to-leading order called two-loop give a remarkable agreement with cosmological N-body simulations at weakly non-linear scales and have been applied to observations. However, this does not ensure the convergence of perturbative calculations including higher order corrections, e.g. three-loop. In fact, a direct calculation at three-loop order suggests a very large ultraviolet contribution to the large-scale modes through non-linear mode coupling (Bernardeau, Taruya & Nishimichi 2014; Blas, Garny & Konstandin 2014) and indicates a breakdown of higher order perturbative expansions even at large scales.

Deficiency of PT calculations has been also highlighted in recent numerical analyses. Nishimichi, Bernardeau & Taruya (2016) directly measured the coupling between different scales in cosmological N-body simulations and found that the actual contribution from small scales to the large-scale modes is suppressed, as opposed to the prediction of PT based on the single-stream approximation. These facts imply that the validity of the single-stream treatment is questionable even at large scales, and higher order perturbative corrections need to be remedied with a proper account of small-scale dynamics, where the multistream flow contribution is important.

One way to account for small-scale dynamics consists of using some ansatz to summarize the main physical effects in the multistream regime, such as Burgers' equation (e.g. Gurbatov, Saichev & Shandarin 1989; Bernardeau & Valageas 2010) and various alternatives (e.g. Sahni & Coles 1995), which roughly amount to adding a source term in the right-hand side of equation (5). In particular, it was proposed recently to use effective fluid equations giving account of the non-vanishing stress tensor arising in equation (5) when calculating the first velocity moment of Vlasov equation in the multistream regime. This effective-field theory approach has attracted a lot of interest and has been studied in detail (e.g. Baumann et al. 2012; Carrasco, Hertzberg & Senatore 2012; Hertzberg 2014; Baldauf, Mercolli & Zaldarriaga 2015). The drawback of this approach is however that the parameters in the stress tensor characterizing the small-scale dynamics need to be calibrated with N-body simulations to keep predictions of perturbative calculations under control. Furthermore, these parameters generally vary with cosmology and with time and no prediction with PT is really possible in this framework independently of N-body simulations.

An alternative approach that we consider in this paper consists in going back to a more fundamental description, i.e. the Vlasov-Poisson system, equations (1) and (2), and trying to follow accurately the phase-space distribution function in the multistream regime, which is essential for describing the formation of dark matter haloes. For this purpose, it is useful to employ Lagrangian PT (e.g. Zel'dovich 1970; Shandarin & Zeldovich 1989; Bouchet

et al. 1992; Buchert 1992; Buchert & Ehlers 1993; Bouchet et al. 1995; Bernardeau et al. 2002, and references therein), where the small parameter is the displacement field. We shall consider one-dimensional gravitational dynamics in standard (three-dimensional) cosmology. In this case, large-scale structure dynamics is described by the gravitational interaction of massive parallel infinite planes moving left and right along a fixed axis, while Hubble expansion is taking place as usual in all the directions following standard Friedman-Lemaître equations. Despite its simplicity, one-dimensional dynamics defined as such still displays a rich physical content which somewhat shares the same features as 3D clustering. This is partly the reason why the 1D model has been studied thoroughly in the literature in various formulations of the cosmological set up (e.g. Rouet, Feix & Navet 1990; Aurell & Fanelli 2001; Fanelli & Aurell 2002; Miller & Rouet 2002; Yano et al. 2004; Manfredi et al. 2016, and references therein) and has recently attracted much attention (e.g. Benhaiem, Joyce & Sicard 2013; McQuinn & White 2016; Vlah et al. 2016; Baldauf et al. 2016). In particular, the Zel'dovich solution provides an exact solution for the dynamics of massive sheets before shell-crossing (Novikov 1969; Zel'dovich 1970; Shandarin & Zeldovich 1989), and thus, starting with Zel'dovich solution, a tractable perturbative treatment of multistream flows is made possible based on a Lagrangian description. The analysis in this paper is an extension of the method developed in Colombi (2015) to the cosmological setup. We will describe perturbatively post-collapse dynamics around the shell-crossing region and apply it to several cases including random initial conditions.

In addition, we shall present a novel regularization scheme reducing the impact of small-scale clustering and improving greatly our post-collapse PT predictions of large-scale structure statistics. The idea is to apply adaptive smoothing to initial density peaks and to better capture the bulk properties of phase-space structures in the post-collapse regions, where the interaction or merger of haloes is supposed to be significant. The idea is similar to the peak-patch treatment proposed by Bond & Myers (1996) and subsequent works (see, e.g. Monaco, Theuns & Taffoni 2002), but we here implement it in the PT prescription in order to better describe the late-time post-collapse dynamics. Indeed, our perturbative approach will not allow us to follow post-collapse dynamics beyond next-crossing time, although an iterative prescription such as proposed by Colombi (2015) could make this possible but is out of the scope of this work. Here, we will show instead that the predictions with adaptive smoothing reproduce simulations remarkably well even at non-linear scales, where Zel'dovich solution significantly deviates from simulations.

Note finally that an analytical study in 1D of course represents only a first step towards a proper description of 6D phase-space dynamics. Apart from a few examples including self-similar solutions (Fillmore & Goldreich 1984a,b; Bertschinger 1985; Ryden 1993; Lithwick & Dalal 2011; Alard 2013), little is indeed known analytically for the phase-space dynamics of Vlasov-Poisson systems, particularly in cosmology. Therefore, further development of analytical treatment in this framework seems indispensable, complementary to simulations, and even helpful to cross-check simulation codes. This is all the more motivated by the fact that simulations in 6D phase space have now become available thanks to recent efforts on the development of pure Vlasov codes (Yoshikawa, Yoshida & Umemura 2013; Sousbie & Colombi 2016; Hahn & Angulo 2016).

This paper is organized as follows. In Section 2, we begin by describing the basic setup of our calculations in one-dimensional cosmology. We then discuss in Section 3 the analytic treatment beyond shell-crossing and develop post-collapse PT. Section 4 introduces
a novel filtering scheme to improve the perturbative description of post-collapse dynamics, suited for a system with interacting clusters. In Section 5, analytic calculations in the post-collapse perturbative framework are tested against controlled N-body experiments. Several cases are considered, including the single initial sine wave and random initial conditions. Finally, Section 6 is devoted to discussion and conclusion.

## 2 1D COSMOLOGY

### 2.1 Basic setup

We consider an ensemble of massive parallel infinite planes moving along x-axis and interacting through gravitational force in the expanding Universe. The Lagrangian equations of motion of the planes are given by

$$ \frac{dx}{dt} = \frac{v}{a}, \qquad (7) $$

$$ \frac{dv}{dt} + H v = -\frac{1}{a} \nabla_x \phi, \qquad (8) $$

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

where $x(t)$ and $v(t)$ are respectively the comoving position and peculiar velocity of each plane, $\phi$ the gravitational potential, $\bar{\rho}_m$ the average matter density, $\delta$ the density contrast and $a$ the expansion factor of the Universe.

To simplify the equations, it is useful to introduce the super-conformal time $\tau$ defined by (e.g. Doroshkevich, Ryaben’kii & Shandarin 1973; Martel & Shapiro 1998),

$$ d\tau = \frac{dt}{a^2}, \qquad (10) $$

the new velocity $u$ and potential $\Phi$:

$$ u \equiv a v, \qquad \Phi \equiv a^2 \phi. \qquad (11) $$

Then, equations (7)–(9) are simplified:

$$ \frac{dx}{d\tau} = u, \qquad (12) $$

$$ \frac{du}{d\tau} = -\nabla_x \Phi, \qquad (13) $$

$$ \nabla_x^2 \Phi = 4\pi G \bar{\rho}_m a^4 \delta = \frac{3}{2} \Omega_{m,0} H_0^2 a \delta, \qquad (14) $$

where $\Omega_{m,0}$ is the matter density parameter of the Universe and $H_0$ is the Hubble constant.

To deal with this system, especially for the dynamics after shell-crossing, it is useful to introduce in the cold case considered here the Lagrangian coordinate $q$ defining the initial position of the planes in the absence of perturbation

$$ q \equiv x(\tau \to 0), \qquad (15) $$

and to express the subsequent position and velocity of each mass element as $x(q, t)$ and $u(q, t)$. In particular, mass conservation implies

$$ dq = [1 + \delta(x)] dx \implies \delta(x) = \left(\frac{\partial x}{\partial q}\right)^{-1} - 1. \qquad (16) $$

This equation is valid in the single-stream regime, i.e. as long as $x(q, \tau)$ remains monotonic as a function of $q$.

With the new expressions above, the solution of the equations of motion can be formally written as:

$$ x(q; \tau) = x(q; \tau_{\text{ini}}) + \int_{\tau_{\text{ini}}}^{\tau} d\tau' u(q; \tau'), \qquad (17) $$

$$ u(q; \tau) = u(q; \tau_{\text{ini}}) - \int_{\tau_{\text{ini}}}^{\tau} d\tau' \nabla_x \Phi[x(q; \tau'); \tau'], \qquad (18) $$

where $x(q; \tau_{\text{ini}})$ and $u(q; \tau_{\text{ini}})$ are the initial positions and velocities given at a starting time $\tau_{\text{ini}}$, which will be specified below.

In what follows, we consider the dynamics of the cosmological system in a finite-size box with periodic boundaries, $0 \le x \le L$. The solution of Poisson equation, equation (14), can be expressed in an integral form as:

$$ \Phi(x) = \frac{3}{2} \Omega_{m,0} H_0^2 a \times \int_0^L dx' \left[ -\frac{L}{2} \left\{ \left( \frac{|x-x'|}{L} - \frac{1}{2} \right)^2 - \frac{1}{12} \right\} \right] \delta(x'). \qquad (19) $$

The derivation of this expression is presented in Appendix A. Then, the force exerted on a mass element at position $x$ is given by:

$$ \begin{aligned} F(x) &\equiv -\nabla_x \Phi(x) \\ &= -\frac{3}{2} \Omega_{m,0} H_0^2 a \times \left[ \int_0^L dx' \frac{\delta(x')}{2} \{ \Theta(x-x') - \Theta(x'-x) \} \right. \\ &\qquad \left. + \frac{1}{L} \int_0^L dx' x' \delta(x') \right], \end{aligned} \qquad (20) $$

where function $\Theta(x)$ represents the Heaviside step function. In the above, we used the fact that fluctuations averaged over space vanish, $\int_0^L dx' \delta(x') = 0$. In the limit $L \to \infty$, equation (20) of course converges to the well-known result in the infinite space.

### 2.2 Pre-collapse dynamics

In one dimension and in the cold case, Zel'dovich approximation (Novikov 1969; Zel'dovich 1970) is known to provide the exact solution for the dynamics of mass elements before shell-crossing. It can be explicitly written as

$$ x(q; \tau) = q + \psi(q) D_{+}(\tau), \qquad u(q; \tau) = \psi(q) \frac{dD_{+}(\tau)}{d\tau}. \qquad (21) $$

Here, function $D_+$ corresponds to the linear growth factor satisfying the following equation:

$$ \left[ \frac{d^2}{d\tau^2} - \frac{3}{2} \Omega_{m,0} H_0^2 a(\tau) \right] D_+(\tau) = 0. \qquad (22) $$

Note that in terms of cosmic time $t$, equation (22) reduces to the standard form of linear evolution equation:

$$ \left[ \frac{d^2}{dt^2} + 2H(t) \frac{d}{dt} - \frac{3}{2} \frac{\Omega_{m,0} H_0^2}{a^3(t)} \right] D_+(t) = 0. \qquad (23) $$

The Zel'dovich solution in equation (21) contains an arbitrary function $\psi(q)$ that we call displacement field. It is related to the linear density field $\delta_L(q)$ given at a very early time ($\tau_{\text{ini}} \to -\infty$ or $t_{\text{ini}} \to 0$) through

$$ \frac{d\psi(q)}{dq} D_+(\tau_{\text{ini}}) = -\delta_L(q; \tau_{\text{ini}}) = -\delta_L(q) D_+(\tau_{\text{ini}}). \qquad (24) $$
Since Zel'dovich solution is exact until shell-crossing, we do not necessarily assume that the evolved density field $\delta(x)$ is small. One may thus consider the situation where in a region around a Lagrangian coordinate $q_0$, the density field becomes large and the region will undergo shell-crossing at time $\tau_0$. Recalling the fact that in 1D, a shell-crossing point corresponds to an inflection point of the mapping from Lagrangian to Eulerian space, the actual conditions for shell-crossing are given by

$$ \left. \frac{\partial x}{\partial q} \right|_{q_0} = 0, \quad \left. \frac{\partial^2 x}{\partial q^2} \right|_{q_0} = 0, \quad \left. \frac{\partial^3 x}{\partial q^3} \right|_{q_0} > 0. \qquad (25) $$

At the time of shell-crossing, $\tau_0$, the solution (21) around the shell-crossing region can be expanded as follows:

$$ x(q; \tau_0) \approx q_0 + \psi(q_0) D_+(\tau_0) + \left\{ 1 + \frac{d\psi(q_0)}{dq_0} D_+(\tau_0) \right\} (q - q_0) + \sum_{n=2} \frac{1}{n!} \frac{d^n \psi(q_0)}{dq_0^n} D_+(\tau_0) (q - q_0)^n. \qquad (26) $$

Using equation (24), the conditions (25) for shell-crossing can be rewritten as

$$ \delta_L(q_0) = \frac{1}{D_+(\tau_0)}, \quad \left. \frac{d\delta_L(q)}{dq} \right|_{q_0} = 0, \quad \left. \frac{d^2\delta_L(q)}{dq^2} \right|_{q_0} < 0. \qquad (27) $$

In other words, shell-crossing takes place at local density peaks. Hence, in the 1D case, the conditions for shell-crossing are equivalent to peak constraints and collapse take places exactly when the Eulerian linear density contrast at the peak positions becomes equal to unity.

### 3 PERTURBATIVE TREATMENT OF POST-COLLAPSE DYNAMICS

We are interested in the dynamics of mass elements after shell-crossing, when the Zel'dovich solution is no longer valid. In this section, extending the work of Colombi (2015), we develop perturbative calculations to deal with the multistream motion around the shell-crossing region.

#### 3.1 Post-collapse PT

The basic formalism to treat post-collapse dynamics can be described as follows. Starting with the cold initial conditions in Section 2.2, we first follow pre-collapse dynamics with the exact Zel'dovich solution. Then, in regions undergoing shell-crossing, we switch to a perturbative treatment and compute the backreaction to the Zel'dovich flow based on an explicit functional form of the displacement field in the shell-crossing region. To be precise, we compute the force from equation (20) exerted at each position using the extrapolation of the Zel'dovich flow. Integrating the force over time, we obtain the correction for the velocity to the Zel'dovich motion from equation (18). Further integrating the corrected velocity over time, a correction for the position is obtained from equation (17). Throughout these calculations, we assume that the collapsing region is small, which allows us to express the phase-space configuration in terms of polynomial forms of low order of the Lagrangian coordinate.

Let us focus on a collapse point centred on Lagrangian position $q_0$ and perform a perturbative description of the multistream flow around $q = q_0$. As we mentioned earlier, the conditions for shell-crossing are given by equation (25). Shortly after shell-crossing

<figure><img src="image_2.png" /><figcaption>Figure 1. Geometrical configuration of phase space around a local density peak after shell-crossing, extrapolating the Zel'dovich solution (see equation 28). The structure of the system around the density peak is symmetric with respect to the Lagrangian coordinate $Q \equiv q - q_0$, where $q_0$ is the Lagrangian position of the shell-crossing point (see also Colombi 2015).</figcaption></figure>

time $\tau_0$, the multistream flow has just started growing in a small Lagrangian region around $q_0$ and the displacement field $\psi(q)$ can be approximately described by a third-order polynomial in $q - q_0$. Since the deviation from Zel'dovich flow is small, the motion of a mass element around $q = q_0$ may be expanded at third order in $q$ for $\tau > \tau_0$ as

$$ x(q; \tau) \approx A(q_0; \tau) - B(q_0; \tau) (q - q_0) + C(q_0; \tau) (q - q_0)^3 + \dots \qquad (28) $$

with the time-dependent coefficients A, B and C defined by

$$ A(q_0; \tau) \equiv x(q_0; \tau) = q_0 + D_+(\tau) \psi(q_0), \qquad (29) $$

$$ \begin{aligned} B(q_0; \tau) &\equiv - \left. \frac{\partial x}{\partial q} \right|_{q_0} = -1 - D_+(\tau) \psi'(q_0) \\ &= \{D_+(\tau) - D_+(\tau_0)\} \delta_L(q_0), \end{aligned} \qquad (30) $$

$$ \begin{aligned} C(q_0; \tau) &\equiv \frac{1}{6} \left. \frac{\partial^3 x}{\partial q^3} \right|_{q_0} = \frac{1}{6} D_+(\tau) \psi'''(q_0) \\ &= -\frac{1}{6} D_+(\tau) \delta''_L(q_0), \end{aligned} \qquad (31) $$

where the prime denotes the derivative with respect to $q$. Here, we used equations (24) and (27). The above expressions imply that the shell-crossing point slightly moves from $x = x(q_0; \tau_0)$ to $x = A(q_0; \tau)$ and the shell-crossing structure develops as an 'S' shape, as shown in Fig. 1.

We now introduce the new Lagrangian coordinate, $Q \equiv q - q_0$, for which the structure of the local density peak is symmetric. Since we expand the displacement field at third order, the equation $x(q, \tau) = x_0$ inside the shell-crossing region has three solutions. An example is shown in Fig. 1 for the equation $x(Q) = x(Q_*)$ with solutions $Q_- < Q_* < Q_+$. In general, the three-ordered roots of the equation $x(Q) = \tilde{x}$, when they exist, are related to each other
through the following expressions:

$$Q_* = \frac{1}{2} \left\{ -Q_{\pm} \pm \sqrt{3(\hat{Q}_c^2 - Q_{\pm}^2)} \right\}, \quad (32)$$

$$ \begin{aligned} Q_{\pm} &= \frac{1}{2} \left\{ -Q_* \mp \sqrt{3(\hat{Q}_c^2 - Q_*^2)} \right\} \\ &= \frac{1}{2} \left\{ -Q_{\mp} \pm \sqrt{3(\hat{Q}_c^2 - Q_{\mp}^2)} \right\}. \end{aligned} \quad (33) $$

Here, the Lagrangian extent $\hat{Q}_c$ of the boundary of the multivalued region, also related to the Lagrangian position $Q_c$ of the left caustic in Fig. 1 through $x(\hat{Q}_c) = x(Q_c)$, is expressed as

$$\hat{Q}_c = 2Q_c = \sqrt{\frac{4B}{3C}}. \quad (34)$$

Because of the time dependence of the coefficients (see equations 30 and 31), the position of the boundary gradually changes in time. The leading-order expression in powers of time ($\tau - \tau_0$) reads

$$\hat{Q}_c \approx \left\{ \frac{8}{\kappa(q_0, \tau_0)} \right\}^{1/2} (\tau - \tau_0)^{1/2} + \mathcal{O}((\tau - \tau_0)^{3/2}), \quad (35)$$

with

$$\kappa(q_0, \tau_0) \equiv \frac{-\delta''_L(q_0) D_+(\tau_0)}{\delta_L(q_0) \frac{dD_+(\tau_0)}{d\tau_0}}. \quad (36)$$

In other words, it also means that a given element of fluid of Lagrangian position $Q$ will enter the multivalued region at some time $\hat{\tau}_c$ defined by $\hat{Q}_c(\hat{\tau}_c) = |Q|$ and then subsequently coincide exactly with the Lagrangian position of the caustic at a time $\tau_c > \hat{\tau}_c$ with $Q_c(\tau_c) = |Q|$. By inverting the relations $|Q| = Q_c(\tau_c)$ and $|Q| = \hat{Q}_c(\hat{\tau}_c)$, we obtain

$$\hat{\tau}_c(Q) - \tau_0 \approx \frac{1}{8} \kappa(q_0, \tau_0) Q^2 \quad (37)$$

and

$$\tau_c(Q) - \tau_0 = \hat{\tau}_c(2Q) - \tau_0 \approx \frac{1}{2} \kappa(q_0, \tau_0) Q^2. \quad (38)$$

The quantities defined above thus play a crucial role to disentangle the single-stream regime from post-collapse dynamics, as well as the inner part of the 'S' shape of Fig. 1 from its outer part, the Lagrangian position of the caustic acting as a separator in the latter case. The expression of the force will indeed be different between the Zel'dovich single-flow regime $|Q| > \hat{Q}_c$ and the multivalued regime composed of an outer caustic region, $Q_c < |Q| < \hat{Q}_c$, and an inner caustic region, $|Q| < Q_c$.

One important remark for the subsequent calculations is that the boundaries coordinates $\hat{Q}_c$ and $Q_c$ are assumed to be small and that the backreaction to the position and velocity inside the multivalued region can be described perturbatively in a polynomial form of $Q = q - q_0$ [partly allowing the fractional power of $(\hat{Q}_c^2 - Q^2)$]. Thus, the resultant expressions are, rigorously speaking, only valid for a short period after collapse time, but we shall see in practice, as also shown by Colombi (2015) in the non-cosmological case, that they remain impressively accurate even up to next-crossing time (see Section 5), which will allow us to set up the framework for a powerful self-adaptive scheme.

## 3.2 Computing the force in the multivalued region

To derive the corrections to the motion, we first compute the force exerted on a mass element inside the multivalued region,

$-\hat{Q}_c \le Q \le \hat{Q}_c$ (see Fig. 1). Note again that the outer regions $x < x(-\hat{Q}_c)$ and $x > x(\hat{Q}_c)$ are described by the Zel'dovich solution.

The force in the multivalued region is computed using equation (20), dividing each integral of the right-hand side into three contributions:

$$\int_0^L dx \rightarrow \left( \int_0^{x(-\hat{Q}_c)} + \int_{x(-\hat{Q}_c)}^{x(\hat{Q}_c)} + \int_{x(\hat{Q}_c)}^L \right) dx. \quad (39)$$

Assuming that the collapse region, $|Q| \le \hat{Q}_c$, is small enough, the contributions to the integrals from each domain can be computed analytically, based on the geometrical setup in Fig. 1. The detailed calculations are presented in Appendix B. Summing up all the contributions given in equations (B4), (B5), (B8) and (B12), the force exerted on the mass element at $x = x(Q)$ inside the multivalued region becomes

$$F(x(Q; \tau)) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau) [\mathcal{J}(Q; q_0, \tau) + \mathcal{F}(q_0, \tau)] \quad (40)$$

with the functions $\mathcal{J}$ and $\mathcal{F}$ respectively defined by

$$\mathcal{J}(Q; q_0, \tau) = \begin{cases} \{1 + B(q_0; \tau)\} Q - C(q_0; \tau) Q^3 \\ - \operatorname{sgn}(Q) \sqrt{3(\hat{Q}_c^2 - Q^2)}; \\ Q_c < |Q| < \hat{Q}_c, \\ \\ \{-2 + B(q_0; \tau)\} Q - C(q_0; \tau) Q^3; \\ \\ |Q| < Q_c, \end{cases} \quad (41)$$

and

$$\mathcal{F}(q_0, \tau) = -\psi(q_0) D_+(\tau), \quad (42)$$

where the quantities A, B and C are defined by equations (29)–(31). Note that in deriving equation (40), we have assumed that the system follows Zel'dovich solution if $|Q| > \hat{Q}_c$. Since the resultant expressions are written in terms of the local quantities characterizing the density peak at position $q_0$ and the shell-crossing time $\tau_0$, equation (40) is in fact still applicable to other shell-crossing regions possibly appearing elsewhere in the region $|Q| > \hat{Q}_c$, allowing one to generalize the result to a smooth random initial field with multiple peaks. One issue, discussed later, is then to treat mergers, i.e. the case when for instance two multistream regions start overlapping.

## 3.3 Corrections to the Zel'dovich flow: basic post-collapse PT results

Given the explicit expression for the force in the multistream region and using the formal solution given by equations (17) and (18), we now compute corrections to the Zel'dovich flow that we write as follows:

$$\Delta u(Q; \tau, \hat{\tau}_c) = \int_{\hat{\tau}_c}^{\tau} d\tau' F(x(Q, \tau')), \quad (43)$$

$$\Delta x(Q; \tau, \hat{\tau}_c) = \int_{\hat{\tau}_c}^{\tau} d\tau' \Delta u(Q; \tau', \hat{\tau}_c). \quad (44)$$

We noticed in previous section that, depending on the Lagrangian position of interest, the expression for the force is different and we have to divide the domain of the integrals in equations (43) and (44) into several pieces:

(i) $\tau_0 \le \tau < \hat{\tau}_c(Q)$: the position $Q$ is located in the single-valued region (i.e. $|Q| > Q_c$) and the motion is still described by the
Zel'dovich solution. We have

$$x(Q; \tau) = x_{\text{Zel}}(Q; \tau) \equiv q + \psi(q) D_{+}(\tau), \quad (45)$$

$$u(Q; \tau) = u_{\text{Zel}}(Q; \tau) \equiv \psi(q) \frac{dD_{+}(\tau)}{d\tau}. \quad (46)$$

(ii) $\hat{\tau}_c(Q) \le \tau < \tau_c(Q)$: the position $Q$ is in the multivalued region and satisfies $Q_c < |Q| \le \hat{Q}_c$, i.e. lies in the outer part of the caustic. In addition to the Zel'dovich flow, the corrections arising from the multistream flow need to be added and we note them as follows:

$$x(Q; \tau) = x_{\text{Zel}}(Q; \hat{\tau}_c(Q)) + \Delta x_{\text{out}}(Q; \tau, \hat{\tau}_c(Q)), \quad (47)$$

$$u(Q; \tau) = u_{\text{Zel}}(Q; \hat{\tau}_c(Q)) + \Delta u_{\text{out}}(Q; \tau, \hat{\tau}_c(Q)). \quad (48)$$

In this region, which correspond to the tails of the $S$ shape in Fig. 1, the system is globally expanding in phase space and gaining energy, at variance with the central part, which on the contrary, contracts. This can be understood from computing the variation of energy of a typical test particle during a fraction of orbit. This phenomenon was studied in the non-cosmological case by Colombi (2015) and the results should not be fundamentally different in the cosmological case studied here.

(iii) $\tau_c(Q) \le \tau$: this corresponds to $|Q| \le Q_c$, i.e. the position $Q$ now lies in the inner part of the multistream region with respect to the caustic. Similarly to the above case, the backreaction to Zel'dovich flow needs to be computed, including at present both the multistream dynamics of the inner part and the incoming flow from the outer part. We may write

$$x(Q; \tau) = x_{\text{Zel}}(Q; \hat{\tau}_c(Q)) + \Delta x_{\text{in}}(Q; \tau, \hat{\tau}_c(Q)), \quad (49)$$

$$u(Q; \tau) = u_{\text{Zel}}(Q; \hat{\tau}_c(Q)) + \Delta u_{\text{in}}(Q; \tau, \hat{\tau}_c(Q)). \quad (50)$$

Note that the corrections $\Delta x_{\text{in}}$ and $\Delta u_{\text{in}}$ partly come from outer part contributions, $\Delta x_{\text{out}}(Q; \tau_c, \hat{\tau}_c)$ and $\Delta u_{\text{out}}(Q; \tau_c, \hat{\tau}_c)$. As mentioned in the previous point (ii), in this region, which corresponds to the central part of the $S$ on Fig. 1, the system is globally contracting and losing energy in favour of the tails of the $S$.

In what follows, we compute the backreaction to Zel'dovich flow and derive the expressions for $\Delta x$ and $\Delta u$ in each domain. The calculation for the corrected motion is rather straightforward but needs several steps. Here, we give a brief sketch of the calculation, deferring details to Appendix C. The final results are given in equations (52) and (54) for the outer part, equations (56) and (58) for the inner part, together with the coefficients in Tables 1 and 2. Note that although our expansion is rigorously valid only at third order in $Q$, higher order contributions will appear in the corrections to provide a continuous solution, up to fifth and seventh orders for the velocity and the position, respectively.

### 3.3.1 Velocity and position in the outer part: $Q_c < |Q| \le \hat{Q}_c$

In the outer part of the multivalued region, (ii), the correction to the velocity can be expressed as

$$\Delta u_{\text{out}}(Q; \tau, \hat{\tau}_c) = -\frac{3}{2} H_0^2 \Omega_{\text{m},0} \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' a(\tau') \times \{\mathcal{J}(Q; q_0, \tau') + \mathcal{F}(q_0, \tau')\}. \quad (51)$$

In the above, while the first integral is performed with the help of formulae in Appendix E, the second integral is computed exactly.

The whole derivation is provided in Appendix C1. The resultant expression can be summarized as

$$\begin{aligned} \Delta u_{\text{out}}(Q; \tau, \hat{\tau}_c) = & -\frac{3}{2} H_0^2 \Omega_{\text{m},0} a(\tau_0) [\tilde{\alpha}_1(\tau) Q + \tilde{\beta}_1(\tau) Q^3 \\ & + \tilde{\gamma}_1(\tau_0) \{\hat{Q}_c^2(\tau) - Q^2\}^{3/2} + \tilde{\delta}_1(\tau_0) Q^5 \\ & + \tilde{\epsilon}_1(\tau, \hat{\tau}_c), \end{aligned} \quad (52)$$

with the time-dependent coefficients given in Table 1. Note that the coefficient $\tilde{\epsilon}_1$ implicitly depends on the Lagrangian position $Q$ through $\hat{\tau}_c(Q) \simeq \tau_0 + (\kappa/8) Q^2$ but is not Taylor expanded with respect to $Q$, for simplicity. A fully analytical theory, in particular to predict the power spectrum of the projected density field, would in principle require such a Taylor expansion. Performing it should not change significantly the performances of post-collapse PT as presented in this article.

Once $\Delta u_{\text{out}}$ is obtained, the expression for the correction $\Delta x_{\text{out}}$ is evaluated by further integrating equation (52) over time:

$$\Delta x_{\text{out}}(Q; \tau, \hat{\tau}_c) = \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \Delta u_{\text{out}}(Q, \tau'). \quad (53)$$

The resultant expression becomes (see Appendix C2 for derivation):

$$\begin{aligned} \Delta x_{\text{out}}(Q; \tau, \hat{\tau}_c) = & -\frac{3}{2} H_0^2 \Omega_{\text{m},0} a(\tau_0) [\tilde{\alpha}_2(\tau) Q + \tilde{\beta}_2(\tau) Q^3 \\ & + \tilde{\gamma}_2(\tau_0) \{\hat{Q}_c^2(\tau) - Q^2\}^{5/2} + \tilde{\delta}_2(\tau) Q^5 \\ & + \tilde{\zeta}_2(\tau) Q^7] + \tilde{\epsilon}_2(\tau, \hat{\tau}_c), \end{aligned} \quad (54)$$

with the time-dependent coefficients presented in Table 1. Note again the dependence on $Q$ of parameter $\tilde{\epsilon}_2(\tau, \hat{\tau}_c)$ through $\hat{\tau}_c$ and the fact that we did not Taylor expand it in polynomials of $Q$, as it would be required for a fully analytical theory.

### 3.3.2 Velocity and position in inner part: $|Q| \le Q_c$

Let us consider now the inner part of the multivalued region (iii). For the correction to the velocity, $\Delta u_{\text{in}}$, the expression to be evaluated is divided into three contributions:

$$\begin{aligned} \Delta u_{\text{in}}(Q; \tau, \hat{\tau}_c) = & -\frac{3}{2} H_0^2 \Omega_{\text{m},0} \left[ \int_{\hat{\tau}_c(Q)}^{\tau_c(Q)} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \right. \\ & \left. + \int_{\tau_c(Q)}^{\tau} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \right. \\ & \left. + \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' a(\tau') \mathcal{F}(q_0, \tau') \right]. \end{aligned} \quad (55)$$

Evaluating each contribution in Appendix C3, we obtain

$$\begin{aligned} \Delta u_{\text{in}}(Q; \tau, \hat{\tau}_c) = & -\frac{3}{2} H_0^2 \Omega_{\text{m},0} a(\tau_0) [\tilde{\alpha}_3(\tau) Q + \tilde{\beta}_3(\tau) Q^3 \\ & + \tilde{\delta}_3(\tau_0) Q^5] + \tilde{\epsilon}_1(\tau, \hat{\tau}_c), \end{aligned} \quad (56)$$

with the time-dependent coefficients given in Table 2. Note that the coefficient $\tilde{\epsilon}_1$ in the above equation is the same one as in equation (52).
<table><caption>**Table 1.** Coefficients for the corrections to Zel'dovich flow in the outer part of the collapse region, $Q_c < |Q| \le \tilde{Q}_c$.</caption><thead><tr><th>Coefficients</th><th>Velocity Δu<sub>out</sub> [equation 52]</th><th>Position Δx<sub>out</sub> [equation 54]</th></tr></thead><tbody><tr><td>α̃</td><td>T(≡τ − τ<sub>0</sub>)</td><td>T<sup>2</sup> / 2</td></tr><tr><td>β̃</td><td>T (δ<sub>L</sub>''(q<sub>0</sub>) / 6) D<sub>+</sub>(τ<sub>0</sub>) − κ / 8</td><td>−κ / 8 T + δ<sub>L</sub>''(q<sub>0</sub>) / 12 D<sub>+</sub>(τ<sub>0</sub>) T<sup>2</sup></td></tr><tr><td>γ̃</td><td>−sgn(Q) κ / (4√3)</td><td>−sgn(Q) κ<sup>2</sup> / (80√3)</td></tr><tr><td>δ̃</td><td>−(δ<sub>L</sub>''(q<sub>0</sub>) / 48) κ D<sub>+</sub>(τ<sub>0</sub>)</td><td>(1/2) (κ/8)<sup>2</sup> − (κ/8) δ<sub>L</sub>''(q<sub>0</sub>) / 6 D<sub>+</sub>(τ<sub>0</sub>) T</td></tr><tr><td>ζ̃</td><td>—</td><td>(1/2) (κ/8)<sup>2</sup> δ<sub>L</sub>''(q<sub>0</sub>) / 6 D<sub>+</sub>(τ<sub>0</sub>)</td></tr><tr><td>ε̃</td><td>ψ(q<sub>0</sub>) [dD<sub>+</sub>(τ') / dτ']<sup>τ</sup> |<sub>τ̃c(Q)</sub></td><td>ψ(q<sub>0</sub>) { D<sub>+</sub>(τ) − D<sub>+</sub>(τ̃c(Q)) − dD<sub>+</sub>/dτ |<sub>τ̃c(Q)</sub> (τ − τ̃c(Q)) }</td></tr></tbody></table>

<table><caption>**Table 2.** Coefficients for the corrections to Zel'dovich flow in the inner part of the collapse region, $|Q| \le Q_c$.</caption><thead><tr><th>Coefficients</th><th>Velocity Δu<sub>in</sub> [equation 56]</th><th>Position Δx<sub>in</sub> [equation 58]</th></tr></thead><tbody><tr><td>α̃</td><td>−2T</td><td>−T<sup>2</sup></td></tr><tr><td>β̃</td><td>T (δ<sub>L</sub>''(q<sub>0</sub>) / 6) D<sub>+</sub>(τ<sub>0</sub>) + κ/8</td><td>(5κ/8) T + δ<sub>L</sub>''(q<sub>0</sub>) / 12 D<sub>+</sub>(τ<sub>0</sub>) T<sup>2</sup></td></tr><tr><td>δ̃</td><td>− (κ/8) δ<sub>L</sub>''(q<sub>0</sub>) / 6 D<sub>+</sub>(τ<sub>0</sub>)</td><td>− (κ/4)<sup>2</sup> (67/40) − (κ/8) δ<sub>D</sub>''(q<sub>0</sub>) / 6 D<sub>+</sub>(τ<sub>0</sub>) T</td></tr><tr><td>ζ̃</td><td>—</td><td>(1/2) (κ/8)<sup>2</sup> δ<sub>L</sub>''(q<sub>0</sub>) / 6 D<sub>+</sub>(τ<sub>0</sub>)</td></tr><tr><td>ε̃</td><td>ψ(q<sub>0</sub>) [dD<sub>+</sub>(τ') / dτ']<sup>τ</sup> |<sub>τ̃c(Q)</sub></td><td>ψ(q<sub>0</sub>) { D<sub>+</sub>(τ) − D<sub>+</sub>(τ̃c(Q)) − dD<sub>+</sub>/dτ |<sub>τ̃c(Q)</sub> (τ − τ̃c(Q)) }</td></tr></tbody></table>

Similarly, the correction to the position, $\Delta x_{in}$, is divided into three pieces:

$$ \Delta x_{in}(Q; \tau, \hat{\tau}_c) = \left( \int_{\hat{\tau}_c(Q)}^{\tau_c(Q)} d\tau' + \int_{\tau_c(Q)}^{\tau} d\tau' \right) \Delta \bar{u}(Q; q_0, \tau') + \tilde{\epsilon}_2(\tau, \hat{\tau}_c), \quad (57) $$

where the function $\Delta \tilde{u}$ represents the correction of the velocity given in either equation (52) or equation (56), but without the contribution corresponding to the Zel'dovich motion, i.e. the last term in their expressions. The term $\tilde{\epsilon}_2$ is the same quantity as in equation (54) (see Table 1).

Based on the calculation presented in Appendix C4, we obtain the expression for $\Delta x_{in}$:

$$ \Delta x_{in}(Q; \tau, \hat{\tau}_c) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) [\tilde{\alpha}_4(\tau) Q + \tilde{\beta}_4(\tau) Q^3 + \tilde{\delta}_4(\tau) Q^5 + \tilde{\zeta}_4(\tau) Q^7] + \tilde{\epsilon}_2(\tau, \hat{\tau}_c), \quad (58) $$

with the time-dependent coefficients given in Table 2.

## 3.4 Higher order corrections

The perturbative description of post-collapse dynamics introduces two kinds of approximations around the singularities (shell-crossing points). First, we assumed the position and velocity to be described by third-order polynomials in $Q$ in the multistream regions. To preserve continuity between the multistream and the single-stream regime, we added higher order terms in $Q$ to the solutions for $x$ and

$v$, but strictly speaking, these solutions are valid locally only at third order in $Q$ around the singularities. Secondly, these solutions are correct only very shortly after collapse time $\tau_0$ up to second order in $\tau - \tau_0$. Our post-collapse dynamics can thus be improved at a twofold level:

(i) Improvement at the spatial level: we could try, starting from initial conditions, to use a description of the $S$ shape of Fig. 1 at higher order in $Q$, e.g. fourth or fifth instead of third order. In this latter case, deriving the solution of equation $x(Q) = x(Q_*)$ comes down to find the roots of a third- or fourth-order polynomial instead of a second-order one, which makes, in addition to the symmetry breaking of the $S$ shape, the calculations much more cumbersome, even though a series expansion of the roots at fourth or fifth order in $Q$ remains tractable and might bring interesting improvements on the description of the tails of the $S$ shape.

(ii) Higher order in time: we can introduce higher order corrections in time to follow more accurately the evolution of the central part of the singularity after collapse. Note that after integration of the equations of motion, such corrections bring out higher order terms in $Q$. Strictly speaking, the right way to proceed would consist of performing 'post-post-collapse' dynamics, i.e. of computing a counter term to the Zel'dovich+leading-order post-collapse solution. This procedure is cumbersome, and seems to make sense only if performed simultaneously with a description at higher order in $Q$, in particular because the tails of the $S$ structure in Fig. 1 also
influence the dynamics of the central part of the system (but this remains to be proved).

Going beyond leading order in post-collapse dynamics thus seems rather challenging even in the simple one-dimensional case considered in this article.

However, we can still try to investigate some potential improvements of the theory while maintaining the complexity of the calculations at an acceptable level. For this purpose, we consider here mainly two approaches. The first one just consists in incorporating next-to-leading-order terms into the expressions for $\hat{Q}_c$ and $\hat{\tau}_c$ (and $Q_c$ and $\tau_c$). These quantities, which determine the boundary between the single- and the multistream regimes, are indeed crucial because they correspond to the main scale/time-scale of the system. The second approach consists, in addition, of improving the evaluation of the time integrals in equations (43) and (44). In Section 3.3, when evaluating these time integrals, a part of the integrands is Taylor-expanded and the results are presented at leading order in time. It is thus possible to include higher order corrections from the Taylor-expanded integrands.

Again, one has to be aware that the two corrections we propose might not improve the results, since a correct higher order time approach, as discussed in point (ii) above, would actually require to compute counter terms to the standard post-collapse result, which we do not do here. In fact, our higher order corrections might even worsen the results since the corrections brought by such additional terms might have opposite sign compared to the iterated counter-term approach.

The derivation of the higher order corrections we propose here and summarize just below is rather straightforward compared to its rigorous alternatives (i) and (ii) above, but still cumbersome. The resultant expressions now include higher order polynomials of $Q$ as well as new contributions, which are given in Appendix D. To be more specific, we consider the following cases:

(i) Higher order for the expressions of critical time $\hat{\tau}_c$ and position $\hat{Q}_c$ (hc): based on the basic results in Section 3.3, we include the contributions arising from the next-to-leading-order expressions for $\hat{Q}_c$ and $\hat{\tau}_c$ (as well as $Q_c$ and $\tau_c$). The calculation to derive higher order terms is almost the same as presented in Section 3.3, but we replace the leading-order expressions for $\hat{Q}_c$ and $\hat{\tau}_c$ in equations (35) and (37) with those including the next-to-leading order, given in equations (D1) and (D3). The higher order results, denoted by $\Delta x^{(hc)}$ and $\Delta u^{(hc)}$, correspond to the expressions summarized in Appendix D, but where the terms involving the quantities $g(\tau_0)$ or $H(\tau_0)$ are all set to zero.

(ii) Higher order in time (ho): in addition to the above treatment, we also include higher order corrections to the time integrals, based on the formulae given in Appendix E. The resultant expressions, $\Delta x^{(ho)}$ and $\Delta u^{(ho)}$, are summarized in Appendix D.

Finally, we also consider a semi-analytic approach in order to smoothly connect the inner part of the multistream region to the Zel'dovich solution:

(iii) Third-order spline (spl): basic post-collapse PT is applied to the inner part, $|Q| \le Q_c$, but we take the higher order expression for $Q_c$. In the outer part, $Q_c < |Q| < \hat{Q}_c$, on the other hand, a third-order spline interpolation is used to smoothly connect the inner post-collapse region to the Zel'dovich solution, with again, the higher order expression for $\hat{Q}_c$.

## 4 IMPROVED TREATMENT WITH ADAPTIVE SMOOTHING

Post-collapse PT as described in Section 3 provides a way to describe the dynamics around an initial density peak shortly beyond collapse time. In practice, the quality of the description is expected to be good until next-crossing time, at least this is what is known in the non-cosmological case (Colombi 2015). Using the leading-order expression given in equation (58), the duration of time spent between collapse and next shell-crossing can be estimated from equation (49),

$$ \tau_{\text{cross}}(q_0) = \frac{\delta_L(q_0) \frac{dD_+(\tau_0)}{d\tau_0}}{\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0)}, \qquad (59) $$

which is derived by solving the equation $(\partial x/\partial q)_{q=q_0} = 0$.

Beyond time $\tau = \tau_0 + \tau_{\text{cross}}$, post-collapse PT is expected to fail, similarly as Zel'dovich approximation fails beyond $\tau = \tau_0$. The only way to fix this problem at the fine grain level would be to use an iterative procedure similar to what was proposed by Colombi (2015) to reconstruct the spiral structure built up during the course of the dynamics. In addition, post-collapse PT assumes the initial density distribution to be sufficiently smooth. While this is certainly expected at very small scales in the standard cold dark matter (CDM) scenario, the issue of mergers arises. Indeed, post-collapse PT is valid only if the forming halo is isolated. If there are two adjacent collapsing regions, they might, at some point, overlap in Lagrangian space, hence merge with each other in Eulerian space. After the merger, a separate description of each halo becomes inadequate. In this case, the prediction of post-collapse PT, which is in principle able to handle multiple but non-overlapping structures in Lagrangian space, will breakdown much earlier than expected.

One simple way to fix both limitation in time and mergers issues, already applied in the past to the Zel'dovich approximation (see, e.g. Coles, Melott & Shandarin 1993; Sahni & Coles 1995, and references therein), consists in trying to describe the dynamics at the coarse-grained level by introducing a smoothing procedure in Lagrangian space, so that post-collapse PT remains applicable, in particular that the condition $\tau \lesssim \tau_0 + \tau_{\text{cross}}$ remains valid for the system with coarse-grained initial conditions. This means that single or composite haloes as well as mergers are locally summarized by a simplified $S$ shape (Fig. 1) corresponding to a single peak of the coarse-grained initial density field and such that $\tau \lesssim \tau_0 + \tau_{\text{cross}}$, while retaining the correct large-scale dynamics. The hope is that this simplified structure still captures the most important features of the non-linear system, e.g. its mass and its typical size.

Of course, the coarsening scale should depend on environment, hence on Lagrangian position $Q$, so for best results, adaptive smoothing of initial conditions has to be introduced. In this sense our method is very similar to the 'peak-patch' treatment by Bond & Myers (1996) (see also Monaco et al. 2002), but we implemented it in the post-collapse PT framework. To be precise, the practical set up of post-collapse dynamics with adaptive smoothing is summarized as follows:

(i) The initial density field is smoothed at various scales, employing the sharp-k filter function in Fourier space:

$$ W(k; k_{\text{cut}}) = \Theta(k_{\text{cut}} - k) \qquad (60) $$

with $\Theta(x)$ being the Heaviside step function. The cut-off wavenumber $k_{\text{cut}}$ ranges from the fundamental mode $k_{\text{min}} = 2\pi/L$ to a certain value $k_{\text{max}} = 2\pi m_{\text{max}}/L$, where $m_{\text{max}}$ is an integer corresponding to the effective resolution of the simulation, as discussed more in
details in Section 5.1. In our investigations, we covered all the possible values of $k_{\text{cut}}$, i.e. $k_{\text{cut}} = 2\pi m/L$ with $m = 1, 2, \dots, m_{\text{max}}$, but this is costly and probably not required, in practice.

(ii) Zel'dovich solutions are constructed from the smoothed density fields. Then, we identify the critical points $q_0$ in Lagrangian space that have undergone shell-crossing at time $\tau_0 < \tau$, where $\tau$ is the time of interest for the analyses. For each critical point, various quantities needed for perturbative calculations, such as $\delta_L(q_0)$ and $\delta''_L(q_0)$, are also computed.

(iii) Starting from the largest smoothing scale (i.e. $k_{\text{min}}$), we compute the next-crossing time $\tau_{\text{cross}}$ associated to each critical point. For the points satisfying $\tau \ge \tau_0 + f_{\text{cross}} \tau_{\text{cross}}$ with $f_{\text{cross}} \le 1$ a constant parameter, the post-collapse correction to the Zel'dovich flow is obtained in the interval $I_0 = [q_0 - \hat{Q}_c(q_0), q_0 + \hat{Q}_c(q_0)]$, which is labelled.

(iv) The procedure (iii) is repeated, but with a smaller smoothing scale, i.e. larger $k_{\text{cut}}$, and only for peaks lying in the unlabelled regions, until the cut-off wavenumber reaches $k_{\text{max}}$. Note that if two Lagrangian regions overlap, e.g. $I_0(k_{\text{cut},1})$ and $I_0(k_{\text{cut},2})$, corresponding to two different values of the smoothing scale, $k_{\text{cut},1} < k_{\text{cut},2}$, the interval $I_0(k_{\text{cut},2})$ corresponding to the smallest smoothing scale, i.e. largest wavenumber $k_{\text{cut},2}$, prevails. This is the way we solve the cloud-in-cloud problem, but not the only possible choice.

(v) At the end of the iterative process, post-collapse dynamics is applied as well for the rest of the critical points corresponding to the highest resolution scale $2\pi/k_{\text{max}}$ that still require treatment, i.e. for those with $\tau_0 < \tau < \tau_0 + f_{\text{cross}} \tau_{\text{cross}}$. The procedure is the same as described in point (iv) above, i.e. only critical points which do not fall in a labelled region are treated, and if it is the case, then their Lagrangian interval $I_0$ prevails for applying post-collapse dynamics. For the rest of the unlabelled regions, standard Zel'dovich dynamics is applied and the mapping uses the initial density corresponding, of course, to the highest resolution scale.

Note that the approximation of the phase-space density distribution constructed this way is generally discontinuous, because we collect perturbative solutions for the displacement field in phase space coming from different coarse-graining scales without imposing smoothness at the transition between these solutions. In this respect, the procedure given above is not entirely consistent and remains to be improved in the future. Nevertheless, the discontinuities only affect the results at small scales and in regions where the projected density is small. Indeed, the discontinuities correspond to transitions between tails of the S shape representing a halo and a non-collapsed region or transitions between tails of two S shapes in the case of a merger. As illustrated below through numerical examples, these regions of phase space have small density contrast and thus contribute weakly to second- or higher order statistics, hence to the power spectrum which is the statistics under focus in this work. With an appropriate choice of $f_{\text{cross}}$, we shall show that the adaptive smoothing procedure provides a substantial improvement on the post-collapse prediction for the overall phase-space structure of the system and the power spectrum of the projected density.

## 5 COMPARISON WITH CONTROLLED N-BODY EXPERIMENTS

We are now in position to compare the predictions of post-collapse PT with simulation results. In Section 5.1, we briefly describe the one-dimensional cosmological N-body code we wrote to perform the numerical analyses of this article and discuss various important parameters settings. We then present the results of the comparison between theory and numerical experiments in various cases with cold initial conditions: formation of an isolated single halo (Section 5.2), merging of two neighbouring haloes (Section 5.3), as well as random initial conditions with CDM-like and power-law power spectra (Sections 5.4 and 5.5). In the later case, we shall discuss more in details highly non-linear dynamics, in particular predictions from stable clustering, the expected outcome on the power spectrum from haloes with power-law density profiles, as well as the statistical properties of the phase-space sheet resulting from a caustic distribution, which will be reminiscent of a similar analysis performed previously by Yano et al. (2004).

### 5.1 Simulations

To solve numerically Vlasov-Poisson equations in one dimension, we wrote a FORTRAN 90 particle-mesh code, Vlafroid. In this code, publicly available through the following web page, www.vlasix.org, Poisson equation is solved on a periodic grid of fixed resolution $N_{\text{grid}}$ using fast Fourier transform, after projecting the $N_{\text{particles}}$ particles on the grid using simple cloud-in-cell interpolation (see, e.g. Hockney & Eastwood 1988).<fnref n="1" /> The particles positions and velocities are updated at each time-step using a standard predictor-corrector integration scheme. After the predictor step,

$$x(\tau + \Delta\tau/2) = x(\tau) + v(\tau)\Delta\tau/2, \quad (61)$$

Poisson equation is solved and the force is estimated on the grid by direct derivation in Fourier space. The acceleration $\gamma(t + \Delta t/2)$ of each particle is computed by linear interpolation of the force. This is followed by the corrector step:

$$v(\tau + \Delta\tau) = v(\tau) + \gamma(\tau + \Delta\tau/2)\Delta\tau, \quad (62)$$

$$x(\tau + \Delta\tau) = x(\tau + \Delta\tau/2) + v(\tau + \Delta\tau)\Delta\tau/2. \quad (63)$$

The value of the time-step $\Delta\tau$ is slowly varying with time and is bounded by three constraints:

(i) a limit on the relative variations of the expansion factor during a time-step to have a good description of the linear growing modes at early stages of the simulation,

$$\Delta\tau \le C_{\text{dloga}} \left( \frac{d \log a}{d \tau} \right)^{-1}, \quad (64)$$

(ii) a dynamical condition which states that the time-step should be a small fraction of the smallest harmonic oscillator time-scale (see, e.g. Colombi & Touma 2014),

$$\Delta\tau \le \frac{C_{\text{dyn}}}{\sqrt{\frac{3}{2} a \Omega_{m,0} \rho_{\text{max}} H_0^2}}, \quad (65)$$

where $\rho_{\text{max}}$ is the maximum value of the projected density (normalized to be unity in average),

<footnote n="1">¹ Note that in one dimension, one can use equation (20) to estimate the force directly as a quantity proportional to the difference between the total mass on the right of the particle of interest and the total mass on the left of it, along with a centre of mass correction to account for periodicity. However, to reduce discreteness effects, we aim to compute a force smoothed with a softening kernel of fixed length (in the particle-mesh code, the size of a cell). Performing such a softening on the exact force is possible but costly, this is why we use the particle-mesh technique, which is less ‘clean’ than the exact calculation in the sense that softening is performed only approximately on a per-cell basis, but more efficient.</footnote>
(iii) the standard Courant-Friedrichs-Lewy (CFL) condition

$$ \Delta\tau \le C_{\text{CFL}} \frac{L}{N_{\text{grid}}v_{\text{max}}}, \qquad (66) $$

with $v_{\text{max}}$ the maximum magnitude of the velocity, stating that a particle should not travel more than a fraction of a grid cell size during a time-step.

For all the runs performed in this article, our practical and safe choice of the time-step parameters was:

$$ C_{\text{dloga}} = 0.1, \quad C_{\text{dyn}} = 0.01, \quad C_{\text{CFL}} = 0.25. \qquad (67) $$

The evolution of single or composite objects as studied in Sections 5.2 and 5.3 does not pose any problem from the numerical point of view. Indeed, we consider here only a small number of dynamical times, so we do not expect significant defects due to the discrete nature of the particle distribution (for instance, such as described in Melott et al. 1997). The case of a random field, as studied in Section 5.4 for CDM and in Section 5.5 for power-law power spectra, is less trivial. Indeed, one needs accurate numerical modelling of the dynamics of each single structure, which has to be sampled with many particles and be resolved with a sufficient number of grid elements. In particular, to avoid contamination by the discrete nature of the particle distribution, we require, in average, many particles per grid element:

$$ N_{\text{particles}} \gg N_{\text{grid}}. \qquad (68) $$

To avoid the formation of small structures that would be insufficiently resolved by the computational grid, there must be a coarse-graining scale below which the initial conditions are smooth. This also facilitates, of course, the calculation of (post-collapse) PT predictions. Here, we impose a cut-off on the power spectrum of initial conditions using the sharp-k filter introduced in equation (60), with

$$ k_{\text{max}} \ll N_{\text{grid}} \frac{2\pi}{L}. \qquad (69) $$

In practice, a choice which works well and adopted in this paper is

$$ N_{\text{particles}} = 10 N_{\text{grid}} = 100 k_{\text{max}} \frac{L}{2\pi}, \qquad (70) $$

which means that $k_{\text{max}}$ is five time smaller than the Nyquist frequency of the grid. Since the cut-off scale on the power spectrum can affect the dynamics compared to the true, unsmooth system, we shall perform some convergence studies to determine the available dynamic range in which the results are not influenced by $k_{\text{max}}$.

## 5.2 Formation of a single structure

Let us first present the simplest case of the formation of a single structure seeded with the following sinusoidal density perturbation:

$$ \delta_L(x) = -A \cos \left( \frac{2\pi}{L} x \right). \qquad (71) $$

We adopt the Einstein-de Sitter cosmology ($\Omega_{m,0} = 1$, $\Omega_{\Lambda} = 0$ with $h \equiv H_0/100 = 0.7$) and the box size of the simulation is arbitrarily set to $L = 1$. The simulation was started at expansion factor $a = 0.01$ with initial amplitude $A = 0.1$ and was run using ($N_{\text{particles}}$, $N_{\text{grid}}$) = (10, 000; 1, 000). PT analyses were performed on the fly as part of routines of the code, when needed.

Fig. 2 shows snapshots of the phase-space structure (upper inserts) and projected density profiles (lower inserts) obtained from the simulation (red) and post-collapse PT (blue) based on the basic prescription described in Section 3.3. Also, for reference, the

Zel'dovich solution is plotted as a green dashed line on each panel. Although the three curves are identical at early time, differences show up, as expected, after collapse, which occurs at $a \simeq 0.1$. A closer look at the central part of the system reveals however that post-collapse PT captures the trends of the N-body simulation quite well and succeeds to reproduce the structure of the inner part of the system until next-crossing time, which occurs at $a \approx 0.25$. This is in marked contrast with Zel'dovich solution, which, without the backreaction contribution present in the post-collapse correction, overshoots the displacement of the phase-space sheet after first crossing time, in particular of the caustics. After next-crossing time, $a \gtrsim 0.25$, the phase-space structure given by post-collapse PT goes away from the N-body result. Nevertheless, the predicted projected density profile still provides a reasonably good approximation of the evolved density structure, which is clearly not the case for the Zel'dovich solution.

Fig. 3 presents results of post-collapse PT when including the higher order corrections discussed in Section 3.4. While the three prescriptions examined here, namely *hc*, *ho* and *spl*, account for higher order terms in different manners, they are all based on the force expression in equation (40), which is valid only shortly after collapse. Indeed, as discussed in Section 3.4, our implementation of higher order corrections is not performed in a fully rigorous way, that would require to compute a backreaction to post-collapse PT, similarly as post-collapse PT introduces a backreaction correction to Zel'dovich solution. It also does not propose an iteration procedure that would allow us, similarly as in Colombi (2015), to follow the system behind next-crossing time. Hence, all of the treatments breakdown at the end of the simulation. Nevertheless, taking the higher order corrections to $\hat{Q}_c$ and $\hat{\tau}_c$ into account improves the prediction of the outer boundary of the multivalued region, and the phase-space description becomes visually better than that of the basic PT prescription in Fig. 2. Even after next-crossing time ($a \gtrsim 0.25$), the predicted phase-space structure still remains consistent, although the size of the halo is prone to be overestimated, leading to an extended density profile, the best result being obtained with the *spl* prescription. Still, one can notice that if only the projected density is under consideration, basic post-collapse PT does better than any of its higher order counterpart, except for *spl* which, not surprisingly, provides nearly equivalent results, since it is exactly the same as post-collapse PT in the centre of the system. As far as statistics depending on the projected density are concerned, such as the power spectrum that will be studied later, we thus expect from the analyses of this single halo that the simplest prescription of post-collapse PT is going to perform best.

## 5.3 Merger

We examine now the second simplest case of two merging structures in Einstein-de Sitter cosmology ($\Omega_{m,0} = 1$, $\Omega_{\Lambda} = 0$ and $h \equiv H_0/100 = 0.7$). The initial density contrast is given by

$$ \delta_L(x) = A \left[ \exp \left\{ - \left( \frac{x - x_1}{\sigma} \right)^2 \right\} + \exp \left\{ - \left( \frac{x - x_2}{\sigma} \right)^2 \right\} - c \right]. \qquad (72) $$

The constant $c$ is determined so as to satisfy the condition $\int_0^L dx \delta_L(x) = 0$. With $x_1 = 0.35L$, $x_2 = 0.65L$ and $\sigma = 0.07L$, we have $c = 0.248L$. The overall amplitude is set to $A = 0.3$ and the simulation, of box size unity ($L = 1$), was started at $a = 0.01$ and was run using ($N_{\text{particles}}$, $N_{\text{grid}}$) = (10, 000; 1, 000).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Snapshots of the phase-space structure (upper insert of each panel) and projected density profile (lower insert of each panel) of a single halo at different times in an Einstein-de Sitter universe. The initial conditions represented in upper left panel correspond to the projected initial density contrast given by equation (71). On each panel, results of *N*-body simulations are depicted as a red curve, the Zel'dovich solution as a green dotted line and the basic post-collapse PT prediction as a blue solid line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as in Fig. 2, but variants of the post-collapse PT calculation that include higher order corrections are compared with the *N*-body simulation, still in red: higher order for critical times (*hc*, cyan dot-dashed), higher order for critical times and integrands of the equations of motion (*ho*, black dotted), and higher order for critical times combined with spline interpolation to connect inner part with Zel'dovich solution (*spl*, dashed magenta).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Phase-space structure (upper inserts in each panel) and density profile (lower insert in each panel) of merging clusters in an Einstein-de Sitter universe. For the initial density contrast given in equation (72), results of N-body simulations are depicted as red lines, the Zel'dovich solution is shown as green dotted lines and the blue solid lines correspond to the prediction of post-collapse PT. Upper and lower panels, respectively, show the results without and with adaptive smoothing.</figcaption></figure>

The results are shown in Fig. 4 for basic post-collapse PT and Fig. 5 for the higher order variants of post-collapse PT. Here, we focus on the merger phase which takes place after collapse of each individual structure. Figs 4 and 5 demonstrate how the adaptive smoothing introduced in Section 4 improves PT prescriptions to describe the overall phase-space structure after the merger. In each figure, the upper and lower panels show the results without and with adaptive smoothing, respectively. The free parameter controlling adaptive smoothing, $f_{cross}$, is set here to 1 for post-collapse PT and 0.5 for the Zel'dovich solution.

Without adaptive smoothing, as soon as the two clusters cross each other ($a = 0.2$), both post-collapse PT and Zel'dovich solution start failing to describe the real dynamics in the N-body simulation. Indeed, while post-collapse PT tries to provide a local perturbative correction for each single cluster, it is clearly unable to account for non-linear couplings taking place at larger scales, where it simply follows Zel'dovich motion. As a result, the location of multivalued regions predicted by post-collapse PT, as well as their shape, largely deviates from the correct one.

On the other hand, implementing adaptive smoothing in the PT calculations significantly improves the predicted phase-space structure, starting already from the beginning of the merger phase ($a = 0.2$). After merging, dynamics of the central part of the system is now described by the displacement field derived from the smoothed initial density field. By construction, it summarizes the composite structure made of two sub-haloes into a single cluster. While this is a rough approximation of the real dynamics, the description of the outer part of the system is substantially

improved as well as the prediction for projected density profiles. Introducing both adaptive smoothing and higher order corrections to post-collapse PT further improves the results from the visual point of view (Fig. 5), although the predicted size of the halo tends to be slightly overestimated. As already noticed in Section 5.2, the best results are obtained for *spl*, which also gives, after adaptive smoothing, comparable projected density profiles to basic post-collapse PT.

The above results demonstrate that post-collapse PT with adaptive smoothing is effective in capturing the main features of phase-space structures. Strictly speaking, it does not give an accurate prescription for the fine structure of high-density regions, but it provides a way of regularizing or mitigating the impact of highly non-linear dynamics, keeping the location and size of haloes reasonably accurate. As we will see next, adaptive smoothing brings a drastic improvement on the prediction of power spectra for systems with random initial conditions. Furthermore, the introduction of adaptive smoothing makes the analytic calculations much less sensitive to the choice of the small-scale cut-off in the initial conditions so that PT predictions become more robust. In these respects, the determination of the extension of multistreaming regions for a given smoothing scale of initial conditions, as performed in step (iii) of Section 4, is essential, making the choice of $f_{cross}$ crucial. Concerning the particular merger experiment studied here, $f_{cross} = 1$ and 0.5 seem to be the best choices for tuning post-collapse PT and Zel'dovich solutions, respectively, and we shall adopt these values in subsequent section. We shall however discuss again about the choice of $f_{cross}$ later.
<figure><img src="image_6.png" /><figcaption>Figure 5. Same as in Fig. 4, but the simulation results are compared with improved treatments of post-collapse PT without (upper panels) and with (lower panels) adaptive smoothing. The coding for line types is the same as in Fig. 3.</figcaption></figure>

## 5.4 Random initial conditions: CDM-like spectrum

Let us now consider random CDM initial conditions. Although there is no realistic setup in 1D, a way of mimicking as well as possible the 3D case may be to consider initial conditions given by a random Gaussian field with the following power spectrum (see, e.g. McQuinn & White 2016):

$$P_{1D}(k) = \frac{k^2}{2\pi} P_{3D}(k), \qquad (73)$$

with $P_{3D}$ being the 3D matter power spectrum of the initial density fluctuations. Here, we use for $P_{3D}$ the transfer function of Eisenstein & Hu (1998) who wrote the corresponding FORTRAN module that we included in Vlafroid. We set the cosmological parameters to those of the base ΛCDM model determined by Planck (Planck Collaboration XIII 2016): total matter density parameter $\Omega_{m,0} = 0.3121$, cosmological constant $\Omega_{\Lambda} = 0.6879$, baryonic matter density parameter $\Omega_{b} = 0.04884$, Hubble constant $H_{0} = 67.51$ km s$^{-1}$ Mpc$^{-1}$, effective slope $n_{s} = 0.9653$ and root mean square deviation $\sigma_{8} = 0.815$.<fnref n="2" /> The simulations were performed with a box size $L = 1000$ Mpc and started at initial redshift $z_{i} = 99$. Convergence of simulations results was tested by varying the

number of particles $N_{particle}$, the resolution $N_{grid}$ of the grid used to compute the force and the cut-off wavenumber of the initial power spectrum, $k_{max}$, introduced in Section 5.1.<fnref n="3" /> Here, we present the results of the analyses for $N_{particle} = 200000$, $N_{grid} = 20000$ and $k_{max} = 2000(2\pi/L) = 12.6$ Mpc$^{-1}$. Our convergence study shows, unless specified otherwise below, that the results can be trusted for $k \lesssim 10$ Mpc$^{-1}$, i.e. in the available dynamic range shown in left-hand panels of Figs 6 and 7. To have sufficient statistics for measuring power spectra, we ran 50 simulations with different random seeds.

Figs 6 and 7 show the measured power spectra at different times (left-hand panel of each figure) as well as snapshots of phase space for a given realization (right-hand panel). Right-hand panel of Fig. 6 is supplemented with Fig. 8, which displays in approximately the same scale interval, Eulerian positions and velocities as functions of Lagrangian coordinates. One can first notice that, in contrast to the 3D case, the amplitude of the power spectrum at small scales is not strongly enhanced in 1D, but damped. For the CDM-like initial conditions considered here, function $kP_{1D}(k)$ becomes nearly flat in the highly non-linear regime, which is probably a consequence of stable clustering as discussed further in details in next section, where the case of scale-free initial conditions is considered.

$$\sigma_{8}^{2} = \int \frac{d^{3}k}{(2\pi)^{3}} P_{3D}(k)\{W_{\text{th}}(kR_{8})\}^{2} \qquad (74)$$

$$= \int_{0}^{\infty} \frac{dk}{\pi} P_{1D}(k)\{W_{\text{th}}(kR_{8})\}^{2} \qquad (75)$$

with $R_{8} = 8 h^{-1} \text{ Mpc}$, $h = H_{0}/100$ and $W_{\text{th}}(x) = 3\{\sin(x) - x \cos(x)\}/x^{3}$.

<footnote n="2">² To be precise, the power spectrum amplitude is normalized by $\sigma_{8}$ through</footnote>

<footnote n="3">³ Namely, we performed three sets of 50 simulations for $(N_{particle}, N_{grid}, k_{max}L/2\pi) = (4 \times 10^{5}, 4 \times 10^{4}, 4 \times 10^{3})$, $(2 \times 10^{5}, 2 \times 10^{4}, 2 \times 10^{3})$ and $(10^{5}, 10^{4}, 10^{3})$. While the highest resolution runs with $(N_{particle}, N_{grid}, k_{max}L/2\pi) = (4 \times 10^{5}, 4 \times 10^{4}, 4 \times 10^{3})$ only allowed us to test the N-body results and Zel'dovich solution without adaptive smoothing, lower resolution runs included full PT predictions with and without adaptive smoothing.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of the power spectrum (left-hand panel) and phase-space structure (right-hand panel) in a CDM-like cosmology. From top to bottom, the results are shown at $z = 5.3$, 1.5 and 0. On the left-hand panel, symbols with error bars give the power spectrum measured in the simulations. The error bars correspond to the variance over all the modes contained in each bin and sampled from the 50 realizations. Predictions based on basic post-collapse PT (blue solid) and Zel'dovich solution (green dotted) are plotted, with and without adaptive smoothing, respectively, for the thick upper and thin lower curves. In addition, the analytic result obtained from equation (76) is shown as a thin black solid curve and linear theory displayed as a thin black dashed curve. On the right-hand panel, the left and right inserts compare, at various times, the simulated phase-space structure to the prediction of post-collapse PT (blue) and Zel'dovich solution (green), without and with adaptive smoothing, respectively. The data are displayed for a particular realization and in a small interval of scales.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as in Fig. 6, but simulations results are now compared to the variants of post-collapse PT predictions, *hc* (cyan dot-dashed), *ho* (black dotted) and *spl* (dashed magenta).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Lagrangian/Eulerian correspondence of phase-space structures shown in right-hand panel of Fig. 6. Left-hand panel shows the velocity as a function of Lagrangian position *q* at various times, while right-hand panel plots the relation between Lagrangian and Eulerian positions. Predictions presented here correspond to Zel'dovich solution (green) and basic post-collapse PT (blue) with adaptive smoothing.</figcaption></figure>

In Fig. 6, the simulations results are compared to basic PT predictions (solid blue) and Zel'dovich solution (dotted green), while, in Fig. 7, they are compared to various higher order variants of post-collapse PT. Note, that except for the solid line that we explain below, all power spectra predicted by post-collapse PT or Zel'dovich solution are computed directly from the particle distribution evolved according to the displacement field predicted by each theoretical model under scrutiny. So it is important to mention that we did not actually compute analytically power spectra in these cases, although it would be in fact possible. The thin black solid lines, on the other hand, provide the analytic results for the Zel'dovich solution without adaptive smoothing, $P_{ZA}(k)$ (e.g. Bond & Couchman 1988; Schneider & Bartelmann 1995; McQuinn & White 2016):

$$ 
\begin{aligned}
P_{ZA}(k; z) &= \int_{0}^{\infty} dq \cos(k q) \left[ e^{-k^2 \{I(0) - I(q)\} D_{+}(z)^2} - 1 \right]; \\
I(q) &= \int_{0}^{\infty} \frac{dp}{\pi} \cos(p q) \frac{P_{1D}(p)}{p^2},
\end{aligned} 
\quad (76) 
$$

where $P_{1D}(p)$ is the un-truncated power spectrum, i.e. without suppression of power at $k \ge k_{max}$. The black solid curves therefore differ slightly from the green dotted ones on left-hand panel of Fig. 6, because of the effect of the cut-off at $k_{max}$ (and not because of discreteness effects or finite grid resolution). While the cut-off wavenumber $k_{max}$ has, as expected, a noticeable effect on the un-smoothed Zel'dovich solution as well as the un-smoothed post-collapse PT prediction and its variants (a property that we do not show but that we verified), its influence becomes negligible in the range of values of $k$ under consideration, both for the simulations and the theoretical models with adaptive smoothing.

As expected, without adaptive smoothing, both post-collapse PT prescriptions and Zel'dovich solution deviate from the simulations

as soon as the non-linear scale becomes larger than the cut-off scale $2\pi/k_{max}$, which is already the case in upper insert of left-hand panel of Figs 6 and 7, i.e. at $z = 5.3$. Still, at this redshift, post-collapse PT (or its variants) without adaptive smoothing is able to capture the main features of phase space better than Zel'dovich solution, although this is not obvious to decipher on upper left insert of right-hand panel of Figs 6 and 7. At lower redshift, haloes relax to a highly non-linear state and can merge together, which cannot be described well by Zel'dovich solution that gives elongated structures in phase space (middle left and lower left inserts of right-hand panels of Fig. 6. As already argued in previous section, post-collapse PT provides only local corrections to Zel'dovich flow and therefore cannot account for non-linear dynamics at large scales, hence presents the same defect. This elongation effect obviously leads to a strongly underestimated power spectrum. Note however that this depends also on the choice of the cut-off scale $2\pi/k_{max}$: if this scale would be kept, for the theoretical predictions, variable so that it remains close to the typical transition scale to the non-linear regime, the results would improve greatly both for Zel'dovich solution and for post-collapse PT, as it is the case for upper insert of left-hand panels of Figs 6 and 7. But here, instead of varying the effective smoothing scale $2\pi/k_{max}$ for the PT predictions, we now focus on adaptive smoothing, which should be even better because it provides locally the 'optimal' softening of initial conditions.

When adaptive smoothing is employed, agreement between PT predictions and simulation improves strikingly both from the visual and statistical points of view, even for Zel'dovich solution which performs comparably well to post-collapse PT at $z = 1.5$, reproducing simulations results in the whole available dynamic range, including the highly non-linear regime, with a slight underestimation of the power spectrum, though. Of course, superiority of post-collapse PT shows up at $z = 0$ but Zel'dovich solution remains quite good. Keep in mind, however, that Zel'dovich solution with adaptive smoothing as presented here uses the calculation of next-crossing time to estimate the local smoothing scale, hence relies as well partly on post-collapse PT.

These results confirm the conclusions of Section 5.3: adaptive smoothing allows one to summarize composite structures into a single S shape like halo with the right size so that two-point statistics matches the true one even in the highly non-linear regime. In fact, PT predictions reproduce, at least partly, the plateau that can be observed at small scales in middle and lower insert of left-hand panels of Figs 6 and 7. In this regime that we will discuss more in details in next section, the system is most probably following stable clustering. Of course, because of the way we implement our adaptive smoothing, continuity of the phase-space sheet is not preserved anymore, as illustrated by right-hand panel of Figs 6 and 7, but the overall description of phase-space structures is nevertheless improved tremendously. This is illustrated even more clearly by Fig. 8, which also shows, as expected, that positions predicted by post-collapse PT or Zel'dovich solution match the simulation significantly better than predicted velocities. This is a trivial consequence of the fact that positions correspond to velocities integrated over time: by construction, velocities are better approximated by theoretical models at early than at late times, hence their integrated counterpart, which corresponds to some averaged behaviour, will compare better to the simulations. Note thus on right-hand panel of Fig. 8, the excellent agreement with the simulation for the predicted position from basic post-collapse PT, even when significant mergers take place: one just needs to compare bottom right insert of right-hand panel of Fig. 6 to bottom insert of right-hand panel of Fig. 8 to be convinced by this state of fact.
<figure><img src="image_10.png" /><figcaption>Figure 9. Power spectrum (left-hand panel) and phase-space diagram (right-hand panel) from scale-free simulations. Left-hand panel shows the matter power spectrum at $z=0$. Upper and lower inserts correspond, respectively, to $n=0$ and 1. The symbols with error bars give simulation results, while predictions with adaptive smoothing are plotted as a blue solid curve for post-collapse PT and as a green dotted curve for Zel'dovich solution. For reference, linear theory and analytic Zel'dovich power spectra are also plotted as black dashed and black solid curves as well as the prediction from stable clustering as thin straight lines. Right-hand panel shows, for the same values of $n$, a snapshot of a small region of phase space at $z=0$. Predictions from basic post-collapse PT (blue) and Zel'dovich solutions (green) are plotted together with simulation results (red). Note that the result shown here is generated with the same random seed as in right-hand panels of Figs 6 and 7.</figcaption></figure>

As a consequence, we obtain an excellent match between basic post-collapse PT with adaptive smoothing and simulations measurements for the power spectrum, which is a two-point statistics not depending directly on the quality of representation of the velocity. These arguments also apply, to a lesser extent, obviously, to the Zel'dovich solution.

When focusing again on the power spectrum (left-hand panel of Figs 6 and 7), a detailed comparison between various variant of post-collapse PT models reveals that the basic prescription remains the best among all of them, although spl performs nearly as well. This confirms the results of Sections 5.2 and 5.3, where we noticed that higher order prescriptions for post-collapse PT tended to overpredict the size of haloes and mergers (although this is not obvious at first sight when examining phase-space diagrams of right-hand panels of Figs 6 and 7), leading to a smaller power spectrum than the basic post-collapse PT prediction at small scales.

Note that, in agreement with intuition, the results of PT with adaptive smoothing are rather insensitive to the choice of cut-off wavenumber $k_{max}$ (as long as it is kept large enough compared to the wavenumber corresponding to transition towards non-linearity), but the quality of the agreement with simulations depends, on the other hand, on the value of the parameter $f_{cross}$ introduced in Section 4. Here, we adopt the intuitive choices $f_{cross} = 1$ for post-collapse PT and 0.5 for Zel'dovich solution, but, as illustrated in Appendix F, it is possible to optimize the value of $f_{cross}$. For instance at $z=0$, we find for the analyses performed in this section that $f_{cross} = 0.6$ and 0.3 would provide better agreement with simulation measurements. However, we also noticed that the optimal value of $f_{cross}$ depends slightly on redshift as well as on initial conditions, this is why we decided for the analyses presented in this article to stick to our fiducial, slightly sub-optimal choice.

## 5.5 Random initial conditions: power-law power spectrum

We now turn to random initial conditions with a power-law power spectrum:

$$P_{1D}(k) = A k^n. \qquad (77)$$

This allows us to test post-collapse PT for different kinds of initial conditions and to understand more in details small-scale dynamics. Assuming Einstein-de Sitter universe ($\Omega_{m,0} = 1$, $\Omega_{\Lambda} = 0$ and $h = H_0/100 = 1$), we performed simulations for $n=0$ and 1, with the same resolution as our main CDM runs, i.e. $N_{particle} = 200000$, $N_{grid} = 20000$ and $k_{max} = 12.6 \text{ Mpc}^{-1}$. Although the initial conditions are scale-free, we set an actual size for the simulation box, $L = 1000 \text{ Mpc}$, which allows us to normalize the power spectrum using $\sigma_8 = 1.0$ (see footnote 2). With such a normalization, the level of non-linearity of $P_{1D}(k)$ at the end of the simulations is roughly the same as in the CDM runs studied in previous section. To measure power spectra with sufficient statistics, we performed sets of 50 simulations with different initial random seed for each value of $n$.

Fig. 9 displays the power spectra (left-hand panel) and phase-space diagrams (right-hand panel) measured at $z=0$ in the scale-free simulations. The numerical results are compared to basic post-collapse PT (blue solid) and to Zel'dovich solution (green dotted) with adaptive smoothing. For reference, linear theory (black dashed) and the analytic prediction given by equation (76) for the Zel'dovich solution without smoothing (black solid) are also plotted in left-hand panel.<fnref n="4" /> In addition, predictions from stable clustering (Joyce

<footnote n="4">⁴ While computing $I(q)$ in equation (76), we introduced the actual cut-offs in the initial conditions of the simulations, i.e. below $k_{min} = 2\pi/L$ and above $k_{max} = 12.6 \text{ Mpc}^{-1}$.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Lagrangian/Eulerian correspondence of phase-space structures shown in right-hand panel of Fig. 9. Upper and lower inserts correspond, respectively, to *n* = 0 and 1, while left- and right-hand panels correspond, respectively, to the velocity and the position as functions of Lagrangian coordinate *q*. Predictions presented here correspond to Zel'dovich solution (green) and basic post-collapse PT (blue) with adaptive smoothing.</figcaption></figure>

& Sicard 2011; Benhaiem, Joyce & Sicard 2013) are displayed as a thin solid line and will be discussed more in details below. To complete visual inspection of phase space, Fig. 10 plots velocities and positions as functions of Lagrangian coordinate *q* in the same region of phase space as in right-hand panel of Fig. 9. Note that the simulations considered in right-hand panel of Figs 9 and 10 have the same random initial seed as for right-hand panel of Figs 6 and 8, so one can compare directly phase-space features between scale-free and CDM-like simulations.

Examination of Figs 9 and 10 demonstrates again the success of post-collapse PT with adaptive smoothing. While the results are nearly as good for *n* = 0 as for the CDM-like case, post-collapse PT does slightly less well for *n* = 1. The same trend can be seen for Zel'dovich solution with adaptive smoothing. This result is not surprising, for two reasons. First, the level of non-linearity in the *n* = 1 simulations is slightly larger than for *n* = 0 and the CDM-like initial conditions. Indeed, deviation from linear theory or pure Zel'dovich solution happens at smaller values of *k* for *n* = 1 than for other cases. Secondly and more importantly, adaptive smoothing using a sharp low-pass filter is less efficient, from the dynamical point of view, when considering a power spectrum with a high slope. The results would look even worse in the case *n* = 2 that we simulated but do not show here. To explain this, one can just take the example of the single *k* = *k*<sub>single</sub> mode system, e.g. initial conditions with a single sine wave as shown in Fig. 2. In this case, with a sharp-*k* filter, the fluctuation at *k*<sub>single</sub> is either kept, hence post-collapse PT becomes invalid after the system reaches next-crossing time, or it is completely suppressed, which means that the procedure with adaptive smoothing is at the end equivalent to no smoothing at all. Increasing the power spectrum slope *n* basically makes us approaching the single-mode case with *k*<sub>single</sub> = *k*<sub>max</sub>, since we always introduce a cut-off at some value *k*<sub>max</sub> to have a sufficiently smooth representation of the phase-space sheet at very small scales.

Of course, this situation is unrealistic in practice, because no such problem occurs when considering CDM-like power spectra, but one has to stay aware of the fact that adaptive smoothing does not always work. A possible way to improve it could consist in changing the filtering window so that all the wavenumbers are affected by the softening procedure, even if it is performed predominantly at a given scale: to this respect, adaptive smoothing with a Gaussian window might represent an interesting alternative to our sharp-*k* filter.

Another interesting result of our measurements is the excellent agreement between the simulations and the stable clustering prediction from Joyce & Sicard (2011) and Benhaiem et al. (2013) in the non-linear regime<fnref n="5" />, who predict that the slope of the power spectrum should be equal to

$$ \alpha = \frac{d \log k P_{1D}(k)}{d \log k} = \frac{n+1}{n+7}. \qquad (78) $$

In fact, even in the CDM-like case, which roughly corresponds, in one dimension, to an asymptotic value of *n* equal to −1 at large *k*, stable clustering prediction works, since it predicts a plateau, *α* = 0, which is indeed what we observe on left-hand panel of Fig. 6. Note that what Joyce & Sicard (2011) and Benhaiem et al. (2013) call ‘stable clustering’ is very specific. In the 1D case considered here, it is equivalent to assuming that relaxed objects become of constant size in the coordinate *r*′ = *a*<sup>1/3</sup>*x* instead the physical coordinate *r* = *ax* as normally supposed in the three-dimensional case (see, e.g. Davis & Peebles 1977; Peebles 1980). It is important to notice here that post-collapse PT seems to be able to reproduce correctly the measured power spectrum in the beginning of the stable clustering regime. This suggests that stable clustering regime is established only in a few dynamical times, since post-collapse PT is only able to follow the dynamics of a single halo until next-crossing time.

On the other hand, if one considers a ‘halo model’ approach (see, e.g. Ma & Fry 2000; Peacock & Smith 2000; Scoccimarro et al. 2001; Cooray & Sheth 2002, and references therein) and assumes that at large *k* the power spectrum is dominated by the ‘one-halo term’, interesting properties are expected, at least when initial density fluctuations are sufficiently smooth at small scale. Indeed, one-dimensional simulations of single cold haloes suggest that their central density profile nearly behaves like a power law

$$ \rho(x) \propto x^{-\gamma}, \qquad (79) $$

with 0.4 ≲ γ ≲ 0.5. The actual value of γ is controversial. While Binney (2004) postulated it to be equal to 0.5, detailed N-body simulations suggest instead γ ≃ 0.47 (Schulz et al. 2013). On the other hand, pure Vlasov simulations results using the waterbag method suggest γ = 0.4 (Colombi & Touma 2014). Other values to keep in mind are γ = 2/3 from the singularity occurring exactly at crossing times and γ = 1/2 from the caustics (see, e.g. Shandarin & Zel'dovich 1989). These singularities imply, as discussed below, some expected behaviour for *P*(*k*) at large *k* (Gouda & Nakamura 1989; Yano et al. 2004).

A slope −γ for the density profile corresponds, after some trivial algebraic calculations, to a power-law slope

$$ \alpha = 2\gamma - 1 \qquad (80) $$

<footnote n="5">⁵ We also note here that the asymptotic slope seen in the simulations is rather different from the one predicted by the so-called adhesion model based on Burgers’ equation, which gives $kP(k) \to k$ at very high-*k*, irrespective of the initial power-law slope (Valageas 2009a,b; Valageas & Bernardeau 2011).</footnote>
for the power spectrum. For instance, one expects the Zel'dovich approximation to be dominated by the caustics at very small scales, i.e. $kP(k) \to$ constant at large $k$, which is indeed the case on left-hand panels of Figs 6 and 9 in both analytic (black) and numerical (green) cases, irrespective of adaptive smoothing. Turning to post-collapse PT with adaptive smoothing, since, by using $f_{\text{cross}} = 1$, we push the (smoothed) system exactly to next-crossing time, we are dominated by the singularity with $\gamma = 2/3$. Hence, in this case one expects $\alpha$ close to $1/3$ when $k$ becomes very large, as can be observed on all the figures. With a smaller value of $f_{\text{cross}}$, the system is dominated, similarly as for Zel'dovich solution, by the caustics at very small scales and then $kP(k)$ reaches a plateau at large $k$ as can be seen in Fig. E1 of Appendix F. Without adaptive smoothing, one also intuitively expect post-collapse PT power spectrum to be dominated by caustics at very small scales, however the objects have a more complex structure than Zel'dovich solution as shown by lower right panel of Fig. 2, hence the plateau at large $k$ is not clearly visible on left-hand panel of Figs 6, 7 and 9 probably because not reached as quickly as for the pure Zel'dovich solution when increasing $k$.

Turning to the actual non-linear evolution of single objects, assuming as found in the literature $\gamma = 0.4$ and $0.5$ for the density profile of a halo would respectively give $\alpha = -0.2$ and $0$ for the 'one halo' power spectrum. This of course supposes that statistical averaging over many haloes of different masses does not affect the slope of the power spectrum predicted by the single halo. Again, we have to be aware of the fact that at very large $k$, one must be dominated by the caustics, since we introduce a cut-off at large $k$. Such a cut-off enforces some level of smoothness of the curve representing the system in phase space and thus actually imposes the existence of well-defined caustic structures in the non-linear regime,<fnref n="6" /> hence a plateau for $kP(k)$ at very large $k$, even for the simulations. However, this plateau is not visible in the simulations, because it lies outside the range of $k$ shown in the figures. On the other hand, there must be, according to our 'halo model' prescription, a regime where $kP(k)$ might decrease ($\gamma = 0.4$) or present another plateau ($\gamma = 0.5$) or at least some kind of inflection at some intermediate values of $k$. This is what can indeed be observed, e.g. on lower insert of left-hand panel of Fig. 6 for $0.6 \lesssim k \lesssim 2 \text{ Mpc}^{-1}$ for the CDM-like case and on lower insert of left-hand panel of Fig. 9 for $1 \lesssim k \lesssim 2 \text{ Mpc}^{-1}$ for the $n = 1$ simulation. There is no such signature obviously visible in the $n = 0$ simulation, although the measurements are too noisy to conclude in this case. Interestingly, when the inflection is visible, its extension in $k$ space is correctly predicted by basic post-collapse PT with adaptive smoothing. This just shows that basic post-collapse PT is able to predict approximately the correct size for the haloes, which we already knew.

Note that here, we did not try to position our single halo term discussion in the context of stable clustering predictions. These latter can be considered as the outcome of some statistical averaging after weighting the one-halo term by the mass function. They are also expected to be valid only in a finite range of values of $k$ if one assumes a cut-off at $k_{\text{max}}$: what we mean here is that our 'one-halo term' discussion is in fact not incompatible with the concept of stable clustering although the link between the halo model and the stable clustering predictions remains to be performed in details in the present case. Clearly, however, the small inflection discussed above marks a small deviation from stable clustering.

# 6 CONCLUSION

In this paper, we have developed a Lagrangian PT method for solving, beyond crossing time, Vlasov-Poisson equations in the expanding universe for initially cold systems in the infinite parallel planes geometric set up. The proposed approach captures post-collapse dynamics by computing at leading order a counter term to the Zel'dovich solution just after first crossing time. It extends earlier work of Colombi (2015) to the cosmological case and, thanks to adaptive smoothing, to random initial conditions instead of a single halo. By performing a local Taylor expansion of the velocity and position of each mass element as functions of Lagrangian coordinate around each initial density peak, we are able to compute the force field in the multivalued region. A correction to the single-stream flow motion, described by the Zel'dovich solution, is systematically obtained. The results are given as polynomial functions of Lagrangian position and the coefficient of each power is expressed in terms of quantities associated to density peaks. A simple analysis of the formation of a single cluster shows that post-collapse PT reproduces well the phase-space structure up to next-crossing time and can describe the projected density profile quite well even at later time.

To deal with the bulk properties of phase-space structures, we propose an adaptive filtering scheme. The idea is very close to the peak-patch treatment by Bond & Myers (1996). In particular, the local smoothing scale of initial conditions is determined by the calculation of a dynamical time, but in post-collapse PT this dynamical time is given by next-crossing time instead of collapse time. This allows us to produce a realistic coarse-grained version of the phase-space structure, even in the presence of mergers. After demonstrating that adaptive smoothing indeed works well in a simple setup with coalescing clusters, we applied it to more general, random initial conditions that include CDM-like as well as scale-free power spectra. The results of our analyses show that post-collapse PT predictions with adaptive smoothing reproduce the power spectrum of our numerical simulations remarkably well even at scales sufficiently non-linear to probe the beginning of the stable clustering regime.

On important thing to keep in mind for further work is that the critical quantity in post-collapse PT is the time of next crossing, which can in fact be used to also create a spectacularly good pure Zel'dovich prescription with adaptive smoothing. While the calculations of post-collapse PT can be cumbersome, the expression (59) of this time is relatively simple and this can be very interesting when trying to generalize results of post-collapse PT to the three-dimensional case. When analysing equation (59), one indeed notes that $\tau_{\text{cross}}$ depends both on collapse time $\tau_0$ and on the height $\delta_L$ of the initial density peak. Since $\tau_0$ depends solely on $\delta_L$, we can see that basically the dynamical time we consider for defining the local smoothing scale is a single function of $\delta_L$, larger than $\tau_0$. In fact, things are very simple: when performing adaptive smoothing one enforces $\tau_0 + \tau_{\text{cross}}$ to be equal to present time, or in other words, adaptive smoothing of initial conditions is performed so that the local peak height given by linear theory at present time, $\delta_L$, is a number larger than unity (instead of exactly unity for the collapse time). For instance, in the Einstein-de Sitter case, $\Omega_{m,0} = 1$ and $\Omega_{\Lambda} = 0$, we have

$$ \delta_L = \frac{1}{(1 - f_{\text{cross}}/3)^2}, \qquad (81) $$

with $f_{\text{cross}}$ chosen in this article to be equal to unity for post-collapse PT and one-half for the Zel'dovich solution. Note that

<footnote n="6">⁶ Things would be more difficult to apprehend if the curve representing the phase-space sheet would be unsmooth.</footnote>
the optimal value of $f_{cross}$ is smaller than these intuitive choices (see Appendix F).

The present post-collapse PT treatment together with the adaptive smoothing scheme seems a very promising tool. A successful extension to the 3D case should lead to a major breakthrough in the precision PT calculations of large-scale structure statistics beyond the single-stream approximation. Our 1D study is however only a first step. Towards a practical application to the 3D case, there are several issues to be addressed. One is the analytic calculation of statistical quantities such as the power spectrum. Indeed, in this work, we measured directly the power spectrum from a particle distribution following the dynamics prescribed by the various PT prescriptions under scrutiny. In fact, analytic calculations of the power spectrum are tractable, at least in the presence of fixed smoothing. Adaptive smoothing itself can in principle be implemented in the analytic framework by inspiring for instance from the theory of excursion sets discussed in e.g. Bond et al. (1991), although we can foresee that actual calculations will probably be very cumbersome.

Another issue, when generalizing to the 3D case, is that the Zel'dovich solution no longer provides a sufficiently accurate description of the dynamics before collapse. At least third-order Lagrangian PT is required to estimate correctly collapse times as well as the power spectrum beyond tree-level approximation, which clearly adds some significant level of complexity. Furthermore, the variety and complexity of topological configurations of singularities appearing in the 3D case might in fact represent a nearly impassable obstacle in three dimensions (see, e.g. Hidding, Shandarin & van de Weygaert 2014). We will tackle these issues in future work.

Finally, while we have focused in this article on the impact of post-collapse dynamics in 1D on the phase-space structure and on the power spectrum of the matter distribution, it would be interesting to see quantitatively how the small-scale modes responsible for the multistream flows couple to large-scale fluctuations described by the single-stream flows. This can be presumably clarified by measuring the response of the evolved power spectrum to a small initial perturbation of high-$k$ modes, through the comparison between a pair of particle distributions modelled with post-collapse PT that slightly differs in their initial conditions (Nishimichi et al. 2016). In the 3D case, a strong damping is found for the mode transfer from small to large scales, but the physical origin of it is still unclear. A detailed analysis in the 1D cosmological setup would provide important clues on the nature of mode transfers.

## ACKNOWLEDGEMENTS

AT thanks Patrick Valageas for references and useful comments on the adhesion model. SC greatly acknowledges hospitality of Yukawa Institute for Theoretical Physics (YITP), where this work was initiated. We also acknowledge the support of YITP in organizing the workshop 'Vlasov-Poisson: towards numerical methods without particles' in Kyoto, funded by grant YITP-T-15-02, ANR grant ANR-13-MONU-0003 and by Institut Lagrange de Paris (ANR-10-LABX-63 and ANR-11-IDEX-0004-02). This work was also supported in part by MEXT/JSPS KAKENHI grant numbers JP15H05889 and JP16H03977 (AT) as well as ANR grant ANR-13-MONU-0003 (SC).

## REFERENCES

* Alard C., 2013, MNRAS, 428, 340

* Aurell E., Fanelli D., 2001, preprint (arXiv:cond-mat/0106444)

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

* Baldauf T., Codis S., Desjacques V., Pichon C., 2016, MNRAS, 456, 3985

* Baumann D., Nicolis A., Senatore L., Zaldarriaga M., 2012, J. Cosmol. Astropart. Phys., 7, 51

* Benhaiem D., Joyce M., Sicard F., 2013, MNRAS, 429, 3423

* Bernardeau F., Valageas P., 2010, Phys. Rev. D, 81, 043516

* Bernardeau F., Colombi S., Gaztanaga E., Scoccimarro R., 2002, Phys. Rep., 367, 1

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

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

* Bertschinger E., 1985, ApJS, 58, 39

* Binney J., 2004, MNRAS, 350, 939

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

* Bond J., Couchman H., 1988, in Kaiser N., Lasenby A. N., eds, The Post-Recombination Universe, Proceedings of the NATO Advanced Study Institute, Kluwer Academic Publishers, p. 263

* Bond J. R., Myers S. T., 1996, ApJS, 103, 1

* Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440

* 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., 9, 82

* Coles P., Melott A. L., Shandarin S. F., 1993, MNRAS, 260, 765

* Colombi S., 2015, MNRAS, 446, 2902

* Colombi S., Touma J., 2014, MNRAS, 441, 2414

* Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

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

* Davis M., Peebles P. J. E., 1977, ApJS, 34, 425

* Doroshkevich A. G., Ryaben'kii V. S., Shandarin S. F., 1973, Astrophysics, 9, 144

* Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

* Fanelli D., Aurell E., 2002, A&A, 395, 399

* Fillmore J. A., Goldreich P., 1984a, ApJ, 281, 1

* Fillmore J. A., Goldreich P., 1984b, ApJ, 281, 9

* Gouda N., Nakamura T., 1989, Prog. Theor. Phys., 81, 633

* Gradshteyn I. S., Ryzhik I. M., Jeffrey A., Zwillinger D., 2007, Table of Integrals, Series and Products, 7th edn. Elsevier Academic Press

* Gurbatov S. N., Saichev A. I., Shandarin S. F., 1989, MNRAS, 236, 385

* Hahn O., Angulo R. E., 2016, MNRAS, 455, 1115

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

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

* Hockney R. W., Eastwood J. W., 1988, Computer Simulation using Particles. Iwanami shoten, Tokyo

* Joyce M., Sicard F., 2011, MNRAS, 413, 1439

* Lithwick Y., Dalal N., 2011, ApJ, 734, 100

* Ma C.-P., Fry J. N., 2000, ApJ, 543, 503

* Manfredi G., Rouet J.-L., Miller B., Shiozawa Y., 2016, Phys. Rev. E, 93, 042211

* Martel H., Shapiro P. R., 1998, MNRAS, 297, 467

* McQuinn M., White M., 2016, J. Cosmol. Astropart. Phys., 1601, 043

* Melott A. L., Shandarin S. F., Splinter R. J., Suto Y., 1997, ApJ, 479, L79

* Miller B. N., Rouet J. L., 2002, Phys. Rev. E, 65, 056121

* Monaco P., Theuns T., Taffoni G., 2002, MNRAS, 331, 587

* Moriguchi S., Udagawa K., Ichimatsu S., 1987, Fourier Series Analysis (Iwanami mathematical formulas 2). Hilger, Bristol

* Nishimichi T., Bernardeau F., Taruya A., 2016, Phys. Lett., B762, 247

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

* Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144

* Peebles P. J. E., 1980, The Large-scale Structure of the Universe. Princeton Univ. Press, Princeton, NJ

* Planck Collaboration XIII, 2016, A&A, 594, A13

* Rouet J. L., Feix M. R., Navet M., 1990, Vistas Astron., 33, 357

* Ryden B. S., 1993, ApJ, 418, 4

* Sahni V., Coles P., 1995, Phys. Rep., 262, 1

* Schneider P., Bartelmann M., 1995, MNRAS, 273, 475
Schulz A. E., Dehnen W., Jungman G., Tremaine S., 2013, MNRAS, 431, 49

Scoccimarro R., Sheth R. K., Hui L., Jain B., 2001, ApJ, 546, 20

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

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

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

Valageas P., 2009a, J. Stat. Phys., 134, 589

Valageas P., 2009b, J. Stat. Phys., 137, 729

Valageas P., Bernardeau F., 2011, Phys. Rev. D, 83, 043508

Vlah Z., Seljak U., Chu M. Y., Feng Y., 2016, J. Cosmol. Astropart. Phys., 1603, 057

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

Yoshikawa K., Yoshida N., Umemura M., 2013, ApJ, 762, 116

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

## APPENDIX A: DERIVATION OF INTEGRAL EXPRESSION FOR THE POTENTIAL Φ

In this appendix, we derive the expression in equation (19). To do this, we first consider the 1D Green's function which satisfies the following equation:

$$ \frac{\partial^2}{\partial x^2} G(x, x') = -\delta_D(x - x'). \qquad (A1) $$

With this Green function, the solution of the Poisson equation is formally expressed as

$$ \Phi(x) = -\frac{3}{2} H_0^2 \Omega_{m,0} a \int_0^L dx' G(x, x') \delta(x'). \qquad (A2) $$

To solve equation (A1) under periodic boundaries condition, we move to harmonic space and express the Green's function as:

$$ G(x, x') = \sum_{n=-\infty}^{+\infty} G_n e^{i 2n \pi x'/L}. \qquad (A3) $$

Recalling that the Fourier transform of Dirac's delta function reads $\delta_D(y) = (1/L) \sum_n \exp[i 2n \pi y/L]$,<fnref n="7" /> substitution of equation (A3) into equation (A1) gives

$$ G_n = \frac{L}{(2n\pi)^2} e^{-i 2n \pi x'/L}; \quad n \neq 0. \qquad (A4) $$

Thus, we obtain

$$ G(x, x') = \sum_{n=-\infty}^{+\infty} \frac{L}{(2n\pi)^2} e^{-i 2n \pi (x-x')/L} = \sum_{n=1}^{+\infty} \frac{2L}{(2n\pi)^2} \cos \left\{ \frac{2n\pi(x-x')}{L} \right\}. \qquad (A5) $$

Note that there is a subtlety in the above equation. In deriving it, we ignore the contribution from $n = 0$ mode in $G_n$, while the $n = 0$ mode does not vanish in Dirac's delta function. The result remains however correct as long as the density field does not contain a $n = 0$ mode, or equivalently, $\int_0^L dx \delta(x) = 0$.

Keeping this point in mind, we can further simplify the expression. We use the following formula (see, e.g. Chap. 1.44 of Gradshteyn et al. 2007; Moriguchi, Udagawa & Ichimatsu 1987):

$$ \sum_{n=1}^{\infty} \frac{\cos(n y)}{n^2} = \frac{1}{4}(y - \pi)^2 - \frac{\pi^2}{12}; \quad 0 \le y \le 2\pi. \qquad (A6) $$

The final expression for the Green's function is then given by

$$ G(x, x') = \frac{L}{2} \left[ \left\{ \frac{|x - x'|}{L} - \frac{1}{2} \right\}^2 - \frac{1}{12} \right], \quad |x - x'| \le L. \qquad (A7) $$

Substituting this into equation (A2), we obtain the expression in equation (19).

## APPENDIX B: DERIVATION OF THE EXPRESSION FOR THE FORCE

In this appendix, we derive the expression for the force in the multistream region, given in equations (40)–(42).

To start with, we divide the integral in equation (20) into several pieces according to equation (39):

$$ F(x) = -\frac{3}{2} \Omega_{m,0} H_0^2 a \left[ \int_{x(-\bar{Q}_c)}^{x(\bar{Q}_c)} dx' \frac{\delta(x')}{2} \left\{ \Theta(x(Q) - x') - \Theta(x' - x(Q)) \right\} + \frac{1}{L} \int_{x(-\bar{Q}_c)}^{x(\bar{Q}_c)} dx' x' \delta(x') + \left( \int_{0}^{x(-\bar{Q}_c)} + \int_{x(\bar{Q}_c)}^{L} \right) dx' \frac{\delta(x')}{2} \right] $$

<footnote n="7">⁷ Here, we assume $0 \le y \le L$.</footnote>
$$ \times \{\Theta(x-x') - \Theta(x'-x)\} + \frac{1}{L} \left( \int_{0}^{x(-\hat{Q}_{c})} + \int_{x(\hat{Q}_{c})}^{L} \right) dx' x' \delta(x') \quad (B1) $$

We shall below evaluate these integrals one by one.

Consider the first integral in right-hand side of equation (B1). We evaluate it at position $x_* = x(Q_*) = x(Q_\pm)$ in Fig. 1. Using equation (16), the integral is rewritten

$$ \int_{x(-\hat{Q}_c)}^{x(\hat{Q}_c)} dx' \frac{\delta(x')}{2} \{\Theta(x_* - x') - \Theta(x' - x_*)\} = \int_{-\hat{Q}_c}^{\hat{Q}_c} \frac{dQ'}{2} \left\{1 - \left(\frac{\partial x}{\partial q}\right)_{Q'}\right\} [\Theta(x(Q_*; \tau) - x(Q'; \tau)) - \Theta(x(Q'; \tau) - x(Q_*; \tau))]. \quad (B2) $$

Exploiting the geometric properties of the configuration in Fig. 1, substitution of the explicit form of $x(Q)$ (equation 28) into the above leads to

$$ \left( \int_{-\hat{Q}_c}^{Q_-} - \int_{Q_-}^{Q_*} + \int_{Q_*}^{Q_+} - \int_{Q_+}^{\hat{Q}_c} \right) \frac{dQ'}{2} \{1 + B(q_0; \tau) - 3C(q_0; \tau) Q'^2\} = \{1 + B(q_0; \tau)\} \times (Q_- + Q_+ - Q_*) - C(q_0; \tau) (Q_-^3 + Q_+^3 - Q_*^3). \quad (B3) $$

Here, we used the fact that $Q_- < Q_* < Q_+$. Recalling that $Q_\pm$ and $Q_*$ are related with each other through equations (32) and (33), equation (B3) can be recast as a single function of either $Q_\pm$ or $Q_*$, which respectively provides the general expression valid for $Q_c < |Q| < \hat{Q}_c$ and $|Q| < Q_c$. Omitting the subscript $\pm$ or $*$ in either cases, we finally obtain

$$ \int_{x(-\hat{Q}_c)}^{x(\hat{Q}_c)} dx' \frac{\delta(x')}{2} \{\Theta(x-x') - \Theta(x'-x)\} = \begin{cases} \{1 + B(q_0; \tau)\} Q - C(q_0; \tau) Q^3 & - \operatorname{sgn}(Q) \sqrt{3(\hat{Q}_c^2 - Q^2)}; \quad Q_c < |Q| < \hat{Q}_c, \\ -2 + B(q_0; \tau)\} Q - C(q_0; \tau) Q^3; & |Q| < Q_c. \end{cases} \quad (B4) $$

Next, we consider the second integral in equation (B1). Rewriting the integrand in terms of the variable $Q$ with the help of equations (16) and (28), it is straightforward to obtain

$$ \frac{1}{L} \int_{x(-\hat{Q}_c)}^{x(\hat{Q}_c)} dx' x' \delta(x') = \frac{1}{L} \int_{-\hat{Q}_c}^{\hat{Q}_c} dx' \{A(q_0; \tau) - B(q_0; \tau) Q' + C(q_0; \tau) Q'^3\} \{1 + B(q_0; \tau) - 3C(q_0; \tau) Q'^2\} \\ = \frac{2}{L} A(q_0; \tau) \{1 + B(q_0; \tau) - C(q_0; \tau) \hat{Q}_c^2\} \hat{Q}_c. \quad (B5) $$

The remaining integrals to be computed, i.e. the third and fourth terms in equation (B1), are performed over intervals outside the multivalued region, which allows us to use Zel'dovich solution to evaluate the integrands. For the third integral in equation (B1), we obtain

$$ \left( \int_{0}^{x(-\hat{Q}_c)} + \int_{x(\hat{Q}_c)}^{L} \right) dx' \frac{\delta(x')}{2} \{\Theta(x-x') - \Theta(x'-x)\} = \left( \int_{Q(x=0)}^{-\hat{Q}_c} - \int_{\hat{Q}_c}^{Q(x=L)} \right) dQ' \frac{1}{2} \left( -\frac{\partial \psi}{\partial Q'} \right) D_+(\tau) \\ = \frac{1}{2} \left\{ \psi(Q(x=0)) + \psi(Q(x=L)) - \psi(-\hat{Q}_c) - \psi(\hat{Q}_c) \right\} D_+(\tau). \quad (B6) $$

Although this expression is exact, our perturbative treatment implicitly assumes that the multistream region is small and the quantities inside this region are written in polynomial forms of $Q$. We may expand the last two terms of equation (B6) as

$$ \psi(-\hat{Q}_c) + \psi(\hat{Q}_c) \approx 2\psi(q_0) + \mathcal{O}(\hat{Q}_c^4), \quad (B7) $$

where we used the fact that $\psi''(q_0) = 0$ (see equation 27). Additionally, periodic boundaries imply $\psi(Q(x=0)) = \psi(Q(x=L))$. Therefore, equation (B6) can be recast as

$$ \left\{ \int_{0}^{x(-\hat{Q}_c)} + \int_{x(\hat{Q}_c)}^{L} \right\} dx' \frac{\delta(x')}{2} \{\Theta(x-x') - \Theta(x'-x)\} = \left[ \psi(q) \right]_{q_0}^{q(x=L)} D_+(\tau). \quad (B8) $$

In similar manner, we evaluate the fourth integral in equation (B1) as follows:

$$ \frac{1}{L} \left( \int_{0}^{x(-\hat{Q}_c)} + \int_{x(\hat{Q}_c)}^{L} \right) dx' x' \delta(x') = \frac{1}{L} \left( \int_{q(x=0)}^{q(x=L)} - \int_{q_0-\hat{Q}_c}^{q_0+\hat{Q}_c} \right) dq' \{q' + \psi(q') D_+(\tau)\} \left( -\frac{\partial \psi(q')}{\partial q'} \right) D_+(\tau). \quad (B9) $$

The first of these two integrals is performed exactly to give

$$ \frac{1}{L} \int_{q(x=0)}^{q(x=L)} dq' \{q' + \psi(q') D_+(\tau)\} \left(-\frac{\partial \psi(q')}{\partial q'}\right) D_+(\tau) \\ = \frac{1}{L} \left\{ -\left[q' \psi(q')\right]_{q(x=0)}^{q(x=L)} - \frac{1}{2} \left[\psi(q')^2\right]_{q(x=0)}^{q(x=L)} + \int_{q(x=0)}^{q(x=L)} dq' \psi(q') \right\} D_+(\tau) = -\psi(q(x=0)) D_+(\tau), \quad (B10) $$
where we exploited properties of periodic boundaries, $q(x=L) - q(x=0) = L$ and $\psi(q(x=0)) = \psi(q(x=L))$, as well as momentum conservation, which implies $\int_{q(x=0)}^{q(x=L)} dq' \psi(q') = 0$. Turning to the second integral in equation (B9), we use again the fact that the domain of integration is assumed to be narrow to Taylor expand the integrand using the formula

$$ \int_{-\Delta s}^{\Delta s} ds' f(s) \approx 2 f(0) \Delta s + \frac{1}{3} f''(0) (\Delta s)^3, $$

to obtain the following expression, valid at third order in $\hat{Q}_c$:

$$ \begin{aligned} & \frac{1}{L} \int_{-\hat{Q}_c + q_0}^{\hat{Q}_c + q_0} dq' \{ q' + \psi(q') D_+(\tau) \} \left( -\frac{\partial \psi(q')}{\partial q'} \right) D_+(\tau) \\ & \approx \frac{2}{L} \{ q_0 + \psi(q_0) D_+(\tau) \} \left\{ \delta_L(q_0) \hat{Q}_c + \frac{1}{6} \delta_L''(q_0) \hat{Q}_c^3 \right\} D_+(\tau), \\ & = \frac{2}{L} A(q_0; \tau) \left\{ 1 + B(q_0; \tau) - C(q_0; \tau) \hat{Q}_c^2 \right\} \hat{Q}_c. \end{aligned} \quad (B11) $$

Here, equations (29)–(31) are used to rewrite the expression in the last line. Summing up the results in equations (B10) and (B11), we obtain, at the end,

$$ \frac{1}{L} \left( \int_{0}^{x(-\hat{Q}_c)} + \int_{x(\hat{Q}_c)}^{L} \right) dx' x' \delta(x') = -\frac{2}{L} A(q_0; \tau) \left\{ 1 + B(q_0; \tau) - C(q_0; \tau) \hat{Q}_c^2 \right\} \hat{Q}_c - \psi(q(x=0)) D_+(\tau). \quad (B12) $$

Plugging equations (B4), (B5), (B8) and (B12) into equation (B1), we finally obtain the expression of the force in the multivalued region (equations 40–42).

# APPENDIX C: DERIVATION OF BASIC POST-COLLAPSE PT RESULTS

In this appendix, we detail the calculations of basic post-collapse PT predictions summarized in equations (52), (54), (56) and (58), together with the coefficients in Tables 1 and 2.

## C1 Derivation of $\Delta u_{\text{out}}$

Our starting point is equation (51), which involves two integrals. To evaluate the first integral, the approximate formula given by equation (E5) is applied. Recalling that we are considering the outer part of the multivalued region, $Q_c < |Q| < \hat{Q}_c$, we have

$$ \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \approx a(\tau_0) \left[ T Q + \left\{ -\frac{\kappa}{8} + \frac{1}{6} \delta''_L(q_0) D_+(\tau_0) T \right\} Q^3 - \operatorname{sgn}(Q) \frac{\kappa}{4\sqrt{3}} (\hat{Q}_c(\tau)^2 - Q^2)^{3/2} - \frac{\kappa}{48} \delta''_L(q_0) D_+(\tau_0) Q^5 \right], \quad (C1) $$

with $T \equiv \tau - \tau_0$. Here, we used the leading-order expression for critical time $\hat{\tau}_c$, $\hat{\tau}_c - \tau_0 \simeq (\kappa/8) Q^2$, and ignore in equation (E5) the higher order terms involving the quantities $g(\tau_0)$ or $H(\tau_0)$.

On the other hand, the second integral in equation (51) is analytically performed. Using the evolution equation in equation (22), we obtain

$$ -\frac{3}{2} H_0^2 \Omega_{m,0} \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' a(\tau') \mathcal{F}(q_0, \tau') = \psi(q_0) \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \frac{d^2 D_+(\tau')}{d\tau'^2} = \psi(q_0) \left[ \frac{dD_+(\tau')}{d\tau'} \right]_{\hat{\tau}_c(Q)}^{\tau}. \quad (C2) $$

Summing up the above two contributions, we obtain the expression for $\Delta u_{\text{out}}$ summarized in equation (52).

## C2 Derivation of $\Delta x_{\text{out}}$

To derive the expression for $\Delta x_{\text{out}}$, we just need to integrate over time equation (52), as written in equation (53):

$$ \begin{aligned} \Delta x_{\text{out}}(Q; \tau, \hat{\tau}_c) = & \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \Delta u_{\text{out}}(Q, \tau') = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \left[ \tilde{\alpha}_1(\tau') Q + \tilde{\beta}_1(\tau') Q^3 \right. \\ & \left. + \tilde{\gamma}_1(\tau_0) \left\{ \hat{Q}_c^2(\tau) - Q^2 \right\}^{3/2} + \tilde{\delta}_1(\tau_0) Q^5 \right] + \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \hat{\epsilon}_1(\tau', \hat{\tau}_c). \end{aligned} \quad (C3) $$

Making use of the expressions summarized in Table 1, we evaluate the first integral. A straightforward calculation leads to

$$ \begin{aligned} & \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \left[ \tilde{\alpha}_1(\tau') Q + \tilde{\beta}_1(\tau') Q^3 + \tilde{\gamma}_1(\tau_0) \left\{ \hat{Q}_c^2(\tau) - Q^2 \right\}^{3/2} + \tilde{\delta}_1(\tau_0) Q^5 \right] \approx \frac{T^2}{2} Q + \left\{ -\frac{\kappa}{8} T + \frac{\delta''_L(q_0)}{6} D_+(\tau_0) \right\} Q^3 \\ & - \operatorname{sgn}(Q) \frac{\kappa^2}{80\sqrt{3}} \left\{ \hat{Q}_c^2(\tau) - Q^2 \right\}^{5/2} + \left[ \frac{1}{2} \left(\frac{\kappa}{8}\right)^2 - \left(\frac{\kappa}{8}\right) \frac{\delta''_L(q_0)}{6} D_+(\tau_0) T \right] Q^5 + \frac{1}{2} \left(\frac{\kappa}{8}\right)^2 \frac{\delta''_L(q_0)}{6} D_+(\tau_0) Q^7. \end{aligned} \quad (C4) $$
Note again that we used the leading-order expression for $\hat{\tau}_c$ (see equation 37).

The second integral in equation (C3) is analytically performed to give

$$ \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \tilde{\epsilon}_1(\tau, \hat{\tau}_c) = \psi(q_0) \int_{\hat{\tau}_c(Q)}^{\tau} d\tau' \left[ \frac{dD_+(\tau'')}{d\tau''} \right]_{\hat{\tau}_c(Q)}^{\tau'} = \psi(q_0) \left\{ D_+(\tau) - D_+(\hat{\tau}_c(Q)) - \left. \frac{dD_+}{d\tau} \right|_{\hat{\tau}_c(Q)} (\tau - \hat{\tau}_c(Q)) \right\}. \quad (C5) $$

Collecting equations (C4) and (C5), we obtain equation (54).

### C3 Derivation of $\Delta u_{in}$

For $\Delta u_{in}$, the expression to be evaluated is given by equation (55), which involves three integrals. The first integral is obtained from equation (C1) by setting $\tau - \tau_0 \equiv \tau_c(Q) - \tau_0 \simeq (\kappa/2) Q^2$ (equation 38). Recalling the fact that $\hat{Q}_c(\tau_c) = 2Q$ (equation 34), we have

$$ \int_{\hat{\tau}_c(Q)}^{\tau_c(Q)} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \approx a(\tau_0) \left[ -\frac{3}{8} \kappa Q^3 + \frac{\kappa}{16} \delta''_L(q_0) D_+(\tau_0) Q^5 \right]. \quad (C6) $$

For the second integral in equation (55), we use the approximate formula given by equation (E5) in the regime $|Q| < Q_c$. Ignoring higher order terms involving the quantities $g(\tau_0)$ or $H(\tau_0)$, we obtain

$$ \int_{\tau_c(Q)}^{\tau} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \approx a(\tau_0) \left[ -2T Q + \left\{ \kappa + \frac{T}{6} \delta''_L(q_0) D_+(\tau_0) \right\} Q^3 - \frac{\kappa}{12} \delta''_L(q_0) D_+(\tau_0) Q^5 \right]. \quad (C7) $$

Finally, the last integral in equation (55) is evaluated using the formula given in equation (C2). Collecting the three contributions above, we obtain the expression for $\Delta u_{in}$ (equation 56).

### C4 Derivation of $\Delta x_{in}$

The starting expression for $\Delta x_{in}$ is given in equation (57), which involves two integrals. To evaluate the first integral in equation (57), we use equation (C4). Setting $T = \tau_c - \tau_0 \simeq (\kappa/2) Q^2$, an expression at seventh order in $Q$ is obtained:

$$ \int_{\hat{\tau}_c(Q)}^{\tau_c(Q)} d\tau' \Delta \tilde{u}(Q; q_0, \tau') \approx -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \left( -\frac{27}{640} \kappa^2 Q^5 + \frac{3}{4} \left(\frac{\kappa}{8}\right)^2 \delta''_L(q_0) D_+(\tau_0) Q^7 \right). \quad (C8) $$

To compute the second integral in equation (57), we substitute equation (56) into the integrand, but without the term $\tilde{\epsilon}_1$. Using the explicit expression for the coefficients in Table 2, a straightforward calculation gives

$$ \int_{\tau_c(Q)}^{\tau} d\tau' \Delta \tilde{u}(Q; q_0, \tau') = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \int_{\tau_c(Q)}^{\tau} d\tau' \left[ \tilde{\alpha}_3(\tau') Q + \tilde{\beta}_3(\tau') Q^3 + \tilde{\delta}_3(\tau_0) Q^5 + \tilde{\zeta}_3(\tau_0) Q^7 \right] = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \times \left( -T^2 Q + \left\{ \frac{5\kappa}{8} T + \frac{\delta''_L(q_0)}{12} D_+(\tau_0) T^2 \right\} Q^3 - \left\{ \left(\frac{\kappa}{4}\right)^2 + \left(\frac{\kappa}{8}\right) \frac{\delta''_D(q_0)}{6} D_+(\tau_0) T \right\} Q^5 - 4 \left(\frac{\kappa}{8}\right)^2 \frac{\delta''_L(q_0)}{6} D_+(\tau_0) Q^7 \right). \quad (C9) $$

Combining the above two results, we obtain equation (58), which gives the final expression for $\Delta x_{in}$ along with the coefficients in Table 2.

## APPENDIX D: HIGHER ORDER CORRECTIONS TO POST-COLLAPSE PT

In this appendix, we detail calculations of the various higher order corrections to post-collapse PT we propose in Section 3.4.

One possible way to improve on basic post-collapse PT is to Taylor expand at next-to-leading order the expression for $\hat{Q}_c$ (as well as $Q_c$) that determines the boundary of the post-collapse region (equations 34 and 35):

$$ \hat{Q}_c = 2Q_c = \sqrt{\frac{4B}{3C}} \approx \left\{ \frac{8}{\kappa(q_0, \tau_0)} \right\}^{1/2} T^{1/2} \left\{ 1 - \frac{\eta(\tau_0)}{2} T + \dots \right\}, \quad (D1) $$

with

$$ \eta(\tau_0) = \frac{1}{D_+(\tau_0)} \frac{dD_+(\tau_0)}{d\tau_0} - \frac{1}{2} \frac{d^2 D_+(\tau_0)}{d\tau_0^2} \bigg/ \frac{dD_+(\tau_0)}{d\tau_0}. \quad (D2) $$

Inverting the relation $Q = \hat{Q}_c(\hat{\tau}_c) [Q = Q_c(\tau_c)]$, we obtain the corresponding critical times at next-to-leading order,

$$ \hat{\tau}_c(Q) - \tau_0 \approx \frac{\kappa(q_0, \tau_0)}{8} Q^2 + \left\{ \frac{\kappa(q_0, \tau_0)}{8} \right\}^2 \eta(\tau_0) Q^4, \quad (D3) $$

and

$$ \tau_c(Q) - \tau_0 = \hat{\tau}_c(2Q) - \tau_0 \approx \frac{\kappa(q_0, \tau_0)}{2} Q^2 + \left\{ \frac{\kappa(q_0, \tau_0)}{2} \right\}^2 \eta(\tau_0) Q^4. \quad (D4) $$
Using these expressions, we repeat below the same calculations as in Section 3.3. Higher order expressions for the correction to Zel'dovich flow are derived on top of basic post-collapse PT results and are denoted by $\Delta x^{(hc)}$ and $\Delta u^{(hc)}$. In Section 3.4, a semi-analytic treatment designed by spl is also discussed. In this case, we adopt the basic post-collapse PT results in the inner part of the multistream region, but the extension of this latter is calculated using the higher order expression (D1) for $Q_c$. Then, a straightforward third-order spline interpolation, that we do not detail here, is used to connect the inner part to the Zel'dovich solution, assuming again that the Lagrangian boundary $\hat{Q}_c$ of the multistream region is given by the higher order expansion (D1).

Another improvement may come from the calculation of the integrals in equations (43) and (44). In Section 3.3, a part of the integrands is Taylor expanded in time and the integration is performed for the leading-order terms. The resultant expressions include time-dependent terms up to $O(T^1)$ for the velocity and $O(T^2)$ for the position. Here, using the approximate integral formulae for the Taylor-expanded integrands in Appendix E, we keep terms up to $O(T^2)$ for the velocity, $O(T^3)$ for the position, and derive the corresponding higher order corrections to the Zel'dovich flow, $\Delta x^{(ho)}$ and $\Delta u^{(ho)}$. These corrections assume, of course, higher order expansions for critical times $\hat{\tau}_c$ and $\tau_c$.

Below, we present the expressions for $\Delta x^{(ho)}$ and $\Delta u^{(ho)}$. The expressions for $\Delta x^{(hc)}$ and $\Delta u^{(hc)}$ are obtained by simply setting $g(\tau_0)$ and $H(\tau_0)$ to zero.

## D1 Velocity and position in the outer part of the multistream region: $Q_c < |Q| \le \hat{Q}_c$

The expression for the higher order correction to the velocity, $\Delta u^{(ho)}$, is

$$ \Delta u_{\text{out}}^{(\text{ho})}(Q; \tau) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \left[ \tilde{\alpha}_1(\tau) Q + \tilde{\beta}_1(\tau) Q^3 + \tilde{\delta}_1(\tau_0) Q^5 + \tilde{\zeta}_1(\tau_0) Q^7 \right], \quad (\text{D5}) $$

with

$$ \tilde{\alpha}_1(\tau) = \frac{g(\tau_0)}{2} T^2, \quad (\text{D6}) $$

$$ \tilde{\beta}_1(\tau) = \frac{g(\tau_0)}{2} T^2 \frac{\delta''_{L}(q_0)}{6} D_+(\tau_0), \quad (\text{D7}) $$

$$ \tilde{\delta}_1(\tau_0) = -\left(\frac{\kappa}{8}\right)^2 \left\{ \frac{g(\tau_0)}{2} + \eta \right\}, \quad (\text{D8}) $$

$$ \tilde{\zeta}_1(\tau_0) = -\left(\frac{\kappa}{8}\right)^2 \left[ \frac{\delta''_{L}(q_0)}{6} D_+(\tau_0) \left\{ \frac{g(\tau_0)}{2} + \eta \right\} + \frac{\kappa}{8} g(\tau_0)\eta \right]. \quad (\text{D9}) $$

Here, we have introduced the following function:

$$ g(\tau_0) = H(\tau_0) \{a(\tau_0)\}^2 + \frac{d \ln D_+(\tau_0)}{d \tau_0}, \quad (\text{D10}) $$

where $H(\tau_0)$ is the value of the Hubble parameter at time $\tau_0$. For the higher order correction to the position, $\Delta x_{\text{out}}^{(ho)}$, we have

$$ \Delta x_{\text{out}}^{(\text{ho})}(Q; \tau) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) \left[ \tilde{\alpha}_2(\tau) Q + \tilde{\beta}_2(\tau) Q^3 + \tilde{\delta}_2(\tau) Q^5 + \tilde{\zeta}_2(\tau) Q^7 + \tilde{\mu}_2(\tau) Q^9 + \tilde{\nu}_2(\tau) Q^{11} \right], \quad (\text{D11}) $$

with

$$ \tilde{\alpha}_2(\tau) = \frac{g(\tau_0)}{6} T^3, \quad (\text{D12}) $$

$$ \tilde{\beta}_2(\tau) = \frac{\delta''_{L}(q_0)}{6} D_+(\tau_0) \left\{ \frac{T^2}{2} + \frac{g(\tau_0)}{6} T^3 \right\}, \quad (\text{D13}) $$

$$ \tilde{\delta}_2(\tau) = -\left(\frac{\kappa}{8}\right)^2 \left(\frac{g(\tau_0)}{2} + \eta\right) T, \quad (\text{D14}) $$

$$ \tilde{\zeta}_2(\tau) = -T \left\{ \frac{g(\tau_0)}{2} + \eta \right\} \left( \frac{\kappa}{8} \right)^2 \frac{\delta''_{L}(q_0)}{6} D_+(\tau_0) + \left( \frac{\kappa}{8} \right)^3 \left\{ \frac{g(\tau_0)}{3} + \eta \left( 1 - g(\tau_0)T \right) \right\}, \quad (\text{D15}) $$

$$ \tilde{\mu}_2(\tau) = \left(\frac{\kappa}{8}\right)^3 \frac{\delta''_{L}(q_0)}{6} D_+(\tau_0) \left\{ \frac{g(\tau_0)}{3} + \eta \left(1 - g(\tau_0)T\right) \right\} + \left(\frac{\kappa}{8}\right)^4 \eta \left\{ g(\tau_0) + \frac{\eta}{2} \left(1 - g(\tau_0)T\right) \right\}, \quad (\text{D16}) $$

$$ \tilde{\nu}_2(\tau) = \frac{1}{49152} \kappa^4 \eta^2 \delta''_{L}(q_0) D_+(\tau_0). \quad (\text{D17}) $$

Note the term proportional to $Q^{11}$, which is necessary to enforce continuity of the solution.
## D2 Velocity and position in the inner part of the multistream region: $|Q| \le Q_c$

The expression for the higher order corrections to the velocity, $\Delta u_{in}^{(ho)}$, becomes

$$ \Delta u_{in}^{(ho)}(Q; \tau) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) [\tilde{\alpha}_3(\tau) Q + \tilde{\beta}_3(\tau) Q^3 + \tilde{\delta}_3(\tau_0) Q^5 + \tilde{\zeta}_3(\tau_0) Q^7], \quad (D18) $$

with

$$ \tilde{\alpha}_3(\tau) = \frac{T^2}{2} [g(\tau_0) - 3 \{a(\tau_0)\}^2 H(\tau_0)], \quad (D19) $$

$$ \tilde{\beta}_3(\tau) = \frac{g(\tau_0)}{2} T^2 \frac{\delta_L''(q_0)}{6} D_+(\tau_0), \quad (D20) $$

$$ \tilde{\delta}_3(\tau_0) = \frac{47}{64} \kappa^2 \eta - \frac{\kappa^2}{128} [g(\tau_0) - 48\{a(\tau_0)\}^2 H(\tau_0)], \quad (D21) $$

$$ \tilde{\zeta}_3(\tau_0) = -\left(\frac{\kappa}{8}\right)^2 \left\{ \frac{g(\tau_0)}{2} + \eta \right\} \frac{\delta_L''(q_0)}{6} D_+(\tau_0) - \left(\frac{\kappa}{8}\right)^3 \eta \left[ g(\tau_0) - 192\{a(\tau_0)\}^2 H(\tau_0) \right]. \quad (D22) $$

Finally, the higher order corrections to the position, $\Delta x_{in}^{(ho)}$, are given by

$$ \Delta x_{in}^{(ho)}(Q; \tau) = -\frac{3}{2} H_0^2 \Omega_{m,0} a(\tau_0) [\tilde{\alpha}_4(\tau)Q + \tilde{\beta}_4(\tau)Q^3 + \tilde{\delta}_4(\tau)Q^5 + \tilde{\zeta}_4(\tau)Q^7 + \tilde{\mu}_4(\tau)Q^9 + \tilde{\nu}_4(\tau)Q^{11}], \quad (D23) $$

with

$$ \tilde{\alpha}_4(\tau) = \frac{T^3}{6} [g(\tau_0) - 3 \{a(\tau_0)\}^2 H(\tau_0)], \quad (D24) $$

$$ \tilde{\beta}_4(\tau) = \frac{\delta_L''(q_0)}{6} D_+(\tau_0) \frac{g(\tau_0)}{6} T^3, \quad (D25) $$

$$ \tilde{\delta}_4(\tau) = -\left(\frac{\kappa}{4}\right)^2 \left(\frac{g(\tau_0)}{8} T - 6\{a(\tau_0)\}^2 H(\tau_0) T\right) + 47\left(\frac{\kappa}{8}\right)^2 \eta T, \quad (D26) $$

$$ \tilde{\zeta}_4(\tau) = -\left(\frac{\kappa}{8}\right)^2 T \left(\frac{g(\tau_0)}{2} + \eta\right) \frac{\delta_L''(q_0)}{6} D_+(\tau_0) + \left(\frac{\kappa}{8}\right)^3 \left\{ \frac{g(\tau_0)}{3} - 64 \{a(\tau_0)\}^2 H(\tau_0) - \eta \left(95 + g(\tau_0)T - 192 \{a(\tau_0)\}^2 H(\tau_0)T\right) \right\}, \quad (D27) $$

$$ \tilde{\mu}_4(\tau) = \kappa^3 \delta_L''(q_0) D_+(\tau_0) \left( \frac{g(\tau_0)}{9216} + \frac{1}{3072}\eta - \frac{g(\tau_0)\eta T}{3072} \right) + \kappa^4 \eta \left\{ \frac{g(\tau_0)}{4096} - \frac{3}{16}\{a(\tau_0)\}^2 H(\tau_0) - \frac{767\eta}{8092} - \frac{g(\tau_0)\eta T}{8192} + \frac{3}{32}\eta \{a(\tau_0)\}^2 H(\tau_0)T \right\}, \quad (D28) $$

$$ \tilde{\nu}_4(\tau) = \frac{1}{49152} \kappa^4 \eta^2 \delta_L''(q_0) D_+(\tau_0). \quad (D29) $$

## APPENDIX E: APPROXIMATE FORMULAE FOR INTEGRALS

In this appendix, we provide useful formulae for primitive integrals over time, which are used to derive the post-collapse PT results in Appendices C and D.

Our aim here is to evaluate the time-integral of the force in the multivalued region, given in equation (40). The time integral of the second term, involving function $F(q_0, \tau)$, is performed exactly (see equation C2). On the other hand, the integral of the first term needs an approximate treatment which will be valid as long as the multivalued region is sufficiently small, or equivalently, shortly after collapse time $\tau_0$. Below, we evaluate the integral by Taylor expanding the integrand. We first give the results term by term:

$$ \int d\tau' a(\tau') \approx a(\tau_0)(\tau - \tau_0) + \frac{1}{2} H(\tau_0) a(\tau_0)^3 (\tau - \tau_0)^2, \quad (E1) $$

$$ \int d\tau' a(\tau') B(q_0; \tau') \approx \frac{\delta_L(q_0)}{2} a(\tau_0) D_+(\tau_0) \left\{ g(\tau_0) - H(\tau_0) \{a(\tau_0)\}^2 \right\} (\tau - \tau_0)^2, \quad (E2) $$
<figure><img src="image_12.png" /><figcaption>**Figure E1.** Dependence of the PT predictions on the choice of parameter $f_{cross}$ controlling adaptive smoothing in the case of a CDM-like power spectrum. Left-hand panel shows the power spectrum, while right-hand panel plots the phase-space portraits. In each panel, the predictions of post-collapse PT and Zel'dovich solution are, respectively, shown in upper and lower inserts, with various values of $f_{cross}$ (indicated by different colours).</figcaption></figure>

$$ \int d\tau' a(\tau') C(q_0; \tau') \approx \left\{ -\frac{\delta''_L(q_0)}{6} \right\} a(\tau_0) D_+(\tau_0) \left\{ (\tau - \tau_0) + \frac{1}{2} g(\tau_0) (\tau - \tau_0)^2 + \dots \right\}, \quad (E3) $$

where $g(\tau_0)$ is given by equation (D10). For the term with fractional power of time, we obtain

$$ \int d\tau' a(\tau') \sqrt{\hat{Q}_c^2(\tau') - Q^2} \approx a(\tau_0) \frac{\kappa(q_0, \tau_0)}{12} \{\hat{Q}_c^2 - Q^2\}^{3/2}, \quad (E4) $$

where we used the leading-order approximation for $\hat{Q}_c$, i.e. $\hat{Q}_c \approx (8/\kappa)^{1/2}(\tau - \tau_0)^{1/2}$.

Summing up the above formulae, we obtain the primitive of the time integrals of the force involving function $\mathcal{J}(Q; q_0, \tau)$. The result including polynomials of time up to $(\tau - \tau_0)^2$ is

$$ \int_{\tau_i}^{\tau_f} d\tau' a(\tau') \mathcal{J}(Q; q_0, \tau') \approx \begin{cases} a(\tau_0) \left\{ \left[ (\tau - \tau_0) + \frac{1}{2}g(\tau_0)(\tau - \tau_0)^2 \right]_{\tau_i}^{\tau_f} Q + \frac{1}{6}\delta''_L(q_0) D_+(\tau_0) \left[ (\tau - \tau_0) + \frac{1}{2}g(\tau_0)(\tau - \tau_0)^2 \right]_{\tau_i}^{\tau_f} Q^3 \right. \\ \left. - \operatorname{sgn}(Q) \frac{\kappa}{4\sqrt{3}} \left[ (\hat{Q}_c(\tau)^2 - Q^2)^{3/2} \right]_{\tau_i}^{\tau_f} \right\}; & Q_c < |Q| < \hat{Q}_c, \\ \\ a(\tau_0) \left\{ \left[ -2(\tau - \tau_0) + \frac{1}{2}\left\{-3H(\tau_0)a(\tau_0)^2 + g(\tau_0)\right\} (\tau - \tau_0)^2 \right]_{\tau_i}^{\tau_f} Q \right. \\ \left. + \frac{1}{6}\delta''_L(q_0) D_+(\tau_0) \left[ (\tau - \tau_0) + \frac{1}{2}g(\tau_0)(\tau - \tau_0)^2 \right]_{\tau_i}^{\tau_f} Q^3 \right\}; & |Q| < Q_c. \end{cases} \quad (E5) $$

## APPENDIX F: ON THE CHOICE OF PARAMETER IN ADAPTIVE SMOOTHING

In this appendix, we study the effect of varying the parameter $f_{cross}$ on the performances of PT predictions with adaptive smoothing (see Section 4 for definition). In the main text, we adopt the intuitive setup $f_{cross} = 1$ for post-collapse PT and 0.5 for Zel'dovich solution. However, PT predictions are prone, after shell-crossing, to get worse over time, so we can naively expect that choosing a smaller value of $f_{cross}$ will provide a better agreement with simulations.

Fig. E1 shows the results obtained in the case of the CDM-like power spectrum studied in Section 5.4 when varying the value of $f_{cross}$. Here, the calculations are performed using $N_{particle} = 10^5$, $N_{grid} = 10^4$ and $k_{cut} = 6.3 \text{ Mpc}^{-1}$ for a box size $L = 1000 \text{ Mpc}$, which is a relatively
low-resolution set up, but this will not have any consequence on the discussion that follows. Left-hand panel shows the power spectrum at $z=0$ and the resultant predictions for different values of $f_{\text{cross}}$ are plotted as solid lines with different colours. As anticipated, decreasing $f_{\text{cross}}$ results in an enhancement of small-scale power and a better agreement between post-collapse PT and the simulation is obtained. But a too small value of $f_{\text{cross}}$ overshoots the simulation and an optimal choice of $f_{\text{cross}}$ roughly corresponds to $f_{\text{cross}} \sim 0.6$. The same trend can be also seen for the Zel'dovich solution, but the improvement of the power spectrum amplitude is rather mild, while $f_{\text{cross}} = 0.3$ seems to provide the best overall behaviour for the Zel'dovich solution, the result is clearly still not as good as what post-collapse PT can provide, in particular in the regime $0.5 \lesssim k \lesssim 1 \text{ Mpc}^{-1}$.

Right-hand panel of Fig. E1 shows the phase-space portraits and compares the simulation results (red) to predictions obtained with our ‘standard’ value of $f_{\text{cross}} = 1$ for post-collapse PT (upper insert) and $f_{\text{cross}} = 0.5$ for Zel’dovich solution (lower insert) and to predictions obtained with respective nearly optimal values of $f_{\text{cross}} = 0.6$ and 0.3. It seems difficult to judge by eye whether the prediction with a smaller value of $f_{\text{cross}}$ really improves the description at small scales, although more structures seems to be captured for post-collapse PT with $f_{\text{cross}} = 1$, which contradict the results obtained for the power spectrum. However, the simple examination of this figure corresponding to a single realization of the random initial conditions is not conclusive, obviously, it does not preclude the fact that a smaller value of $f_{\text{cross}}$ can give a better result for the power spectrum obtained from averaging over many realizations.

Finally, to conclude this section, although we do not show the results here, we noticed as well that the best choice of $f_{\text{cross}}$ could vary slightly according to redshift or initial condition, for example for the CDM cosmology considered here, a better choice of $f_{\text{cross}}$ at redshift $z = 1.5$ is $f_{\text{cross}} = 0.7$ and 0.4, respectively, for post-collapse PT and Zel'dovich solution. This finally explains why we decided to keep, for simplicity, the generic values $f_{\text{cross}} = 1$ and 0.5, respectively, for post-collapse PT and Zel'dovich solution, even if they are sub-optimal.

This paper has been typeset from a TeX/LA TeX file prepared by the author.