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

HAL
open science

# Shell-crossing in quasi-one-dimensional flow

Cornelius Rampf, Uriel Frisch

► To cite this version:

Cornelius Rampf, Uriel Frisch. Shell-crossing in quasi-one-dimensional flow. Monthly Notices of the Royal Astronomical Society, 2017, 471 (1), pp.671-679. 10.1093/mnras/stx1613. hal-01582412

HAL Id: hal-01582412

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

Submitted on 29 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Shell-crossing in quasi-one-dimensional flow

Cornelius Rampf<fnref n="1" />,<fnref n="2" />* and Uriel Frisch<fnref n="3" />

Accepted 2017 June 21. Received 2017 May 24

## ABSTRACT

Blow-up of solutions for the cosmological fluid equations, often dubbed shell-crossing or orbit crossing, denotes the breakdown of the single-stream regime of the cold-dark-matter fluid. At this instant, the velocity becomes multi-valued and the density singular. Shell-crossing is well understood in one dimension (1D), but not in higher dimensions. This paper is about quasi-one-dimensional (Q1D) flow that depends on all three coordinates but differs only slightly from a strictly 1D flow, thereby allowing a perturbative treatment of shell-crossing using the Euler-Poisson equations written in Lagrangian coordinates. The signature of shell-crossing is then just the vanishing of the Jacobian of the Lagrangian map, a regular perturbation problem. In essence, the problem of the first shell-crossing, which is highly singular in Eulerian coordinates, has been desingularized by switching to Lagrangian coordinates, and can then be handled by perturbation theory. Here, all-order recursion relations are obtained for the time-Taylor coefficients of the displacement field, and it is shown that the Taylor series has an infinite radius of convergence. This allows the determination of the time and location of the first shell-crossing, which is generically shown to be taking place earlier than for the unperturbed 1D flow. The time variable used for these statements is not the cosmic time $t$ but the linear growth time $\tau \sim t^{2/3}$. For simplicity, calculations are restricted to an Einstein-de Sitter universe in the Newtonian approximation, and tailored initial data are used. However it is straightforward to relax these limitations, if needed.

**Key words:** dark matter – large-scale structure of Universe – cosmology: theory.

## 1 INTRODUCTION

It is widely known that exact analytic solutions to the cosmological fluid equations exist for initial data that only depend on one space variable. These play an important role in cosmology, not only because they are simple but also because the breakdown of smooth three-dimensional (3D) solution through the development of infinite-density caustics begins generically as an almost 1D phenomenon with the formation of pancakes (Melott & Shandarin 1989, 1993). When the problem is exactly 1D, the fluid equations become linear when expressed in Lagrangian coordinates. As a consequence, the linear Lagrangian solution solves the fully non-linear problem in 1D (Novikov 1970; Zentsova & Chernin 1980). Expressed in comoving coordinates and enabling the linear growth time $\tau$ (see below for details), the resulting Lagrangian map is in one spatial dimension exactly

$$x(q; \tau) = q + \tau v^{(\text{init})}(q), \quad (1)$$

where $q$ and $v^{(\text{init})}$ are, respectively, the initial position and velocity of the fluid particle. The last term on the right-hand side in equation (1) is the 1D Lagrangian displacement field. It is linear in the time variable, and, evidently, the 1D displacement could be viewed as the first and only non-zero term of an infinite time-Taylor series around $\tau = 0$. Obviously, the 1D Lagrangian map (1) is analytic in $\tau$, has no singularities and thus an infinite radius of convergence. Singularities appear, however, when reverting back to Eulerian coordinates, since the Lagrangian map is not invertible anymore when its Jacobian vanishes for the first time. At this instant, commonly referred to as shell-crossing, the fluid enters the multi-stream regime, which implies that the single-stream fluid description breaks down (in both Lagrangian and Eulerian coordinates).

Generally, departing from 1D leads to a non-zero population of higher-order time-Taylor coefficients of the displacement that should be taken into account. For generic 3D initial data, lower-order solutions of the displacement are well known; see for example the first-order solution, which is called the Zel'dovich approximation (ZA; Zel'dovich 1970); for a generalization of this approximation, see Buchert (1992). Explicit solutions to the second (Bouchet et al. 1992; Buchert & Ehlers 1993), to the third

<footnote n="1">¹Institute for Theoretical Physics (ITP), Philosophenweg 16, University of Heidelberg, D-69120 Heidelberg, Germany</footnote>

<footnote n="2">²Department of Physics, Israel Institute of Technology – Technion, Haifa 32000, Israel</footnote>

<footnote n="3">³Lab. Lagrange, UCA, OCA, CNRS, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote>* E-mail: rampf@physics.technion.ac.il</footnote>
(Buchert 1994; Bouchet et al. 1995), and to the fourth order (Rampf & Buchert 2012) are known as well. Furthermore, truncated approximations for the 3D displacement up to the third order have been applied to numerically extrapolate for the particle trajectories [see e.g. Buchert & Bartelmann (1991), Buchert, Melott & Weiss (1994), Melott, Buchert & Weiss (1995), Buchert et al. (1997), Tassev & Zaldarriaga (2012)].

For generic 3D initial data, the radius of convergence for the time-Taylor series of the displacement field in Lagrangian coordinates is, most likely, not infinite anymore but determined by complex-time singularities, not related to shell-crossing. Analytical bounds on the radius of convergence can be obtained by investigating the large-order behaviour of the infinite time-Taylor series. This, of course, requires explicit all-order recursion relations, as obtained by Zheligovsky & Frisch (2014), who used their recursion relations for the displacement to obtain a lower bound on the radius of convergence. Such a lower bound amounts to finding a time T, as large as possible, such that the time-Taylor series around $\tau = 0$ is guaranteed to converge for $0 \le \tau \le T$. Furthermore, Rampf, Villone & Frisch (2015) have shown that shell-crossing is ruled out in that time-domain. Obtaining the time of first shell-crossing with generic 3D initial data can probably be investigated only by numerical means, for example by employing the multi-time stepping algorithm called the Cauchy-Lagrangian method, an algorithm that has so far been implemented only for incompressible flow (Podvigina, Zheligovsky & Frisch 2016).

In the present paper, we show that much more can be handled analytically when restricting the initial data to being close to 1D. We prove that the relevant Lagrangian map has infinitely many non-vanishing terms in its time-Taylor series, but that the series is entire in time, that is, it has an infinite radius of convergence. Thus, particle trajectories can be evaluated in a single time-step from initial time all the way up to the first shell-crossing (but not beyond, for reasons that will be discussed later). To unravel the above, we make use of novel recursion relations which are, most importantly and crucially, tailored to initial data that are perturbatively close to 1D. By contrast, all-order recursive solutions for generic 3D initial data (see also Goroff et al. 1986; Ehlers & Buchert 1997; Matsubara 2015; Rampf et al. 2015) are not suitable for the considered problem. Indeed, the usage of Q1D initial data introduces another expansion parameter, namely a parameter which parametrizes the perturbative departure from 1D. Thus, the power counting in the perturbative expansion in Q1D is formally different from the generic 3D case. For this reason, our results are not contained in the commonly used Lagrangian perturbation solutions.

This paper is organized as follows. In Section 2 we review the 3D Euler-Poisson equations, first in the Eulerian- and then in the Lagrangian-coordinates approach. The latter approach serves as our starting point for the present paper. In Section 3, we show how to embed the Q1D problem into three-dimensional space, and particularly discuss the used initial conditions (ICs) and perturbation Ansatz. The resulting equations can be easily solved to a given order in the book-keeping perturbation variable $\epsilon$. The zeroth-order solution in $\epsilon$, which we call the solution of the *unperturbed problem*, is the one for which the initial data depend only on one space variable (i.e. the 1D case). The first-order equations to order $\epsilon$, which resemble the perturbed equations with respect to the unperturbed problem, are given in Section 3.4. We solve this *perturbed problem* by using a time-Taylor series in Section 4. The proof of the absence of singularities in the perturbed Lagrangian equations is given in Section 5. In Section 6 we show how to obtain the time and location of the first shell-crossing. In Section 7 we give a concrete example

involving a three-sine wave Q1D IC. Concluding remarks are presented in Section 8.

## 2 EULER-POISSON EQUATIONS IN 3D

### 2.1 Basic equations in Eulerian coordinates

The Euler-Poisson equations are usually formulated in comoving coordinates $x = r/a$, where $r$ is the proper space coordinate and $a$ the cosmic scale factor. The latter parametrizes the global background/Hubble expansion, and its evolution is given by the usual Friedmann equations. In the present work, we restrict our analysis, for simplicity, to an Einstein-de Sitter (EdS) cosmology, where the universe is filled only with a cold-dark-matter (CDM) fluid; the cosmological constant, usually denoted by $\Lambda$, is set to zero. This and many other approximations are, however, easily rectified if needed [see e.g. Rampf et al. (2015) for an analysis within the $\Lambda$CDM model and beyond].

We denote by $v$ the peculiar velocity with respect to the Hubble flow, by $\delta = (\rho - \bar{\rho})/\bar{\rho}$ the density contrast with background density $\bar{\rho} \sim a^{-3}$, by $\phi_g$ the cosmological potential, and by $\tau$ the linear growth time [often denoted with $D$ or $D(t)$]. For an EdS universe, $\tau$ is related to the cosmic time $t$ via $\tau \sim t^{2/3}$. The Euler-Poisson equations for an EdS universe are (Brenier et al. 2003; for more general cosmologies see Shandarin 1992, 1994)

$$ \partial_{\tau} v + (v \cdot \nabla) v = -\frac{3}{2\tau} (v + \nabla \phi_g), \qquad (2a) $$

$$ \partial_{\tau} \delta + \nabla \cdot [(1 + \delta)v] = 0, \qquad (2b) $$

$$ \nabla^2 \phi_g = \frac{\delta}{\tau}, \qquad (2c) $$

(here $\nabla^2$ is the Eulerian Laplacian). Enabling $\tau$ as the time variable is convenient when studying the well-posedness of the fluid equations at short times, as explained hereafter. It is actually essential when investigating the time-analyticity of the Lagrangian map (Zheligovsky & Frisch 2014; Rampf et al. 2015).

Formally linearizing around the steady state $v = 0$ and $\delta = 0$, the above equations can be written in terms of a single differential equation for the density contrast (Peebles 1980),

$$ \partial_{\tau\tau}^2 \delta = -\frac{3}{2\tau} \left( \partial_{\tau} \delta - \frac{\delta}{\tau} \right). \qquad (3) $$

This equation has two power-law solutions. One is called the decaying mode which behaves as $\tau^{-3/2}$ and thus blows up when $\tau \to 0$, thereby invalidating linearization. The other one is linear in $\tau$ and thus trivially analytic; furthermore it stays analytic in the presence of a non-vanishing cosmological constant (Rampf et al. 2015). Therefore, $\tau$ is the physically appropriate variable for describing the growth of density fluctuations at short times.

Before investigating the fully non-linear theory in Lagrangian coordinates, let us briefly discuss an important feature of the Euler-Poisson equations. Observe the presence of the linear growth time in the denominators of the right-hand side of equations (2a) and (2c). This indicates that the solution at $\tau = 0$ is singular unless the following *slaving conditions* are satisfied [denoted by the superscript (init)]:

$$ \delta^{(\text{init})} = 0, \qquad v^{(\text{init})} = -\nabla \phi_g^{(\text{init})}. \qquad (4) $$
As argued by Rampf et al. (2015), the second condition immediately implies that the velocity is potential at initial time, a feature that, because of (2a), persists also at later times in Eulerian coordinates,

$$ \nabla \times \mathbf{v} = 0. \qquad (5) $$

(In Lagrangian coordinates, the potential character of the velocity is generally lost.) Finally, we note that in our mathematical description, we allow the time variable $\tau$ to become arbitrarily small. From the point of view of boundary layer analysis, this amounts to subsuming the whole primordial physics into the slaving conditions [see Zheligovsky & Frisch (2014); Rampf et al. (2015) for detailed discussions]. Nevertheless, initialization times with $\tau^{(\text{init})} > 0$ can be incorporated in our formalism which, however, requires some additional calculational steps; we shall come back to this issue in Section 8.

## 2.2 Basic equations in Lagrangian coordinates

We denote by $\mathbf{q}$ the Lagrangian coordinates with components $q_i$ ($i = 1,2,3$); a partial derivative with respect to $q_i$ acting on a given function $f$ is denoted by $f_i$, and, occasionally, by $\partial_i^L f$. Summation over repeated indices is implied, and, for simplicity, since we work in the Newtonian limit, we do not distinguish between contra- and covariant coordinate indices. Let $\mathbf{q} \mapsto \mathbf{x}(\mathbf{q}; \tau)$ be the Lagrangian map from the initial ($\tau = 0$) position $\mathbf{q}$ to the Eulerian position $\mathbf{x}$ at time $\tau$. The map satisfies $\mathbf{v}(\mathbf{x}(\mathbf{q}; \tau); \tau) = \partial_{\tau}^L \mathbf{x}(\mathbf{q}; \tau)$, where $\partial_{\tau}^L$ is the Lagrangian time derivative – the latter also denoted with an overdot in the following. At initial time, $\tau = 0$, the velocity is

$$ \mathbf{v}^{(\text{init})}(\mathbf{q}) = \mathbf{v}(\mathbf{x}(\mathbf{q}; 0); 0), \qquad (6) $$

which agrees with the initial Eulerian velocity. Mass conservation is, until the first shell-crossing, given by

$$ \delta = 1/J - 1, \qquad (7) $$

where $J = \det(x_{i,j})$, the determinant of the Jacobian matrix with entries $x_{i,j}$ is called the Jacobian (as long as it is non-negative). With these definitions, the Euler-Poisson equations can be written in Lagrangian coordinates in the compact form

$$ \varepsilon_{ikl} \varepsilon_{jmn} x_{k,m} x_{l,n} \Re \tau x_{i,j} = 3(J-1), \qquad (8a) $$

$$ \varepsilon_{ijk} \dot{x}_{l,j} x_{l,k} = 0, \qquad (8b) $$

where we have defined the operator $\Re_{\tau} \equiv \tau^2 (\partial_{\tau}^L)^2 + (3\tau/2)\partial_{\tau}^L$, and $\varepsilon_{ijk}$ is the fundamental antisymmetric tensor. Equation (8a) is a scalar equation that results from combining equations (2a) and (2c), and by taking mass conservation (7) into account (for a derivation see e.g. Rampf et al. 2015, and in there set $\Lambda = 0$). Equations (8b) are the Cauchy invariants; these are Lagrangian (kinematical) constraints on the Lagrangian map that must be satisfied in order to maintain the curlfree motion in Eulerian space. The Cauchy invariants can thus be understood as the corresponding Lagrangian counterpart of equation (5). See e.g. Rampf et al. (2016) for a detailed derivation of the Cauchy invariants.

Equations (8) constitute the well-known closed system of Lagrangian equations for CDM. As can be easily checked, power-series solutions for these equations are non-singular at $\tau = 0$, provided one makes use of the slaving conditions (4). Actually, such expansions in powers of the linear growth time $\tau$ are very common in the Lagrangian perturbation theory [see e.g. Buchert (1994); Matsubara (2008); Rampf & Buchert (2012); Zheligovsky & Frisch (2014)].

# 3 THE QUASI-ONE-DIMENSIONAL PROBLEM EMBEDDED IN 3D

The aim of this paper is to analyse three-dimensional shell-crossing with ICs that are close to 1D, i.e. the ICs depend, to the zeroth order in a perturbation parameter $\epsilon$, only on one space variable, and, to first order in $\epsilon$, in general on all space variables. Appropriate ICs and our perturbation Ansatz are introduced in the following two sections. Equations to order $\epsilon^0$ and $\epsilon^1$ are then given in Sections 3.3 and 3.4, respectively.

## 3.1 Initial conditions

Quasi-1D initial conditions can be formulated in terms of a superposition of two contributions for the initial gravitational potential. The first is an arbitrary function in the space variable $q_1$ and the second one a small perturbation, proportional to $\epsilon$, which depends generally on all space variables. Although the former, which characterizes the ICs for the purely 1D problem, could be taken quite arbitrarily (within a class of function guaranteeing well-posedness for at least a finite time), it is advantageous to choose this function wisely: We know that in the 1D case, the occurrence of the first shell-crossing will appear downstream at the spatial position $q_1$ where the initial velocity gradient achieves its most negative value. By a suitable spatial translation, we can take this location to be $q_1 = 0$. Then, by a suitable Galilean transformation, we can take the velocity at this location to be zero. An instance is to take for this 1D gravitational potential the function $-\cos q_1$. Considerations of normal-form reduction indicate that we can actually make this choice of ICs without loss of generality (within the class of $2\pi$-periodic functions).

As to the perturbation, it must be taken fairly general. We thus use the initial data

$$ \varphi_{g}^{(\text{init})}(q_1, q_2, q_3) = -\cos q_1 + \epsilon \varphi^{(\text{init})}(q_1, q_2, q_3), \qquad (9) $$

where $\epsilon > 1$ is a small perturbation parameter and $\varphi^{(\text{init})}$ an arbitrary $2\pi$-periodic function of $(q_1, q_2, q_3)$. Without loss of generality, we can take the transverse average to be zero, namely

$$ \langle \varphi^{(\text{init})} \rangle = \int_{0}^{2\pi} \frac{dq_2}{2\pi} \int_{0}^{2\pi} \frac{dq_3}{2\pi} \varphi^{(\text{init})}(q_1, q_2, q_3) = 0. \qquad (10) $$

Indeed, if this average is a non-trivial function of $q_1$, we can incorporate it into the unperturbed flow.

Taking into account the equality of the initial gravitational and velocity potentials (imposed by slaving), the initial velocity is, in index notation,

$$ \mathbf{v}_{i}^{(\text{init})}(\mathbf{q}) = -\delta_{i1} \sin q_1 - \epsilon \partial_{\tau}^{L} \varphi^{(\text{init})}(\mathbf{q}), \qquad (11) $$

where $\delta_{ij}$ is the Kronecker delta.

## 3.2 The Lagrangian perturbation Ansatz

We are going to use a perturbation method in which the solutions to the Lagrangian equations (8) are expanded in powers of the small parameter $\epsilon$. Namely, we look for a solution in which the Lagrangian map is given by the perturbation Ansatz

$$ \mathbf{x}(\mathbf{q}; \tau) = \mathbf{q} + \xi^{(0)}(\mathbf{q}; \tau) + \epsilon \xi^{(1)}(\mathbf{q}; \tau) + \epsilon^2 \xi^{(2)}(\mathbf{q}; \tau) + \dots, \qquad (12) $$

where $\xi^{(n)}(\mathbf{q})$ is the coefficient of $\epsilon^n$ in the expansion of the displacement $\mathbf{x} - \mathbf{q}$. For $\epsilon = 0$, we are back to exactly 1D; the displacement
depends only on $q_1$ and is in the direction of the first coordinate axis. We can thus write

$$ \xi_i^{(0)}(\boldsymbol{q}) = \delta_{i1} F(q_1; \tau). \qquad (13) $$

In this paper, we do not expand beyond first order in $\epsilon$ and, for brevity, we write $\xi^{(1)}(\boldsymbol{q}; \tau) = \xi(\boldsymbol{q}; \tau)$. Henceforth all calculations are extended only to first order in $\epsilon$. For example, from equations (12) and (13), we have thus

$$ x_i(\boldsymbol{q}; \tau) = q_i + \delta_{i1} F(q_1; \tau) + \epsilon \xi_i(\boldsymbol{q}; \tau). \qquad (14) $$

From equation (14) it follows that the Jacobian matrix is given by

$$ x_{i,j} = \delta_{i,j} + \delta_{i1}\delta_{j1}F_{,1} + \epsilon\xi_{i,j}, \qquad (15) $$

so that, ignoring $O(\epsilon^2)$ terms, its determinant, the Jacobian, is given by

$$ J = 1 + F_{,1} + \epsilon (\xi_{1,1} + \xi_{2,2} + \xi_{3,3}) + \epsilon F_{,1} (\xi_{2,2} + \xi_{3,3}). \qquad (16) $$

The vanishing of the Jacobian is evidence of shell-crossing. In the following we solve for the displacement and evaluate the Jacobian at a given order in the $\epsilon$-perturbation. To zeroth order the problem is effectively 1D, with the unperturbed Jacobian $J^{(0)} = 1 + F_{,1}$. Since the unperturbed problem is embedded in the three-dimensional space (see Buchert & Goetz 1987), the first shell-crossing will occur at a given value $q_{*}^{(0)} = 0$ (modulo $2\pi$) and in the two-dimensional plane that is spanned by the 2-tuple $(q_2, q_3)$. Switching on the $\epsilon$-perturbation, the two-dimensional plane that marked the instance of shell-crossing in the unperturbed problem will collapse to a given point $(q_{1*}, q_{2*}, q_{3*})$ (modulo $2\pi$). Also, the time of the first shell-crossing will generally change when switching from the unperturbed to the perturbed problem. Analytical and calculational details to the respective values, $\tau_{*}^{(0)}$ and $\tau_{*}$, are discussed in the following sections.

### 3.3 Lagrangian equations to zeroth order in $\epsilon$

To lowest (zeroth) order in $\epsilon$, the problem is exactly 1D and has the IC $\varphi_g^{(\text{init})}(q_1, q_2, q_3) = -\cos q_1$. In this instance, as is well known, the ZA is exact until shell-crossing (Novikov 1970; Zel'dovich 1970). In Lagrangian coordinates, this is particularly obvious. Indeed, substituting Ansatz (14) into (8b) gives a trivial identity, whereas for equation (8a), we obtain

$$ 2\Re_{\tau} F_{,1} = 3F_{,1}. \qquad (17) $$

This linear equation has, on the one hand, the decaying solution $\sim \tau^{-2/3}$, which is unbounded as $\tau \to 0$ and not acceptable and, on the other hand the growing solution

$$ F_{,1}(q_1; \tau) = \tau \chi_{,1}(q_1), \qquad (18) $$

where $\chi(q_1)$ can be taken arbitrarily, but by (11) is actually equal to the initial velocity

$$ \chi(q_1) = -\sin q_1. \qquad (19) $$

To zeroth order in $\epsilon$, the particle trajectory is thus

$$ x_i = q_i - \delta_{i1} \tau \sin q_1, \qquad (20) $$

and the unperturbed Jacobian is

$$ J^{(0)} = 1 - \tau \cos q_1. \qquad (21) $$

The instance of first shell-crossing is determined by the lowest time-value for which $J^{(0)}(q_{1*}^{(0)}, \tau_{*}^{(0)}) = 0$. It is easily checked that this happens at $\tau_{*}^{(0)} = 1$ and within the two-dimensional $(q_2, q_3)$-plane at $q_{1*}^{(0)} = 0$ (modulo $2\pi$).

Note that the particle trajectory (20) is a linear function in time $\tau$. Hence its time-Taylor series has just two terms and, trivially, has an infinite radius of convergence, i.e. it is an entire function. A singularity for the 1D case appears only when reverting to Eulerian coordinates. Indeed, when the Jacobian vanishes, the Lagrangian map ceases to be smoothly and uniquely invertible. In a purely Eulerian formulation, by contrast, infinities appear explicitly when the fluid enters the multi-stream regime. Actually, the Eulerian time-Taylor series for both the density and the velocity have a singularity at shell-crossing. As we shall now see, this picture remains basically unchanged, when we switch on the perturbation – which makes the problem fully three-dimensional.

### 3.4 Lagrangian equations to first order in $\epsilon$

Collecting all terms $O(\epsilon)$, we obtain from equation (8a)

$$ \Re_{\tau} \xi_{1,1} + [1 - \tau \cos q_1] \Re_{\tau} (\xi_{2,2} + \xi_{3,3}) = \frac{3}{2} \xi_{1,1}, \qquad (22a) $$

and from the three components of (8b), i.e. $i = 1, 2, 3$, respectively,

$$ \dot{\xi}_{2,3} - \dot{\xi}_{3,2} = 0, \qquad (22b) $$

$$ [1 - \tau \cos q_1] \dot{\xi}_{1,3} + \xi_{1,3} \cos q_1 = \dot{\xi}_{3,1}, \qquad (22c) $$

$$ [1 - \tau \cos q_1] \dot{\xi}_{1,2} + \xi_{1,2} \cos q_1 = \dot{\xi}_{2,1}. \qquad (22d) $$

Observe that equations (22) are, by construction, the Lagrangian Euler-Poisson equations (8), linearized around the exact 1D solution (20). Since the latter depends explicitly on $q_1$, so do the linearized equations. But there is no explicit dependence on $q_2$ and $q_3$. As a consequence, if we are able to solve equations (22) for the case where the initial potential perturbation is $\varphi(q_1, q_2, q_3) = \bar{\varphi}(q_1) \exp[i(k_2q_2 + k_3q_3)]$, i.e. with a single (transverse) Fourier harmonic, then we can handle the general case by linear superposition. In the single-harmonic case, derivatives with respect to $q_2$ and $q_3$ can be replaced with $ik_2$ and $ik_3$, respectively. Using this and defining

$$ \chi \equiv \xi_1, \quad \zeta \equiv ik_2\xi_2 + ik_3\xi_3, \quad k_{\perp}^2 \equiv k_2^2 + k_3^2, \qquad (23) $$

we are left with only two unknowns, and equations (22) reduce to just two equations, which are

$$ \Re_{\tau} \chi_{,1} + [1 - \tau \cos q_1] \Re_{\tau} \zeta = \frac{3}{2} (\chi_{,1} + \zeta), \qquad (24a) $$

$$ k_{\perp}^2 \{[1 - \tau \cos q_1] \dot{\chi} + \chi \cos q_1\} = -\dot{\zeta}_{,1}. \qquad (24b) $$

This is the basic set of linearized Lagrangian equations, which we solve in the following section.

## 4 TAYLOR EXPANSION AND RECURSION RELATIONS

We observe that the linearized equations (24) constitute a system of two linear partial differential equations in the variables $\tau$ (second order) and $q_1$ (first order), in which the transverse coordinates appear only parametrically through the wavenumber $k_{\perp}$. Our method of solution will use time-Taylor expansions to arbitrary high order, based on novel recursion relations for the Taylor coefficients.

For this, we seek a solution to (24) in the form of a Taylor series in the $\tau$-time for the displacement components,

$$ \chi(\boldsymbol{q}; \tau) = \sum_{n=1}^{\infty} \chi^{(n)}(\boldsymbol{q}) \tau^n, \quad \zeta(\boldsymbol{q}; \tau) = \sum_{n=1}^{\infty} \zeta^{(n)}(\boldsymbol{q}) \tau^n. \qquad (25) $$
Substituting this Ansatz into equations (24), and collecting all the terms containing a given power in $\tau^n$, yields the following relations between the time-Taylor coefficients,

$$ \left[n^2 + \frac{n}{2} - \frac{3}{2}\right] (\chi_{,1}^{(n)} + \xi^{(n)}) = \left[n^2 - \frac{3n}{2} + \frac{1}{2}\right] \xi^{(n-1)} \cos q_1, \quad (26a) $$

$$ n \chi^{(n)} + k_{\perp}^{2} n \xi_{,1}^{(n)} = (n-2) \chi^{(n-1)} \cos q_1. \quad (26b) $$

Here and in the following, by construction, coefficients vanish if their index is zero or negative. Equations (26) can be simplified by Fourier transforming also in the $q_1$ variable. For this we define

$$ \chi^{(n)} = \hat{\chi}_{k_1}^{(n)} e^{ik_1 q_1}, \quad \xi^{(n)} = \hat{\xi}_{k_1}^{(n)} e^{ik_1 q_1}, \quad (27) $$

and then, making use of Euler's formula $\cos q_1 = (\exp\{iq_1\} + \exp\{-iq_1\})/2$, we obtain for equations (26)

$$ \left[n^2 + \frac{n}{2} - \frac{3}{2}\right] (ik_1 \hat{\chi}_{k_1}^{(n)} + \hat{\xi}_{k_1}^{(n)}) = \left[n^2 - \frac{3n}{2} + \frac{1}{2}\right] \times \frac{1}{2} \left( \hat{\xi}_{k_1+1}^{(n-1)} + \hat{\xi}_{k_1-1}^{(n-1)} \right), \quad (28a) $$

$$ n \hat{\chi}_{k_1}^{(n)} + k_{\perp}^{2} n i k_1 \hat{\xi}_{k_1}^{(n)} = \frac{n-2}{2} \left( \hat{\chi}_{k_1+1}^{(n-1)} + \hat{\chi}_{k_1-1}^{(n-1)} \right). \quad (28b) $$

For $n=1$, the first of these equations amounts to an identity, and the last equation gives

$$ \hat{\chi}_{k_1}^{(1)} = -ik_1 k_{\perp}^{-2} \hat{\xi}_{k_1}^{(1)}, \quad (29) $$

and, together with the definition of the Lagrangian map and (11), thus

$$ \hat{\chi}_{k_1}^{(1)} = -ik_1 \hat{\phi}^{(\text{init})}, \quad \hat{\xi}_{k_1}^{(1)} = k_{\perp}^2 \hat{\phi}^{(\text{init})}. \quad (30) $$

For $n > 1$, we obtain from equations (28) the following explicit recursion relations:

$$ \hat{\xi}_{k_1}^{(n)} = \left(1 + \frac{k_1^2}{k_{\perp}^2}\right)^{-1} \left[ \frac{n-1/2}{2n+3} \left( \hat{\xi}_{k_1+1}^{(n-1)} + \hat{\xi}_{k_1-1}^{(n-1)} \right) - ik_1 \frac{n-2}{2n} \left( \hat{\chi}_{k_1+1}^{(n-1)} + \hat{\chi}_{k_1-1}^{(n-1)} \right) \right], \quad (31a) $$

$$ \hat{\chi}_{k_1}^{(n)} = \left(1 + \frac{k_1^2}{k_{\perp}^2}\right)^{-1} \left[ -ik_1 k_{\perp}^{-2} \frac{n-1/2}{2n+3} \left( \hat{\xi}_{k_1+1}^{(n-1)} + \hat{\xi}_{k_1-1}^{(n-1)} \right) + \frac{n-2}{2n} \left( \hat{\chi}_{k_1+1}^{(n-1)} + \hat{\chi}_{k_1-1}^{(n-1)} \right) \right]. \quad (31b) $$

We then construct the nth-order time-Taylor coefficient of the displacement field using a Helmholtz-Hodge decomposition. The latter reads in Fourier space ($n \ge 1$):

$$ \hat{\xi}_{k_1}^{(n)} = -(k_1^2 + k_2^2 + k_3^2)^{-2} (ik \cdot [ik \cdot \hat{\xi}_{k_1}^{(n)}] - ik \times \hat{T}_{k_1}^{(n)}), \quad (32) $$

with

$$ i \mathbf{k} \cdot \hat{\xi}_{k_1}^{(n)} = ik_1 \hat{\chi}_{k_1}^{(n)} + \hat{\xi}_{k_1}^{(n)}, \quad (33) $$

$$ \hat{T}_{k_1}^{(n)} = i(0, k_3, -k_2)^T (\hat{\chi}_{k_1}^{(n)} + ik_1 k_{\perp}^{-2} \hat{\xi}_{k_1}^{(n)}). \quad (34) $$

The respective right-hand sides of the two last equations are combinations of time-Taylor coefficients to order $n$. By virtue of equations (28), however, these time-Taylor coefficients can be written in terms of the lower-order time-Taylor coefficients $n-1$. We

thus can construct the time-Taylor coefficients of the displacement in a recursive way. We find in Fourier space ($n \ge 1$)

$$ \hat{\xi}_{k_1}^{(n)} = -ik \hat{\phi}^{(\text{init})} \delta_{1n} - (k_1^2 + k_2^2 + k_3^2)^{-2} \left( ik \frac{n-1/2}{2n+3} \left[ \hat{\xi}_{k_1+1}^{(n-1)} + \hat{\xi}_{k_1-1}^{(n-1)} \right] + \frac{n-2}{2n} (-k_{\perp}^2, k_1 k_2, k_1 k_3)^{\text{T}} \left[ \hat{\chi}_{k_1+1}^{(n-1)} + \hat{\chi}_{k_1-1}^{(n-1)} \right] \right), \quad (35) $$

and in real space

$$ \hat{\xi}^{(n)} = -\nabla^{\text{L}} \phi^{(\text{init})} \delta_{1n} + \nabla_{\text{L}}^{-2} \left( \frac{2n-1}{2n+3} \nabla^{\text{L}} [\xi^{(n-1)} \cos q_1] \right) + \frac{n-2}{n} \left( \partial_{\perp}^{\text{L}} \partial_{\perp}^{\text{L}} + \partial_{\perp}^{\text{L}} \partial_{\perp}^{\text{L}}, -\partial_{\perp}^{\text{L}} \partial_{\perp}^{\text{L}}, -\partial_{\perp}^{\text{L}} \partial_{\perp}^{\text{L}} \right)^{\text{T}} [\chi^{(n-1)} \cos q_1], \quad (36) $$

where $\nabla_{\text{L}}^{-2}$ is the inverse Laplacian in Lagrangian coordinates.

From this we obtain, to order $\epsilon$ and to all orders in time, respectively, the particle trajectory and then, using (16), the Jacobian,

$$ x_i(\mathbf{q}; \tau) = q_i - \delta_{1i} \tau \sin q_1 + \epsilon \sum_{n=1}^{\infty} \xi_i^{(n)}(\mathbf{q}) \tau^n, \quad (37) $$

$$ J = 1 - \tau \cos q_1 + \epsilon \sum_{n=1}^{\infty} \left( \chi_{,1}^{(n)} + [1 - \tau \cos q_1] \xi^{(n)} \right) \tau^n, \quad (38) $$

where $\chi_{,1}^{(n)} = \xi_{1,1}^{(n)}$ and $\xi^{(n)} = \xi_{2,2}^{(n)} + \xi_{2,3}^{(n)}$. Equations (37) and (38) constitute the main technical results of this paper.

In the following section, we show that these formal solutions are actually convergent series and free of any singularities. Then, in Section 6, we comment on $\tau_*$, the time of first shell-crossing in the perturbed problem.

## 5 NO SINGULARITIES IN LAGRANGIAN SOLUTIONS

Observe that in equation (24a), there is a term $[1 - \tau \cos q_1]\Re_\tau \zeta$, involving the second-order time derivative of one of the unknowns, $\zeta$, whose coefficient is the Jacobian of the unperturbed problem, i.e. $J^{(0)} = 1 - \tau \cos q_1$. This term vanishes at the first shell-crossing for the unperturbed (1D) problem. As it is known from e.g. the Fuchsian theory (Moser 1959), the vanishing of the coefficient in front of the highest-derivative term may easily lead to a singularity (at least for ODEs). Should this happen here, we would have to face a singular perturbation problem. Fortunately, this is not the case, and the Lagrangian map determined by the first-order perturbation equations (24) is an entire function of time, as we now show. Here, we must stress that after shell-crossing, because of multi-streaming, the true Lagrangian map ceases to be governed by the Euler-Poisson equations (2), but this does not matter for the determination of the first shell-crossing.

It is easily shown that the entire character of the time-Taylor series is related to the behaviour at large orders $n$ of the Taylor coefficients $\chi^{(n)}(\mathbf{q})$ and $\xi^{(n)}(\mathbf{q})$. These satisfy the recursion relations (31), which can be simplified for large $n$ and approximated by their asymptotic form:

$$ \chi_{,1}^{(n)} + \xi^{(n)} = \zeta^{(n-1)} \cos q_1, \quad (39a) $$

$$ \chi^{(n)} + k_{\perp}^{-2} \xi_{,1}^{(n)} = \chi^{(n-1)} \cos q_1. \quad (39b) $$
The large-$n$ recursion relations (39) can actually be solved explicitly. Paradoxically, to achieve this, we shall return to $\tau$ space rather than working directly with the time-Taylor coefficients and their recursion relations. However, we shall not work with the full Taylor series but only with lower-truncated series, i.e.

$$ \chi_N(\mathbf{q}; \tau) = \sum_{n=N}^{\infty} \chi^{(n)}(\mathbf{q}) \tau^n, \quad \zeta_N(\mathbf{q}; \tau) = \sum_{n=N}^{\infty} \zeta^{(n)}(\mathbf{q}) \tau^n, \quad (40) $$

where $N$ is taken large enough to be able to use the asymptotic form (39) of the recursion relations. Comparison of equations (25) and (40) shows that $\chi_N(\mathbf{q}; \tau)$ and $\chi(\mathbf{q}; \tau)$ differ by a polynomial in $\tau$ (with $\mathbf{q}$-dependent coefficients) of degree $N-1$. The same statement holds for $\chi_N(\mathbf{q}; \tau)$ and $\chi(\mathbf{q}; \tau)$. As a consequence, it is equivalent to show the entire character in $\tau$ of the pair $(\chi(\tau), \zeta(\tau))$ or of the pair $(\chi_N(\tau), \zeta_N(\tau))$.

At this point, to simplify the calculations, and without loss of generality, we can set $q_\perp = 1$ (if not, rescale the transverse variables $q_\perp \equiv (q_2, q_3)$ and the amplitude of the first component of the perturbation $\chi \equiv \xi_1$ suitably).

We now multiply (39a) and (39b) by $\tau^n$ and sum on $n$ from $N$ to infinity, to obtain the following equations:

$$ \chi_{N,1} + (1 - \tau \cos q_1) \zeta_N = \tau^N \cos q_1 \zeta^{(N-1)}, \quad (41a) $$

$$ (1 - \tau \cos q_1) \chi_N + \zeta_{N,1} = \tau^N \cos q_1 \chi^{(N-1)}. \quad (41b) $$

The functions $\chi_N$ and $\zeta_N$ are coupled by a system of two equations. However, by simply introducing the sum and the difference

$$ Z_N^{\pm} \equiv \chi_N \pm \zeta_N, \quad (42) $$

we obtain two decoupled equations:

$$ Z_{N,1}^{\pm} = \mp[1 - \tau \cos q_1]Z_N^{\pm} + F_N^{\pm}(q_1; \tau), \quad (43) $$

where

$$ F_N^{\pm}(q_1; \tau) = \tau^N \cos q_1 (\zeta^{(N-1)} \pm \chi^{(N-1)}). \quad (44) $$

Observe that equations (43) are first-order ordinary differential equations in $q_1$ in which the time appears just as a parameter; indeed, thanks to the large-$n$ asymptotics, $\tau$-derivatives have dropped out. Also observe that $F_N^{\pm}(q_1; \tau)$ are polynomials in $\tau$ and trigonometric polynomials in $q_1$, and thus entire functions of $\tau$ and $q_1$.

We now show that equations (43) have unique solutions within the class of functions of $q_1$ that are $2\pi$-periodic. We begin with the case of $Z_N^+(q_1; \tau)$. For this, we first consider the (left spatial) initial value problem for which we prescribe the value of $Z_N^+$ for some initial value $q_{in}$ of $q_1$ and look for the solution to its right, i.e. for $q \ge q_{in}$. This has the explicit solution

$$ Z_N^+(q_1; \tau) = G^+(q_1, q_{in}; \tau) Z_N^+(q_{in}) + \int_{q_{in}}^{q_1} dq' G^+(q_1, q'; \tau) F_N^+(q'; \tau), \quad (45) $$

in terms of Green's function $G^+$ of the associated linear differential equation (without the $F_N^+$ term), given, for $q_1 \ge q'$, by

$$ \begin{aligned} G^+(q_1, q'; \tau) &= \exp \left\{ -\int_{q'}^{q_1} dq'' \left[1 - \tau \cos q''\right] \right\} \\ &= \exp \left\{ \tau \left( \sin q_1 - \sin q' \right) - (q_1 - q') \right\}. \end{aligned} \quad (46) $$

The solution (45) is in general not periodic in $q_1$, but a periodic solution can be constructed by letting $q_{in} \to -\infty$, because far to the right of the point $q_{in}$, the solution relaxes to spatial periodicity. This is proved by decomposing the interval $[q_{in}, q_1]$ into adjacent

intervals all-but-the-first (which may be smaller) of length $2\pi$. The number $M$ of intervals of length $2\pi$ is the integer part of $|q_1 - q_{in}|/(2\pi)$ and thus tends to infinity when $q_{in} \to -\infty$. We observe that, when the argument $q'$ of the integrand on the right-hand side of (45) is shifted from one interval to the neighbouring left interval by subtracting $2\pi$, the integrand is multiplied by a factor $\exp(-2\pi)$. Indeed, the only $q'$-dependent term which is not periodic in $q'$ is $(q_1 - q')$ in the exponential, which generates the stated factor. Hence, the sum over all the $M$ intervals produces a geometric series of ratio $\exp(-2\pi)$. As $q_{in} \to -\infty$ at fixed $Z_N^+(q_{in})$, the first term on the right-hand side of (45) tends to zero, and the second term is given by the sum of an infinite convergent geometric series, namely

$$ Z_N^+(q_1; \tau) = \frac{1}{1 - e^{-2\pi}} \int_{q_1 - 2\pi}^{q_1} dq' F_N^+(q'; \tau) \times \exp\{\tau(\sin q_1 - \sin q') - (q_1 - q')\}. \quad (47) $$

It is easily checked that (47) is a $2\pi$-periodic solution of (43), and, furthermore, is the unique one. Indeed, let $Z_N^+(q_1; \tau)$ and $\tilde{Z}_N^+(q_1; \tau)$ be two such solutions. Their difference $\Delta(q_1; \tau)$ satisfies the homogeneous equation

$$ \Delta_{,1} = (\tau \cos q_1 - 1)\Delta, \quad (48) $$

which implies that

$$ \Delta(q_1 + 2\pi; \tau) = e^{-2\pi} \Delta(q_1; \tau). \quad (49) $$

Since $Z_N^+(q_1; \tau)$ and $\tilde{Z}_N^+(q_1; \tau)$ are $2\pi$-periodic, so is their difference $\Delta$, which by (49) vanishes, hence we confirm the uniqueness of the solution. Furthermore, using the fact that $F_N^+(q'; \tau)$ is polynomial in $\tau$, we easily check that, for any real $q_1$, $Z_N^+(q_1; \tau)$, given by (47), is an entire function of $\tau$.

The case of $Z_N^-(q_1; \tau)$, which satisfies (42) with the minus sign, is handled similarly, except that we must replace the left spatial initial value problem by a right spatial initial value problem, where we seek the solution for $q \le q_{in}$ (or, equivalently, we can just change $q_1$ into $-q_1$). Hence $Z_N^-(q_1; \tau)$ is also an entire function of $\tau$. As a consequence $\chi(q_1; \tau)$ and $\zeta(q_1; \tau)$ are, for any real $q_1$, entire functions of $q_1$. It is also easily shown that they are also entire functions of $\tau$.

## 6 THE TIME OF PERTURBED SHELL-CROSSING

In the absence of perturbations, when the flow is exactly 1D and with our IC, shell-crossing happens at the time $\tau_*^{(0)} = 1$ and location $q_{1*} = 0$, and for arbitrary $q_2$ and $q_3$. Thus, the whole plane $q_1 = 0$ shell-crosses at $\tau_*^{(0)} = 1$. When the perturbation is switched on, translation invariance in the directions of $q_2$ and $q_3$ is broken and, generically, shell-crossing takes place at a time $\tau_* \ne 1$ and at a single location $(q_1, q_2, q_3)$. We shall now show that $\tau_*^*$ is generically happening earlier than $\tau_*^{(0)} = 1$ and explain how the precise time and location can be obtained.

For this purpose we will assume that the initial perturbation is a finite-order trigonometric polynomial

$$ \epsilon \phi^{(\text{init})} = \epsilon \sum_{m=1}^{M} \sum_{n=1}^{N} [a_{mn}^{(\text{init})}(q_1) \cos(mq_2 + nq_3) + b_{mn}^{(\text{init})}(q_1) \sin(mq_2 + nq_3)]. \quad (50) $$

Infinite Fourier series can also be handled, but this requires some functional analysis which we would rather avoid here.

Because of the linearity of (22a)–(22d) and of their autonomous character in $q_2$ and $q_3$, it is enough to know how to solve the
linearized equations with an IC given by a single term in the sum (50). As we shall see, the solution is needed only at $\tau = 1$; it can be obtained by summing the time-Taylor series (25) to a suitable order, depending on the desired accuracy (sixth to eighth order is usually more than enough). In this way, one obtains the following expression for the first-order perturbation of the displacement:

$$ \epsilon \xi = \epsilon \sum_{m=1}^{M} \sum_{n=1}^{N} \left[ \mathbf{a}_{mn}(q_1; \tau) \cos(mq_2 + nq_3) + \mathbf{b}_{mn}(q_1; \tau) \sin(mq_2 + nq_3) \right]. \quad (51) $$

From this, one can calculate the Jacobian up to first order in $\epsilon$ to obtain

$$ J = 1 - \tau \cos q_1 + \epsilon[(1 - \tau \cos q_1)(\xi_{2,2} + \xi_{3,3}) + \xi_{1,1}]. \quad (52) $$

As we know, for $\epsilon = 0$ the Jacobian vanishes for the first time at $\tau_*^{(0)} = 1$ and $q_{1,*}^{(0)} = 0$ (the coordinate system was actually chosen to ensure this, without loss of generality for 1D flow.). For small $\epsilon$, by continuity, the perturbed Jacobian will vanish at a time and place close to $\tau = 1$ and $q_1 = 0$. For such values, we have

$$ 1 - \tau \cos q_1 \approx (1 - \tau) + \frac{q_1^2}{2} + \text{h.o.t.,} \quad (53) $$

where h.o.t. stands for higher-order terms in $1 - \tau$ and $q_1$. To determine the leading order of the perturbed first shell-crossing we may thus discard in (52) the higher-order term involving the factors $\epsilon(1 - \tau \cos q_1)$ and use the following approximation:

$$ J \approx (1 - \tau) + \frac{q_1^2}{2} + \epsilon \xi_{1,1}(0, q_2, q_3; 1) + \text{h.o.t.,} \quad (54) $$

where $\xi_{1,1}(0, q_2, q_3; 1)$ denotes the $q_1$ derivative of the first component of the displacement $\xi$, evaluated at time $\tau = 1$, at $q_1 = 0$ and arbitrary $q_2$ and $q_3$ (so far). Using (51), the Jacobian takes the following form

$$ J \approx (1 - \tau) + \frac{q_1^2}{2} + \epsilon \sum_{m=1}^{M} \sum_{n=1}^{N} [\alpha_{mn} \cos(mq_2 + nq_3) + \beta_{mn} \sin(mq_2 + nq_3)] + \text{h.o.t.,} \quad (55) $$

where the coefficients $\alpha_{mn}$ and $\beta_{mn}$ are easily expressed in terms of the coefficients $a_{mn}$ and $b_{mn}$ that appear in (51).

Although the determination of $a_{mn}$ and $b_{mn}$ cannot be done by purely analytic means, the very form of (51) allows us to conclude that the first shell-crossing takes place at time $\tau_* = 1 + \epsilon C$, where

$$ C = \min_{q_2, q_3} \sum_{m=1}^{M} \sum_{n=1}^{N} [\alpha_{mn} \cos(mq_2 + nq_3) + \beta_{mn} \sin(mq_2 + nq_3)]. \quad (56) $$

The first shell-crossing is at $q_{1,*} = 0$ and at those values of $q_{2,*}$ and $q_{3,*}$ for which the minimum in (56) is achieved. That this minimum is negative, and thus that the perturbed shell-crossing takes place slightly before the unperturbed one, follows from the observation that a trigonometric polynomial (in one or several variables) without a constant term necessarily takes both positive and negative values, because it is continuous and its space average over the period(s) is zero.

## 7 A CONCRETE EXAMPLE: THE THREE-SINE WAVES MODEL

Let us apply the developed tools to a concrete example, for which we determine, to first order in the perturbation, the Jacobian of

the Lagrangian map and the time of first shell-crossing. We set the initial gravitational potential to

$$ \varphi^{(\text{init})} = - \cos q_1 + \epsilon_2 \sin q_2 + \epsilon_3 \sin q_3. \quad (57) $$

Here, $\epsilon_2 = \epsilon C_2$ and $\epsilon_3 = \epsilon C_3$ with $\epsilon > 0$, $C_2 > 0$, $C_3 > 0$, so that the relative amplitudes of the $q_2$-dependent perturbation and of the $q_3$-dependent perturbation can be taken arbitrarily.

From the analysis of the previous section, obtaining the time of first shell-crossing requires the knowledge of the perturbed Jacobian at $\tau = 1$ and $q_1 = 0$. The latter is given by (38) which involves a time-Taylor series to all orders in $\tau$, which is guaranteed to converge at $\tau = 1$ because we are dealing, as we have seen in Section 5, with an entire function of $\tau$. In practice, to obtain numerical approximations for the time of first shell-crossing, we can truncate this time-Taylor series to a finite order $N$, using instead of (38), the truncated Jacobian

$$ J_N = 1 - \tau \cos q_1 + \epsilon \sum_{n=1}^{N} (\chi_{,1}^{(n)} + [1 - \tau \cos q_1] \zeta^{(n)}) \tau^n. \quad (58) $$

Given that the initial perturbation (57) is composed of two sine waves in the transverse coordinates $q_2$ and $q_3$, the Jacobian and its truncations will also have this property and we can write

$$ J_N = 1 - \tau \cos q_1 + \sum_{n=1}^{N} \kappa^{(n)}(q_1; \tau) [\epsilon_2 \sin q_2 + \epsilon_3 \sin q_3] \tau^n, \quad (59) $$

where the coefficients $\kappa^{(n)}(q_1; \tau)$ are easily computed by using our recursion relations (36) and symbolic algebra tools. The first few coefficients to order $N$ read

$$ \kappa^{(1)} = 1 - \tau \cos q_1, \quad (60) $$

$$ \kappa^{(2)} = \frac{3}{14} (2 - \tau \cos q_1) \cos q_1, \quad (61) $$

$$ \kappa^{(3)} = \frac{\cos q_1}{420} [50 \cos q_1 + \tau \{\cos(2q_1) - 25\}], \quad (62) $$

$$ \kappa^{(4)} = \frac{\cos q_1}{184800} [7000 - 280 \cos(2q_1) - 2715 \tau \cos q_1 + 443 \tau \cos(3q_1)]. \quad (63) $$

Higher orders are somewhat bulky and are given in Appendix A (for $5 \le N \le 8$).

To obtain the perturbed time of first shell-crossing to leading order in $\epsilon$ we can set $\tau = 1$ and $q_1 = 0$ in the coefficients $\kappa^{(n)}$, because the discrepancies will only contribute to higher orders. In this way, we finally obtain the following Taylor-truncated approximations

$$ \tau_* \approx 1 - c_N (\epsilon_2 + \epsilon_3), \quad (64) $$

where

$$ c_4 = 0.3003, \quad c_6 = 0.3073, \quad c_7 = 0.3076, \quad c_8 = 0.3077. \quad (65) $$

As to the spatial location of this first shell-crossing, it is found to be at $q_{1,*} = 0$ and $q_{2,*} = q_{3,*} = -\pi/2$ (modulo $2\pi$).

## 8 CONCLUDING REMARKS

For any time-Taylor series, the radius of convergence is the distance between the expansion point and the closest singularity in the
complex-time plane. Applying this statement to the Euler-Poisson equations in an Eulerian formulation, where the density (and velocity) is expanded in a time-Taylor series, it is evident that the radius of convergence cannot be infinite, neither in 1D nor beyond, because of the explicit appearance of real-space density singularities at shell-crossing. In a Lagrangian formulation, by contrast, the use of the Lagrangian map acts as the desingularization transformation of the problem, and thus, shell-crossing can be investigated in a rigorous way. In particular, we have shown that in Q1D, the time-Taylor coefficients of the linearized displacement field are all non-vanishing, but we also found that its Taylor series has an infinite radius of convergence. Without linearization this is unlikely to remain true, and we expect that there will be complex-time singularities and a finite radius of convergence of the Taylor series, which will however be large when the perturbation parameter $\epsilon$ is small.

The time and location of the first shell-crossing – which in Lagrangian coordinates is not a singularity but just a vanishing of the Jacobian, can be found by perturbation theory. In practice, one needs to calculate a sufficient number of time-Taylor coefficients $\xi^{(n)}$ for the displacement field $\xi = \sum_{n=1}^{\infty} \xi^{(n)} \tau^n$, which are easily generated by the use of our novel recursion relations (36). As a general rule of thumb – verified for various ICs of the type (50) – the time of shell-crossing can be determined to four digit accuracy when truncating the time-Taylor expansion of the displacement up to the order of $n = 7$.

In our Lagrangian formulation, the classical ZA is achieved by setting $n = 1$ in the time-Taylor series, and discarding all higher-order time-Taylor coefficients. Traditionally, the ZA has not only been applied to the purely 1D but also to the 2D and 3D cases. Most famous in this context is the prediction of the so-called Zel'dovich pancake that originates from the gravitational collapse of an ellipsoidal distribution of pressureless matter (see e.g. Arnold, Shandarin & Zeldovich 1982). In 1D the ZA is exact until shell-crossing, so one could hope that in a quasi-1D situation, the ZA would still give meaningful results. This, however, need not be the case, as is now demonstrated with a simple counterexample. Setting $\varphi_{g}^{(\text{init})} = -\cos q_1 + \epsilon \sin q_2$, it is straightforward to determine the particle trajectory within the classical ZA; it reads $x_i^{ZA} = q_i - \tau \delta_{i1} \sin q_1 - \epsilon \tau \delta_{i2} \cos q_2$. Analysing the entries of the Jacobian matrix $x_{i,1}^{ZA}$, it is seen that, within the ZA and to any order in $\epsilon$, shell-crossing happens at $\tau_{*}^{ZA} = 1$ and at $q_{1*}^{ZA} = 0$ (modulo $2\pi$) and arbitrary values of $q_2$ and $q_3$. However, using our tools and determining the displacement up to the eighth order in the time-Taylor series, we find that shell-crossing happens already at $\tau_{*} = 1 - 0.3077\epsilon$, at the position $q_{1*} = 0$ and $q_{2*} = -\pi/2$ (and arbitrary $q_3$). Thus, shell-crossing happens earlier than predicted by the ZA, and furthermore it occurs at a specific value of $q_{2*}$, and not for arbitrary values of $q_2$.

We note that in our model we have assumed that quasi-one-dimensionality holds already at initial time $\tau = 0$. An improved model would have quasi-one-dimensionality holding after a pancake has formed at some $\tau^{(\text{init})} > 0$. To handle this one should use time-Taylor expansions around $\tau^{(\text{init})}$ (and not around $\tau = 0$). Since the perturbed Euler-Poisson equations (22) are non-autonomous in the time variable, it follows that the resulting time-Taylor coefficients around the 'shifted' expansion point will differ from the ones we have obtained. Developing recursion relations for the time-shifted expansion is in principle fairly straightforward, but will be left for future work.

What happens after shell-crossing? Although this is an important question, in the present paper we have focused on the time before and at shell-crossing. We thus leave this issue to follow-up

studies. Qualitatively, it is expected that for sufficient short times after shell-crossing, the fluid description should still deliver physically meaningful results, provided that the Poisson equation and Lagrangian mass conservation are appropriately generalized – to take into account the multiple branches of the Lagrangian map. Deep into the multi-stream regime, however, a phase-space description becomes eventually mandatory. Such a phase-space description, in Q1D and beyond, is still missing in the literature (for the 1D case, see Colombi & Touma 2014; Colombi 2015; Taruya & Colombi 2017, and for approximative models beyond 1D, see Buchert & Domínguez 1998). In principle, cosmological N-body simulations aim to solve the said phase-space dynamics to high accuracy, however, being a brute-force method by nature and furthermore relying on a particle description, it is quite a challenge to gain mathematical insight about shell-crossing and the time after. In this context, we note the novel cosmological simulations of (Hahn, Abel & Kaehler 2013; Hahn & Angulo 2016) where N-body particles are used as tracers of (adaptively refineable) phase-space elements. Here, a smooth representation of the gravitational field is obtained, which improves the force computation especially near caustics, and delivers the phase-space dynamics to a good approximation.

Finally, let us comment on the possibility of applying our methodology to the relativistic shell-crossing in Q1D, an outstanding problem within the field of General Relativity. For reasons similar to those given in this paper, we expect that a Lagrangian-coordinates approach would be more fruitful than an Eulerian one. In particular, the use of the Lagrangian map could possibly desingularize the relativistic problem as well. For an irrotational and pressureless matter fluid, it is known that a synchronous-comoving coordinate system resembles the relativistic Lagrangian frame of reference (see e.g. Rampf & Wiegand 2014; Rampf, Rigopoulos & Valkenburg 2014). A promising starting point for such an investigation could be the relativistic Lagrangian equations (24)–(27) of Alles et al. (2015), which are closely related to our starting point, the Newtonian Lagrangian equations (8).

## ACKNOWLEDGEMENTS

We thank Tom Abel, Thomas Buchert, Stéphane Colombi, Vincent Desjacques, Adi Nusser, Fabian Schmidt, Sergei Shandarin and Barbara Villone for useful discussions and/or comments on the manuscript. CR is supported by the DFG through the SFB-Transregio TRR33 'The Dark Universe'. CR thanks the Lagrange Laboratory of the Observatoire de la Côte d'Azur and UF thanks the Technion Department of Physics for their hospitality and support.

## REFERENCES

* Alles A., Buchert T., Al Roumi F., Wiegand A., 2015, Phys. Rev. D, 92, 023512

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

* 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

* Brenier Y., Frisch U., Henon M., Loeper G., Matarrese S., Mohayaee R., Sobolevskii A., 2003, MNRAS, 346, 501

* Buchert T., 1992, MNRAS, 254, 729

* Buchert T., 1994, MNRAS, 267, 811

* Buchert T., Bartelmann M., 1991, A&A, 251, 389

* Buchert T., Domínguez A., 1998, A&A, 335, 395

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

* Buchert T., Goetz G., 1987, J. Math. Phys., 28, 2714
Buchert T., Melott A. L., Weiss A. G., 1994, A&A, 288, 349

Buchert T., Karakatsanis G., Klafl R., Schiller P., 1997, A&A, 318, 1

Colombi S., 2015, MNRAS, 446, 2902

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

Ehlers J., Buchert T., 1997, Gen. Relativ. Grav., 29, 733

Goroff M. H., Grinstein B., Rey S. J., Wise M. B., 1986, ApJ, 311, 6

Hahn O., Abel T., Kaehler R., 2013, MNRAS, 434, 1171

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

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

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

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

Melott A. L., Shandarin S. F., 1993, ApJ, 410, 469

Melott A. L., Buchert T., Weiss A. G., 1995, A&A, 294, 345

Moser J., 1959, J. Math. Z., 72, 379

Novikov E. A., 1970, Sov. Phys. JETP, 30, 512

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

Podvigina O., Zheligovsky V., Frisch U., 2016, J. Comput. Phys., 306, 320

Rampf C., Buchert T., 2012, J. Cosmol. Astropart. Phys., 1206, 021

Rampf C., Wiegand A., 2014, Phys. Rev. D, 90, 123503

Rampf C., Rigopoulos G., Valkenburg W., 2014, Class. Quant. Grav., 31, 234004

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

Rampf C., Villa E., Bertacca D., Bruni M., 2016, Phys. Rev. D, 94, 083515

Shandarin S. F., 1992, NASA STI/Recon Technical Report N, 95

Shandarin S. F., 1994, Phys. D Nonlinear Phenomena, 77, 342

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

Tassev S., Zaldarriaga M., 2012, J. Cosmol. Astropart. Phys., 1212, 011

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

Zentsova A. S., Chernin A. D., 1980, Astrophysics, 16, 108

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

## APPENDIX A: HIGHER-ORDER TAYLOR COEFFICIENTS FOR THE THREE-SINE WAVES MODEL

In Section 7 we have applied our tools to an explicit example, the so-called three-sine waves model, for which we have generated so-

lutions for the Jacobian up to the order of $N = 8$ (see equation 59). For brevity, we have skipped in the main text the terms beyond the order of $N > 4$; for the higher-order terms we find

$$ \kappa^{(5)} = \frac{\cos q_1}{680680000} [17\{76316\tau \cos(2q_1) + 407250 \cos q_1 - 75(886 \cos(3q_1) + 2715\tau)\} + 59485\tau \cos(4q_1)], \quad (A1) $$

$$ \kappa^{(6)} = -\frac{\cos q_1}{2654652000000} [\tau\{2659396870 \cos(q_1) - 962623441 \cos(3q_1) + 28312075 \cos(5q_1)\} + 2860(1297372 \cos(2q_1) + 59485 \cos(4q_1) - 3461625)], \quad (A2) $$

$$ \kappa^{(7)} = \frac{10^{-7} \cos q_1}{198703356852} [37\{3\tau(3999222278389 \cos(2q_1) - 55818626690 \cos(4q_1)) + 15470(-962623441 \cos(3q_1) + 28312075 \cos(5q_1) - 1329698435\tau) + 41140869578900 \cos q_1\} - 1201960568875\tau \cos(6q_1)], \quad (A3) $$

$$ \kappa^{(8)} = \frac{10^{-8} \cos q_1}{130879277713184} [\tau\{407511682759452297 \cos(3q_1) - 827808374798014805 \cos q_1 + 335214948722225 \cos(7q_1) - 32682412580801125 \cos(5q_1)\} + 5200(761106087209650 - 443913672901179 \cos(2q_1) + 6195867562590 \cos(4q_1) + 1201960568875 \cos(6q_1))]. \quad (A4) $$

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