Semiclassical path to cosmic large-scale structure

TL;DR

This paper introduces a semiclassical propagator approach to model nonlinear gravitational dynamics of cold dark matter, regularizing shell-crossing and avoiding spurious vorticity, with potential for improved numerical simulations.

Contribution

It develops a semiclassical framework using a Schrödinger equation for large-scale structure, connecting to Lagrangian perturbation theory and preserving Hamiltonian structure.

Findings

Propagator solutions resemble LPT at small ar.

Spurious vorticity is avoided in the semiclassical approach.

The method regularizes dynamics at shell-crossing.

Abstract

We chart a path toward solving for the nonlinear gravitational dynamics of cold dark matter by relying on a semiclassical description using the propagator. The evolution of the propagator is given by a Schr\"odinger equation, where the small parameter $\hbar$ acts as a softening scale that regulates singularities at shell-crossing. The leading-order propagator, called free propagator, is the semiclassical equivalent of the Zel'dovich approximation (ZA), that describes inertial particle motion along straight trajectories. At next-to-leading order, we solve for the propagator perturbatively and obtain, in the classical limit the displacement field from second-order Lagrangian perturbation theory (LPT). The associated velocity naturally includes an additional term that would be considered as third order in LPT. We show that this term is actually needed to preserve the underlying Hamiltonian structure, and ignoring it could lead to the spurious excitation of vorticity in certain implementations of second-order LPT. We show that for sufficiently small $\hbar$ the corresponding propagator solutions closely resemble LPT, with the additions that spurious vorticity is avoided and the dynamics at shell-crossing is regularised. Our analytical results possess a symplectic structure that allows us to advance numerical schemes for the large-scale structure. For times shortly after shell-crossing, we explore the generation of vorticity, which in our method does not involve any explicit multi-stream averaging, but instead arises naturally as a conserved topological charge.

Figures

Equations

$$S_0(\mathbf{x},\mathbf{q};a)=\frac{1}{2}(\mathbf{x}-\mathbf{q})\cdot \frac{\mathbf{x}-\mathbf{q}}{a},$$

$$K_0(\mathbf{x},\mathbf{q};a)=N\exp \left\{\frac{\mathrm{i}}{\hslash }{S}_0(\mathbf{x},\mathbf{q};a)\right\},$$

$${\psi }_0(\mathbf{x};a)=\int {\mathrm{d}}^{3}q{K}_0(\mathbf{x},\mathbf{q};a){\psi }_0^{(\mathrm{ini})}(\mathbf{q}),$$

$$\mathrm{i}\hslash {\partial }_a{K}_0$$

$$\mathrm{i}\hslash {\partial }_a{\psi }_0$$

$$\mathrm{i}\hslash {\partial }_{a}K$$

$$\hat{H}\equiv -\frac{{\hslash }^2}{2}{\nabla }_x^2+V_{\mathrm{eff}}(\mathbf{x};a),$$

$${\partial }_a{\varphi }_{\mathrm{v}}-\frac{1}{2}|\nabla {\varphi }_{\mathrm{v}}{|}^2$$

$${\partial }_a\delta -\nabla \cdot [(1+\delta )\nabla {\varphi }_{\mathrm{v}}]$$

$${\nabla }^2{\phi }_{\mathrm{g}}$$

$$V_{\mathrm{eff}}\equiv \frac{3}{2a}\left({\phi }_{\mathrm{g}}-{\varphi }_{\mathrm{v}}\right),$$

$${\delta }^{(\mathrm{ini})}=0,{\phi }_{\mathrm{g}}^{(\mathrm{ini})}={\varphi }_{\mathrm{v}}^{(\mathrm{ini})}.$$

$$\delta (\mathbf{x};a)$$

$$V_{\mathrm{eff}}(\mathbf{x};a)$$

$$V_{\mathrm{eff}}^{(2)}=\frac{3}{7}{\nabla }^{-2}\left[{\left({\nabla }^2{\phi }_{\mathrm{g}}^{(\mathrm{ini})}\right)}^2-{\left({\nabla }_i{\nabla }_j{\phi }_{\mathrm{g}}^{(\mathrm{ini})}\right)}^2\right],$$

$$\mathbf{q}\mapsto \mathbf{x}(\mathbf{q};a)=\mathbf{q}+\mathbf{\xi }(\mathbf{q};a)$$

$$\mathbf{v}(\mathbf{x}(\mathbf{q};a);a)={\partial }_a^{\mathrm{L}}\mathbf{x}(\mathbf{q};a)\equiv \dot{\mathbf{x}}(\mathbf{q};a).$$

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

$$\ddot{\mathbf{x}}(\mathbf{q};a)+\frac{3}{2a}\dot{\mathbf{x}}(\mathbf{q};a)=-\frac{3}{2a}{\nabla }_x{\phi }_{\mathrm{g}}(\mathbf{x}(\mathbf{q};a);a),$$

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

$${\xi }_i(\mathbf{q};a)=\sum _{n=1}^{\infty }{\xi }_i^{(n)}(\mathbf{q})a^n,$$

$${\xi }_i^{(1)}$$

$${\xi }_i^{(2)}$$

$${\mu }_2(\mathbf{q})=\frac{1}{2}\left[{\left({\nabla }_{\mathrm{L}}^2{\phi }_{\mathrm{g}}^{(\mathrm{ini})}\right)}^2-{\left({\nabla }_i^{\mathrm{L}}{\nabla }_j^{\mathrm{L}}{\phi }_{\mathrm{g}}^{(\mathrm{ini})}\right)}^2\right].$$

$${\mathcal{C}}_i^{[n]}:={\epsilon }_{ijk}\left({\delta }_{lj}+{\xi }_{l,j}^{[n]}\right){\dot{\xi }}_{l,k}^{[n]},$$

$${\mathcal{C}}_i^{[2]}$$

$$=\frac{a^2}{2}{\epsilon }_{ijk}{\phi }_{\mathrm{g},lj}{V}_{\mathrm{eff},lk}^{(2)}+\mathcal{O}(a^3),$$

$$\Gamma \equiv {\oint }_{C(a)}\mathbf{v}\cdot \mathrm{d}\mathbf{x}={\int }_{S^{(\mathrm{ini})}}({\nabla }^{\mathrm{L}}\times {\mathbf{v}}^{(\mathrm{ini})})\cdot {\mathrm{d}\mathbf{S}}^{(\mathrm{ini})}$$

$$\psi (\mathbf{x};a)=\int {\mathrm{d}}^{3}qK(\mathbf{x},\mathbf{q};a){\psi }^{(\mathrm{ini})}(\mathbf{q}),$$

$$f_{\mathrm{W}}(\mathbf{x},\mathbf{p})=\int \frac{{\mathrm{d}}^3{x}^{\prime }}{(2\pi {)}^3}\exp \left[\frac{-\mathrm{i}\mathbf{p}\cdot {\mathbf{x}}^{\prime }}{a^{3/2}}\right]\psi (\mathbf{x}+\frac{\hslash }{2}{\mathbf{x}}^{\prime })\bar{\psi }(\mathbf{x}-\frac{\hslash }{2}{\mathbf{x}}^{\prime }),$$

$$\varrho (\mathbf{x})$$

$$\mathbf{j}(\mathbf{x})$$

$$\underset{\hslash \to 0}{\lim }{f}_{\mathrm{W}}(\mathbf{x},\mathbf{p})=\varrho (\mathbf{x}){\delta }_{\mathrm{D}}^{(3)}\left(\frac{\mathbf{p}}{a^{3/2}}-\mathbf{v}(\mathbf{x})\right):=f_{\mathrm{fl}}(\mathbf{x},\mathbf{p}),$$

$$f_{\mathrm{fl}}(\mathbf{x},\mathbf{p})=\int {\mathrm{d}}^{3}q{\delta }_{\mathrm{D}}^{(3)}\left[\mathbf{x}-\mathbf{q}-\mathbf{\xi }(\mathbf{q})\right]{\delta }_{\mathrm{D}}^{(3)}\left[\frac{\mathbf{p}}{a^{3/2}}-{\mathbf{v}}^{\mathrm{L}}(\mathbf{q})\right],$$

Full text

I Introduction

Analytical methods for the evolution of cold dark matter (CDM), relevant for investigating the large-scale structure of our Universe, are valuable for gaining theoretical understanding and efficiently parametrising the cosmology-dependence of observables. Furthermore, analytical methods can also assist in improving numerical computations, for example by setting up initial conditions for cosmological $N$-body simulations Hahn and Abel (2011); Crocce et al. (2006), or as input for multi-time-stepping algorithms that are used to generate fast mock catalogues Tassev et al. (2013); Feng et al. (2016).

At sufficiently early times, the distribution for standard cold dark matter is characterised by the (perfect) fluid variables, a density and a single-valued velocity field. These fields are the dynamical quantities in the Eulerian description. In the corresponding perturbative framework, Eulerian perturbation theory Bernardeau et al. (2002), one assumes a perturbative smallness of the density and velocity from their background values. As a consequence, Eulerian perturbation theory struggles to capture large densities that naturally arise from gravitational dynamics. Moreover, due to the collisionless nature of dark matter, the infall induced by the gravitational potential leads to the crossing of fluid trajectories, which is often called ‘shell-crossing’. The instance of shell-crossing leads to the formation of caustics with a (formally) infinite density. Such non-analytic behaviour prevents a smooth continuation into the regime beyond shell-crossing, in which regions of multiple fluid streams with distinct velocities arise.

Some problems associated with shell-crossing are alleviated in a phase-space description, where the dark matter evolution corresponds to describing the embedding of an initially flat, 3-dimensional thin sheet into 6-dimensional phase-space, the ‘Lagrangian submanifold’. In this case, the singular behaviour of the density at shell-crossing does not pose a problem per se, because it is merely induced by projecting the folded phase-space sheet into 3-dimensional position space. A convenient way of describing the embedding of the sheet in phase-space is to use Lagrangian coordinates Villone and Rampf (2017); Zeldovich (1970); Buchert (1989); Bouchet et al. (1992); Buchert and Ehlers (1993); Rampf and Buchert (2012). The central quantity is the Lagrangian displacement field that encodes how fluid elements are displaced as a function of time and initial (Lagrangian) position. The corresponding perturbative framework is usually called Lagrangian perturbation theory (LPT), and, although a challenge, allows to investigate the instance of shell-crossing by analytic means Zheligovsky and Frisch (2014); Rampf et al. (2015); Taruya and Colombi (2017); Rampf and Frisch (2017); Rampf (2019); Saga et al. (2018). However, to update the gravitational potential that is responsible for displacing the fluid elements, one still requires, effectively, the Eulerian density. Since an irregular density carries over to irregularities in the tidal or force field Pietroni (2018), a theory that regularises the dynamics at shell-crossing and provides a smooth transition into the multi-stream regime is highly desirable.

In this paper, we provide a novel method to evolve CDM by relying on semiclassical dynamics using the propagator – which we motivate in more detail in Sec. II. Our primary motivation is to formulate an analytical framework that allows to go beyond shell-crossing while admitting perturbative solutions for the mildly nonlinear dynamics. To realise this goal, we adopt a formalism that respects key dynamical invariants, which additionally offers interesting perspectives for advancing numerical schemes for gravitational dynamics.

In essence, we introduce a suitable perturbative framework for the propagator and show that, at the leading order, our approach amounts to solving a ‘free-particle Schrödinger equation’, with a solution that is closely related to the classical ZA dynamics. At the leading order, our approach constitutes a heuristic derivation for a free-particle Schrödinger equation, which has been introduced in Refs. Coles and Spencer (2003); Short and Coles (2006a). At higher orders, our approach includes gravitational effects beyond the ZA, which are encoded in an effective potential. The resulting evolution equation for the propagator is then a Schrödinger equation that includes the said effective potential. Our formalism naturally respects the underlying Hamiltonian structure and allows us to determine the regularised fluid variables in the classical limit. Additionally, our formalism does not suffer from singularities at shell-crossing, thereby allowing it to enter the multi-stream regime smoothly.

The outline of the paper is as follows. In Sec. II, we begin by motivating our semiclassical description, which at the leading order returns a free propagator that reduces to the ZA in the classical limit. We then generalise the related Schrödinger equation for the propagator to include an effective potential. In Sec. III, we relate this effective potential to the cosmological fluid equations and describe standard perturbative solutions. In Sec. IV, we briefly review LPT since it will be helpful to analyse the results of the propagator method. Specifically, we utilise the so-called Cauchy invariants, that in our case encode the conservation of vanishing vorticity, as a diagnostic tool to investigate the Hamiltonian structure of the propagator method. In Sec. V, we show how the classical limits for the propagator relate to the classical fluid variables. As a demonstrative example, we perform those limits on the free propagator, from which we reproduce the ZA. Furthermore, we discuss how vorticity is generated by shell-crossing in the semiclassical picture. In Sec. VI, we go to the next-to-leading order in the propagator method and show that we recover exactly the second-order Lagrangian displacement. We show that the associated velocity receives a higher-order correction term compared to LPT, restoring compatibility with the Cauchy invariants thanks to the underlying symplectic structure of our method. In Sec. VII, we present concrete examples comparing numerical implementations of our propagator method to LPT, and show qualitative results for the density and vorticity after shell-crossing. We conclude in Sec. VIII, where we also provide an outlook on potential applications of our method.

II Motivating a propagator for gravitational interactions

Despite its simplicity, the classical Zel’dovich approximation (ZA) Zeldovich (1970) reproduces the cosmic web remarkably well. The ZA is based on a ballistic motion of fluid particles with prescribed velocity. While the ZA is the exact solution of the underlying fluid equations before shell-crossing in one spatial dimension, it receives corrections due to tidal gravitational fields in two and three dimensions. Furthermore, after the intersection of particle trajectories at shell-crossing, the ZA is doomed to fail, partially because it does not take the time evolution of the gravitational potential into account, which leads to nonzero particle acceleration and secondary infall.

In the following, we develop another perspective for the ZA, and outline how this model can be translated into the language of transition probabilities – with the aim to investigate novel avenues to surpass some of the shortcomings of the ZA. We shall make use of the dimensionless temporal coordinate $a=a(t)$, which is the cosmic scale factor with evolution governed by the usual Friedmann equations. For simplicity, we shall limit the analysis to a CDM dominated universe, the so-called Einstein–de Sitter (EdS) model. A generalisation to the $\Lambda$CDM model, where $\Lambda$ is the cosmological constant, is however straightforward.

II.1 Free propagator from Zel’dovich approximation

In a Cartesian coordinate system which comoves with the particles, the ZA amounts to particles moving along straight lines with constant velocity. For fluid particles that are initially ($a=0$) at position $\boldsymbol{q}$, and at time $a$, at position $\boldsymbol{x}$, the classical action is

[TABLE]

which, essentially, is the product of particle displacements and their velocities. Following standard methods inspired by Feynman Feynman and Hibbs (1965) the classical action yields a transition amplitude $K_{0}$ for a particle from being at time [math] at position $\boldsymbol{q}$ to being at time $a$ at position $\boldsymbol{x}$,

[TABLE]

where $N=(2\pi{\rm{i}}\hbar a)^{-3/2}$ is a normalisation coefficient chosen such that the transition amplitude returns initially the Dirac delta, K_{0}(\boldsymbol{x},\boldsymbol{q};\text{{\smalla!=!0}})=\delta_{\rm D}^{(3)}(\boldsymbol{x}-\boldsymbol{q}).

Now, as is well known in the context of quantum mechanics, the transition amplitude $K_{0}$ (and any other transition amplitude too) propagates a wave function $\psi_{0}$ from some initial to some final state,

[TABLE]

where \psi_{0}^{\rm(ini)}(\boldsymbol{q})\equiv\psi_{0}(\boldsymbol{q};\text{\smalla=0}). It is then elementary to show that, the evolution of the transition amplitude $K_{0}$ is governed by a potential-free Schrödinger equation,

[TABLE]

which is a free Schrödinger equation for the ZA wave function. This provides a heuristic derivation for the ‘free particle approximation’, which was introduced as an ad-hoc approximation of gravitational dynamics in Coles and Spencer (2003)111It appears that the work of Coles and Spencer (2003) was inspired by the Schrödinger-Poisson description of Widrow and Kaiser (1993), however, as we explain in App. B, the two approaches are formally not equivalent., and compared to the linearised fluid dynamics in Short and Coles (2006a).

Even at the leading order, our approach goes beyond the so-called adhesion model Gurbatov et al. (1989); Weinberg and Gunn (1990), which avoids shell-crossing through the introduction of an artificial viscosity, and has also been formulated in terms of propagators Jones (1999); Matarrese and Mohayaee (2002); Chavanis (2011); Rigopoulos (2015). In contrast to the adhesion approximation, our semiclassical propagator method possesses a conserved current that arises from the unitary evolution of the wave function and allows us to propagate through shell-crossing into the multi-stream regime.

The above equation is easily solved either by numerical or analytical means (see the following sections for details). As a demonstrative example, Fig. 1 shows the evolution of the wave function $\psi_{0}$ with the 1D initial data $\psi_{0}^{\rm(ini)}(q)\!=\!\exp\left(-{\rm{i}}\phi_{\rm v}^{\rm(ini)}/\hbar\right)$, with $\phi_{\rm v}^{\rm(ini)}\!=-\cos q$. The graph uses the domain coloring technique to assign a unique color to each point of the complex plane (cf. Wegert (2012)), so that the amplitude of $\psi_{0}$ is mapped to the brightness and the phase to the color hue for each point in the space-time plane $(x,a)$.

From Fig. 1 it becomes evident that, from the time of shell-crossing $a_{\rm sc}$ onwards, interference patterns arise. These patterns can be understood from a combined perspective based on fluid dynamics and wave mechanics: waves resemble fluid particles, and when individual fluid trajectories cross, the waves create interference patterns. The crossing of particle trajectories implies the transition from single to multi-stream regime. In cosmological fluid dynamics, this instance is often denoted with shell-crossing with the appearance of infinite densities Arnold et al. (1982); Hidding et al. (2014); Feldbrugge et al. (2018). In the wave-function approach, by contrast, infinite densities are naturally regularised by a nonzero $\hbar$ which acts as a softening scale. The particular semiclassical diffraction pattern emerging here from a classical caustic (called a ‘fold’ or ‘cusp’) is also known from wave optics, in particular in the framework of catastrophe optics Berry and Upstill (1980).

Furthermore, Fig. 1 displays the dispersion of the outward propagating caustics, emanating around $a=a_{\rm sc}$ and $x=0$. This smearing effect is expected from the quantum dispersion relation, and leads to an uncertainty in the space location and time of shell-crossing for $\hbar>0$.

In Fig. 2 we show the density and velocity associated with the free wave function (3), and compare the results against the classical ZA. We make use of the same 1D initial data and $\hbar$-values as before. From the time-evolution of the wave function, written as $\psi=\sqrt{1+\delta}\exp(-{\rm{i}}\phi_{\rm v}/\hbar)$, we extract the normalised density $1+\delta$, the velocity $v=-\nabla\phi_{\rm v}$ and phase $\phi_{\rm v}$. Before shell-crossing ($a\!=\!a_{\rm sc}/2$), the wave approach agrees well with the ZA, while at the time of (classical) shell-crossing, when the density in the ZA becomes infinite, the wave-density remains finite thanks to the taming scale $\hbar$. After shell-crossing ($a\!=\!2a_{\rm sc}$), the ZA leads to the known overshooting in the multi-stream regime—which is most clearly seen in the velocity plot—while the wave approach leads to interference patterns. Since the semiclassical method relies on a finite phase-space resolution associated with the small parameter $\hbar$, the density and velocity associated with the wave-function should be interpreted in a coarse-grained sense, which washes out small-scale oscillations; see Sec. V for further details.

II.2 Propagators with an effective potential

Beyond the ZA, which is the exact solution of the one-dimensional collapse problem before shell-crossing, particles will not move along straight trajectories but will be influenced by gravitational tidal effects. To transmit such tidal interactions, an effective (gravitational) potential should be included in the Hamiltonian operator. Consequently, the associated propagator $K(\boldsymbol{x},\boldsymbol{q};a)$ and wave function $\psi(\boldsymbol{x};a)$ satisfy respectively the following Schrödinger equations ($a>0$)

[TABLE]

with the Hamiltonian operator

[TABLE]

where $V_{\rm eff}$ can be considered as an external potential. An explicit expression for the effective potential is given in Eq. (8). We remark that while this formulation indeed can capture both the effects from going beyond the ZA and beyond shell-crossing, in the present paper we focus on the dynamics before (and at) shell-crossing. For a discussion of the relation of our semiclassical method to other quantum-inspired formalisms relying on a Schrödinger equation we refer to App. B.

II.3 Strategy for solving the propagator equation

The aim of this paper is to investigate perturbative solutions for the Schrödinger equation (5). For this, we discuss the following aspects of our propagator formalism in the context of perturbative schemes.

As an external potential for the Hamiltonian operator $\hat{H}$ in Eq. (6), we employ an effective potential that is determined by using standard perturbative techniques for the classical fluid equations, see Sec. III. 2. 2.

As a benchmark for our perturbative results, we will use LPT and the Cauchy invariants as diagnostic tool for detecting spurious effects in perturbative truncations, which we introduce in Sec. IV. 3. 3.

We show how quantities constructed from the wave function are related to the classical fluid variables, by employing classical limits; see Sec. V. 4. 4.

Finally, in Sec. VI, we provide a perturbative expansion for the propagator, solve the problem up to the next-to-leading order (NLO), and determine the associated fluid variables in the classical limit.

III The cosmological fluid equations and the effective potential

In this section, we relate the external potential in the Schrödinger equation (5) to an effective potential within the cosmological fluid equations. We show that, in the fluid description, this effective potential is a combination of the gravitational potential and a term due to the overall expansion of the Universe.

On sufficiently large scales and before shell-crossing, dark matter can be treated as a perfect fluid described in terms of a single-valued velocity and density. The corresponding equations are usually formulated in comoving coordinates $\boldsymbol{x}=\boldsymbol{r}/a$, where $\boldsymbol{r}$ is the physical space coordinate and $a$ the cosmic scale factor. For convenience, we decompose the fluid density $\rho(\boldsymbol{x};a)$ into a background part $\bar{\rho}(a)$ and a density contrast $\delta$, which are related to each other via $\delta=(\rho-\bar{\rho})/\bar{\rho}$. Furthermore, we make use of a peculiar velocity $\boldsymbol{v}={\rm{d}}\boldsymbol{x}/{\rm{d}}a$ which is related to the total velocity via $\boldsymbol{U}=H\boldsymbol{r}+Ha^{2}\boldsymbol{v}$, where $H$ is the Hubble parameter. In the present work, we restrict our analysis to a spatially flat universe solely filled with cold dark matter, the so-called Einstein–de Sitter (EdS) universe. A generalisation to a $\Lambda$CDM Universe (and beyond) is however straightforward.

For the case of a potential velocity $\boldsymbol{v}\equiv-\boldsymbol{\nabla}\phi_{\rm v}$, the fluid equations in an Eulerian coordinate system consist of the Bernoulli, continuity and Poisson equation, which are (see e.g. Rampf et al. (2015))

[TABLE]

where $\varphi_{\rm g}$ is the gravitational potential, and we have defined the effective potential

[TABLE]

which is a combined potential that encapsulates the effects from the cosmological potential and the Hubble friction.

Formally linearising the fluid variables around its background and evaluating the linearised fluid equations at arbitrary early times $a\to 0$, it is found that analytic solutions at $a=0$ exist only provided one makes use of the following boundary conditions

[TABLE]

These boundary conditions select the growing-mode solutions, and are in accordance with our requirement of a potential velocity.

Equipped with these boundary conditions, it is straightforward to investigate the standard perturbation theory (SPT) for the fluid equations (7). The following power series Ansätze lead to simple recursion relations that consistently solve the fluid equations order by order,

[TABLE]

All-order results for the density and velocity potential are well known (see e.g. Bernardeau et al. (2002)), from which one can easily construct all-order perturbative results for the effective potential. Explicitly, at first order we have simply $V_{\rm eff}^{(1)}=0$, whereas at second order we find

[TABLE]

with the inverse Laplacian $\nabla^{-2}$, and implied summation over dummy indices. For calculational details, including explicit all-order expressions for $V_{\rm eff}^{(n)}$, see App. A. Observe that $V_{\rm eff}^{(2)}$ is exactly zero for 1D initial conditions, in which case $\varphi_{\rm g}^{\rm(ini)}$ only depends on one spatial coordinate. In the Lagrangian-coordinates formulation, this is just the statement that the ZA is exact only in 1D; see the following section. Beyond 1D, $V_{\rm eff}^{(2)}$ is generally nonzero and constitutes the most important tidal correction to the leading order, and thus should be included for realistic modelling of the gravitational instability.

Before shell-crossing, one can relate the fluid picture to the Schrödinger equation, by using the Madelung polar form for the wave function $\psi=\sqrt{1+\delta}\exp(-{\rm{i}}\phi_{\rm v}/\hbar)$ in Eq. (5). This way, one reproduces the fluid-type equations (7), with a Bernoulli equation (7a) that receives a new term $\sim\hbar^{2}$ Madelung (1927); Short and Coles (2006b). Solving these fluid-type equations perturbatively with the Ansätze (III) reveals that $V_{\rm eff}$ agrees up to second order with the one obtained from the pure fluid approach, see App. A. Thus, up to second order in perturbation theory, which is our focus here, the effective potential can be identically obtained from either of the two approaches.

Equipped with an expression for the effective potential appearing in the Schrödinger equation (5), we have the ingredients to perturbatively compute solutions for the propagator and hence the wave function. Before proceeding with the perturbative treatment of the propagator in Sec. VI, let us provide a motivation based on Lagrangian perturbation theory in Sec. IV and discuss how our semiclassical formalism can be used to infer classical observables in Sec. V.

IV Lagrangian perturbation theory and Cauchy invariants

Perturbative solutions to the fluid equations can also be obtained in Lagrangian coordinates, and these solutions are determined within the framework of Lagrangian perturbation theory (LPT). Utilising a Lagrangian-coordinates approach has several advantages over the Eulerian approach. From the theory side, there exists a mathematical proof of time-analyticity only for the Lagrangian-coordinates approach. LPT is therefore a convergent perturbation theory for the fluid equations and can solve them in the single-stream regime to arbitrary high accuracy Zheligovsky and Frisch (2014); Rampf et al. (2015).

From the numerical side, $N$-body simulations of cosmic structure formation naturally initialise and evolve particles by employing coordinates that comove with the fluid. In a combined effort to close the gap between theory and numerics, a Lagrangian-coordinates approach appears thus to be most suitable.

We denote by $\boldsymbol{q}$ the Lagrangian coordinates, and partial derivatives with respect to the component $q_{i}$ on a given function $f$ as $f_{,i}$ (or by $\nabla_{i}^{\rm L}f$). As before, summation over repeated indices is implied. Let

[TABLE]

be the Lagrangian map from initial ($a\!=\!0$) position $\boldsymbol{q}$ to final Eulerian position $\boldsymbol{x}$ at time $a$, where $\boldsymbol{\xi}(\boldsymbol{q};a)$ is the Lagrangian displacement field. The map is defined in such a way that its Lagrangian (convective) $a$-time derivative $\partial_{a}^{\rm L}$, henceforth denoted by an overdot, returns the Lagrangian representation of the velocity, i.e.,

[TABLE]

At initial time, the velocity is $\boldsymbol{v}^{\rm(ini)}(\boldsymbol{q})=\boldsymbol{v}(\boldsymbol{x}(\boldsymbol{q};0);0)$ which agrees with the initial Eulerian velocity. Until the first shell-crossing, mass conservation $\bar{\rho}[1+\delta(\boldsymbol{x})]{\rm{d}}^{3}x=\bar{\rho}\,{\rm{d}}^{3}q$ can be exactly integrated to give

[TABLE]

where $J=\det x_{i,j}$ is the determinant of the Jacobian matrix $x_{i,j}$. Momentum conservation for the Lagrangian map can be conveniently written as

[TABLE]

The considered flow is potential in Eulerian coordinates, which implies that the Eulerian curl of the Eulerian velocity vanishes, i.e., $\boldsymbol{\nabla}_{\!x\,}\!\!\times\boldsymbol{v}\!=\!\boldsymbol{0}$. By employing the Lagrangian map, the statement of zero vorticity translates in Lagrangian coordinates to the so-called Cauchy invariants, with components ($i=1,2,3$) Rampf (2012); Zheligovsky and Frisch (2014)

[TABLE]

where $\varepsilon_{ijk}$ is the Levi-Civita symbol. See Matsubara (2015) for a simple derivation of the Cauchy invariants which however holds only for zero vorticity; see Rampf et al. (2016) where initial vorticity is included. We note that in 1D, the Cauchy invariants are trivially satisfied, since in 1D there is no vorticity.

Equations (15) constitute a closed set of Lagrangian evolution equations, with all the dynamical information being encoded in the displacement field $\boldsymbol{\xi}$ (from which one can recover the density and velocity). To obtain perturbative solutions, one then writes for the displacement components

[TABLE]

where the time-Taylor coefficients are given by simple all-order recursion relations Zheligovsky and Frisch (2014). For example, the first coefficients are

[TABLE]

with the kernel

[TABLE]

Observe the structural similarity of $\mu_{2}$ and $\nabla^{2}V_{\rm eff}^{(2)}$ (equation (11)). At second order, the only difference between these two expressions is their Lagrangian vs. Eulerian coordinate dependences and derivatives.

The expressions $\xi_{i}^{\rm ZA}(\boldsymbol{q};a)$ and $\xi_{i}^{\rm 2LPT}(\boldsymbol{q};a)$ are the well-known Zel’dovich Zeldovich (1970); Buchert and Ehlers (1993) and 2LPT Buchert (1994); Bouchet et al. (1995) approximations, respectively. The ZA, which is exact in one dimension, states that fluid elements move along straight trajectories with a velocity set by the gradient of the initial gravitational potential $\varphi_{\rm g}^{\rm(ini)}$. The 2LPT extension includes gravitational tidal effects, thereby taking into account the fact that gravity is non-local (beyond 1D).

Observe that both the Zel’dovich and 2LPT solutions give a displacement that is a gradient field and hence possesses a scalar potential. This property ceases to be true beyond second order, where one needs to include transverse modes, which is precisely the stage where one is forced to evaluate the Cauchy invariants (15b).

IV.1 Spurious vorticity in perturbation theory

While the Cauchy invariants are non-perturbative expressions, in LPT they translate into perturbative relations. At fixed order $n$ in LPT, we define the corresponding “truncated” Cauchy invariants by

[TABLE]

where $\xi_{l}^{[n]}(\boldsymbol{q};a):=\sum_{i=1}^{n}\xi_{l}^{(i)}(\boldsymbol{q})\,a^{i}$ and similarly for $\dot{\xi}_{l}^{[n]}$. Even for vanishing vorticity, Eqs. (19) are generically nonzero, as a consequence of the said truncation. Indeed, only at first order it is easily seen that ${\cal C}_{i}^{[1]}$ is still zero, but at the next truncated order we have

[TABLE]

where for later convenience we have expressed the Lagrangian kernel $\mu_{2}(\boldsymbol{q})$ in terms of $V_{\rm eff}^{(2)}(\boldsymbol{q})$. This shows that the truncated Cauchy invariants are evidently nonzero at $\mathcal{O}(a^{2})$, and hence consistent with the Cauchy invariants in a perturbative sense. However, those higher-order terms can excite unwanted spurious vorticity effects that are inconvenient when setting up initial conditions for $N$-body simulations. Such unwanted spurious effects could be artificially amplified when a 2LPT-scheme with insufficient time-stepping is employed to generate fast mock simulations. Furthermore, from theoretical grounds it is expected that such spurious effects could be enhanced at late times. Indeed, the convergence radius of LPT, which should set the maximal step-size for such algorithms, depends on the inverse of the norm of the velocity gradients Zheligovsky and Frisch (2014); Rampf et al. (2015). As a consequence, at late times, when velocity gradients become large, the convergence radius will naturally shrink. If a fixed time step is used in these algorithms that is larger than the radius of convergence, then the error in the velocity will grow as $\sim a^{2}$. In the central panel of Fig. 6 we provide numerical evidence of spurious vorticity generation by using 2LPT.

In the subsequent Sec. V, we employ classical limits to relate the propagators and wave functions introduced in Sec. II to the Lagrangian fluid variables discussed here. As we will demonstrate explicitly in Sec. VI, our propagator formalism avoids leading-order spurious vorticity effects. This is a consequence of the underlying Hamiltonian structure of the Schrödinger equation (5), which ensures the conservation of certain integral invariants associated with vorticity, which we review next.

IV.2 Invariants associated with vorticity

Even for an initially potential flow, in two or higher dimensions, an effective vorticity will be generated beyond shell-crossing. In classical collisionless dynamics, this effective vorticity arises through an averaging of the density weighted multi-stream velocities. The fact that vorticity can only arise due to multi-stream averaging is owed to the Helmholtz theorem for the conservation of vorticity flux (which can be linked to the more global Poincaré invariant). This theorem is closely related to the circulation theorem which is usually attributed to Lord Kelvin. Using the Cauchy invariants and Stoke’s theorem, we obtain the (combined) Kelvin-Helmholtz theorem that states that Frisch et al. (2017)

[TABLE]

is an integral invariant, which is restricted to closed integral curves $C(a)$ in configuration space that follow the inviscid flow. Those curves describe the boundary of a surface, which initially is $S^{\rm(ini)}$ and has a corresponding oriented surface element $\boldsymbol{{\rm{d}}S}^{\rm(ini)}$. The RHS of (21) is evaluated at initial time, which highlights that the integral invariant is a constant of motion. In the present case where there is no initial vorticity, $\Gamma=0$.

The Cauchy invariants are essentially a local form of this integral invariant. The conservation of those invariants associated with vorticity, especially beyond shell-crossing, can also be understood from a fluid perspective. For an inviscid fluid that is at some earlier time in the single-stream regime and irrotational, the evolution equation for the displacement field is sourced by a gradient of the gravitational potential, which thus implies that there is no source of vorticity (see Eq. (15a)). Now, after shell-crossing, when transiting from a single-fluid to a multi-fluid description, the evolution of each fluid stream is still governed by the same fluid-like evolution equation, with the only addition that the individual streams are now coupled gravitationally to the other streams through a common gravitational potential. However, crucially, the evolution equation for a given fluid stream is still sourced by a gradient of a gravitational potential. Thus, even in the multi-fluid regime, there is no source of vorticity, and each fluid stream remains potential at all times. However, as outlined above, there is generally an effective vorticity that results from averaging over the multiple streams.

V Relating the propagator to fluid observables

Since our goal is to use the propagator formalism to obtain improved perturbative solutions for the fluid equations, we need a dictionary for relating the propagator and wave function to the classical fluid variables. This is most conveniently done by relying on a phase-space formulation of quantum mechanics, which translates propagator solutions for nonzero $\hbar$ to observables and is also ideally suited to determine the classical limits ($\hbar\to 0$).

As we have introduced in Sec. II, the propagator evolves the wave function from its initial to the final state

[TABLE]

with the initial condition \psi(\boldsymbol{x};\text{\footnotesizea!=!0})=\psi^{\rm(ini)}(\boldsymbol{q}). We will consider wave functions constructed from propagators that satisfy the Schrödinger equation (5). To establish a connection between the wave function (22) and the classical fluid variables, we employ a method from quantum mechanics for studying quantum corrections to classical statistical mechanics — which is a closely related problem. Following those ideas, we construct a phase-space distribution function $f(\boldsymbol{x},\boldsymbol{p})$ from the wave function $\psi(\boldsymbol{x})$ using the Wigner function Wigner (1932),222We note that the Wigner distribution function is technically not a proper phase-space distribution when resolved on phase-space scales smaller than $\hbar$, since it can be negative and thus is a quasi-probability distribution that escapes the simple interpretation as a probability density. Hence, one should interpret equation (23) in a coarse-grained sense avoiding violation of uncertainty relations, which can be formalised by using the Husimi distribution Husimi (1940). Since we will be interested in the classical limit, this coarse-graining scale will ultimately become superfluous.

which depends explicitly on a phase-space coarse-graining scale $\hbar$,

[TABLE]

where both $\psi$ and $f_{\rm W}$ are functions of time $a$, and $\bar{\psi}$ indicates the complex conjugated wave function. For convenience we have absorbed the particle mass $m$ in the parameter $\hbar$, and have included the factor $a^{-3/2}$ in front of the momentum $\boldsymbol{p}$. This factor stems from the fact that our wave function is defined in terms of a peculiar velocity that is related to the conjugate momentum via $a^{3/2}$ in EdS; see App. B for details.

The way the Wigner distribution $f_{\rm W}(\boldsymbol{x},\boldsymbol{p})$ is constructed guarantees that all phase-space information is encoded in the wave function. It is built in such a way, that the normalised density $\varrho\equiv 1+\delta$ and the mean peculiar momentum $\boldsymbol{j}=(1+\delta)\boldsymbol{v}$ are obtained as the first two kinetic moments333Within the single particle probabilistic Copenhagen interpretation of quantum mechanics, $\rho$ is usually called the ‘probability density’ and $\boldsymbol{j}$ the (conserved) ‘probability flux’.

[TABLE]

Note that the velocity can be written as a gradient of the phase, $\boldsymbol{v}(\boldsymbol{x})=-\boldsymbol{\nabla}\phi_{\rm v}(\boldsymbol{x})$, of the wave function $\psi=\sqrt{1+\delta}\exp(-{\rm{i}}\phi_{\rm v}/\hbar)$, if and only if the amplitude and phase are sufficiently smooth. Shell-crossing however causes strongly oscillatory behaviour, see Fig. 2, which also generates vorticity as we shall discuss later in Secs. V.2 and VII.3. This vorticity can be nonetheless extracted from the velocity field $\boldsymbol{v}=\boldsymbol{j}/\rho$ using the just introduced kinetic moments.

While in principle we could work with explicit expressions for density and momentum, the Wigner function provides a concise and elegant way of simultaneously encoding density and velocity information, which allows us to infer the Lagrangian displacement and corresponding velocity from its classical limit.

V.1 Fluid variables from the classical limit

Taking the classical limit $\hbar\rightarrow 0$, after having obtained the solutions for the Wigner distribution (23) for nonzero $\hbar$, we obtain the phase-space distribution of a perfect fluid

[TABLE]

with a velocity $\boldsymbol{v}(\boldsymbol{x})$ that is single-valued before shell-crossing. Note that the wave-function $\psi$ itself depends on $\hbar$, as illustrated by the split in amplitude and phase, $\psi=\sqrt{\varrho}\exp(-{\rm{i}}\phi_{\rm v}/\hbar)$. Hence, the limit $\hbar\rightarrow 0$ needs to be taken with care and gives a nonzero peculiar velocity despite the $\hbar$ prefactor in Eq. (25). Using mass conservation (14), we can formulate the distribution function of the perfect fluid in Lagrangian coordinates

[TABLE]

where $\boldsymbol{\xi}(\boldsymbol{q})$ is the Lagrangian displacement (12) and $\boldsymbol{v}^{\rm L}(\boldsymbol{q})=\boldsymbol{v}(\boldsymbol{x}(\boldsymbol{q};a);a)$ is the Lagrangian representation of the velocity evaluated at the Eulerian position $\boldsymbol{x}(\boldsymbol{q};a)$. Hence, by performing the classical limit of the Wigner phase-space distribution (23) for a given wave function, we can straightforwardly read off the corresponding Lagrangian displacement and velocity.

Let us demonstrate the outlined technique for obtaining the fluid variables, by using the free theory as an instructive example. The corresponding wave function $\psi_{0}\!=\!\int{\rm{d}}^{3}q\,K_{0}(\boldsymbol{x},\boldsymbol{q},a)\,\psi^{\rm(ini)}(\boldsymbol{q})$, obtained from the free theory propagator $K_{0}$ from Eq. (2), reads

[TABLE]

where the part $\exp[-{\rm{i}}\varphi^{\rm(ini)}_{\rm g}(\boldsymbol{q})/\hbar]\equiv\psi^{\rm(ini)}(\boldsymbol{q})$ reflects the initial condition for the wave function, in accordance with the used boundary conditions (9). Plugging $\psi_{0}$ into the Wigner distribution (23), we have three integrals over $\boldsymbol{x}^{\prime}$, $\boldsymbol{q}$ and $\boldsymbol{q}^{\prime}$. The latter two integrals can be simplified with a change of variables, using center of mass $\boldsymbol{q}_{+}=(\boldsymbol{q}+\boldsymbol{q}^{\prime})/2$ and difference coordinates $\boldsymbol{q}_{-}=\boldsymbol{q}-\boldsymbol{q}^{\prime}$. We obtain

[TABLE]

where we have defined

[TABLE]

In App. C we show that this classical limit is closely related to the so-called stationary phase approximation. Returning to the integrand, performing the integrations over $\boldsymbol{x}^{\prime}$ and $\boldsymbol{q}_{-}$, we then obtain

[TABLE]

where we have renamed the integration variable according to $\boldsymbol{q}_{+}\to\boldsymbol{q}$ for convenience. Comparing this to the fluid distribution function in Lagrangian coordinates, Eq. (27), we can read off the displacement and velocity

[TABLE]

These solutions agree with those obtained from the ZA (cf. Eq. (17a)), and thus, in the classical limit and to leading order in perturbation theory, we reproduce results from classical fluid dynamics. (See Eq. (46) for the classical limit at second order of our propagator method.)

Equipped with a method to relate semiclassical propagators and wave functions to fluid observables, we will proceed to perturbatively solve the Schrödinger equations for the propagator and translate our solutions to the Lagrangian displacement and velocity field in Sec. VI.

V.2 The appearance of vorticity after shell-crossing

Under certain circumstances, the Kelvin-Helmholtz invariant $\Gamma$, given in Eq. (21), also persists for quantum and semiclassical systems. In particular, for sufficiently smooth initial conditions and by using a Madelung transformation, in Ref. Damski and Sacha (2003) it has been shown that $\Gamma$ is also an invariant under evolution with a quantum Hamiltonian (i.e., under the Schrödinger equation), if one simply replaces $\boldsymbol{v}$ with $\boldsymbol{j}/\rho$ and ensures that the integral contour goes only through regions where the velocity is well defined in the course of the evolution.

In quantum systems, vorticity is quantised Gross (1961); Feynman (1958). Since the wave function is always single-valued, quantised vorticity can only arise from topological defects where the phase factor, $\phi_{\rm v}/\hbar$, undergoes a localised phase jump of integer multiples of $2\pi$. Since $\Gamma$ has to vanish for initially irrotational systems, it is topologically required that vortices can only be produced in pairs Savchenko and Zel’dovich (1999), called rotons, i.e.,

[TABLE]

and thus, the sum of negative $n_{-}$ and positive $n_{+}$ topological charges is conserved.

Later in Sec. VII we show that at times shortly after shell-crossing, where vorticity is generated, we indeed observe the appearance of such rotons.

VI Perturbative treatment of the propagator

In the following, we will use SPT results for $V_{\rm eff}$, as discussed in Sec. III as an input to the propagator equation (5) and solve it in a perturbative fashion. From this perturbative solution, we will extract the Lagrangian displacement and velocity using the method described in Sec. V.

Since we already know the solution in the absence of the effective potential, we split the nonlinear propagator into the free propagator and an exponential term

[TABLE]

When combining the exponentials, one recognises the total action $S=S_{0}+S_{\rm tid}$ as a sum of the free particle contribution $S_{0}$ from Eq. (1) and the tidal interaction terms encoded in $S_{\rm tid}$. Plugging the Ansatz (35) for the total propagator $K$ into the evolution Eq. (5), and using Eq. (4a), one obtains a differential equation for the interaction term $S_{\rm tid}=S_{\rm tid}(\boldsymbol{q},\boldsymbol{x};a)$

[TABLE]

where $\hat{D}_{a}(\boldsymbol{q},\boldsymbol{x})\equiv\partial_{a}+(1/a)[\boldsymbol{x}-\boldsymbol{q}]\cdot\boldsymbol{\nabla}_{\!x\,}$. The source term is given in terms of a time-Taylor series for the effective potential $V_{\rm eff}(\boldsymbol{x};a)\equiv\sum_{n=2}^{\infty}V_{\rm eff}^{(n)}(\boldsymbol{x})\,a^{n-2}$, and for the interaction part of the action one can impose a PT Ansatz

[TABLE]

where $S_{n}$ is ideally ${\cal O}(a^{n})$. In the following section, we explicitly solve for the NLO part $S_{1}$, which is the leading-order contribution to $S_{\rm tid}$.

VI.1 Next-to-leading order propagator

At next-to-leading order (NLO), the effective potential $V_{\rm eff}^{(2)}(\boldsymbol{x})$ is time-independent and given by Eq. (11). Hence, the NLO contribution to $S_{\rm tid}$, called $S_{1}$, is expected to be of order $a$. Since the derivative operator, $\hat{D}_{a}$, decreases the power of $a$ by one, the other two terms on the LHS of Eq. (36) are of higher order and do not contribute to $S_{1}$. The evolution equation (36) thus simplifies to

[TABLE]

The solution of this equation is given by

[TABLE]

which is an integral over the time-independent effective potential along a straight line connecting the initial position $\boldsymbol{q}$ and final position $\boldsymbol{x}$. Since it is impracticable to evaluate this integral at every point, let us make further approximations. Preferably, we want to preserve the symmetry between initial and final positions $\boldsymbol{x}\leftrightarrow\boldsymbol{q}$ when exchanging $a\leftrightarrow-a$, which implies a time-reversible propagation. To this end, we use a two-endpoint approximation

[TABLE]

in accordance with the employed PT Ansatz (37). This approximation corresponds to a (numerical) kick-drift-kick scheme with ${\cal O}(a^{3})$ accuracy, in which the potential $V_{\rm eff}^{(2)}$ is evaluated at the initial position $\boldsymbol{q}$, the particle propagated to its final position $\boldsymbol{x}$ and then the potential evaluated there.444Instead of employing an approximation that amounts to a kick-drift-kick scheme, we could have also used a midpoint approximation that would correspond to a drift-kick-drift scheme, where the particle is first propagated from $\boldsymbol{q}$ to the midpoint $(\boldsymbol{q}+\boldsymbol{x})/2$, the potential evaluated there, and then the particle is propagated from the midpoint to the endpoint $\boldsymbol{x}$. In the present case, we prefer the kick-drift-kick procedure, because it updates the effective potential only once. Now, we can combine the solution $S_{1}$ from Eq. (40) with the free kernel $K_{0}$ from Eq. (2) into the NLO kernel $K_{\rm NLO}$ from Eq. (35) and hence the NLO wave function using Eq. (22) to get

[TABLE]

where $V_{\rm eff}^{(2)}$ is given by Eq. (11). The propagator has the structure of a numerical kick-drift-kick scheme, which will simplify the numerical implementation later on. Equations (41) constitute one of the main technical results of this paper. We note that alternatively to the above derivation, one could also solve the Schrödinger equation (5) in operator notation, as shown in App. D.

VI.2 Next-to-leading order observables in the classical limit

To extract the Lagrangian displacement and velocity from the NLO wave function, we proceed along the lines as described in Sec. V. In analogy to the free wave function (V.1), we get for $f_{\rm W,1}=f_{\rm W,1}(\boldsymbol{x},\boldsymbol{p};a)$

[TABLE]

where, in addition to $\delta\varphi$ from Eq. (30), we have defined

[TABLE]

As argued before, in the classical limit, the dominant contributions to the integral come from the terms (43) for which $\boldsymbol{q}_{-}$ is small. Thus we can write

[TABLE]

Using this, the integrations over $\boldsymbol{x}^{\prime}$ and $\boldsymbol{q}_{-}$ are easily performed. After some straightforward calculations, we obtain

[TABLE]

with the NLO displacement and (Lagrangian) velocity

[TABLE]

To arrive at the final expression for the velocity, we expanded $V_{{\rm eff},i}^{(2)}(\boldsymbol{x})\simeq V_{{\rm eff},i}^{(2)}(\boldsymbol{q})-a\varphi_{,l}(\boldsymbol{q})V_{{\rm eff},il}^{(2)}(\boldsymbol{q})$ to the leading order in the displacement, in accordance with the expansion scheme employed for the classical limit. Evidently, the NLO displacement (46a) agrees with its classical counterpart (Eq. (17)). The NLO velocity, by contrast, differs from the second-order velocity in LPT which is $v_{i}^{\rm L,LPT}\simeq-\varphi_{{\rm g},i}^{\rm(ini)}-aV_{{\rm eff},i}^{(2)}$, since the NLO velocity includes an additional term that within the LPT expansion would be considered as higher order.

We now prove that this additional term in the NLO velocity is crucial for maintaining the Hamiltonian structure of the system, and thereby not to excite spurious vorticity. We do this by using the Cauchy invariants as a diagnostic tool which we have introduced in Eq. (19). For this we need the Lagrangian map and its time derivative from our NLO formalism, which are easily obtained from the above results; they read respectively

[TABLE]

Plugging these expressions into the Cauchy invariants formula, we obtain

[TABLE]

and hence there is no vorticity generated at order $a^{2}$, in contrast to the 2LPT case, see Eq. (IV.1). For a numerical analysis related to vorticity, see Sec. VII.3.

In summary, we have thus established a direct correspondence between our semiclassical propagator method and classical fluid mechanics. In particular, the NLO displacement field coincides exactly with its LPT counterpart up to second order, whereas the NLO velocity receives naturally an additional term that is missing in LPT at this fixed perturbation order. However, as we claim, this additional term is crucial to respect the underlying Hamiltonian structure; as outlined in Sec. IV.1, ignoring this term could lead to the artificial generation of vorticity.

VII Results beyond 1D collapse

In the field of cosmological fluid dynamics, it is known that the effective gravitational potential, here dubbed $V_{\rm eff}$, is exactly zero if the initial conditions depend only on one space variable. As a consequence, in 1D, the Zel’dovich solution becomes exact in the single-stream regime. Deviating from 1D, we have generically $V_{\rm eff}\neq 0$, and, as a result, the ZA performs rather poorly as gravitational tidal effects become non-negligible. In particular, this reflects in inaccurate predictions of the ZA for the shell-crossing time that worsen successively when deviating more and more from the 1D collapse Rampf and Frisch (2017); Rampf (2019); Saga et al. (2018). Very similar performance issues are expected for the propagator method. Specifically, the free propagator (2) is only accurate in the absence of tidal forces. Based on such considerations, we have determined the NLO propagator (41b) in the presence of a nonzero effective potential by using perturbation theory.

In the following, we present numerical implementations for quasi-one-dimensional (Sec. VII.1), as well as two-dimensional collapse problems (Sec. VII.2), both with an appropriate choice of tailored initial conditions (ICs). In both cases, we provide a quantitative comparison between the ZA and 2LPT predictions versus the results from the free and NLO propagators. Numerical results for cosmological ICs will be investigated in a forthcoming paper.

Beyond 1D, the phenomenology of the flow is much richer and symmetry principles of the dynamics, such as the conservation of vorticity flux, become manifest. Therefore, in Sec. VII.3, we will analyse essential features of the vortical flow, which arise after shell-crossing and reflect the conservation of certain invariants.

VII.1 Quasi-one-dimensional collapse

For ICs that go beyond 1D, higher-order effects for both the classical and propagator formalism become non-vanishing. As outlined above, this is due to the fact that for departures from 1D the effective potential $V_{\rm eff}$ is generally nonzero. For an initial potential that is (perturbatively) close to 1D, the bulk part stemming from $V_{\rm eff}$ is captured accurately by 2LPT and the NLO propagator. This is so, since the perturbation series converges very fast for such so-called quasi-1D ICs Rampf and Frisch (2017); Saga et al. (2018). For quasi-1D ICs, we thus expect that the NLO propagator should deliver a very accurate description of the collapse. To test the performance of the NLO propagator compared to 2LPT, we investigate quasi-1D ICs for the initial wave function $\psi^{\rm(ini)}(\boldsymbol{q})\!=\!\exp\left(-{\rm{i}}\phi_{\rm v}^{\rm(ini)}(\boldsymbol{q})/\hbar\right)$, with the initial velocity potential

[TABLE]

where we have chosen $\epsilon=1/4$. We show the resulting density profiles shortly before shell-crossing in Fig. 3. In the propagator case, we first evaluate the effective potential at initial position in real space, propagate the particles in Fourier space, and then re-evaluate the potential at the final real-space position. As mentioned above, this routine corresponds to a symplectic kick-drift-kick scheme (see also App. D). To obtain the realisation for the propagator, we have used a $512^{2}$ mesh to evaluate the solution with $\hbar=3\times 10^{-2}$, while in the classical cases, we have used $256^{2}$ particles and evaluate the density field using the tessellation method of Hahn et al. (2015) on a $512^{2}$ grid. We also show the relative differences between the predictions from LPT and the propagator method in the lower panels. Examining first the density profile along the main collapse direction $x$, one sees that all methods agree reasonably well with one another. This is expected since the perturbation along the perturbatively small second dimension is not large so that the evolution is still close to being one-dimensional, for which the free and ZA solutions are good approximations. The relative differences between free and ZA are at the accuracy level of our numerical experiments here at $\sim 10^{-4}$, except at the location of future shell-crossing, where quantum corrections lead to a very localised large deviation at $\lesssim 10^{-1}$ which is a direct result of the reduction of the peak density due to finite $\hbar$. The differences between NLO and 2LPT, by contrast, are much more prominent, especially along the direction of the density ridge (right panel in Fig. 3). For both classical and propagator methods, the peak densities at leading order and next-to-leading order are significantly different. Apart from these differences at $x=-\pi/2$, the error in lower density regions is in fact mainly due to the linear interpolation used in the sheet reconstruction from the tessellation.

VII.2 2D collapse – the density field

For reasons explained above, it is expected that the impact of $V_{\rm eff}$ becomes more prominent when departing greatly from 1D initial conditions. To demonstrate this we use in the present section the 2D initial potential

[TABLE]

which represents a (strongly) phased plane wave (cf. Hahn et al. (2013)). In Fig. 4, we show the resulting two-dimensional density field for the four approaches at $a=1$, i.e., well after shell-crossing. For the Lagrangian perturbation theory results, we have used $1024^{2}$ fluid particles that have been evolved under the ZA or 2LPT in a single time step. In order to compute the multi-stream density field accurately, we have projected the tessellated dark matter sheet using the technique of Abel et al. (2012) onto a uniform mesh of $1024^{2}$. For the propagator method, we have evolved the free and NLO propagator directly on a $1024^{2}$ grid with $\hbar\!=\!5\times 10^{-3}$. Overall the results show excellent agreement between the ZA and the free propagator, as well as between 2LPT and NLO in terms of global shape of the caustics in the various multi-stream regions. Naturally, after shell-crossing, the propagator solutions show rapid oscillations that encode the multi-stream behaviour. This is related to the appearance of higher-dimensional caustics that translate into more complex diffraction patterns than in the one-dimensional case (cf. Berry and Upstill (1980)). When interpreted in a coarse-grained sense, the rapid oscillations disappear from the physical density and velocity, but encode the properties beyond the perfect pressureless fluid, in particular the vorticity that is induced by shell-crossing.

VII.3 2D collapse – the vorticity

As a final aspect of this paper, we investigate how vorticity arises in the classical and semiclassical picture, respectively. For our propagator method, we determine the vorticity $\boldsymbol{w}\equiv\boldsymbol{\nabla}\times\boldsymbol{v}$ in Fourier space by first obtaining an expression for the velocity $\boldsymbol{v}=\boldsymbol{j}/(1+\delta)$, where the RHS is evaluated by using Eqs. (24)–(25). We thus calculate

[TABLE]

where $\mathcal{F}\{\cdot\}$ is a fast Fourier transform (FFT). Since the vortices are point-like, the inverse FFT produces heavy ringing, so that we have to additionally filter the vorticity fields. In order to avoid convolving transversal and longitudinal velocity components, we multiply with a Gaussian filter $\exp\left[-k^{2}/k_{s}^{2}\right]$, where $k_{s}$ is a filter scale, directly in Fourier space when also taking the cross product with $\boldsymbol{k}$. For the LPT prediction of vorticity generation, by contrast, we use the method of Hahn et al. (2015) to explicitly carry out the multi-stream average.

In Fig. 5, we show the wave function using domain colouring in the left panel, along with the semiclassical vorticity for two different smoothing scales $k_{s}$ ($1/4$ and $1/16$ of the Nyquist wave number) in the two middle panels. When using the smaller smoothing scale (second panel from left), one can clearly see that the vortices are indeed point-like objects in two dimensions, which have a positive (red) or negative (blue) sign, and are concentrated around the caustics. In comparison with the full wave function (leftmost panel), one sees that the vortices are always associated with dark regions, where the amplitude of the wave function vanishes. As discussed in Sec. V.2, vorticity is conserved and has thus to be pair-produced with opposite topological charge. These neighbouring positive and negative vortices are clearly visible in Fig. 5. For the larger softening (third panel from left), it becomes obvious that, when averaging over multiple such quantum vortices, one obtains a large-scale limit which is very similar to the vorticity pattern obtained for the ZA (rightmost panel). The agreement of the properties of the two-dimensional flow between classical and quantum dynamics after filtering has been discussed in detail for the cosmological Schrödinger equation by Kopp et al. (2017).

Remarkably, the propagator method allows us to explicitly predict the generation of vorticity without requiring a numerical solution of the Schrödinger equation. Furthermore, the propagator method does not require multi-stream averaging. This should be contrasted to the classical (multi)-fluid picture, where multi-stream averaging is mandatory Hahn et al. (2015) and computationally involved, even for simple cases like the ZA Pichon and Bernardeau (1999).

Finally, in Fig. 6, we display the vorticity at times shortly after the first shell-crossing for the smoothing scale $k_{s}=k_{\rm Ny}/8$. For reasons of comparison we have also performed an $N$-body simulation (leftmost panel) which has been initialised at $a_{\rm ini}=1/30$. In all panels the generation of vorticity through multi-streaming is visible by the appearance of “vorticity islands”, while 2LPT predicts the spurious generation of vorticity in single-stream regions (see Sec. IV.1 for a related discussion). We also note that the white lines in the 2LPT panel indicate where the spurious vorticity changes sign. Our NLO results, by contrast, are free of spurious vorticity; see e.g. Eq. (48) for the provided theoretical argument.

VIII Conclusions

A. Summary. We have introduced a novel semiclassical method for evolving CDM that, for nonzero $\hbar$, is free from any singular behaviour near the crossing of particle trajectories (shell-crossing). The key quantity of our method is the propagator, which fulfils a Schrödinger equation and encodes the transition amplitude of the wave function for the particle trajectories. This propagator can be determined by using perturbation theory. In a suitable coordinate system, the leading-order propagator dictates a ballistic motion of the particles with prescribed velocity, which amounts to the classical Zel’dovich approximation.

To incorporate gravitational interactions in the propagator, we have motivated the inclusion of an effective potential in the Schrödinger equation, which is also present in the cosmological fluid equations [Eq. (8)]. Using standard perturbative techniques for the effective potential as an input for the Schrödinger equation, we have solved the associated propagator equation at next to leading order (NLO). The NLO solution for the propagator delivers the associated wave function, Eq. (41), from which the NLO density and velocity are easily obtained [Eqs. (24)–(25)].

By performing the classical limit, we have shown that our NLO result returns a displacement field (46a) and a corresponding density that are in agreement with the ones from Lagrangian perturbation theory (LPT) up to second order. The associated NLO velocity (46b) receives naturally an additional term ${\cal O}(a^{2})$, that in LPT would be considered as third order, but is actually important to preserve the underlying Hamiltonian structure of the system and to avoid a spurious excitation of vorticity in certain implementations of second-order LPT (see Sec. IV.1).

We compared our NLO results from the propagator method against LPT for two types of initial conditions, and find overall good quantitative agreement, see Figs. 3–4. The propagator method regularises classical caustics, which leaves subtle imprints while preserving an overall good agreement in their global shapes and positions. Furthermore, based on our analytical solutions for the propagator, we have demonstrated that our method is free of spurious vorticity generation (Fig. 6), as well as is capable to predict the generation of vorticity after shell-crossing (Fig. 5). Although the latter vorticity indeed arises solely through multi-stream dynamics, in our propagator method no explicit multi-stream averaging is required (see Sec. VII.3 for a related discussion). In our formalism, vorticity manifests by the pair creation of topological defects, usually called rotons, along classical caustics.

B. Outlook. For simplicity, in the present work we have applied our propagator method to two-dimensional collapse problems only, however all provided tools are ready to use for full 3D calculations. Surely, the phenomenology for cosmological ICs will be much richer than explored for the present case studies, and therefore will be investigated in a forthcoming paper. Another interesting avenue would be to include higher-order corrections in the propagator. For this one requires the external potential to third order, which we discuss in App. A. Finally, our propagator method already shed some light on the highly complicated regime shortly after shell-crossing (including the generation of vorticity, velocity dispersion, etc.). For computations well beyond shell-crossing, however, the effective potential should be updated in order to grasp the full-fledged multi-stream regime, which is required to approach virialisation.

ACKNOWLEDGMENTS

We thank Diego Blas, Lam Hui, Michael Kopp, Giorgio Krstulovic and Aseem Paranjape for many useful discussions and comments on the draft. We also thank the Wolfgang Pauli Institute, and in particular Norbert Mauser, for hospitality during the initial stages of this work and the participants of the workshop on ‘Numerics in cosmology: around the Schrödinger method’ for interesting discussions. CU kindly acknowledges funding by the STFC grant no. RG84196 ‘Revealing the Structure of the Universe’. CR thanks the Lagrange Laboratory of the Observatoire de la Côte d’Azur for financial support. CR acknowledges funding from the People Programme (Marie Curie Actions) of the European Union H2020 Programme under Grant Agreement No. 795707 (COSMO-BLOW-UP). MG acknowledges support from the Marsden Fund of the Royal Society of New Zealand and from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 616170. OH acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 679145, project ’COSMO-SIMS’).

Appendix A Effective potential from standard perturbation theory

Here we will show how the standard recursion relations from SPT can be used to determine the effective potential (8), appearing in the Bernoulli equation (7a), in an $a$-time expansion (III). For notational simplicity, we use the following shorthand notations $\varphi_{\rm g}^{\rm(ini)}\rightarrow\varphi^{\rm(ini)}$ for the initial gravitational potential, and $\phi_{\rm v}\rightarrow\phi$ for the velocity potential.

The recursion relations for the perturbative expansion of the density contrast $\delta$ and velocity potential $\phi$ from Eq. (III) give the following first two terms Bernardeau et al. (2002)

[TABLE]

where the derivatives and dependences are w.r.t. Eulerian coordinates. To get an expression for $V_{\rm eff}^{(n)}$, we simply plug $\phi^{(n)}$ into the perturbed Bernoulli equation (7a), which can be written as

[TABLE]

The first three solutions are

[TABLE]

Note that $V_{\rm eff}$ in the expansion (III) is only time-independent up to second order, while the third-order term is proportional to $a$-time $aV_{\rm eff}^{(3)}$.

As mentioned earlier, before shell-crossing, one can transform the Schrödinger equation (5) into a fluid-like system using a polar form for the wave function $\psi=\sqrt{1+\delta}\exp(-{\rm{i}}\phi/\hbar)$ Madelung (1927). One obtains a system equivalent to (7), but with an extra term in the Bernoulli equation (7a) that adds to the effective potential

[TABLE]

Note that, if one perturbatively solves the corresponding wave-mechanical fluid equations, in analogy to SPT for the fluid equations, one obtains identical solutions up to second order. The reason is that the leading order correction term in (60) is proportional to $\boldsymbol{\nabla}^{2}\delta\propto a$ and hence only enters in the effective potential at third order, where it modifies (59) to

[TABLE]

Appendix B Relation between the propagator formalism and Schrödinger-Poisson

While the Schrödinger equation (5) used here might look similar to the Schrödinger method Widrow and Kaiser (1993); Uhlemann et al. (2014) for approximating classical dynamics through the quantum-classical correspondence, they are physically distinct and rely on different assumptions, as we show now. The cosmological Schrödinger-Poisson equation in the formulation of Widrow & Kaiser Widrow and Kaiser (1993) reads

[TABLE]

where for simplicity we have absorbed the mass in $\hbar$ and set $4\pi G\bar{\rho}_{0}=3/2$. To distinguish the wave function in the different formulations, we attach a tilde to the wave function in the Widrow & Kaiser formulation. The Schrödinger-Poisson equation can be regarded as an approximate treatment of the phase-space dynamics described by the Vlasov-Poisson equation Uhlemann et al. (2014); Uhlemann (2018). When solved numerically, it provides a field-based method that complements particle-based $N$-body simulations and whose accuracy is controlled by the phase-space resolution $\hbar$ Kopp et al. (2017); Garny and Konstandin (2018); Mocz et al. (2018). Additionally, the Schrödinger-Poisson equation is a physical model for scalar field (or wavelike) dark matter such as ultralight axions, where it has attracted considerable attention Schive et al. (2014); Veltmaat and Niemeyer (2016); Hui et al. (2017); Li et al. (2018); Veltmaat et al. (2018); Edwards et al. (2018) recently.

Evidently, the physical time Schrödinger equation (62) does contain a gravitational potential $V$, but not a difference of a gravitational and velocity potential as our effective potential $V_{\rm eff}$. This is related to the fact that the associated wave function encodes the conjugate momentum. Employing the Madelung split, we can write the wave function as $\tilde{\psi}=\sqrt{1+\delta}\exp(-{\rm{i}}\tilde{\phi}_{\rm p}/\hbar)$, where the phase $\tilde{\phi}_{\rm p}$ is the potential of the conjugate momentum. Since the conjugate momentum is defined as $\boldsymbol{p}/m=a^{2}{\rm{d}}\boldsymbol{x}/{\rm{d}}t$, which is related to our peculiar velocity $\boldsymbol{v}={\rm{d}}\boldsymbol{x}/{\rm{d}}a$ via $\boldsymbol{p}/m=\boldsymbol{v}a^{2}{\rm{d}}a/{\rm{d}}t=a^{3/2}\boldsymbol{v}$ in an EdS universe, the conjugate momentum potential is related to our peculiar velocity potential as $\tilde{\phi}_{\rm p}=a^{3/2}\phi_{\rm v}$. Substituting this relation in the wave function one has $\tilde{\psi}=\sqrt{1+\delta}\exp(-{\rm{i}}a^{3/2}\phi_{\rm v}/\hbar)$.

Similarly as done before, we can convert the Schrödinger equation for $\tilde{\psi}$ into fluid-type equations. Separating real and imaginary parts of Eq. (62), leads to the fluid equations in physical time (see Eqs. (14) in Uhlemann et al. (2014)). Those fluid equations can be rewritten in $a$-time to obtain a fluid-like system (7) with a modified effective potential that sources the Bernoulli equation (7a) for $\phi_{\rm v}=a^{-3/2}\tilde{\phi}_{\rm p}$ according to

[TABLE]

which carries an extra $a^{-3}$ dependence in the $\hbar^{2}$-dependent ‘quantum’ term in the effective potential, as compared to our expression (60). Using the peculiar velocity potential, we obtain an effective potential in an $a$-time Schrödinger equation for $\psi=\sqrt{1+\delta}\exp(-{\rm{i}}\phi_{\rm v}/\hbar)$, but since the quantum potential term has a different time-dependence, it cannot be absorbed in the Laplacian, as was the case for our Schrödinger equation (5). One can view the difference in the fluid equations, or the effective potential, as a time-dependent phase-space coarse-graining scale $\hbar$ (or a time-dependent mass). This property of a time-dependent coarse-graining persists in the full-fledged multi-stream regime, where the evolution is governed by the Vlasov-Poisson equations, see Uhlemann (2018).

Appendix C Stationary-phase approximation

The stationary phase approximation (SPA) Hörmander (2015) can be used to estimate integrals of the following type

[TABLE]

In the present context of the paper, evaluating such integrals in the vicinity of $\lambda\to 0$ is relevant when performing the classical limit.

For $\lambda\to 0$, the SPA states that the dominant contribution from the above integral comes from

[TABLE]

where one sums over all critical points $\boldsymbol{q}_{c}$ for which the first-order Taylor coefficient $\partial g(\boldsymbol{q})/\partial q_{i}|_{\boldsymbol{q}=\boldsymbol{q}_{c}}$ vanishes and $H_{ij}=\partial^{2}g(\boldsymbol{q})/(\partial q_{i}\partial q_{j})|_{\boldsymbol{q}=\boldsymbol{q}_{c}}$ is the Hessian. The signature of the Hessian, which is the difference between the number of negative and positive eigenvalues, determines the prefactor in the last exponential of (65). Note that this is just a Wick rotated version of the formula for the method of steepest descent, or so-called saddle-point approximation.

Application to free theory observables

Let us apply the SPA to the free wave function (28). This wave function has the form (64) with the constants $\lambda=a\hbar$, $h_{0}^{\psi}=\exp(-3{\rm{i}}\pi/4)$ and the exponential

[TABLE]

Considering the classical limit $\hbar\rightarrow 0$, we obtain the condition for a point of stationary phase

[TABLE]

which implicitly determines the critical point $\boldsymbol{q}_{c}$. Before shell-crossing, for every $\boldsymbol{x}$ there exists only one such critical point $\boldsymbol{q}_{c}$. To perform the integral over $\boldsymbol{q}$ in (28), we expand $g_{0}(\boldsymbol{x},\boldsymbol{q})$ in a Taylor series up to quadratic order around the critical point $\boldsymbol{q}_{c}$

[TABLE]

Then, one can shift the integration $\boldsymbol{q}\rightarrow\tilde{\boldsymbol{q}}=\boldsymbol{q}-\boldsymbol{q}_{c}$ and perform the following Gaussian integral

[TABLE]

The corresponding wave function is

[TABLE]

where $\boldsymbol{q}_{c}=\boldsymbol{q}_{c}(\boldsymbol{x})$ according to the stationary phase condition (65). From this, we can easily compute the density according to the definition (24), which reads

[TABLE]

from which one evidently recovers the Zel’dovich displacement field (33). By exactly the same arguments, the SPA for the NLO wave function delivers the 2LPT displacement field.

To derive the velocity, one needs a SPA for the spatial gradient of the wave function. In this case one has the prefactor $h_{0}^{\nabla\psi}={\rm{i}}(\boldsymbol{x}-\boldsymbol{q})\exp(-3{\rm{i}}\pi/4)/(a\hbar)$, which is evaluated using the stationary phase condition (65) to obtain

[TABLE]

From this expression, we can now compute the velocity according to Eq. (25) to get

[TABLE]

which agrees with the previous result from Eq. (33) for the velocity in Lagrangian coordinates.

Appendix D Perturbative propagator in operator notation

A particularly concise perturbative expansion can be obtained from an operator expansion. To solve the Schrödinger equation (5) in operator notation, we write the Hamiltonian operator (6) at leading order, where it is time-independent

[TABLE]

For simplicity, we denote the operators associated to kinetic and potential energy as $\hat{T}$ and $\hat{V}$, respectively. Since neither of those operators have explicit time dependence, the Schrödinger equation (5) can be integrated

[TABLE]

The Baker-Campbell-Hausdorff formula (or the equivalent reverse Zassenhaus formula) allows to express exponentials of sums of operators as products of exponentials according to

[TABLE]

where $\epsilon=-{\rm{i}}a/\hbar$, and higher-order terms come from nested commutators which are denoted with $[\cdot,\cdot]$. Using a threefold decomposition of the Hamiltonian $\hat{H}$ in $VTV:=\hat{V}/2+\hat{T}+\hat{V}/2$ allows to both arrive at a time-symmetric formula and cancel all even-order correction terms (and thus the leading order is proportional to $\epsilon^{2}\sim a^{2}$)

[TABLE]

This provides another motivation for the VTV approximation of (40), which leads to the NLO propagator (41b). The approach outlined above is of course equivalent also to the usual operator-split approach followed by symplectic schemes when numerically integrating classical Hamiltonian systems, cf. Yoshida (1990).