Cosmological Vlasov–Poisson equations for dark matter

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.


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. Aghanim et al. 2018;Arcadi et al. 2018;Smorra et al. 2019;Kunz et al. 2016), 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 dark-matter distribution f = f (x, p, t) evolves in the six-dimensional phase-space, where is the matter density contrast defined through the matter density Here, a = a(t) is the cosmic scale factor that encodes the mean expansion of the Universe, and is subject to the Friedmann equation (̇a∕a) 2 = 8 Ḡ∕3 − Kc 2 ∕a 2 + Λc 2 ∕3 , where K is a uniform curvature of space and Λ is a cosmological constant associated with dark energy, while ̄∼ a −3 is the mean matter density which is diluting in time, due to the so-called Hubble expansion of the spatial volume of the Universe. The Friedmann equation determines the functional relationship between a and t for given matter/energy content, and pins down the speed of the Hubble expansion. For example, for the Einstein-de Sitter Universe (O'Raifeartaigh et al. 1933;Bernardeau et al. 2002), which is a simplified cosmological model with K = 0 = Λ and a common starting point for theoretical derivations, we have a ∝ t 2∕3 . For the observationally preferred " Λ cold dark matter" ( ΛCDM) model (Peebles 1993;Hinshaw et al. 2013) which is spatially flat ( K = 0 ) and has a cosmological constant of Λ ≃ 1.1 × 10 −52 m −2 in Planck units (Aghanim et al. 2018), we have a ∝ Λ −1∕3 sinh 2∕3 (3 √ Λct∕2) ∝ t 2∕3 + O(t 8∕3 ). The distribution function f depends on the "co-moving" variables x and p (co-moving with the Hubble expansion), which are linked through a canonical transformation to the physical position X and physical momentum P with (2) =̄(t) [1 + (x, t)], = ∫ f(x, p, t) d 3 p.

Fig. 1
Simulation result showing the logarithm of the over-density + 1 = ∕̄ of the dark-matter distribution on cosmic scales (thickness of slice is l ∼ 10 21 m); figure adapted from Stücker et al. (2018). Here, is the total density while ̄(t) is an isotropic dilution factor stemming from the volume expansion of the Universe 2 Topology of the dark-matter phase-space and shell-crossings

The dark-matter sheet and suitable parametrisation
Cosmological observations indicate that the initial dark-matter distribution is extremely cold (e.g. Aghanim et al. 2018;Kunz et al. 2016;Thomas et al. 2016;Gilman et al. 2020;Ilić et al. 2021). Furthermore, any residual weak thermalisation of the initial dark-matter distribution, if present, should have negligible impact on the gravitational dynamics on the considered macroscopic scales. Therefore, it is customary in cosmology to assume the limiting case of initially perfect coldness, which is 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 v ∶= p∕(ma 2 t a) . 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 setup 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.
(4) H(x, p, t) = |p| 2 2ma 2 + m (x, t), 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 q 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).
Indeed, introducing the Lagrangian map q ↦ x(q, ) from initial ( = 0 ) position q to the current Eulerian position x at a refined temporal variable ( ∝ t 2∕3 ; see below), the Vlasov-Poisson equations for dark matter take the form (see, e.g. Bernardeau et al. 2002;Matsubara 2008;Uhlemann et al. 2019) where, for convenience, we employ a rescaled Poisson potential = ∕(4 Ḡa 2 ) . Here and in the following, over-dots denote convective (i.e. total) derivatives with respect to the dimensionless temporal variable = a ; 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 Sect. 3 for details. The Lagrangian map is defined such that v =∶̇x(q, ) is the dark-matter velocity field (in units of lengths since = a (6) Fig. 2 Evolution of the dark-matter sheet in the phase-space, employing a simplified one-dimensional setup on the torus with x ∈ [− , ) in arbitrary units in a (fictitious) Universe filled only with dark matter. 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 (available from https:// bitbu cket. org/ ohahn/ cosmo_ sim_ 1d), initialised at = 0 with a sine-wave velocity profile and vanishing density contrast ) ( ∝ t 2∕3 where t is standard time). The second panel shows the critical instant when the dark-matter sheet folds for the first time (here at = 1 and x = 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 is dimensionless), where v relates to the momentum (non-canonically) according to p = vma 2 t a . We note that Eq. (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. Rampf et al. 2015;Matsubara 2015;. 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 Eq. (6) can be exactly expressed through a mass conservation law using the Dirac delta (3) D (Matsubara 2008;Taylor and Hamilton 1996;McDonald and Vlah 2018;Morrison et al. 2020), here for the case when → 0 , initially, which, roughly speaking, may be interpreted as the continuum limit N → ∞ of a density contrast N induced through a (hypothetical) collection of N discrete parti- , where x i ( ) denotes the current position of the ith dark-matter particle (see, e.g. Buchert and Dominguez 2005 for a more precise argument). Equation (7) is the result of linking the law of mass conservation ̄(1 + (x, )) d 3 x =̄d 3 q 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 q x = q ⊗ x T is the direct (dyadic) product, while q n denotes the nth root of the equation x − x(q, ) = 0 . Intuitively, the sum in (8) appears since n fluid streams occupy the same current position x , 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 q 1 = q , implying that mass conservation simplifies to (x(q, )) + 1 = 1∕J in Lagrangian coordinates, where J = det[ q x(q, )] 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∕ )̇x , 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 Eqs. (6) and (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 Sect. 6.
Finally, observe that the Vlasov-Poisson equations (6) maintain their regularity for → 0 , provided one imposes the following boundary conditions on the initial conditions (Brenier et al. 2003;Zheligovsky and Frisch 2014) at time = 0 (denoted by the superscript "ini"): where v ini =̇x(q, = 0) . Thus, with these boundary conditions, only the initial gravitational potential ini at = 0 needs to be prescribed; see, e.g. Michaux et al. (2020) 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 (Rampf et al. 2015). In addition, note that the second condition in (10) immediately implies that x × v ini = 0 which, due to Helmholtz's theorem for the conservation of vorticity flux, translates to x ×ẋ = 0 at all times along each matter stream (but note that effective vorticity is generated in Eulerian coordinates through multistreaming Uhlemann et al. 2019;Pueblas and Scoccimarro 2009;Cusin et al. 2017). In the following, we will mostly focus on solutions using the conditions (10), since The Lagrangian map q ↦ x(q, ) from initial position q to current position x (left panel), and the density contrast = ( −̄)∕̄ (right panel), employing the analytical results of  in the same simplified setup as for Fig. 2. Note specifically that a torus geometry is used with q, x ∈ [− , ) , 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 darkmatter evolution arises when ∇ q x vanishes (corresponding to det[ q x] = 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 (centre 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 those are the only one currently available that have mathematical proofs of analyticity until some finite time; see Sects. 4.2 and 4.3 for details.

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 q x vanishes, which by virtue of Eq. (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 "shell-crossing" (since the first collapse models involved the study of "matter shells" in 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 Zel'dovich (1970). 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. Buchert 1992;Valageas 2011;White 2014;McQuinn and White 2016). We review this model in Sect. 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 (Novikov 1969;Zentsova and Chernin 1980). 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 one-dimensional 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 one-dimensional phenomena with the formation of pancake-like structures (e.g. Melott and Shandarin 1993).
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 shellcrossing (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. Fillmore and Goldreich 1984), it was not predicted by theoretical means until recently, when S. Colombi was able to accurately estimate the involved gravitational forces in multi-streaming regions from first principles Colombi (2015) (see also Taruya and Colombi 2017;Pietroni 2018). 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 singlestream regions is only true in the one-dimensional case).
Starting at shell-crossing (here at = 1 ), some particle trajectories in multistreaming regions are bent back into the centre 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 , Colombi (2015), Taruya and Colombi (2017), 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 Taruya and Colombi 2017;Pietroni 2018 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. Schaye et al. (2015), Vogelsberger et al. (2014), Springel et al. (2018), Aricò et al. (2020), Lewandowski et al. (2015).
In the following sections, we will focus on the mathematical analysis of the dark-matter evolution before and (shortly) after the first shell-crossing-a regime that until very recently was not accessible by analytical means. Nonetheless, in Sects. 5.1 and 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 x , which in Lagrangian coordinates is not an independent variable. Indeed, earlier, we have introduced the Lagrangian map q ↦ x(q, ) which implies for the conversion of derivatives that where, from here on, Latin indices i, j, … = 1, 2, 3 denote the three spatial components, and (Einstein) summation over repeated indices is assumed.
To proceed, one introduces the Lagrangian displacement field (q, ), 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. Buehlmann and Hahn 2019). Introducing the operator T n ∶= (2 2 ∕3) 2 + − n 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 are which, once solved, respectively, for q ⋅ and q × , yield the formal solution for the displacement using a Helmholtz decomposition ( ∇ −2 q is the inverse Laplacian), We have defined two non-linear source terms (see below for physical interpretations): where "i" denotes a partial derivative with respect to Lagrangian component q i , while ijk is the fundamental anti-symmetric tensor and ij the Kronecker delta.
Here and in the following, displacements and current positions x depend on q 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 D and the replacement (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  (12), one first multiplies (6) by 2 2 ∕3 , leading to T 0 x l = − ∇ x l . Multiplying the last equation by the matrix element x l ∕ q j then yields x l,jT0 x l = − ∇ q j . Finally, taking the Lagrangian curl of this leads to ( q x l ) ×T 0 q x l = 0 (Matsubara 2015; Zheligovsky and Frisch 2014), which, once written for the displacement = x − q , leads directly to the second equation in (12).

Derivation of equations
Physical interpretation of source terms. M is only non-zero 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 (monokinetic). Indeed, for single-stream flows, the Lagrangian mapping q ↦ 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) becomes T 1 q ⋅ ≃ 0 at sufficiently early times in the single-stream regime. This is an ordinary differential equation in time for q ⋅ 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 depend 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 w ∶= x × v 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 Zheligovsky and Frisch (2014) for details, and e.g. Matsubara (2015), , Rampf (2012) 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 (q, ) = x(q, ) − q , one obtains where from here on ∶= q , and we remind the reader that T n = (2 2 ∕3) 2 + − n . 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. Bernardeau et al. 2002;Peebles 1980). If the boundary conditions (10) are employed, then only the growing solution ∼ is selected and one arrives at the classical Zel'dovich solution (Zel'dovich 1970;Buchert 1989) where ini is the initial gravitational potential. Of course, other solutions (e.g. with terms that decay in ) exist as well (Buchert 1992;Bouchet et al. 1995;Nadkarni-Ghosh and Chernoff 2013), but require a more sophisticated analysis of the boundary layer (see, e.g. Fidler et al. 2016;Adamek et al. 2017).
(2) Shell-crossing solutions (W ≠ 0, M = 0). As long as the flow is still single stream, the evolution equations for the particle displacements are exactly If we impose the boundary conditions (10), then these equations remain analytic for → 0 (Brenier et al. 2003;Zheligovsky and Frisch 2014), thereby suggesting that a series representation for the displacement in powers of (the linear growth time), around = 0 , is amenable: Here, (1) (q) = − 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 . Plugging (20a) into (19), one obtains explicit all-order recursion relations (Rampf et al. 2015;Matsubara 2015;Zheligovsky and Frisch 2014;Rampf 2012;Schmidt 2021) from follows, and subsequently = ∑ n (n) n . 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 (Rampf et al. 2015;Zheligovsky and Frisch 2014), as well as numerical evidence of convergence at shell-crossing and even shortly after ; see Sect. 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 q = q ⋆ for which the density blows up for the first time (cf. Eq. 8), i.e.
Of course, not only the particle with label q ⋆ but all particles can be evolved until ⋆ ; in the following we denote the family of displacements at ⋆ with ⋆ ∶= ⋆ (q, = ⋆ ) = ∑ n max n=1 (n) (q) n ⋆ , truncated at sufficiently large order n max .
(3) Post-shell-crossing solutions (W ≠ 0, M ≠ 0). Instantly after the fluid bifurcates at shell-crossing time = ⋆ to three streams at location q = q ⋆ , the multistream source term M is non-zero in the neighbourhood of q ⋆ (blue shaded region in Fig. 3). In this case, the evolution equation for ⋅ at ≥ ⋆ can be formally solved using the method of variation of constants, i.e. , Colombi (2015), Taruya and Colombi (2017), Pietroni (2018) for complementary methods applied to the one-dimensional case. Note that the re-occurring terms W and M in equations (22a), (22b) and (22c) are a priori unknowns for times > ⋆ , since they are functions of the unknown post-shell-crossing displacement (cf. equations 14 and 15). However, the flow directions of the matter stream are known until shell-crossing exactly-and at least approximately shortly after, essentially by an argument of momentum conservation. Thus, the unknowns W and M in equations (22a), (22b) and (22c) can be approximately determined by extrapolating the shell-crossing solutions to times shortly after shell-crossing. Specifically, for > ⋆ , we approximate (n) (q) n is the truncated shell-crossing displacement until n = n max with Taylor coefficients (n) , the latter determined through (20b) and (20c). We remark that the computation of the multi-streaming source M(x • ) is still non-trivial and requires normal-form considerations inherited from catastrophe theory (Berry 1977;Arnol'd 1980); see Sect. 5.3 for details.
As first observed by Colombi (2015), 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 Taruya and Colombi 2017;Pietroni 2018); so are applications beyond 1D (e.g. Chen and Pietroni 2020;Zimmermann et al. 2021). 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 Sects. 4.2 and 4.3, we investigate the mathematics of shell-crossings, respectively, for simplified and random initial conditions. (22b) ( ) = 3

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. Bernardeau et al. 2002 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 (a) 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 (b) 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 x ; 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) (Peebles 1980;Fry 1984;Goroff et al. 1986;Blas et al. 2014), and there exist explicit all-order recursion relations (Goroff et al. 1986) as well as a variety of SPT extensions (e.g. Pietroni 2008;Bernardeau et al. 2008;Anselmi et al. 2011;Crocce et al. 2012;Carlson et al. 2013;Kitaura and Hess 2013;Vlah et al. 2015); see also Bernardeau et al. (2002) for an extensive review. A straightforward application of SPT is to determine so-called loop corrections to n-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 (Suto and Sasaki 1991;Nishimichi et al. 2009;Carlson et al. 2009;Taruya et al. 2012;Simonović et al. 2018). Furthermore, certain effective approaches based on SPT are capable of incorporating shell-crossing effects to some extent; see Sect. 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. Matsubara (2011) We remark that, so far, the mathematical convergence of the SPT series until shell-crossing has not been explicitly demonstrated, except for the simplified cases of one-dimensional (McQuinn and White 2016) and spherical collapse (Rampf 2019). Nonetheless, convergence for cosmological initial conditions until-but excluding-shell-crossing is likely, mainly as a consequence of the analyticity results in Lagrangian coordinates Schmidt 2021). 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 n = 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 q ↦ x , 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 (Zel'dovich 1970) in the fundamental coordinate system spanned by q x , 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 (Doroshkevich 1970).
The Zel'dovich approximation has several key advantages over the corresponding linear solution in Eulerian coordinates, where the latter predicts that (x, ) = ∇ 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 the chosen spatial resolution; see Sect. 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 (Klypin and Shandarin 1983;Efstathiou et al. 1985;Scoccimarro 1998;Valageas 2002); see also Michaux et al. (2020), ,  for higher order refinements within the framework of Lagrangian perturbation theory (Buchert 1992;Buchert 1989;Bouchet et al. 1995 al. 1991;Bouchet et al. 1992;Buchert and Ehlers 1993;Buchert 1994;Rampf and Buchert 2012) (LPT), which are predictions based on perturbative truncations of equation (20a).
-Predicting the change of photon momenta emitted from observed luminous tracers of dark matter, due to peculiar motion and the expansion of space ("redshiftspace distortions") (Matsubara 2008;Taylor and Hamilton 1996;Valogiannis et al. 2020;Jackson 1972;Kaiser 1987;Fisher and Nusser 1996).

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. Vogelsberger et al. 2009;Abel et al. 2012;Shandarin et al. 2012;Hahn and Angulo 2016;Sousbie and Colombi 2016;Stucker et al. 2020).
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 Eq. 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 q 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 (equation 20); however, then no 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 = x − q (Rampf and Frisch 2017): where F 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 F = − sin q 1 , while at order 1 , one finds refined recursion relations for the Taylor coefficients appearing in (q, ) = ∑ ∞ n=1 (n) (q) n (see equation (36) in Rampf and Frisch 2017). Most importantly, the refined recursion relation for (n) 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-onedimensional 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 Saga et al. 2018), while the red-dashed line is obtained using (20a) 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 Saga et al. 2018 for the specific setup). 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) (Peebles 1967). 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 (Bernardeau 1992; Munshi et al. 1994;Tatekawa 2005).
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 shellcrossing (Rampf 2019). In that case, the initial gravitational potential can be taken to be where K is a positive [negative] constant amplitude of a spherical over-density [under-density/void], and ini is an arbitrary function that reflects the departure from spherical symmetry. The mathematical form of ini suggests the solution Ansatz where S is a purely time-dependent function that models the purely spherical collapse, while reflects a perturbative departure therefrom. In the present case, both S and can be represented individually by convergent time-Taylor series (Rampf 2019). 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 O( ) 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 Saga et al. 2018 for the highly related triaxial symmetric case, or Nadkarni-Ghosh and Chernoff (2011) for a complementary analysis with similar conclusions). Nonetheless, we should stress that such singularities are basically man-made: 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 (Doroshkevich 1970).

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 (20a) on a three-dimensional numerical grid with N collocation points. However, there are also ways to prove, by entirely theoretical means, that the solutions (20a) are time-analytic at least until some finite time, which goes back to a seminal paper by Zheligovsky and Frisch (2014). 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.
Time-analyticity of dark-matter displacements. Instead of proving the convergence of = ∑ ∞ n=1 (n) n , it is more instructive to first prove the convergence of the tensorial gradient therefrom, i.e. = ∑ ∞ n=1 (n) n , with Taylor coefficients where L (n) and T (n) are given by the recursion relations (20b) and (20c), respectively. Before proceeding let us briefly introduce the 1 norm denoted with ‖ ‖ ∶= ∑ n � k � < ∞ , for any periodic tensor, vector or scalar function (q) = ∑ k k exp(ik ⋅ q) , where k are its Fourier coefficients. The corresponding function space has the algebra property ‖ ‖ ≤ ‖ ‖ ‖ ‖ , and furthermore bounds the operator ∇ −2 to unity from below. Using these properties as well as the explicit forms of L (n) and T (n) (equations 20b and 20c), the 1 norm of (28) turns into the polynomial inequality (Zheligovsky and Frisch 2014) In deriving this inequality, one considers the large-n limit (relevant for investigating convergence), where the coefficients appearing in the various terms in L (n) and T (n) are bounded (rational) functions. Now, introducing the generating function ( ) ∶= ∑ ∞ n=1 ‖ (n) ‖ n , as well as multiplying (29) with n and summing over n from one to infinity, one gets ≤ ‖ ini ‖ + 12 2 + 6 3 , which is equivalently (see Rampf et al. 2015 for a graphical representation as well as for the ΛCDM case). Initially for = 0 , this polynomial p behaves as − for small and asymptotically as 3 ; thus, p( = 0) has three roots with one at phys = 0 . For small positive times, p gets shifted upwards, thereby drifting the root phys to the right. This branch from = 0 until phys marks the physical regime where p > 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 p( ) to zero, leading to Zheligovsky and Frisch (2014) 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 (Zheligovsky and Frisch 2014; Michaux et al. 2020), especially when assumptions about the regularity of ini are imposed. Shell-crossing and radius of convergence. To determine the actual radius of convergence of (q, ) = ∑ ∞ s=1 (s) (q) s 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 n of the displacement series. That is one searches for the location where the perturbatively truncated Jacobian vanishes for the first time (cf. Eq. 21) . For each truncation order n, one then obtains perturbative estimates of the shell-crossing time and location, denoted, respectively, with (n) ⋆ and q (n) ⋆ . By monitoring the trend of these estimates at successively high orders n, one then obtains evidence for or against (weak) convergence.
In , the recursive displacement solutions (20a), (20b) and (20c) as well as J {n} have been realised on a spatial grid with N collocation points for cosmological random initial conditions. The corresponding code 1 is parallelised (MPI+threads), and de-aliases cubic functions in the displacement. Most analysis in  was performed with a grid resolution of N = 256 3 , for which convergence of the shell-crossing time ⋆ accurate to three significant digits was typically reached between orders n = 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 ΔJ {n} ∶= J {n} − J {n−1} from all Lagrangian grid points, evaluated at the time of first shell-crossing (here at ⋆ = 0.0813 for a standard ΛCDM cosmology with N = 256 3 ). Evidently, ΔJ {n} peaks sharper at the origin for higher orders, indicating that truncation errors decrease rapidly at increasing order n, 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 R is the radius of convergence-when that limit exists. A standard way to estimate R is the Domb-Sykes plot (Domb and Sykes 1957), where one draws ‖ (n) ‖∕‖ (n−1) ‖ versus 1/n, from which 1∕ R follows by extrapolation to the y-intercept (see also Michaux et al. 2020;Podvigina et al. 2016 for related applications). In the right panel of Fig. 6 we show the Domb-Sykes plot for the spatial coefficients ‖ (n) ‖ evaluated at the Lagrangian locations of first shell-crossing for three different resolutions, N = 512 3 , 256 3 , 128 3 (respectively, in blue, orange and green; top to bottom), while the length of the grid is fixed at L = 125Mpch −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 ( n ≳ 7 ) the ratios settle into a linear behaviour, thereby justifying a linear extrapolation to the y-intercept (dotted lines), which, using (34), leads to the radius of convergence R . In all cases considered the (temporal) value of R is roughly 30 − 40 % larger than the shell-crossing time (see table 1 in  for specific values), thus indicating that the mathematical validity of (20a) surpasses vastly its physical validity (since equation (20a) 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 Sect. 4.2). Obviously, R 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. Bernardeau et al. (2002), Michaux et al. (2020), ; see, however, Schmidt (2021) 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  where is a singularity exponent. In this case, the ratio of coefficients satisfies the linear relationship (cf. Podvigina et al. 2016;van Dyke 1974) 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 N = 128 3 , 256 3 , 512 3 . Since these exponents are positive non-integers and smaller than unity, it follows that the first time derivative of ‖ (q, )‖ blows up at = R .
Concluding this section, the recursive solutions (20) for = ∑ ∞ n=1 (n) n appear to deliver physically and mathematically meaningful results until the time of shell-crossing ⋆ , while the various bounds and limits are related through 0 < c < ⋆ < R , where c is an entirely analytical estimate of the radius of convergence, while R 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 N-body simulations, where (very massive) macroparticles are employed to coarse sample the dark-matter distribution in phase-space (see, e.g. Bertschinger and Gelb 1991; Dehnen and Read 2011 for reviews). However, these simulations can be prone to errors introduced through the discrete representation of the phase-space distribution (e.g. Abel et al. 2012;Melott et al. 1997;Angulo et al. 2013;Colombi and Touma 2014). 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 (Sect. 5.3), in the following two sections, we attempt to summarise important findings related to semi-analytical/effective avenues (Sect. 5.1) as well as to novel simulation techniques (Sect. 5.2).

Semi-analytical/effective methods
Simple (time-)Taylor expansions, such as (20a) in Lagrangian coordinates, cease to be valid at the time of first shell-crossing (Pietroni 2018;. Of course, the same can be said when seeking theoretical solutions to Vlasov-Poisson in Eulerian coordinates x , where the density and velocity can be Taylor expanded until shell-crossing with equivalent methods as reviewed in Sect. 3.2 (dubbed SPT; see also, e.g. Bernardeau et al. 2002;Buchert and Dominguez 1998). However, there also exist methods that seek to incorporate multi-streaming effects in a rather effective way (Buchert and Dominguez 2005;Pueblas and Scoccimarro 2009), notably the so-called coarse-grained perturbation theory Pietroni et al. (2012) as well as the effective theory of large-scale structure (Baumann et al. 2012;Carrasco et al. 2012;Porto et al. 2014). 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 Colas et al. 2020).
At the technical level, this can be achieved by introducing a spatially smoothed where f is the distribution function that solves equation (1), while W Λ is a Window function which is usually taken to be a normalised Gaussian, i.e. W Λ ∝ exp(−Λ 2 |x| 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 (Baumann et al. 2012;Carrasco et al. 2012) where , u and are, respectively, the zeroth-, first-and second-kinetic moments of the distribution function (see, e.g. Bernardeau et al. 2002), while [ (x)] Λ ∶= ∫ W Λ (x − x � ) (x � ) d 3 x � 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. Bartelmann et al. 2016;. The central building block is the generating functional involving the N-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 shot-noise 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 cannot be made at this stage.
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 N-body techniques.
Phase-space reconstruction. With the standard N-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 N-body technique.
One of these new approaches is dubbed GDE which is short for geodesic deviation equation (Vogelsberger et al. 2009;Vogelsberger and White 2011). 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 (Wang and White 2007;Melott 2007).
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 Hahn and Angulo (2016) (see also Abel et al. 2012;Shandarin et al. 2012;. 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 (Stucker et al. 2020).
Pure phase-space tessellation methods. Instead using tracer particles, the numerical code Coldice employs simplices (tetrahedra) to tessellate the continuous darkmatter sheet (Sousbie and Colombi 2016). Fairly similarly as in the above sheet method, the forces are obtained from continuous density estimates and subsequent interpolations to an undistorted grid. Furthermore, Coldice allows for anisotropic mesh refinements, thereby being able to resolve the phase-space to high precision. Unfortunately, the complexity of the dark-matter 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 (Colombi 2021) for explicit comparisons of Coldice against standard N-body codes during the early violent relaxation phase of halo formation (see also Colombi and Touma 2014).
Vlasov-Poisson simulations in 6D. There exist also full Vlasov-Poisson simulations that are not constrained in evolving the dark-matter sheet (Yoshikawa et al. 2013;Tanaka et al. 2017). 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) (Yoshikawa et al. 2020;Colombi and Alard 2017).
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 (Widrow and Kaiser 1993). Indeed in Zhang et al. (2002) 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. Mocz et al. (2017), Kopp et al. (2017), Garny et al. (2020), Eberhardt et al. (2020).

Mathematical post-shell-crossing analysis
Turning back to the theoretical solution scheme for the Lagrangian master equation (22a), the leading-order behaviour of the involved source terms shortly after shell-crossing time ⋆ is given by W(x(q, )) ≃ W(x • (q, )) and M(x(q, )) ≃ M(x • (q, )) , which is, strictly speaking, exact at = ⋆ , and approximately thereafter, since x • (q, ) = q + ∑ n max n=1 (n) (q) n 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 q ′ 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 x • (q, ) around the location q ⋆ of first shell-crossing, combined with exploiting the expected topology of x • (q, ) 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 q ⋆ = 0 ), the map has an inflection point with local behaviour ∼ q 3 around the shell-crossing location; shortly later (here shown at = 1.3 in blue) the local behaviour of the map in 1D is , Colombi (2015), Taruya and Colombi (2017) where c 1,3 are positive Taylor coefficients. Here, N (q) corresponds to the said normal form in the one-dimensional case, and note specifically that a term ∼ (q − q ⋆ ) 2 is absent (it can be removed by a suitable Galilean transformation). Employing (40) to determine the roots of x • (q, ) − x • (q � , ) = 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 Colombi 2015;Taruya and Colombi 2017). 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(x • (q, )) is determined, the post-shell-crossing displacement is readily obtained using equations (); see, e.g. equation (18) in  for explicit solutions in the 1D case.
These purely analytical solutions are shown in Fig. 7 and compared against numerical simulations at high resolution (available from https:// bitbu cket. org/ ohahn/ cosmo_ sim_ 1d). Panel (a) shows the phase-space in one dimension at time = 1 = 1.001 featuring single-and multi-streaming regions that are spatially separated by infinite densities (panel b). At the same time, panel (c) shows the corresponding acceleration of particles, displaying four non-differentiable sharp features : in single-stream regions the acceleration is exactly zero (very left and very right regions in panel c), but, at the depicted time, then jumps to non-zero values with local behaviour ∼ ( − 1.001) 5∕2 and ∼ (q ± 0.17) 5∕2 at time 1 , thereby indicating that the third derivatives in space and time of the displacement blow up. Furthermore, the third space derivative of the displacement flips sign (here around q = q 2 = ±0.08 ) thereby marking a singularity of ∼ (q − q 2 ) 3 Θ(q 2 − q) and similarly in the temporal dependence. Finally, a non-trivial singularity appears when the initial velocity is not exactly point-symmetric, here for a simplified model with v ini = − sin q + 0.1 sin 4 q − 0.12 sin 6 q : the particle that is initially at q = 0 will remain there until shell-crossing ( = 1 ), at which time an asymmetric multistreaming force arises, essentially as a consequence of momentum conservation, that leads the q = 0 particle to begin moving with (q = 0) ∼ 3 ) (panel d in Fig. 7; see also Pietroni (2018) for similar conclusions in a slightly different setup). 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. . The panels (a)-(c) are evaluated at = 1.001 (shell-crossing is at = 1 ), and show, respectively, the dark-matter phase-space, the projected density contrast, and the particle acceleration. Panel (d) shows the temporal evolution of the particle at q = 0 (i.e. the shell-crossing location), which suddenly begins moving at shell-crossing time = 1 due to an asymmetry in the initial data. Panels (c) and (d) display non-analytic behaviour where derivatives in the phase-space blow up. Figures obtained using the methods of  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 N electrons with charge e and mass m that interact only through electrostatic forces. To warrant charge neutrality, these electrons are immersed in a rigid uniform background of opposite charge, which in the present case is −̄= −Ne∕V where V 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. Baus and Hansen 1980;Koester and Chanmugam 1990). Adopting somewhat the above notation, the equation of motion for the jth electron, formulated in physical coordinates in 3D and employing the standard time t, reads (see, e.g. Baus and Hansen 1980;Escande et al. 2018) where x = E 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 employing the Lagrangian map q ↦ x(q, t) from initial position q to the current Eulerian position x . The plasma velocity can be represented as v(x(q, t)) =̇x(q, t) , where the over-dot denotes now the convective time derivative with respect to standard time, i.e. d∕dt = ∕ t | q = ∕ t| x + (dx∕dt) ⋅ x . Then, equations (41) may be written as a continuous Vlasov-Poisson description 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 = a (since in the cosmic case, time-analyticity holds in ∼ t 2∕3 but not in t). In addition, equations (6) feature a Hubble drag term ∼̇x 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 dt = a 2 d (Buchert 1989;(41) (3) Shandarin 1980). Furthermore, equations (6) have a minus sign in front of x 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 x → x + n(t) , where n(t) 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. Diamond et al. 2010;Escande et al. 2016), however, we are unaware of explicit use of Lagrangian coordinates (though they are implicitly assumed for Eq, 41 in the discretised sense). In this context note that for a single electron particle at current position x(t) exposed to its own electric field, the Poisson potential is in Fourier space for k ≠ 0 (cf. equation 1 in Escande et al. 2016). By contrast, for a set of electrons parametrised in Lagrangian coordinates q , the Fourier transform of the Poisson potential is easily obtained from (42) 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 x and time t. See, e.g. Dawson (1960), Anderson et al. (2001) 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 (Escande et al. 2018;O'Neil et al. 1971;Lesur and Diamond 2013). 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. Carlevaro et al. 2014Carlevaro et al. , 2020. Similarly as in the cosmological case, also here the Vlasov-Poisson equations are directly related to a Hamiltonian principle (Tennyson et al. 1994;Antoniazzi et al. 2006), and thus can be reduced to Newtonian-type equation of motion coupled to a Poisson equation. This set of equations has been provided by O 'Neil et al. (1971) in the case of N discretised electrons in 1D, which, adopting the same notation as above, can be written as where x j denotes the current position of the jth 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 (Carlevaro et al. 2020), thereby justifying linear temporal expansions applied to the dielectric. This allows to incorporate the effect of the plasma density in the Poisson equation by a linear dielectric, and at the same time, cast the Poisson equation into a simple evolution equation (O'Neil and Malmberg 1968). By contrast, the density of the electric beam is just a superposition of the N charged particles, (see, e.g. Tennyson et al. 1994). 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 time-dependent dielectric function (see, e.g. equation 14 in O' Neil et al. 1971). 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 (O'Neil et al. 1971). 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 q ↦ x (q, t) , where =p,b (see also Firpo and Elskens 1998;Antoni et al. 1998 for alternative paths to a continuum description). The beam  Neil et al. (1971); note that the filamentary tails of the distribution have been cut off for reasons of better visibility and plasma velocities are then defined with v (q, t) =̇x (q, t) , and the Vlasov-Poisson equations are for the components =p,b, where we have defined (∇ x ) = ∇ x (x (q, t)) , and mass conservation implies b (x (q, t)) = ∫ D (x (q, t) − x b (q � , t))dq � , 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 , Chen et al. (2019b) (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. Escande et al. 2018), as well as to the case of multiple cold beams, by, e.g. adopting the formalism laid out in Carlevaro et al. (2015). Specifically, for M cold beams, one essentially needs M equations of motion (44) in the continuum limit coupled to a Poisson equation, where the latter then involves the superposition of M charge densities (45).
The weak warm beam instability. Finally, another related plasma problem concerns the wave-particle description of Langmuir waves (Escande et al. 2018;Carlevaro et al. 2014;Besse et al. 2011), 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 so-called 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 (Escande et al. 2018;Antoni et al. 1998;Elskens and Pardoux 2010), which is obtained by setting for the trajectories x j = q j + v j t + x ini j +̇x ini j t in our notation, where x ini j and ̇x ini j are the initial mismatches in position and velocities within the multi-beam-multiarray (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 higher order 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 Sect. 3.2 and specifically Sect. 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 (Sect. 4.1), and/or using effective methods to Vlasov-Poisson (Sect. 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 (Sect. 5.2), that uncover much more details of the dark-matter phase-space as compared to standard particle-incell/N-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 N-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 (Sect. 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 setup 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. In addition, in the long term, 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.