Shell-crossing in quasi-one-dimensional flow

TL;DR

This paper develops a perturbative approach in Lagrangian coordinates to analyze shell-crossing in quasi-one-dimensional cosmological flows, providing a way to predict the timing and location of singularities more accurately than in purely 1D models.

Contribution

It introduces an all-order recursion relation for displacement field coefficients, demonstrating the convergence of the series and enabling precise prediction of shell-crossing in Q1D flows.

Findings

All-order recursion relations derived for displacement field coefficients.

The Taylor series for shell-crossing has an infinite radius of convergence.

Shell-crossing occurs earlier in Q1D flows than in pure 1D models.

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.

Equations

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

$${\partial }_{\tau }\mathbf{v}+(\mathbf{v}\cdot \nabla )\mathbf{v}=-\frac{3}{2\tau }\left(\mathbf{v}+\nabla {\phi }_{\mathrm{g}}\right),$$

$${\partial }_{\tau }\delta +\nabla \cdot \left[(1+\delta )\mathbf{v}\right]=0,$$

$${\nabla }^2{\phi }_{\mathrm{g}}=\frac{\delta }{\tau },$$

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

$${\delta }^{(\mathrm{init})}=0,{\mathbf{v}}^{(\mathrm{init})}=-\nabla {\phi }_{\mathrm{g}}^{(\mathrm{init})}.$$

$$\nabla \times \mathbf{v}=0.$$

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

$$\delta =1/J-1,$$

$${\epsilon }_{ikl}{\epsilon }_{jmn}{x}_{k,m}{x}_{l,n}{\mathfrak{R}}_{\tau }{x}_{i,j}=3\left(J-1\right),$$

$${\epsilon }_{ijk}{\dot{x}}_{l,j}{x}_{l,k}=0,$$

$${\phi }_{\mathrm{g}}^{(\mathrm{init})}(q_1,q_2,q_3)=-\cos q_1+ϵ{\varphi }^{(\mathrm{init})}(q_1,q_2,q_3),$$

$$⟨{\varphi }^{(\mathrm{init})}⟩\equiv {\int }_0^{2\text{uppi}}\frac{\mathrm{d}{q}_2}{2\text{uppi}}{\int }_0^{2\text{uppi}}\frac{\mathrm{d}{q}_3}{2\text{uppi}}{\varphi }^{(\mathrm{init})}(q_1,q_2,q_3)=0.$$

$$v_i^{(\mathrm{init})}(\mathbf{q})=-{\delta }_{i1}\sin q_1-ϵ{\partial }_i^{\mathrm{L}}{\varphi }^{(\mathrm{init})}(\mathbf{q}),$$

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

$${\xi }_i^{(0)}(\mathbf{q})={\delta }_{i1}F(q_1;\tau ).$$

$$x_i(\mathbf{q};\tau )=q_i+{\delta }_{i1}F(q_1;\tau )+ϵ{\xi }_i(\mathbf{q};\tau ).$$

$$x_{i,j}={\delta }_{ij}+{\delta }_{i1}{\delta }_{j1}{F}_{,1}+ϵ{\xi }_{i,j},$$

$$J=1+F_{,1}+ϵ\left({\xi }_{1,1}+{\xi }_{2,2}+{\xi }_{3,3}\right)+ϵF_{,1}\left({\xi }_{2,2}+{\xi }_{3,3}\right).$$

$$2{\mathfrak{R}}_{\tau }{F}_{,1}=3F_{,1}.$$

$$F_{,1}(q_1;\tau )=\tau {\chi }_{,1}(q_1),$$

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

$$x_i=q_i-{\delta }_{i1}\tau \sin q_1,$$

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

$${\mathfrak{R}}_{\tau }{\xi }_{1,1}+\left[1-\tau \cos q_1\right]{\mathfrak{R}}_{\tau }\left({\xi }_{2,2}+{\xi }_{3,3}\right)=\frac{3}{2}{\xi }_{l,l},$$

$${\dot{\xi }}_{2,3}-{\dot{\xi }}_{3,2}=0,$$

$$\left[1-\tau \cos q_1\right]{\dot{\xi }}_{1,3}+{\xi }_{1,3}\cos q_1={\dot{\xi }}_{3,1},$$

$$\left[1-\tau \cos q_1\right]{\dot{\xi }}_{1,2}+{\xi }_{1,2}\cos q_1={\dot{\xi }}_{2,1}.$$

$$\chi \equiv {\xi }_1,\zeta \equiv \mathrm{i}{k}_2{\xi }_2+\mathrm{i}{k}_3{\xi }_3,k_{\perp }^2\equiv k_2^2+k_3^2,$$

$${\mathfrak{R}}_{\tau }{\chi }_{,1}+\left[1-\tau \cos q_1\right]{\mathfrak{R}}_{\tau }\zeta =\frac{3}{2}\left({\chi }_{,1}+\zeta \right),$$

$$k_{\perp }^2\left\{\left[1-\tau \cos q_1\right]\dot{\chi }+\chi \cos q_1\right\}=-{\dot{\zeta }}_{,1}.$$

$$\chi (\mathbf{q};\tau )=\sum _{n=1}^{\infty }{\chi }^{(n)}(\mathbf{q}){\tau }^n,\zeta (\mathbf{q};\tau )=\sum _{n=1}^{\infty }{\zeta }^{(n)}(\mathbf{q}){\tau }^n.$$

$$\left[n^2+\frac{n}{2}-\frac{3}{2}\right]\left({\chi }_{,1}^{(n)}+{\zeta }^{(n)}\right)=\left[n^2-\frac{3n}{2}+\frac{1}{2}\right]{\zeta }^{(n-1)}\cos q_1,$$

$$n{\chi }^{(n)}+k_{\perp }^{-2}n{\zeta }_{,1}^{(n)}=(n-2){\chi }^{(n-1)}\cos q_1.$$

Full text

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.

††pubyear: 2017††pagerange: Shell-crossing in quasi-one-dimensional flow–A

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

[TABLE]

where ${q}$ and ${v}^{\rm(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 one-dimensional 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, low-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 (Buchert & Ehlers, 1993; Bouchet et al., 1992), to the third (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 et al. (1994); Melott et al. (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\leq\tau\leq T$. Furthermore, Rampf et al. (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 et al., 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 one-dimensional. By contrast, all-order recursive solutions for generic 3D initial data (see also Goroff et al., 1986; Ehlers & Buchert, 1997; Rampf et al., 2015; Matsubara, 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 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 depends 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 initial condition. 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 $\boldsymbol{x}=\boldsymbol{r}/a$, where $\boldsymbol{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 $\boldsymbol{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 $\varphi_{\rm 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)

[TABLE]

(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 $\boldsymbol{v}=\boldsymbol{0}$ and $\delta=0$, the above equations can be written in terms of a single differential equation for the density contrast (Peebles, 1980),

[TABLE]

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)]:

[TABLE]

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,

[TABLE]

(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^{\rm(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 $\boldsymbol{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}^{\rm 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 $\boldsymbol{q}\mapsto{\boldsymbol{x}}(\boldsymbol{q};\tau)$ be the Lagrangian map from the initial ($\tau\!=\!0$) position $\boldsymbol{q}$ to the Eulerian position ${\boldsymbol{x}}$ at time $\tau$. The map satisfies $\boldsymbol{v}(\boldsymbol{x}(\boldsymbol{q};\tau);\tau)=\partial_{\tau}^{\rm L}\boldsymbol{x}(\boldsymbol{q};\tau)$, where $\partial_{\tau}^{\rm L}$ is the Lagrangian time derivative – the latter also denoted with an overdot in the following. At initial time, $\tau=0$, the velocity is

[TABLE]

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

[TABLE]

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

[TABLE]

where we have defined the operator $\mathfrak{R}_{\tau}\equiv\tau^{2}\left(\partial_{\tau}^{\rm L}\right)^{2}+(3\tau/2)\partial_{\tau}^{\rm 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 initial conditions (ICs) that are close to one dimension, 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 one-dimensional 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 initial conditions for the purely one-dimensional problem, could be taken quite arbitrary (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 one-dimensional 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 one-dimensional gravitational potential the function $-\cos q_{1}$. Considerations of normal-form reduction indicate that we can actually make this choice of initial conditions without loss of generality (within the class of $2\uppi$-periodic functions).

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

[TABLE]

where $\epsilon>1$ is a small perturbation parameter and $\phi^{\rm(init)}$ an arbitrary $2\uppi$-periodic function of $(q_{1},q_{2},q_{3})$. Without loss of generality, we can take the transverse average to be zero, namely

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 one-dimensional, 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_{1\star}^{(0)}=0$ (modulo $2\uppi$) 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\star},q_{2\star},q_{3\star})$ (modulo $2\uppi$). 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_{\star}^{(0)}$ and $\tau_{\star}$, are discussed in the following sections.

3.3 Lagrangian equations to zeroth order in $\boldsymbol{\epsilon}$

To lowest (zeroth) order in $\epsilon$, the problem is exactly one-dimensional and has the initial condition $\varphi_{\rm g}^{(\rm init)}(q_{1},q_{2},q_{3})=-\cos q_{1}$. In this instance, as is well known, the Zel’dovich approximation 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

[TABLE]

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

[TABLE]

where $\chi(q_{1})$ can be taken arbitrary, but by (11) is actually equal to the initial velocity

[TABLE]

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

[TABLE]

and the unperturbed Jacobian is

[TABLE]

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

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 one-dimensional 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 $\boldsymbol{\epsilon}$

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

[TABLE]

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})=\tilde{\varphi}(q_{1})\exp\{{\rm{i}}(k_{2}q_{2}+k_{3}q_{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 ${\rm{i}}k_{2}$ and ${\rm{i}}k_{3}$, respectively. Using this and defining

[TABLE]

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

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

and then, making use of Euler’s formula $\cos q_{1}=(\exp\{{\rm{i}}q_{1}\}+\exp\{-{\rm{i}}q_{1}\})/2$, we obtain for equations (26)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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\geq 1$)

[TABLE]

and in real space

[TABLE]

where $\nabla_{\rm 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,

[TABLE]

where $\chi_{,1}^{(n)}=\xi_{1,1}^{(n)}$ and $\zeta^{(n)}=\xi_{2,2}^{(n)}+\xi_{2,2}^{(n)}$. Equations (37)–(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_{\star}$, the time of first shell-crossing in the perturbed problem.

5 No Singularities in Lagrangian solutions

Observe that in (24a), there is a term $[1-\tau\cos q_{1}]\mathfrak{R}_{\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. 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)}(\boldsymbol{q})$ and $\zeta^{(n)}(\boldsymbol{q})$. These satisfy the recursion relations (31), which can be simplified for large $n$ and approximated by their asymptotic form:

[TABLE]

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

[TABLE]

where $N$ is taken large enough to be able to use the asymptotic form (39) of the recursion relations. Comparison of (25) and (40) shows that $\chi_{{}_{N}}(\boldsymbol{q};\tau)$ and $\chi(\boldsymbol{q};\tau)$ differ by a polynomial in $\tau$ (with $\boldsymbol{q}$-dependent coefficients) of degree $N-1$. The same statement holds for $\chi_{{}_{N}}(\boldsymbol{q};\tau)$ and $\chi(\boldsymbol{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 $k_{\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:

[TABLE]

The functions $\chi_{{}_{N}}$ and $\zeta_{{}_{N}}$ are coupled by a system of two equations. However, by simply introducing the sum and the difference

[TABLE]

we obtain two decoupled equations:

[TABLE]

where

[TABLE]

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\uppi$-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_{\rm in}$ of $q_{1}$ and look for the solution to its right, i.e. for $q\geq q_{\rm in}$. This has the explicit solution

[TABLE]

in terms of the Green’s function $G^{+}$ of the associated linear differential equation (without the $F_{N}^{+}$ term), given, for $q_{1}\geq q^{\prime}$, by

[TABLE]

The solution (45) is in general not periodic in $q_{1}$, but a periodic solution can be constructed by letting $q_{\rm in}\to-\infty$, because far to the right of the point $q_{\rm in}$, the solution relaxes to spatial periodicity. This is proved by decomposing the interval $[q_{\rm in},q_{1}]$ into adjacent intervals all-but-the-first (which may be smaller) of length $2\uppi$. The number $M$ of intervals of length $2\uppi$ is the integer part of $|q_{1}-q_{\rm in}|/(2\uppi)$ and thus tends to infinity when $q_{\rm in}\to-\infty$. We observe that, when the argument $q^{\prime}$ of the integrand on the right-hand side of (45) is shifted from one interval to the neighbouring left interval by subtracting $2\uppi$, the integrand is multiplied by a factor $\exp(-2\uppi)$. Indeed, the only $q^{\prime}$-dependent term which is not periodic in $q^{\prime}$ is $(q_{1}-q^{\prime})$ in the exponential, which generates the stated factor. Hence, the sum over all the $M$ intervals produces a geometric series of ratio $\exp(-2\uppi)$. As $q_{\rm in}\to-\infty$ at fixed $Z_{N}^{+}(q_{\rm 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

[TABLE]

It is easily checked that (47) is a $2\uppi$-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

[TABLE]

which implies that

[TABLE]

Since $Z_{N}^{+}(q_{1};\tau)$ and $\tilde{Z}_{N}^{+}(q_{1};\tau)$ are $2\uppi$-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^{\prime};\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\leq q_{\rm 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 one-dimensional and with our initial condition, shell-crossing happens at the time $\tau_{\star}^{(0)}=1$ and location $q_{1\star}=0$, and for arbitrary $q_{2}$ and $q_{3}$. Thus, the whole plane $q_{1}=0$ shell-crosses at $\tau_{\star}^{(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_{\star}\neq 1$ and at a single location $(q_{1\star},\,q_{2\star},\,q_{3\star})$. We shall now show that $\tau_{\star}$ is generically happening earlier than $\tau_{\star}^{(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

[TABLE]

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 initial condition given by a single term in the sum (6). 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 (6–8th order is usually more than enough). In this way, one obtains the following expression for the first-order perturbation of the displacement:

[TABLE]

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

[TABLE]

As we know, for $\epsilon=0$ the Jacobian vanishes for the first time at $\tau_{\star}^{(0)}=1$ and $q_{1\star}^{(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

[TABLE]

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:

[TABLE]

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

[TABLE]

where the coefficients $\alpha_{mn}$ and $\beta_{mn}$ are easily expressed in terms of the coefficients $\mathfrak{a}_{mn}$ and $\mathfrak{b}_{mn}$ that appear in (51).

Although the determination of $\mathfrak{a}_{mn}$ and $\mathfrak{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_{\star}=1+\epsilon C$, where

[TABLE]

The first shell-crossing is at $q_{1\star}=0$ and at those values of $q_{2\star}$ and $q_{3\star}$ 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

[TABLE]

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

From the analysis of the previous section, obtaining the time of first shell-crosssing 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 ${\cal N}$, using instead of (38), the truncated Jacobian

[TABLE]

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

[TABLE]

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 ${\cal N}$ read

[TABLE]

Higher orders are somewhat bulky and are given in Appendix A (for $5\leq{\cal N}\leq 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

[TABLE]

where

[TABLE]

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

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

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 $\boldsymbol{\xi}^{(n)}$ for the displacement field $\boldsymbol{\xi}=\sum_{n=1}^{\infty}\boldsymbol{\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 initial conditions of the type (6) – the time of shell-crossing can be determined to four digit accuracy when truncating the time-Taylor expansion of the displacement up to order $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. Tradionally, the ZA has not only been applied to the purely 1D but also to the 2D and 3D case. 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 et al. 1982). In one dimension the ZA is exact until shell-crossing, so one could hope that in a quasi-one-dimensional 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_{\rm g}^{(\rm init)}=-\cos q_{1}+\epsilon\sin q_{2}$, it is straightforward to determine the particle trajectory within the classical ZA, it reads $x_{i}^{\rm 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,j}^{\rm ZA}$, it is seen that, within the ZA and to any order in $\epsilon$, shell-crossing happens at $\tau_{\star}^{\rm ZA}=1$ and at $q_{1\star}^{\rm ZA}=0$ (modulo $2\uppi$) and arbitrary values of $q_{2}$ and $q_{3}$. However, using our tools and determining the displacement up to the 8th order in the time-Taylor series, we find that shell-crossing happens already at $\tau_{\star}=1-0.3077\epsilon$, at the position ${q}_{1\star}=0$ and $q_{2\star}=-\uppi/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\star}$, 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^{\rm(init)}>0$. To handle this one should use time-Taylor expansions around $\tau^{\rm(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 & Dominguez 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 et al. (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 et al., 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.

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 solutions for the Jacobian up to order ${\cal N}=8$ (see equation (59)). For brevity, we have skipped in the main text the terms beyond order ${\cal N}>4$; for the higher-order terms we find

[TABLE]

[TABLE]