Cosmological Vlasov-Poisson equations for dark matter: Recent developments and connections to selected plasma problems

The cosmic large-scale structures of the Universe are mainly the result of the gravitational instability of initially small density fluctuations in the dark-matter distribution. Dark matter appears to be initially cold and behaves as a continuous and collisionless medium on cosmological scales, with evolution governed by the gravitational Vlasov--Poisson equations. Cold dark matter can accumulate very efficiently at focused locations, leading to a highly non-linear filamentary network with extreme matter densities. Traditionally, investigating the non-linear Vlasov--Poisson equations was typically reserved for massively parallelised numerical simulations. Recently, theoretical progress has allowed us to analyse the mathematical structure of the first infinite densities in the dark-matter distribution by elementary means. We review related advances, as well as provide intriguing connections to classical plasma problems, such as the beam-plasma instability.

. Here, is the total density while¯( ) is an isotropic dilution factor stemming from the volume expansion of the Universe.

Basic problem
What is the nature of dark matter and dark energy, which together make up about 95% of the Universe's energy content? One way to address such a question is to investigate the evolution of cosmic structures on the largest observable length scales (Fig. 1).
Dark matter constitutes the bulk part of the overall matter distribution. Apart from gravitational interactions, dark matter appears to be extremely weakly interacting (e.g. [2,3,4,5]), thereby justifying the validity of the collisionless limit on cosmological scales (∼ 10 22 m − 10 26 m). The gravitational evolution of such a collisionless medium is governed by the cosmic Vlasov-Poisson equations, which describe how the darkmatter distribution = ( , , ) evolves in the six-dimensional phase-space, where is the matter density contrast defined through the matter density For convenience we have added individual colour markings to the dark-matter fluid particles (i.e., the Lagrangian coordinate), thereby highlighting the fluid motion during the temporal evolution. Shown results were obtained using a particle in cell simulation with 8192 fluid particles , initialised at = 0 with a sine-wave velocity profile and vanishing density contrast [14] ( ∝ 2/3 where is standard time). The second panel shows the critical instant when the dark-matter sheet folds for the first time (here at = 1 and = 0), which is usually denoted with shell-crossing or stream crossing. The third and fourth panels show multi-streaming regions where multiple matter streams occupy the same current position. also assumed in the present paper. A perfectly cold distribution has vanishing (thermal) velocity dispersion, implying that at sufficiently early times, the whole kinetic information is encoded into a single velocity field := /( 2 ). Nevertheless, non-zero velocity dispersion (and thus an effective non-zero temperature) is generated during the later stages of the gravitational dynamics, which is discussed next.
The topology of the phase-space distribution is such that it is confined on an infinitesimally thin sheet along the momentum direction and thus occupies, at all times, only a three-dimensional hypersurface in the full 6D phase-space (the Lagrangian submanifold). The resulting dark-matter sheet is illustrated in Fig. 2 in a simplified set-up in 1+1 dimensions. During the course of the gravitational evolution, the dark-matter sheet never tears and remains connected, due to the Hamiltonian nature of the system. Nevertheless, gravitational interactions can generate folds in the dark-matter sheet. At those folds, the number of matter streams that occupy the same current position changes. An example of this is shown in Fig. 2, which indicates the emergence of a proliferation of multiple fluid streams.
Obviously, solving the Vlasov-Poisson equations (1) in "full 6D" is vastly excessive, since the dark-matter distribution is confined on a three-dimensional hypersurface. Therefore, it is customary, in both theory and numerics, to employ an efficient parametrisation of the dark-matter sheet, that essentially tracks how the positions and momenta of dark-matter particles change on the three-dimensional hypersurface. Precisely such a "map" is realised when introducing Lagrangian coordinates, which in cosmology are usually denoted with and can be thought of as a continuous set of particle labels (denoted by the colour marking of the dark matter sheet in Fig. 2).
Here and in the following, over-dots denote convective (i.e., total) derivatives with respect to the dimensionless temporal variable = ; using the cosmic scale factor (or more generally, the growth function of linear matter density fluctuations) as the temporal variable is another ingredient for establishing time-analytic solutions for sufficiently short times; see section 3 for details. The Lagrangian map is defined such that =: ( , ) is the dark-matter velocity field (in units of lengths since = is dimensionless), where relates to the momentum (non-canonically) according to = 2 . We note that equations (6) are derived for the case when the Universe is spatially flat and solely composed of dark matter (i.e., the simplified Einstein-de Sitter model), but this could be easily generalised to the more realistic ΛCDM model if needed (see e.g. [17,18,19]). Nonetheless, from here on we will mostly work with this simplified case of an Einstein-de Sitter model to avoid unnecessary cluttering of the expressions.
The matter density contrast appearing in the equations (6) can be exactly expressed through a mass conservation law using the Dirac-delta (3) D [15,20,21,22], here for the case when → 0 initially, which, roughly speaking, may be interpreted as the continuum limit → ∞ of a density contrast induced through a (hypothetical) collection of discrete particles with ( , ) + 1 ∼ =1 (3) D ( − ( )), where ( ) denotes the current position of the th dark-matter particle (see e.g. [28] for a more precise argument). Equation (7) is the result of linking the law of mass conservation¯(1 + ( , )) d 3 =¯d 3 to the determinant of the Lagrangian mapping. Specifically, employing the composition property of the Dirac delta, the mass conservation (7) can be equivalently written as where ∇ = ∇ ⊗ T is the direct (dyadic) product, while denotes the th root of the equation − ( , ) = 0. Intuitively, the sum in (8) appears since fluid streams occupy the same current position , and all fluid streams obviously contribute to the density (see also the discussion further below and Fig. 3). Note that as long as the fluid is single stream, there is only a single root at 1 = , implying that mass conservation simplifies to ( ( , )) +1 = 1/ in Lagrangian coordinates, where = det[∇ ( , )] is the Jacobian.
The Vlasov-Poisson equations (6) have the simple physical interpretation of being essentially a Newtonian equation of motion, with the addition of a "drag" term ∝ (1/ ) , which stems from the present choice of co-moving coordinates (co-moving with respect to a fixed position on a spatial grid that follows the overall Hubble expansion of the Universe). Note the structural similarity of equations (6)-(7) as compared to the Vlasov-Poisson equations for the beam-plasma instability, a classical plasma problem. There, a beam of (discretised) charged particles moves in a background neutral plasma. We will come back to such structural similarities in section 6.
Note that the Vlasov-Poisson equations (6) are invariant under the transformation for arbitrary time-dependent boosts ( ), provided that the gravitational potential transforms as ( ) = ( ) − ( + [2 /3] ) · + ( ), for arbitrary ( ). This non-Galilean invariance has been first analysed by Heckmann & Schücking in 1955 [23], and recently drew renewed attention in cosmology (e.g., [14,24,25,26,27]). This invariance is of particular relevance when determining critical phenomena in the dark-matter phase-space (see section 5.3, in particular panel d in Fig. 7). Finally, observe that the Vlasov-Poisson equations (6) maintain their regularity for → 0, provided one imposes the following boundary conditions on the initial conditions [29,30] at time = 0 (denoted by the superscript "ini"): where ini = ( , = 0). Thus, with these boundary conditions, only the initial gravitational potential ini at = 0 needs to be prescribed; see e.g. [31] for explicit instructions how this can be achieved in practice. Mathematically, initialising the Vlasov-Poisson system at time "zero" (a.k.a. the big bang) using (10) implies that the inflationary physics has been reduced to an infinitely thin boundary layer [17]. Also note that the second condition in (10) immediately implies that ∇ × ini = 0 which, due to Helmholtz's theorem for the conservation of vorticity flux, translates to ∇ × = 0 at all times along each matter stream (but note that effective vorticity is generated in Eulerian coordinates through multi-streaming [16,32,33]). In the following we will mostly focus on solutions using the conditions (10), since those are the only one currently available that have mathematical proofs of analyticity until some finite time; see sections 4.2 and 4.3 for details.
2.2 Theoretical challenge of resolving folds of the dark-matter sheet: Shell-crossings The instant when the dark-matter sheet begins to fold (central panel in Fig. 2), called the first shell-crossing, implies a critical moment in the phase-space evolution: at folds, streams of dark-matter particles accumulate at focused locations leading to formally infinite densities. This can also be seen in Fig. 3 where we show, in the left panel, the current positions of matter streams (the particle trajectories for certain times) and, in the right panel, the corresponding matter density. Whenever the number of matter streams changes, the determinant of ∇ vanishes, which, by virtue of equation (8) implies the said infinite density. The critical instant, in time and space, when the dark-matter continuum changes from single stream to triple stream, is denoted for traditional reasons by "shellcrossing" (since the first collapse models involved the study of "matter shells" in Fig. 3 Shown is the Lagrangian map ↦ → ( , ) from initial position to current position (left panel), and the density contrast = ( −¯)/¯(right panel), employing the analytical results of [14] in the same simplified set-up as for Fig. 2. Note specifically that a torus geometry is used with , ∈ [− , ), while shown is only the physically interesting regime. Initially ( = 0, green lines) the particle positions are perfectly uniformly distributed in space, which is reflected by a linear regression in the left panel, and by = 0 in the right panel resembling an initial homogeneous state. A critical instant in the dark-matter evolution arises when ∇ vanishes (corresponding to det[∇ ] = 0 in three space dimensions; cf. Eq. 8), which appears in the present case at = 1 for the first time (orange line in left panel). This instant is usually called shell-crossing and is accompanied with an infinite density (center region in right plot). Furthermore, it marks the starting point when multiple streams of particles can occupy the same current position (blue shaded region shown at = 1.3), which comes with non-zero velocity dispersion.
spherical symmetry). The first approximative theoretical model of the cosmic Vlasov-Poisson equation that predicts the appearance of the first shell-crossing traces back to the work of Ya. B. Zel'dovich in 1970 [34]. Nowadays the approach is called Zel'dovich solution (in 1D) or Zel'dovich approximation (beyond 1D), and, despite its age, comprises one of the most successful and insightful theoretical models in cosmic structure formation (e.g. [35,36,37,38]). We review this model in section 4.1, but note that the Zel'dovich approximation constitutes an exact solution of the Vlasov-Poisson equations in a fictitious one-dimensional setup until the first shell-crossing [39,40]. Here and in the following, a "one-dimensional setup" can be established by providing 3D initial conditions that only depend on one space coordinate (i.e., a onedimensional embedding in 3D space). Such solutions play an important role in the realistic modelling of cosmic structure formation, since, in the fundamental coordinate system of a given collapsing patch, shell-crossings generically arise as an almost onedimensional phenomena with the formation of pancake-like structures (e.g. [41]).
Secondary gravitational infall. Once matter streams have crossed for the first time, nearby particles experience non-trivial accelerations due to the secondary gravitational infall into the potential well. Subsequently, a self-trapped quasi-stationary structure emerges (e.g., a dark-matter halo in the long-term), forcing distinct matter streams to cross again -the second shell-crossing. For topological reasons, the number of streams generically bifurcates from three to five at the second shell-crossing (and then from five to seven at the third shell-crossing; etc.). While numerical simulations revealed the second shell-crossing several decades ago (e.g. [42,43]), it was not predicted by theoretical means until recently, when S. Colombi was able to accurately estimate the involved gravitational forces in multi- Fig. 4 Evolution of selected dark-matter trajectories, using the same 1D setup as in Fig. 3. Left panel: The analytical Zel'dovich solution [34] predicts the first shell-crossing (here at = 1 and = 0), however not the second one. Instead, non-linear structures are incorrectly washed out after the first shell-crossing. Central panel: By contrast, modern post-shell-crossing theories (here the purely analytical result of [14]) predict the second shell-crossing (here at 2.2 and = 0), however would require higher refinements to predict the third shell-crossing (not shown). streaming regions from first principles [44] (see also [14,45,46]). This is shown in Fig. 4, where we compare the classical Zel'dovich solution (left panel) against the refined theoretical prediction of particle trajectories (central panel) and a numerical simulation (right panel). First notice that in regions with no multi-streaming, the trajectories in all panels are identical and simply straight lines, indicating that particles follow their initially prescribed constant velocity with zero acceleration (note this feature of exact ballistic motion in single-stream regions is only true in the onedimensional case).
Starting at shell-crossing (here at = 1), some particle trajectories in multistreaming regions are bent back into the center of the potential well, a feature that is not encapsulated in Zel'dovich's theory, but is captured in refined theoretical predictions as well as in the numerical simulation (central and right panel of Fig. 4).
Predicting the third and even higher shell-crossings by analytical means is in principle feasible with straightforward extensions of the analytical techniques of [14,44,45], but have not been performed so far; this is also the reason why the shown refined theoretical prediction in Fig. 4 does not agree with the numerical prediction at later times.
Of course, the actual Universe with random initial conditions is more complicated than the depicted scenario in Fig. 4. In particular, due to the random nature, collapsing structures are generally not isolated objects. Instead, structures are formed on smaller spatial scales that eventually merge to combined structures on larger scales (e.g., dark-matter halos). Resolving such merging events is currently beyond the reach of analytical considerations, and only amenable by employing numerical simulations (but see [45,46] for related semi-analytical avenues using adaptive smoothing). On top of that, the process of halo or galaxy formation requires accurate physical modelling that goes well beyond the cosmic Vlasov-Poisson equations, in particular physics from visible (baryonic) matter; see e.g. [47,48,49,50,51].
In the following sections, we will focus on the mathematical analysis of the darkmatter evolution before and (shortly) after the first shell-crossing -a regime that until very recently was not accessible by analytical means. Nonetheless, in sections 5.1-5.2 we will summarise complementary theoretical (effective) approaches as well as numerical simulations that consider the dark-matter evolution at substantially later times.

Lagrangian evolution equations
The Vlasov-Poisson equations (6) are not yet in a form that allows for straightforward analytical considerations. One reason for this is the appearance of the Eulerian space derivative with respect to the coordinate , which in Lagrangian coordinates is not an independent variable. Indeed, earlier we have introduced the Lagrangian map ↦ → ( , ) which implies for the conversion of derivatives that ∇ = ( / )∇ , where, from here on, Latin indices , , . . . = 1, 2, 3 denote the three spatial components, and (Einstein) summation over repeated indices is assumed.
To proceed, one introduces the Lagrangian displacement field ( , ), which encodes the whole dynamical information of the Vlasov-Poisson equations for cold dark matter -at all times. Indeed, once is determined, one can retrieve the matter density (using Eq. 7), the velocity (possibly involving multi-stream averaging), and higher kinetic moments of the dark-matter distribution function, such as the velocity dispersion tensor (see e.g. [52]). Introducing the operatorT := (2 2 /3) 2 + − and applying suitable (derivative) operations to the Vlasov-Poisson equations (6), summarised below, the Lagrangian evolution equations for the divergence and curl of the displacement arê which, once solved respectively for ∇ · and ∇ × , yield the formal solution for the displacement using a Helmholtz decomposition (∇ −2 is the inverse Laplacian), We have defined two non-linear source terms (see below for physical interpretations) where " , " denotes a partial derivative with respect to Lagrangian component , while is the fundamental anti-symmetric tensor and the Kronecker delta. Here and in the following, displacements and current positions depend on and if not otherwise stated. Equations (12)- (15) are derived for an Einstein-de Sitter cosmological model, but are equally valid for a ΛCDM Universe upon identifying with the standard ΛCDM linear growth function and the replacementT →T Λ = Derivation of equations (12). To derive the divergence equation, one first takes the Eulerian space derivative of the Vlasov equation (6), thereby allowing one to express its right-hand side by the Poisson equation and the mass conservation (7). Next, one converts the remaining Eulerian space derivative according to ∇ = ( , , one then obtains the desired result (see e.g. [18,53] for equivalent derivations for the case M = 0). To derive the equation on the right in (12), one first multiplies (6) by 2 2 /3, leading toT 0 = − ∇ . Multiplying the last equation by the matrix element / then yields ,T0 = − ∇ . Finally, taking the Lagrangian curl of this leads to (∇ ) ×T 0 ∇ = 0 [18,30], which, once written for the displacement = − , leads directly to the second equation in (12).
Physical interpretation of source terms. M is only nonzero once there are locations with an overlap of multiple matter streams or, stated the other way around, M is exactly zero everywhere, as long as the flow is still single stream (mono-kinetic). Indeed, for single-stream flows, the Lagrangian mapping ↦ → is injective with det[∇ ] > 0 (corresponding to ∇ > 0 in 1D; cf. By contrast, W is the result of converting Eulerian space derivatives to Lagrangian space, thereby accumulating quadratic and cubic terms in the displacement. With regard to the physical interpretation: for generic cosmological initial conditions in three space dimensions, W constitutes a sub-dominant (perturbatively suppressed) term during the very early gravitational evolution. Thus, the first evolution equation in (12) becomesT 1 ∇ · 0 at sufficiently early times in the single-stream regime. This is an ordinary differential equation in time for ∇ · and, importantly, is a linear equation in the displacement. Supplemented with suitable initial conditions, the solution of this differential equation yields the aforementioned Zel'dovich approximation (see the following section). Furthermore, W is exactly zero at all times, iff the initial data depends only on a single-space coordinate (since derivatives in the transverse flow direction vanish in that case): This is the technical reason why the Zel'dovich theory, a linear description in Lagrangian coordinates, is exact in 1D until shell-crossing.
Finally, if the vorticity := ∇ × vanishes initially, as assumed implicitly when using the boundary conditions (10), then the second equation in (12) amounts to a conservation equation for the zero-vorticity condition. This equation and its variants have a long history tracing all the way back to Cauchy in 1815; see [30] for details and e.g. [18,24,54] for contemporary formulations in cosmology.

Theoretical solution strategy
Obtaining theoretical solutions for the Vlasov-Poisson equations essentially require three ingredients as summarised below.
1 Linear analysis (W = 0, M = 0). Linearising the evolution equations (12) around the steady state of the displacement ( , ) = ( , ) − , one obtainŝ where from here on ∇ := ∇ , and we remind the reader thatT = (2 2 /3) 2 + − . The general solutions of these equations are respectively where , , and are spatial constants to be determined by the boundary conditions (the first of the solutions in 17 is formally identical with the one for the linearised density contrast; see e.g. [7,55]). If the boundary conditions (10) are employed, then only the growing solution ∼ is selected and one arrives at the classical Zel'dovich solution [34,56] where ini is the initial gravitational potential. Of course, other solutions (e.g., with terms that decay in ) exist as well [35,57,58], but require a more sophisticated analysis of the boundary layer (see e.g. [59,60]).
2 Shell-crossing solutions (W ≠ 0, M = 0). As long as the flow is still singlestream, the evolution equations for the particle displacements are exactlŷ If we impose the boundary conditions (10), then these equations remain analytic for → 0 [29,30], thereby suggesting that a series representation for the displacement in powers of (the linear growth time), around = 0, is amenable: Here, (1) ( ) = −∇ ini is the first of infinitely many, purely space-dependent Taylor coefficients. We note that this factorisation into spatial and temporal parts persists also for a ΛCDM Universe [24]. Plugging (20a) into (19), one obtains explicit all-order recursion relations [17,18,30,54,61,62] from which ( ) = ∇ −2 (∇ ( ) −∇× ( ) ) follows, and subsequently = ( ) . Note that for sufficiently smooth ini , as it is the case for cosmological initial conditions, there exists mathematical proofs of time-analytic displacements until some finite time [17,30], as well as numerical evidence of convergence at shell-crossing and even shortly after [61]; see section 4.3 for details. Once the displacement is determined to sufficient accuracy, the time of first shell-crossing, denoted with ★ , can be determined by searching for the particle = ★ for which the density blows up for the first time (cf. equation 8), i.e., Of course, not only the particle with label ★ but all particles can be evolved until ★ ; in the following we denote the family of displacements at ★ with ( ) ( ) ★ , truncated at sufficiently large order max .
3 Post-shell-crossing solutions (W ≠ 0, M ≠ 0). Instantly after the fluid bifurcates at shell-crossing time = ★ to three streams at location = ★ , the multistream source term M is nonzero in the neighbourhood of ★ (blue shaded region in Fig. 3). In this case, the evolution equation for ∇ · at ≥ ★ can be formally solved by using the method of variation of constants, i.e., where ★ = ★ ( , )| = ★ ; see [14,44,45,46] for complementary methods applied to the one-dimensional case. Note that the re-occurring terms W and M in (22) are a priori unknowns for times > ★ , since they are functions of the unknown post-shell-crossing displacement (cf. equations [14][15]. However, the flow directions of the matter stream are known until shell-crossing exactly -and at least approximatively shortly after, essentially by an argument of momentum conservation. Thus, the unknowns W and M in (22) can be approximatively determined by extrapolating the shell-crossing solutions to times shortly after shell-crossing. Specifically, for > ★ , we approximate in the first iteration, where • ( , ) := + max =1 ( ) ( ) is the truncated shellcrossing displacement until = max with Taylor coefficients ( ) , the latter determined through (20b)-(20c). We remark that the computation of the multistreaming source M ( • ) is still non-trivial and requires normal-form considerations inherited from catastrophe theory [63,64]; see section 5.3 for details.
Cosmological Vlasov-Poisson equations 13 As first observed by S. Colombi [44], the resulting dark-matter displacements subject to the leading-order approximation for M are fairly accurate, and already overcome the major obstacle of theoretically predicting the second shell-crossing. Suitable higher-order refinements for resolving the regime between the first and second shell-crossing are straightforward [14,45,46]; so are applications beyond 1D (e.g., [65,66]). Finally, predicting the third shell-crossing by theoretical means is in principle possible, but increasingly challenging. For this one would first need to provide boundary conditions at the second shell-crossing, followed by approximating W and M along the lines as summarised above.

Results before and until shell-crossing (M = 0)
We first provide an overview of approximative solution techniques, valid as long as multi-streaming has not yet occurred. Then, in section 4.2 and 4.3, we investigate the mathematics of shell-crossings respectively for simplified and random initial conditions.

Perturbative methods and overview of applications
Approximative techniques to Vlasov-Poisson have a long history in cosmology. Most of them rely on the single-stream description of the Vlasov-Poisson equations, which, by applying kinetic moments to (1), reduce to a system of fluid equations with vanishing velocity dispersion (see e.g. [7] for details). The basic motivation for using such perturbative techniques is, that the larger and/or the earlier the cosmological scales considered, the more these techniques deliver physically meaningful results. This is so for at least two reasons, namely that (1) on sufficiently large scales and/or early times, the matter evolution follows mostly the overall Hubble flow due to the expanding Universe, which justifies generic expansions around linearised field variables, and that (2) on such spatio-temporal scales, shell-crossing effects are yet negligible. The overall literature is vast, therefore we attempt here to only summarise some of the rather recent applications.
Eulerian methods. Instead of seeking a power series representation for the displacement in Lagrangian coordinates, one can as well perturbatively expand the fluid variables in Eulerian coordinates ; for example the density may be represented as a power series in the linear growth time, i.e., ( , The general framework for this is called Eulerian or standard perturbation theory (SPT) [55,67,68,69,70], and there exists explicit all-order recursion relations [68] as well as a variety of SPT extensions (e.g. [71,72,73,74,75,76,77,78,79]); see also [7] for an extensive review. A straightforward application of SPT is to determine so-called loop corrections to -point statistics of the matter density contrast, e.g., to the matter power-and bispectrum which are the Fourier transforms of the 2-point and 3-point correlation functions of , respectively [80,81,82,83,84,85,86]. Furthermore, certain effective approaches based on SPT are capable of incorporating shell-crossing effects to some extent; see section 5.1 for details.
Finally, in Eulerian perturbative approaches, physics beyond Vlasov-Poisson can be straightforwardly implemented. One important example for this relates to the biasing problem, where one seeks functional relationships between the spatial distributions of the dark-matter and visible (baryonic) density; see e.g. [87,88,89,90]. We remark that, so far, the mathematical convergence of the SPT series has not been explicitly demonstrated, except for the simplified cases of one-dimensional [38] and spherical collapse [91]. Nonetheless, convergence for cosmological initial conditions until -but excluding -shell-crossing is likely, mainly as a consequence of the analyticity results in Lagrangian coordinates [61,62]. However, the speed of convergence of the power series for the density is expected to be substantially slower as compared to the series for the displacement in the Lagrangian case, due to the presence of the convergence limiting density singularity at shell-crossing.
Zel'dovich approximation. One of the most important yet simple theoretical models of cosmic structure formation is the Zel'dovich approximation, which is the linear solution of the Vlasov-Poisson equations in Lagrangian coordinates with 3D initial conditions (taking the boundary conditions 10 into account). Specifically, the Zel'dovich approximation is achieved by keeping only the = 1 contribution in the infinite Taylor series (20a) for the displacement. As we have elucidated above, the Zel'dovich approximation essentially boils down to a ballistic description for darkmatter trajectories, which we repeat here for convenience, Although being a linear approximation to Vlasov-Poisson in Lagrangian coordinates, it leads to a highly non-linear prediction for the density in Eulerian coordinates, essentially due to the inherent non-linearity in the mapping ↦ → , thereby providing a surprisingly accurate description of the non-linear gravitational collapse. Indeed, since for (24) , the mass conservation law (8) can be written as [34] ( ( , )) + 1 = (1 − 1 ( )) (1 − 2 ( )) (1 − 3 ( )) −1 (25) in the fundamental coordinate system spanned by ∇ , where 1,2,3 are the eigenvalues of the Hessian of the initial gravitational potential. From the relation (25) it is clear that shell-crossing is reached when any one of the three eigenvalues goes to zero. We remark that the instantaneous vanishing of two or even three eigenvalues is essentially excluded for random initial conditions [94].
The Zel'dovich approximation has several key advantages over the corresponding linear solution in Eulerian coordinates, where the latter predicts that ( , ) = ∇ 2 ini at first order in SPT (see above). Evidently, first-order SPT provides a completely unrealistic prediction for the time of first shell-crossing which is achieved only at → ∞ for sufficiently smooth initial conditions. By contrast, we know by now that the first shell-crossing occurs at times = ★ 1 for cosmological initial conditions (the precise time depends on spatial resolution; see section 4.3), while the Zel'dovich approximation predicts ★ to about 20% accuracy. Another advantage of the Zel'dovich approximation is that it correctly predicts the emergence of a three-stream regime in the phase-space after shell-crossing (see e.g. Fig. 4), which in SPT is never achieved, even at arbitrary high orders. Important (and rather recent) applications of the Zel'dovich approximation and higher-order refinements include: -Generating initial conditions for cosmological simulations [95,96,97,98]; see also [31,99,100] for higher-order refinements within the framework of Lagrangian perturbation theory [35,24,56,57,101,102,103,104,105] (LPT), which are predictions based on perturbative truncations of equation (20). -Predicting the change of photon momenta emitted from observed luminous tracers of dark matter, due to peculiar motion and the expansion of space ("redshift-space distortions") [15,20,92,93,106,107,108]. -Theoretically modelling the spatial distribution of cosmic neutral hydrogen through its 21cm emission line [109,110,111,112,113]. -Investigating topological features of the cosmic large-scale structure, such as clusters, sheets, and filaments (cf. Fig. 1) [114,115,116,117,118,119,120]. -Modelling the baryonic acoustic oscillation feature imprinted in statistics of cosmic structures [121,122,123,124,125,126,127]; and improving statistical predictions of the large-scale structure by pairing/emulating techniques with numerical simulations [128,129,130,131,132,133,134].

Analytical shell-crossing solutions for simplified initial conditions
Even without multi-streaming (M = 0 in equation 12), the Vlasov-Poisson equations are highly non-linear, especially near shell-crossings where → ∞. Nonetheless, there are a few purely analytical shell-crossing solutions. Analytical solutions play an important role in theoretical cosmology, and can be used as test problems for state-of-the-art simulation techniques (e.g., [135,136,137,138,139,140,141]).
Quasi-one dimensional collapse. The Zel'dovich solution (24) becomes exact until shell-crossing, if the initial data depends only on one space coordinate, since in that case the displacement can only depend on the same space coordinate, but not on the others; hence W = 0 (see discussion after equation 14). Figuratively, the one-dimensional problem embedded in 3D can be thought of as stacked mass sheets/curtains along the coordinate direction. Suppose now that the initial data depends mostly on one coordinate, say 1 , but also depends weakly on the other coordinates. Mathematically and actually without loss of generality, the corresponding initial gravitational potential can be set to where > 0 is a perturbation parameter that controls the smallness of the deviation from the one-dimensional problem, while ini is an arbitrary 2 -periodic function. This is thus an intrinsically three-dimensional problem; within the above picture, the mass sheets have now a weak functional dependence in all coordinate directions. The depicted collapse problem with initial data (26a) could also be solved directly with the recursion relations for the generic case (equations 20 In the ϵ y;z =ϵ x ¼ 1 case, the phase-space structure is highly stretched along the velocity axis, which reflects a noticeable acceleration of the inward mass flow near x ¼ 0, similarly as in spherical collapse [14,15]. This highly contrasted dynamical behavior slows down LPT convergence near ϵ y;z =ϵ x ¼ 1 and even the tenth-order calculation is insufficient (see also Ref. [46] for a discussion about the spherically symmetric case). nd approach, that we denote by Q1D, assumes e-dimensional dynamics, following in the footsteps 1], i.e., jϵ x j ≫ jϵ y;z j. In this case, one takes the exact of the one-dimensional problem along the x axis Zel'dovich approximation as the unperturbed rder state: Ψ ð0Þ A ðq; tÞ ¼ Ψ ini x ðq x ; tÞδ A;x , A ¼ x, y, z, being the Kronecker delta. Transverse fluctuations dered as small first-order perturbations to this setup, The perturbansion is then performed by keeping terms proporϵ i y and ϵ j z up to second order, i þ j ¼ 2 (so in this go one order beyond [41]), while keeping terms nal to ϵ k x up to tenth order, k ¼ 10. v simulations and phase-space structure.-To ively investigate the dynamics of our system up crossing, we perform high resolution simulations state-of-the-art Vlasov-Poisson code COLDICE is code uses a tessellation, i.e., tetrahedra, to t the phase-space sheet. The vertices of the on form initially a homogeneous mesh of resos (which corresponds to 6n 3 s simplices) and then agrangian equations of motion, Eqs. (1) and (2). son equation is solved by fast-Fourier transform on Cartesian grid of resolution n g , after deposit of the ace sheet on the grid using the method of Franklin kanhalli generalized to linear order [43][44][45]. A of simulations, as listed in Table I, were performed Einstein-de Sitter cosmology. 1 shows representative examples of phase-space at shell-crossing time. As the ratios ϵ y;z =ϵ x the Zel'dovich approximation, which is exact rictly one-dimensional case ϵ y;z =ϵ x ¼ 0, starts to ignificantly from the simulation, as expected. The diction provides a substantial improvement, with lent match of the simulation measurements for 1 (top panel). Still, it cannot catch up to the shellstructure when ratios ϵ y;z =ϵ x approach unity and bottom panels). However, with a systematic on of all the contributions up to tenth order, the ediction improves considerably and accurately es the shell-crossing structure seen in the simuor most of the parameter space, except in the of ðϵ y =ϵ x ; ϵ z =ϵ x Þ ¼ ð1; 1Þ (bottom panel).
In the ϵ y;z =ϵ x ¼ 1 case, the phase-space structure is highly stretched along the velocity axis, which reflects a noticeable acceleration of the inward mass flow near x ¼ 0, similarly as in spherical collapse [14,15]. This highly contrasted dynamical behavior slows down LPT convergence near ϵ y;z =ϵ x ¼ 1 and even the tenth-order calculation is insufficient (see also Ref. [46] for a discussion about the spherically symmetric case).  analytical proofs of convergence are available. Instead, it turns out to be convenient to employ simple multi-scaling techniques, which in the present case amounts to impose a refined Ansatz for the components of the displacement = − [142] where is the displacement part of the one-dimensional problem, while is a leading-order correction to the displacement. Supplemented with the boundary conditions (10), it is straightforward to solve the master equations (19) to successive powers in . For example, at order 0 , one recovers the Zel'dovich solution = − sin 1 , while at order 1 , one finds refined recursion relations for the Taylor coefficients appearing in ( , ) = ∞ =1 ( ) ( ) (see equation 36 in [142]). Most importantly, the refined recursion relation for ( ) can be used to investigate, by analytical means, the asymptotics of the Taylor series at large orders, from which it follows that is an entire function of time . Thus, the quasi-one-dimensional problem is analytically solvable to arbitrary high precision. See Fig. 5 for comparisons of various analytical solutions against numerical simulations (black dotted line): the blue (dot-dashed) line denotes the result using the Ansatz (26b) up to second order in (derived in [143]), while the red-dashed line is obtained using (20) up to the 10th order in the series. The left panel in Fig. 5 denotes the stated quasi-one-dimensional problem, while the right panel shows the collapse for tri-axial sine-wave initial conditions (see [143] for the specific set-up). Note that due to the use of trigonometric functions in the initial data, the theoretical solutions for the displacement are combinations of trigonometric functions as well -thus, no numerics are required for the shown theoretical solutions.
Quasi-spherical collapse. There exists an exact parametric solution to the collapse of a homogeneous, spherically symmetric over-density (sometimes called the top-hat model) [144]. This solution is however not obtained from the Vlasov-Poisson equations, but instead by a special solution of the Friedmann equations. Nonetheless, it can be shown that this simplified model can also be realised within the context of a Cartesian formulation of the Vlasov-Poisson equations [145,146,147].
By applying equivalent multi-scaling techniques in Lagrangian coordinates as reviewed above, it has been shown that the matter collapse of arbitrary small departures from spherical symmetry also constitutes an exact solution of (19) until shell-crossing [91]. In that case, the initial gravitational potential can be taken to be where is a positive [negative] constant amplitude of a spherical over-density [underdensity/void], and ini is an arbitrary function that reflects the departure from spherical symmetry. The mathematical form of ini suggests the solution Ansatz where is a purely time-dependent function that models the purely spherical collapse, while reflects a perturbative departure therefrom. In the present case, both and can be represented individually by convergent time-Taylor series [91]. However, the speed of convergence is significantly slower as compared to the quasi-one-dimensional case (cf. the LPT predictions in left versus right panel of Fig. 5). The reason for the slow convergence is the presence of nearby singularities at an ( ) amount of time after shell-crossing. Actually, in the spherical case with → 0, the speed of convergence is even worse and furthermore comes with a singular velocity at shell-crossing (see also Fig. 1 in [143] for the highly related triaxial symmetric case, or [148] for a complementary analysis with similar conclusions). Nonetheless, we should stress that such singularities are basically manmade: a spherical top-hat over-density is a vastly simplified collapse model that has zero probability to occur in a Universe with random initial conditions [94].

Shell-crossing solutions for cosmological initial conditions
Investigating solutions until shell-crossings for random initial conditions is in general only feasible by semi-analytical avenues. This is so, since the Fourier transforms needed to solve for the Helmholtz decomposition (13) cannot be performed explicitly. Instead one typically resorts to Fast Fourier transforms, thereby determining the Taylor coefficients of (20) on a three-dimensional numerical grid with collocation points.
However, there are also ways to prove, by entirely theoretical means, that the solutions (20) are time-analytic at least until some finite time, which goes back to a seminal paper by Zheligovsky & Frisch [30]. In the following we sketch the essence of a normed-space proof, from which a lower bound on the radius of convergence of the series is obtained. Then, in a subsequent paragraph, the actual radius of convergence is determined by employing numerical extrapolation techniques. 18

Time-analyticity of dark-matter displacements. Instead of proving the convergence of
, it is more instructive to first prove the convergence of the tensorial gradient therefrom, i.e., ∇ = ∞ =1 ∇ ( ) , with Taylor coefficients where ( ) and ( ) are given by the recursion relations (20b) and (20c), respectively. Before proceeding let us briefly introduce the ℓ 1 norm denoted with f := |f | < ∞, for any periodic tensor, vector or scalar function f ( ) = f exp(i · ), where f are its Fourier coefficients. The corresponding function space has the algebra property f g ≤ f g , and furthermore bounds the operator ∇∇∇ −2 to unity from below. Using these properties as well as the explicit forms of ( ) and ( ) (equations 20b and 20c), the ℓ 1 norm of (28) turns into the polynomial inequality [30] ∇ ( ) ≤ ∇∇ ini 1 + 12 , as well as multiplying (29) with and summing over from one to infinity, one gets ≤ ∇∇ ini + 12 2 + 6 3 , which is equivalently (see [17] for a graphical representation as well as for the ΛCDM case). Initially for = 0, this polynomial behaves as − for small and asymptotically as 3 ; thus, ( = 0) has three roots with one at phys = 0. For small positive times, gets shifted upwards, thereby drifting the root phys to the right. This branch from = 0 until phys marks the physical regime where > 0 is bounded (i.e., the Taylor coefficients of do not blow up). However, at some critical time, denoted with c , this boundness property will be lost as the asymptotic behaviour 3 kicks in. It is found that this critical time is achieved precisely when the root phys merges with another root; therefore c can be determined by setting the discriminant of ( ) to zero, leading to [30] Thus, the displacement is guaranteed to be time-analytic between 0 ≤ ≤ c , while the value of the lower bound c crucially depends on the Hessian of the initial gravitational potential. We note that the above bound can be slightly improved [30,31], especially when assumptions about the regularity of ini are imposed.

Shell-crossing and radius of convergence.
To determine the actual radius of convergence of ( , ) = ∞ =1 ( ) ( ) for random initial conditions, one has to resort to numerical tests or extrapolation methods. One of such tests is to verify whether theoretical predictions, such as time and location of the first shell-crossing, converges for successively higher orders in the perturbative truncations of the displacement series. That is one searches for the location where the perturbatively truncated Jacobian vanishes for the first time (cf. equation 21) [61]. For each truncation order , one then obtains perturbative estimates of the shell-crossing time and location, denoted respectively with ( ) ★ and ( ) ★ . By monitoring the trend of these estimates at successively high orders , one then obtains evidence for or against (weak) convergence.
In [61], the recursive displacement solutions (20) as well as { } have been realised on a spatial grid with collocation points for cosmological random initial conditions. The corresponding code is parallelised (MPI+threads), and de-aliases cubic functions in the displacement. Most analysis in [61] was performed with a grid resolution of = 256 3 , for which convergence of the shell-crossing time ★ accurate to three significant digits was typically reached between orders = 6 − 10.
As another convergence test, we show in the left panel of Fig. 6 the one-point probability density function (PDF) of the residual Δ { } := { } − { −1} from all Lagrangian grid points, evaluated at the time of first shell-crossing (here at ★ = 0.0813 for a standard ΛCDM cosmology with = 256 3 ). Evidently, Δ { } peaks sharper at the origin for higher orders, indicating that truncation errors decrease rapidly at increasing order , which is an expected feature of a convergent series.
Next we determine numerically the radius of convergence, for which it is useful to consider the closely related series of the ℓ 2 norms of displacement coefficients, i.e., To investigate its convergence, one can perform d'Alembert's ratio test which states where is the radius of convergence -when that limit exists. A standard way to estimate is the Domb-Sykes plot [149], where one draws ( ) / ( −1) versus 1/ , from which 1/ follows by extrapolation to the -intercept (see also [31,150] for related applications). In the right panel of Fig. 6 we show the Domb-Sykes plot for the spatial coefficients ( ) evaluated at the Lagrangian locations of first shell-crossing for three different resolutions, = 512 3 , 256 3 , 128 3 (respectively in blue, orange and green; top to bottom), while the length of the grid is fixed at = 125 Mpcℎ −1 5.51 × 10 24 m. The shaded areas denote its 32 and 68 percentiles obtained from 5 random realisations, while the solid lines denote its median. For sufficiently large Taylor orders (  7) the ratios settle into a linear behaviour, thereby justifying a linear extrapolation to the -intercepts (dotted lines), which, using (34), leads to the radius of convergence . In all cases considered the (temporal) value of is roughly 30 − 40% larger than the shell-crossing time (see table 1 in [61] for specific values), thus indicating that the mathematical validity of (20) surpasses vastly its physical validity (since 20 is only valid until shell-crossing). This is an understood phenomenon, for example the one-dimensional and quasi-one-dimensional shell-crossing solutions have an infinite radius of convergence (see section 4.2).
Obviously, depends on the chosen resolution: since we keep the grid length fixed, this mostly reflects the known dependency of theoretical solutions on UV physics beyond the Nyquist frequency; see e.g. [7,31,151]; see however [62] for a recent investigation of residual aliasing effects.
Finally, the linear behaviour of the ratios of coefficients at large orders suggests that the convergence limiting singularities have a local behaviour of [61] ( , ) ∝ ( − ) , where is a singularity exponent. In this case, the ratio of coefficients satisfies the linear relationship (cf. [150,152]) Thus, the slope of the linear extrapolations reveal the singularity exponents, which in the considered cases are = 0.61, 0.38, 0.56 respectively for the resolutions = 128 3 , 256 3 , 512 3 [61]. Since these exponents are positive non-integers and smaller than unity, it follows that the first time derivative of ( , ) blows up at = .
Concluding this section, the recursive solutions (20) for = ∞

=1
( ) appear to deliver physically and mathematically meaningful results until the time of shellcrossing ★ , while the various bounds and limits are related through 0 < c < ★ < , where c is an entirely analytical estimate of the radius of convergence, while is the actual radius of convergence determined through numerical extrapolation techniques.

Numerical and analytical post-shell-crossing methods (M ≠ 0)
Solving the Vlasov-Poisson equations beyond shell-crossing is a challenge that is usually attempted by numerical -body simulations, where (very massive) macroparticles are employed to coarse sample the dark-matter distribution in phase-space (see e.g. [184,185] for reviews). However, these simulations can be prone to errors introduced through the discrete representation of the phase-space distribution (e.g. [136,186,187,188]).
Recently, however, there have been encouraging avenues at various fronts using elementary methods as well as (semi-)analytical and numerical techniques. While this review predominantly focuses on the progress related to elementary methods for Vlasov-Poisson (section 5.3), in the following two sections, we attempt to summarise important findings related to semi-analytical/effective avenues (section 5.1) as well as to novel simulation techniques (section 5.2).

Semi-analytical/effective methods
Simple (time-)Taylor expansions, such as (20) in Lagrangian coordinates, cease to be valid at the time of first shell-crossing [46,14]. Of course, the same can be said when seeking theoretical solutions to Vlasov-Poisson in Eulerian coordinates , where the density and velocity can be Taylor expanded until shell-crossing with equivalent methods as reviewed in section 3.2 (dubbed SPT; see also e.g. [7,191]).
However, there also exist methods that seek to incorporate multi-streaming effects in a rather effective way [28,32], notably the so-called coarse-grained perturbation theory [157] as well as the effective theory of large-scale structure [158,159,160]. These approaches have a few properties in common, in particular that they are typically not designed to solve for the dark-matter evolution at the deterministic level; instead predictions are of statistical nature (e.g., for the power spectrum of the matter density).
Another important characteristic of effective approaches is that they attempt to circumvent the obvious shortcomings by first "smoothing out" the small-scale physics -which can be strongly influenced by shell-crossing/multi-streaming dynamics, beyond a spatial cut-off scale 1/Λ. Subsequently, the large-scale physics is handled with Taylor/perturbative expansions, while the small-scale/multi-streaming physics is incorporated through certain outputs from numerical simulations (e.g., velocity dispersion), or even through cosmological observations [192,193].
At the technical level, this can be achieved by introducing a spatially smoothed (large-scale) distribution function Λ ( , , where is the distribution function that solves equation (1), while Λ is a Window function 22 Cornelius Rampf which is usually taken to be a normalised Gaussian, i.e., Λ ∝ exp(−Λ 2 | | 2 /2). Then, by taking kinetic moments of the resulting large-scale Vlasov-Poisson equations, one obtains fluid-type equations for an effective large-scale fluid [158,159] where , and are respectively the zeroth, first and second kinetic moments of the distribution function (see e.g. [7]), while [X( )] Λ := ∫ Λ ( − )X( ) d 3 for any scalar, vector or tensor function X( ).
Equations (37) state respectively mass and momentum conservation of a fluid that is subject to a source term due to small-scale physics, here introduced through the velocity dispersion tensor . Note that generally also includes a non-zero pressure contribution which physically arises due to multi-streaming. How precisely is incorporated, depends on specific model details for which we kindly refer to the original references.
Finally, an entirely different, semi-analytical approach is the so-called kinetic field theory, which is a microscopic statistical field theory that is intimately linked to the classical Hamiltonian approach (e.g. [161,162,163]). The central building block is the generating functional involving the -particle action, which physically encodes the classical transition probability from the initial and final field state. To incorporate particle interactions (relevant for multi-streaming), perturbative expansions and/or approximations are necessary. Furthermore, applying suitable limits, particle shotnoise is removed and, therefore, the theory should deliver meaningful approximations to statistical predictions of Vlasov-Poisson. So far, the corresponding methods have been only worked out at leading order, therefore statements about convergence can not be made at this stage.

Numerical simulations techniques
Standard -body methods solve the Hamiltonian equations of motion (5) for a set of macro-particles, where the continuous phase-space distribution is effectively coarse-sampled. There are various types of -body simulations that essentially differ in the way how the Poisson equation is solved, e.g., using standard particle-in-cell / particle mesh [164,165] or brute force methods [166,167], codes with adaptive mesh refinements [168,169,170,171,172,173], tree algorithms [174,175,176] and, finally, combinations of various algorithms [96,165,177,178,179,180]; see e.g. [181,182,183] for extensive reviews about -body methods.
In the following, we provide an overview of recent numerical approaches that aim to retrieve fine-grained details of the dark-matter phase-space, and thus going, in one way or another, beyond standard -body techniques.
Phase-space reconstruction. With the standard -body technique, quasi-continuous estimates of observables, such as the matter density or (mean) velocity, can be obtained by applying suitable smoothing operations on the set of discrete tracer particles (e.g., the cloud-in-cell interpolation that assigns a density cloud around particles). Recently, however, there have been novel approaches that are able to reconstruct the dark-matter phase-space to much higher accuracy, while at the heart still relying on the -body technique.
One of these new approaches is dubbed GDE which is short for geodesic deviation equation [135,195]. There, on top of the standard equations of motion, one solves for the evolution of the tangent space in the 6D phase-space around each tracer particle, thereby locally recovering some details of the dark-matter sheet. This approach has turned out to be fruitful at the very final stages during halo formation, where the number of streams can easily go into the several hundreds. On the other hand, the approach can suffer artificial fragmentation in certain (warm) dark-matter scenarios, which arises during the rather early gravitational evolution [196,197].
Another approach is the sheet method, which employs interpolation/tessellation techniques to recover the dark-matter phase-space sheet from the set of discrete tracer particles [139] (see also [136,137,138]). The sheet reconstruction (with isotropic mesh refinements) provides continuous density estimates and thus an improved computation of particle forces during the subsequent stages of the gravitational evolution. This also alleviates the aforementioned problem of artificial fragmentation. However, the dark-matter sheet reconstruction fails at the very final stages, especially inside halos, due to the increased complexity of phase-space dynamics. Recently, the sheet and GDE approaches have been married, thereby allowing to efficiently exploit the strengths of both methods [141].
Pure phase-space tessellation methods. Instead using tracer particles, the numerical code C DICE employs simplices (tetrahedra) to tessellate the continuous dark-matter sheet [140]. Fairly similarly as in the above sheet method, the forces are obtained from continuous density estimates and subsequent interpolations to an undistorted grid. Furthermore, C DICE allows for anisotropic mesh refinements, thereby being able to resolve the phase space to high precision. Unfortunately, the complexity of the darkmatter sheet becomes computationally too demanding at the late stages during halo formation. Nonetheless, important lessons can be drawn from using such advanced approaches to Vlasov-Poisson; see in particular [194] for explicit comparisons of C DICE against standard -body codes during the early violent relaxation phase of halo formation (see also [188]).

Vlasov-Poisson simulations in 6D.
There exist also full Vlasov-Poisson simulations that are not constrained in evolving the dark-matter sheet [198,199]. Such avenues are especially important when e.g. investigating the phase-space in the presence of warm components (e.g., neutrinos, or warm dark matter) [200,201].
Vlasov-Poisson through Schrödinger-Poisson. Lastly, an effective Vlasov simulation is realised by exploiting the correspondence to the quantum mechanical Schrödinger-Poisson equations [202]. Indeed in [203] it has been proven for the 1D case that this correspondence is exact in the limit ℏ → 0. By contrast, keeping ℏ small but finite in Schrödinger-Poisson delivers meaningful approximations, where ℏ acts now as an effective coarse-graining scale in the phase-space; see e.g. [153,154,155,156]. 24

Mathematical post-shell-crossing analysis
Turning back to the theoretical solution scheme for the Lagrangian master equations (22), the leading-order behaviour of the involved source terms shortly after shell-crossing time ★ is given by W ( ( , )) W ( • ( , )) and M ( ( , )) M ( • ( , )), which is, strictly speaking, exact at = ★ , and approximatively thereafter, since constitutes only a solution of Vlasov-Poisson until ★ .
Let us show how M can be determined, which can also be written as where ì Θ would correspond to the standard Heaviside function in one space dimensions, while in 3D, it may be considered as a gravitoelectric field. Although not strictly necessary, for the proceeding calculations it turns out to be convenient to employ ì Θ, although similar (but substantially longer) arguments hold when using the Dirac delta instead. Evaluating the integral in (38) boils down to determining the length of the branches when the argument of ì Θ is positive. Obviously, for this one needs to determine the roots of which is a transcendental equation for which exact solutions are generally not available.
To proceed, one considers normal forms, which are spatial Taylor expansions of the map • ( , ) around the location ★ of first shell-crossing, combined with exploiting the expected topology of • ( , ) for times shortly after shell-crossing. What do we mean by the last statement? Consider the left panel in Fig. 3, which shows the map in 1D: at shell-crossing time (here ★ = 1 at ★ = 0), the map has an inflection point with local behaviour ∼ 3 around the shell-crossing location; shortly later (here shown at = 1.3 in blue) the local behaviour of the map in 1D is [14,44,45] • ( , ) where 1,3 are positive Taylor coefficients. Here, N ( ) corresponds to the said normal form in the one-dimensional case, and note specifically that a term ∼ ( − ★ ) 2 is absent (it can be removed by a suitable Galilean transformation). Employing (40) to determine the roots of • ( , ) − • ( , ) = 0 is straightforward, and in fact a good approximation for times shortly after the first shell-crossing (at later times higher refinements become necessary; see [14,44,45]). Generalisations of the above to the three-dimensional case are in principle straightforward yet tedious and, so far, have not been reported in the literature. Once M ( • ( , )) is determined, the post-shellcrossing displacement is readily obtained using equations (22); see e.g. equation (18) in [14] for explicit solutions in the 1D case. These purely analytical solutions are shown in Fig. 7 and compared against numerical simulations at high resolution (see footnote 1). Panel (a) shows the phase-space Finally, a non-trivial singularity appears when the initial velocity is not exactly point-symmetric, here for a simplified model with ini = − sin + 0.1 sin 4 − 0.12 sin 6 : the particle that is initially at = 0 will remain there until shell-crossing ( = 1), at which time an asymmetric multi-streaming force arises, essentially as a consequence of momentum conservation, that leads the = 0 particle to begin moving with ( = 0) ∼ 3 [14] (panel d in Fig. 7; see also [46] for similar conclusions in a slightly different set-up). Of course, the above initial data just serves as a simplified model to investigate realistic collapse scenarios: there, the random nature of initial perturbations will always distort this point-symmetry and thus, the appearance of the presently described singularity is a generic phenomenon, expected for realistic (cosmological) initial conditions. We remark that to determine this highly non-trivial feature depicted in panel (d) of Fig. 7, the invariance (9) of Vlasov-Poisson under non-Galilean transformations can be used.
Concluding this section, standard Taylor expansions of the displacement cease to be valid after shell-crossing due to the emergence of non-analytic behaviour in the phase-space, however by now there are promising avenues that allow to enter into the post-shell-crossing regime while unveiling intrinsic properties of Vlasov-Poisson.

Related problems in plasma physics
Vlasov-Poisson is also relevant in plasma physics, in particular when a collisionless, electrically charged medium is exposed to electrostatic interactions. In what follows we discuss various plasma problems and lay out connections to the cosmological case.
One component plasma model (OCP). Possibly one of the simplest plasma models consists of a set of electrons with charge and mass that interact only through electrostatic forces. In order to warrant charge neutrality, these electrons are immersed in a rigid uniform background of opposite charge, which in the present case is −¯= − / where is the spatial volume under consideration. Possible applications of the OCP includes the (rough) modelling of the interior of stars, in particular of white dwarfs (see e.g. [204,205]). Adopting somewhat the above notation, the equation of motion for the th electron, formulated in physical coordinates in 3D and employing the standard time , reads (see e.g. [204,206] where ∇ = corresponds to the electrostatic field, and we note that the vacuum permittivity is 0 = 1/(4 ) in Planck units. By applying similar arguments as in the cosmological case, one may introduce a continuum description for the OCP by In the following we discuss the similarities and differences of these equations as compared to their cosmic counterparts. While the former makes use of the standard time, the cosmic Vlasov-Poisson equations (6) are formulated in the scale-factor time = (since in the cosmic case, time-analyticity holds in ∼ 2/3 but not in ). Also, equations (6) feature a Hubble drag term ∼ stemming from the present choice of spatio-temporal coordinates (note that the Hubble drag term formally disappears when using a new time variable defined through d = 2 d [56,207]). Furthermore, Cosmological Vlasov-Poisson equations 27 equations (6) have a minus sign in front of ∇ while in (42) there is a plus sign; this flip in sign just reflects the change of charge in the gravitational attractive and electrically repulsive case. Finally, observe that equations (42) are invariant under the non-Galilean transformation → + ( ), where ( ) is an arbitrary function of time; this invariance is in close resemblance with the one of the cosmological case that we discussed around equation (9).
Perturbative techniques applied to variants of equations (42) are frequently encountered in plasma physics (e.g. [208,209]), however, we are unaware of explicit use of Lagrangian coordinates (though they are implicitly assumed for equation 41 in the discretised sense). In this context note that for a single electron particle at current position ( ) exposed to its own electric field, the Poisson potential is in Fourier space for ≠ 0 (cf. equation 1 in [209]). By contrast, for a set of electrons parametrised in Lagrangian coordinates , the Fourier transform of the Poisson potential is easily obtained from (42); it reads Finally, note that as long as the plasma is monokinetic (single-beam), equations (42) can equivalently be formulated in Eulerian coordinates as a set of closed fluid equations where now all fields depend on the Eulerian coordinate and time . See e.g. [210,211] for the application of perturbation techniques to equations (43).
Beam-plasma (two-stream) instability. A classical problem in plasma physics is when a beam of low-density, mono-energetic electrons is exposed to a cold thermalised, large-density plasma [206,212,213]. In the so-called single-wave model thereof, certain assumptions about the involved velocities are employed such that the plasma response is effectively non-resonant; therefore, the plasma can be treated as a bulk with respect to the propagation of the electron beam, while the plasma response can be incorporated in a form of a real dielectric function (see e.g. [214,215]). Similarly as in the cosmological case, also here the Vlasov-Poisson equations are directly related to a Hamiltonian principle [217,218], and thus can be reduced to Newtonian-type equation of motion coupled to a Poisson equation. This set of equations has been provided by [212] in the case of discretised electrons in 1D, which, adopting the same notation as above, can be written as where denotes the current position of the th electron, while p and b are respectively the charge densities of the plasma and the beam. It is usually assumed that := b / p 1 is a perturbatively small control parameter of the kinetic wave-particle interaction. Furthermore, it is assumed that the electron beam frequency is (roughly) equal to the plasma frequency, which ensures that the dielectric of the plasma is nearly vanishing [215], thereby justifying linear temporal expansions applied to the dielectric. This allows to incorporate the effect of the plasma density in the Poisson 28 Cornelius Rampf   Fig. 8 Phase-space for the beam-plasma instability, here in suitably rescaled coordinates. Figure from [212]; note that the filamentary tails of the distribution have been cut off for reasons of better visibility. equation by a linear dielectric, and, at the same time, cast the Poisson equation into a simple evolution equation [216]. By contrast, the density of the electric beam is just a superposition of the charged particles, (see e.g. [217]). Thus, the equations of motion for the beam-plasma instability are formally identical with the one of the OCP, equations (41), however except with the addition that in the present case the Poisson equation receives a space and timedependent dielectric function (see e.g. equation 14 in [212]). We show in Fig. 8 the corresponding phase-space of the beam-plasma instability: At early times, the wave amplitude changes exponentially with the linear growth rate; this is thus a regime which can be determined using standard perturbative techniques [212]. At later times, the electrons are getting trapped, which in the figure is exemplified by the appearance of multi-beam regions that are sloshed back and forth. Similarly as in the OCP case, one may employ a continuum description for investigating the beam-plasma instability. For this one can introduce the Lagrangian maps for the beam and plasma with ↦ → ( , ), where =p,b (see also [220,221] for alternative paths to a continuum description). The beam and plasma velocities are then defined with ( , ) = ( , ), and the Vlasov-Poisson equations are = (∇ ) , (∇ 2 ) = +4 p ( ( , )) + 4 b ( ( , )) for the components =p,b, where we have defined (∇ ) = ∇ ( ( , )), and mass conservation implies b ( ( , )) = ∫ D ( ( , ) − b ( , ))d , and similarly for the plasma. Interestingly, equations (46) have a direct cosmological counterpart (in 3D) where they resemble the gravitationally coupled Vlasov-Poisson system of dark matter and visible matter/baryons [19,222] (but only in the large-scale and late-time limit where baryons are effectively collisionless and pressureless).
The mathematical methodology could also be extended to three space dimensions (see e.g. [206]), as well as to the case of multiple cold beams, by e.g. adopting the formalism laid out in [223]. Specifically, for cold beams, one essentially needs equations of motion (44) in the continuum limit coupled to a Poisson equation, where the latter then involves the superposition of charge densities (45).
The weak warm beam instability. Finally, another related plasma problem concerns the wave-particle description of Langmuir waves [206,214,219], where the latter is the response of a thermalised plasma when an electron beam is injected. The corresponding description is achieved when splitting the charged particles into socalled bulk and tail parts, which are subsequently dealt with differently. Here, the tail [bulk] part is the set of particles that is [not] resonant with Langmuir waves. The beam appears as a bump on the tail of the overall velocity distribution function, and the so-called bump-on tail instability is excited through a region near the bump where the velocity distribution has positive slope. The tail particles are dealt with similarly as above, involving a source term as in the Poisson equation (45). By contrast, for the bulk particles, a quasi-ballistic approximation of particle trajectories is employed [206,221,224], which is obtained by setting for the trajectories = + + ini + ini in our notation, where ini and ini are the initial mismatches in position and velocities within the multi-beam-multi-array (a collection of monokinetic beams that provides a multi-valued velocity distribution). Observe here the structural similarities of this ballistic approximation compared against the Zel'dovich solution (18) or its higherorder generalisation (20) in the cosmological case. We deem it possible to investigate the nature of the singularities (such as in the dielectric function of the plasma), possibly also involving methods adopted from catastrophe theory; see section 3.2 and specifically section 5.3.

Summary and outlook
Numerics should not be more than a few steps ahead of theoretical approaches. While this statement can be easily justified, it can be quite a challenge to accommodate in practice. This is particularly the case for Vlasov-Poisson descriptions that, due to the absence of collisions, have the strong tendency to form instabilities and singularities. Solution methods for the cosmological Vlasov-Poisson equations are typically reserved for exploiting physically distinct regimes: theoretical calculations for the very early times of the gravitational evolution -and brute-force numerical simulations that shed light into the highly non-linear regime of cosmic structures. Recently, however, the gap between theoretical and numerical methods decreased by exploiting at least three complementary avenues.
The first to mention are perturbative techniques, either by exploiting low-order truncations especially in Lagrangian space (section 4.1), and/or by using effective methods to Vlasov-Poisson (section 5.1) that circumvent some of the obstacles, albeit in a rather pragmatic way. Still, these methods turn out to be extremely useful when interpreting data from current and forthcoming cosmological surveys.
The second avenue is novel simulation methods (section 5.2), that uncover much more details of the dark-matter phase-space as compared to standard particle-incell/ -body simulations. One of the simulation methods is solving for the tangent space of tracer particles (governed by a geodesic equation). Other methods employ phase-space tessellation techniques to obtain the dark-matter sheet, either constructed from an -particle distribution, or by following the evolution of the vertices of simplices (tetrahedra).
The third avenue addresses directly the mathematical skeleton of Vlasov-Poisson, which is made possible by exploiting the time-analyticity of dark-matter trajectories at sufficiently early times (section 4.3). Straightforward (perturbative) expansion techniques lead to the first non-trivial shell-crossing solutions, which are the phasespace locations when the number of dark-matter streams (beams) changes due to gravitational interactions. These expansion techniques are employed in Lagrangian coordinates, and allow one to analyse the dark-matter phase-space in a controlled set-up at spatio-temporal locations when the density is infinite (see Figs. 3, 5 and 6).
In general, analyticity in the dark-matter phase-space is lost directly at the first shell-crossing, due to the appearance of non-differentiable features in the particle accelerations (panels c and d in Fig. 7). However, thanks to a combination of an iteration technique (in the spirit of Picard's counterpart) as well as adapting Arnold's catastrophe theory, particle trajectories can be followed through shell-crossings by elementary means.
These theoretical methods are very accurate for times shortly after the first shellcrossing, essentially agreeing with independently performed high-resolution simulations. This agreement degrades at later times (Fig. 4), which could be rectified by including higher-order refinements. Nonetheless, resolving the phase-space in darkmatter halos at late times, by elementary methods, remains a major challenge.
Future directions on the cosmological side include the pairing of theoretical methods with simulations or machine-learning techniques. Also, in the longterm, the gained theoretical understanding will necessarily lead to improved (numerical) predictions for cosmic structure formation.
With regard to the plasma case, we have laid out several intriguing connections to the cosmological case of Vlasov-Poisson. In particular, one may apply the continuous methods to the beam-plasma instability in the cold limit. Applications to the (warm) case of multi-beams are in principle straightforward, with the potential in contributing to the fundamental understanding of instabilities and chaos.
In the past, there has been a continuous transfer of knowledge between plasma physics and cosmology. In particular, many numerical simulation techniques in cosmology have their roots in plasma physics. It is in the interest of scientific advancement that the multi-disciplinary transfer of knowledge thrives now and in the future.