Cosmological fluid dynamics in the Schr\"odinger formalism

We investigate the dynamics of a cosmological dark matter fluid in the Schr\"odinger formulation, seeking to evaluate the approach as a potential tool for theorists. We find simple wave-mechanical solutions of the equations for the cosmological homogeneous background evolution of the dark matter field, and use them to obtain a piecewise analytic solution for the evolution of a compensated spherical overdensity. We analyse this solution from a `quantum mechanical' viewpoint, and establish the correct boundary conditions satisfied by the wavefunction. Using techniques from multi-particle quantum mechanics, we establish the equations governing the evolution of multiple fluids and then solve them numerically in such a system. Our results establish the viability of the Schr\"odinger formulation as a genuine alternative to standard methods in certain contexts, and a novel way to model multiple fluids.


INTRODUCTION
On scales smaller than the Hubble distance, the evolution of structure in the universe can be analysed using Newtonian equations for gravitational instability (Jeans 1928). In particular, the phase space evolution of the dominant cold dark matter (CDM) component is governed by a collisionless Boltzmann, or Vlasov, equation, coupled to a Poisson equation for gravity. When velocity dispersion is negligible, a fluid approximation holds, which leads to a much simpler representation, but one that fails to accommodate multiple streams in the fluid, most notably at the onset of so-called 'shellcrossing' in the non-linear evolution of density perturbations.
In the linear regime, one can (by definition) describe structure growth accurately using first-order perturbation theory (Peebles 1980) applied to the standard Newtonian equations for a selfgravitating fluid. In the strongly nonlinear regime, only numerical simulations are capable of modelling general fluid evolution and reproducing the structures observed in galaxy surveys; for a recent review, see Springel et al. 2006. However, the evolution of density fluctuations in the so-called 'weakly nonlinear regime', those characteristically represented on the largest scales in galaxy clustering surveys, have been well described by a variety of semi-analytic methods. These include nonlinear perturbation theory, and related schemes such as the Zeldovich and adhesion approximations (see Bernardeau et al. 2002 for a review). These methods typically extrapolate results from the linear regime into the quasi-linear regime.
In a separate line of research, the correspondence between the equations describing a collisionless fluid and the Schrödinger equation was first identified by Madelung in 1926. Outside cosmology, this correspondence is most commonly used in drawing a hydrodynamic analogy for quantum mechanics (Spiegel ⋆ e-mail: rrj20@mrao.cam.ac.uk 1980;Skodje et al. 1989). Widrow & Kaiser (1993) were the first to apply the ideas to the problem of cosmological structure formation, using a Schrödinger field representation for collisionless matter, but omitting the 'quantum pressure' term to obtain a linear wave-equation. They developed a numerical technique to follow the nonlinear evolution of the field, offering an alternative to N-body particle mesh codes. A fuller explanation and exploration of their method is given by Davies & Widrow (1997). Coles & Spencer (2002) take a different approach, using the so-called Madelung transformation to establish a simpler correspondence between the self-gravitating fluid equations and a Schrödinger-Poisson system. This transformation means that their representation of the wavefunction can only be used for single streaming fluids, with their aim to develop a semi-analytic technique for tracking the development of structure into the quasilinear regime. Short & Coles (2006a,b) develop this foundational work into what they term the free-particle approximation for precisely this purpose, in which the gravitational potential and the quantum pressure are neglected. They show that this approximation is useful into the mildly nonlinear regime and the resulting method is essentially equivalent to the well-known adhesion approximation. Szapudi & Kaiser (2003) develop a Schrödinger perturbation theory very similar in nature to the well-known perturbation theory of the standard Eulerian fluid equations.
With the exception of Szapudi & Kaiser, all previous work in this area has been aimed chiefly at developing better computational techniques for tracking the evolution of structure, based almost incidentally on a wave mechanical formulation, which may confer certain advantages over rival methods; the full wave-mechanical system of equations that govern structure growth has not been a feature of interest. In this paper, our aim is to explore the Schrödinger formulation of the evolution of structure from a fully self-consistent point of view, without recourse to approximations. c 2009 RAS The plan of the paper is as follows. In Section 2, we outline the Schrödinger formalism for fluid dynamics, summarizing the previous work in the field and outlining our own direct approach. In Section 3, we obtain solutions to the Schrödinger-Poisson system for spatially-flat and closed homogeneous cosmological models. These solutions are then used in Section 4 to derive a wave-mechanical solution for a compensated spherical overdensity model. In Section 5, we describe an extension to the Schrödinger formalism, based on techniques from multi-particle quantum mechanics, that allows one to accommodate multiple fluids. We apply this technique to a modified spherical overdensity model that exhibits multi-streaming by integrating the evolution equations numerically. Finally, our conclusions are presented in Section 6

SCHRÖDINGER FORMALISM FOR FLUIDS
We take as our starting point the usual Newtonian equations for a self-gravitating perfect fluid. The continuity, Euler and Poisson equations comprise the first and second momentum moments of the full Vlasov-Poisson system, and are exact provided the velocity dispersion is zero, i.e. for single streaming fluids. These equations are: Here v, ρ and p are the fluid velocity, density and pressure respectively, and V is the gravitational potential. In (3), we have modified the Poisson equation to include a term representing the repulsive force provided by a non-zero cosmological constant Λ.
Adopting the usual description of cold dark matter as a pressureless fluid, and restricting our analysis to irrotational flow, which is normally considered for cosmological evolution in homogeneous and isotropic universe models, one can replace the Euler equation (2) with the Bernoulli equation with v = ∇φ where φ is the velocity potential. Since our equations describe a single streaming fluid, one may introduce the Madelung transformation In this transformation, ρ = |ψ| 2 = α 2 , so the fluid density is always nonnegative, and ν is an adjustable parameter which has the same dimension as φ, namely [L 2 T −1 ]. It is then easy to show that the equations governing the evolution of ψ and V are Thus, we see that the dynamics of a dark matter fluid may be described by the (modified) Poisson equation coupled to a nonlinear wave equation. This permits a direct formal correspondence to be made between the concepts of fluid dynamics and those of quantum mechanics, leading to potential conceptual advantages such as fluid quantities being associated with Hermitian operators. There are also some possible computational simplifications. In particular, the WKB limit for the wave-equation corresponds to the asymptotic study of the limit ν → 0 (Spiegel 1980). Nonetheless, the nonlinear nature of the wave-equation (6) does present some conceptual difficulties. Nonlinear Schrodinger equations and their applications in (non-cosmological) physics are treated thoroughly in Sulem & Sulem (1999) and Pang & Feng (2005). The theory of nonlinear quantum mechanics is used to describe macroscopic quantum effects and quasi-particles such as solitons, which have a Hamiltonian which depends nonlinearly on their wavefunction. Typically used in the modelling of quantum semiconductors and superfluids, the theory differs from linear quantum mechanics in many important ways. For example, the square of the wavefunction is no longer the probability of finding the macroscopic particle at a given point in spacetime, but instead gives the mass density of its constituent microscopic particles at that point. Operators are no longer linear, and hence wavefunctions no longer obey the principle of superposition; the unitary structure of linear quantum mechanics does not apply.
In addition to these conceptual differences with standard quantum mechanics, the nonlinear nature of the coupled Schrödinger-Poisson equations (6) and (7) make them difficult to solve. The wave-equation has two sources of nonlinearity: the first from the coupling to the Poisson equation via the gravitational potential, and the second from the nonlinear 'quantum pressure' term (this term resembles a pressure gradient when interpreting quantum phenomena in terms of classical fluid behavior) which acts like an effective potential. It is useful to insert a multiplicative parameter β in the nonlinear quantum pressure term in the wave-equation (6), such that β = 0 corresponds to leaving the term out, and β = 1 to retaining it. One then finds that the corresponding fluid equations remain unaltered, except for the Bernoulli equation (4), which now takes the form ∂φ ∂t The previous studies mentioned thus far, which apply the Schrödinger formalism to the dynamics of a cosmological dark matter fluid, typically discard the troublesome nonlinear quantum pressure term P . From (8), we see that omitting P in the modelling of a dark matter fluid can be justified only if |ψ| = √ ρ varies very slowly on the scales of interest. Alternatively, for fluids with pressure (such as a baryon component), one could adjust (or omit) the P term in the wave equation to model genuine effects of pressure (Coles & Spencer 2002), but this is clearly not relevant to modelling a dark matter fluid. We note from (8), however, that in the limit ν → 0, equation (9) reduces to the correct Bernoulli equation (4). This has led several authors to discard the quantum pressure term P , thereby making the resulting wave-equation linear, and considering the resulting theory only in the 'correspondence limit' ν → 0. Although, in this limit, the fluid equations to which this modified Schrödinger-Poisson system is equivalent are those that hold physically, the limit ν → 0 is nonetheless meaningless in the context of the Schrödinger formulation itself; moreover, this assumption constrains these previous authors to considering the study of phenomena in a certain limit of a parameter in their theory. Short & Coles, for example, tackle this difficulty by performing simulations for the lowest value of ν possible numerically and, unsurprisingly, find their simulations agree best with N -body methods for smallest possible ν. Szapudi & Kaiser's perturbation theory takes a more rigor-ously analytic approach, but it is nonetheless based on a linearised Schrödinger equation, and thus is only valid in the correspondence limit. Furthermore, in developing their free particle approximation, Short & Coles rely heavily on results from linear perturbation theory, which is an entirely different mathematical representation of the physical process of interest. In fact, no study yet exists of the mathematical formulation of structure formation in the Schrödinger form which is not approximative in some way.
In this paper, we therefore take a different approach and study the full nonlinear Schrödinger-Poisson system of equations (6-7). This general type of system has been studied before, but usually in the context of quantum mechanics, and always with Λ = 0. Spiegel (1980) was chiefly interested in discussing its correspondence to the fluid equations. Its convergence limits, well posedness, and semiclassical limits have been studied by Jungel, Mariani & Rial (2001), Jungel & Wang (2002) and Li and Lin (2003). Similar properties of simpler linear Schrödinger-Poisson systems have also been studied by Castella (1991) and by Abdallah, Mehats & Pinaud (2005). The nature of the coupling between the equations (6) and (7) means that the Schrödinger equation has what is effectively a time-dependent potential; solution methods for these have been studied by e.g. Park (2001) and Devoto & Pomorisac (1992). More recently, Erhardt and Zisowsky (2006) have explored numerical solution techniques for a spherically-symmetric Schrödinger-Poisson system describing the time evolution of an electron in a classical polar crystal.
Adopting our more direct approach leads inevitably to more mathematical complication, but we believe this is more than offset by a clearer physical interpretation. In particular, in our approach the value of the parameter ν is unimportant, since it does not appear in our analytic solutions for the wavefunction ψ, which is contrary to previous approaches. Also contrary to previous authors, we will not recast our equations in comoving coordinates, since we feel that this allows for a more straightforward physical interpretation of our results. Finally, it is important to stress that, in our cosmological context, the coupled Schrödinger-Poisson equations (6-7) describe a fully classical system. Although we wish to utilise any possible links between quantum mechanics and the Schrodinger description of cosmological structure formation, caution must be used when drawing parallels with familiar quantum mechanical ideas, such as the uncertainty principle.

COSMOLOGICAL SOLUTIONS
In this section, we discuss the solutions to the Schrödinger-Poisson (SP) system (6-7) in the simple case of a homogeneous and isotropic cosmological dark matter fluid. In this case, the P term in the SP system is in fact identically zero, since |ψ| 2 = α 2 = ρ, as a physical observable, must also be homogeneous. Hence the wave-equation becomes linear without any need for approximations. On the other hand, the velocity potential φ and the gravitational potential V are not physical observables, and so may be nonhomogeneous; we may only assume that they are functions of time t and radial distance r only. Thus, we take the ansatz α = α(t), φ = φ(t, r) in the wavefunction ψ and V = V (t, r). Consequently the Laplacian operator in the SP system is simply Although the nonlinear quantum pressure term drops out of the SP system, the system itself remains nonlinear due to the coupling of the gravitational potential V to the wavefunction which is provided by the Poisson equation. As mentioned above, we will work in 'real' (i.e. non-comoving) coordinates. Ideally, we would like to obtain results for the most general case, i.e. that of a universe with arbitrary spatial curvature and a non-zero cosmological constant, from which all other cases can be derived. It is well-known, however, that such a solution would include elliptic functions, with which it is difficult to work analytically. Consequently, we will not consider the open universe case. Indeed, to maintain algebraic simplicity, we consider only a spatially-flat universe with a non-zero cosmological constant and a closed universe with zero cosmological constant.

Spatially-flat universe
Rather than solving the nonlinear SP system directly, in this case we can obtain the required solution very quickly by simply assuming the well-known cosmological result that the fluid density in this case evolves with (a globally defined) cosmic time t as This can then be substituted into the modified Poisson equation (7) to obtain an equation for the gravitational potential V . The solution of this equation is easily found and can then be substituted into the (in this case linear) wave-equation (6) to obtain an equation for the velocity potential φ, which is also easily solved. The resulting wavefunction and gravitational potential are found to be where λ ≡ 3 2 p Λc 2 /3. The solutions in the special case of a spatially-flat universe with zero cosmological constant are easily obtained from these results by performing a series expansion of the functions ρ = α 2 and φ/r 2 around λ = 0, which read Setting Λ to zero, we thus obtain the following solution to the SP system: This result is easily checked by deriving it afresh by assuming the well-known result ρ ∝ t −2 in a spatially-flat universe with Λ = 0 and following an analogous procedure to that used above. In any case, we note that, at at fixed t, the potential V has the form of a harmonic oscillator potential in r. Also of interest is the fact that the wavefunction itself is similar to the Green's function for the diffusion equation.

Closed universe
To obtain solutions for a closed universe, it is first useful to make a change of coordinate from cosmic time t to a 'conformal time'. In the usual cosmological approach, the dimensionless conformal time is defined as where R is the scale factor for the closed universe and c is the speed of light. In the Newtonian context in which we are working, however, c itself has no special meaning, and so we need to define a more appropriate conformal time variable. Introducing a constant, a, with dimensions of velocity, we can define a new dimensionless conformal time η as where the constant M ≡ (4π/3)ρR 3 . For notational brevity, we also introduce the constant µ ≡ 3M/(4π).
In an analogous manner to our analysis of the spatially-flat case, rather than solving the SP system directly, we begin by assuming the well-known form of the density evolution in a closed FRW model, as a function of conformal time η, namely where A is a constant. Following the same procedure as that used in the spatially-flat case, one finds that the wavefunction and gravitational potential that solve the SP system are together with the following condition on the constants: Indeed, one can recover the spatially-flat case (with zero cosmological constant) by taking the limits a → 0 and A → 0 in such a way that A/a 3 remains finite. It will be convenient to write the solutions (19) and (20) in a more transparently cosmological form. To assist in this endeavour, it is first useful to obtain parametric equations for R and t in terms of the conformal time η. Using the expressions (17) and (18), and the relationship (21) between the constants, one quickly finds the well-known results where we have assumed the initial condition that η = 0 when t = 0. Indeed, from the first result we see that our velocity parameter a = p 2GM/Rmax, where Rmax is the maximum value of R attained during the evolution.
Using the results (22) and (23), together with the (17), (18) and (21), one can easily find the expressions for the Hubble and density parameters to be Both parameters have the correct limits at the 'big bang', η = 0, and at 'turnaround', or the maximal expansion point, η = π.
Using these results, it is straightforward now to express the wavefunction (19) and gravitational potential (20) in terms of standard cosmological variables as It is simple to verify that these results reduce to the flat-space solutions (14) and (15) in the case where Ω = 1 and H = 2/(3t).

General case in terms of cosmological parameters
It is, in fact, possible to obtain expressions for φ and V , in terms of standard cosmological parameters, that are valid for any pressureless cosmological fluid.
For φ, we have the obvious result which follows since ∇φ gives the velocity, which should be Hr in the cosmological case. For V , we can use the Bernoulli equation (4). Using the above expression for φ, we can immediately deduce Since H ≡Ṙ/R, we havė where q = −RR/Ṙ 2 is the usual deceleration parameter. Thus, the expression for V becomes which can be seen to agree with all the individual expressions given above, and is also valid in the Λ = 0 case, for arbitrary spatial curvature. In that case, where Ωm = 8πGρ/(3H 2 ) is the matter density, and ΩΛ = Λc 2 /(3H 2 ).
In the same spirit of expressing results in terms of cosmological parameters, the general result for α is although this is less explicit than the results for V and φ. The overall wavefunction in the general case is thus

MODELLING A SPHERICAL OVERDENSITY
We now use the two cosmological solutions found in Section 3 as the basis for constructing an analytical wave-mechanical description of the so called 'top-hat' spherical overdensity model, commonly used to study nonlinear gravitational collapse (e.g. Peebles 1980; Padmanabhan 1993). The reader should keep in mind, however, that in all that follows we are simply describing ordinary Newtonian fluid evolution, albeit in a novel way. The model system comprises: an infinite outer fluid region of homogenous density; an inner fluid sphere, also of homogenous density, but overdense relative to the outer region; and a spherical shell of vacuum which separates the two.fluids. The inner fluid sphere is usually 'compensated', that is, we can imagine forming the inner fluid region by removing the mass from the vacuum (gap) region and adding it (homogeneously) to the mass within the inner radius of the gap. As a consequence of Gauss' theorem (or equivalently of Birkhoff's theorem in general relativity), and because both the inner and outer fluid regions in this model are homogeneous, each region of fluid effectively behaves as a distinct 'universe', evolving independently of the other according to cosmological solutions such as those found in Section 3. The inner 'island universe' region remains isolated from the outer universe in which it is embedded via the compensation condition, which effectively ensures that the outer radius of the gap region evolves in such a way as to keep the two fluid regions distinct at all times.
In the cosmological context, the outer region is generally taken to be a spatially-flat Friedmann universe, with compensation ensuring that the mean density of the system as a whole is kept at the critical value. The overdense inner region will evolve as a closed universe. It is also usual to demand that both the inner and outer universes have the same origin in time, i.e. that their 'big bangs' are synchronized. This simple spherical overdensity forms the basis of the 'Swiss cheese' model widely used in studies of cosmological structure formation. For simplicity, we will assume throughout that Λ = 0.

Specification of the spherical overdensity model
Our aim is to obtain a global wavefunction and gravitational potential that solve the SP equations (6-7) for this physical system; we do this by constructing a piecewise solution from the cosmological solutions derived in Section 3 and a vacuum solution. An illustration of the physical system is shown in Fig. 1, where (ψc, Vc), (ψg, Vg) and (ψ f , V f ) represent the wave-mechanical solutions to the Schrödinger-Poisson system in the inner (closed), vacuum (gap) and outer (flat) regions respectively. The radii Rc and R f label the inner and outer boundaries of the vacuum region respectively. We will denote by (ψ, V ) the global wavefunction and gravitational potential that we construct, in a piecewise fashion, from the solutions for each region.
Our first task is to establish a global time coordinate for the system. This is, in fact, straightforward, since equation (23) already links the 'ordinary' (cosmic) time coordinate t appearing in the solution for the spatially-flat outer region with the conformal time η in the solution for the closed inner region. This condition establishes the correct phasing between the two regions by setting the conformal time equal to zero at the zero of ordinary time. Since the overdense inner region is compensated precisely by the surrounding vacuum region, our system evolves correctly in that the volume of the vacuum region vanishes as t → 0 (or η → 0). This is easily verified as follows. The density of the flat universe varies as Figure 1. Illustration of the spherical overdensity. The wavefunction and gravitational potential that solve the SP system in the inner, vacuum, and outer regions, respectively, are (ψc, Vc), (ψg, Vg), and (ψ f , V f ). The radii Rc and R f label the boundaries of the vacuum region.
where M = (4π/3)ρcR 3 c is the total mass of the inner fluid region. This expression for R f can be compared to the equation (22) for the evolution of outer radius, Rc, of the closed inner region. One finds that R f /Rc is a function of conformal time η only, and for small η is given by which correctly tends to unity as η → 0, i.e. as the big bang is approached.
In the context of the finite inner fluid region of total mass M evolving as an isolated closed universe, it is worth noting that the role of the velocity parameter a = p GM/Rmax, introduced in (17), now becomes clear: it is the Newtonian escape velocity at the surface of the inner fluid when it reaches its maximum radius. The limit in which the inner fluid region becomes spatially-flat corresponds to the escape velocity a at this point tending to zero. (We also note that a = c would correspond to the maximum radius Rmax being the Schwarzschild radius for a mass M ).

Boundary conditions
To construct piecewise global solutions for the wavefunction ψ and gravitational potential V describing the full system, it is necessary to determine the appropriate boundary conditions on these functions at the two moving surfaces Rc and R f .
We first consider the boundary conditions on the gravitational potential. At each boundary in this highly idealized system, the fluid 'boundary layer' is infinitessimally thin, and effectively massless. Consequently, it cannot support a force acting on it, which in turn implies that ∂V /∂r, and hence also V itself, must be continuous at each boundary. The second derivative of V may, however, have a discontinuity, which will correspond to the gradient of the force being discontinuous at the boundary. This is indeed physically correct, since in the vacuum region there will be tidal forces, but clearly none in the homogeneous regions. The gravitational potential V is determined by the Poisson equation (7), with Λ = 0, which for our spherical geometry becomes A discontinuity in the second derivative of V thus implies that α should also possess a discontinuity at each boundary, as expected.
We turn now to the boundary conditions on the wavefunction ψ. Griffiths (1992) has previously considered the boundary conditions on a time-independent wavefunction in the presence of a potential possessing a discontinuity (modelled as a Heaviside function, i.e. the derivative of a Dirac delta function). In our time dependent case, however, Griffiths' approach will not yield the appropriate boundary conditions. To determine the consequences of the discontinuity in square-root density α at each boundary, we must turn to the equation governing the time development of α, which is For our solution (ψ, V ) for the full system to be consistent, the form for α that we know to be correct at the boundaries must satisfy this equation at all times. Consider, for example, the boundary r = R f (t), where we make a transition from vacuum to homogeneous fluid. The form of α, as r increases, is where f (t) describes the time evolution of α within the outer (flat) region of homogeneous fluid and H is the Heaviside step function. Substituting this form into (38), one finds In order for the postulated form (39) of α to be consistent, the coefficients of the step and delta functions need to vanish separately. This follows immediately for the coefficient of H, since on closer inspection this is seen to be the equation for the time development of α in the homogenous region, where α = f . The coefficient of the delta function needs to be evaluated at the moving boundary r = R f (t), which means that one requires dR f /dt = dφ/dr there. This is, however, exactly the statement that the boundary is moving at the speed which is given by the velocity potential at that point, and so this coefficient vanishes also. A similar analysis shows that this form for α is also acceptable at the other boundary. It now remains only to consider what, if any, conditions must be satisfied by the phase of the wavefunction, the velocity potential φ, at the boundaries. The equation determining the evolution of φ is simply the Bernoulli equation (4), which in this case reads Examining this equation, it becomes clear first of all that it is not necessary to consider the possibility of singularities in φ at the boundaries: since V has continuous zeroth and first derivatives, then φ will also remain singularity free. Additionally, in regions where ψ = 0, i.e. the vacuum region, φ may take any value at all, since it is not defined there. To summarize, we have shown that the only physical boundary conditions which must be satisfied for the solution ψ, V for this system are those of continuity of the gravitational potential V and its first derivative. The wavefunction ψ itself need not be continuous at the boundaries. We have also established that our wave mechanical description of the system is thus far self consistent, by demonstrating that a step discontinuity in the modulus of the wavefunction at both boundaries satisfies the time dependent Schrödinger-Poisson system, as expected.

Forming the piecewise solution
Having established the boundary conditions on the wave mechanical description (ψ, V ) of the spherical overdensity system, we can now consider how we can form this wavefunction from those which separately describe each region.
We first consider the wavefunction ψ. The modulus of the wavefunction, the root-density α, is already specified everywhere, and thus we have no further work to do here. The phase of the wavefunction, the velocity potential φ, is fully specified in the inner and outer fluid regions by (26) and (14) respectively, but it is undefined in the vacuum region since there is no fluid there. Moreover, our analysis above has shown that we do not require continuity of phase of the wavefunction at either boundary.
We turn now to the gravitational potential V . We first note that, since α = 0 in the vacuum region, the gravitational potential in this region must satisfy the Laplace equation, so that where B and C are functions of time only and must be fixed by appropriate boundary conditions. We showed above that, for a physically meaningful solution, we require V and ∂V /∂r to be continuous at each boundary. Consequently, we must match the vacuum solution (42) and its radial derivative at each boundary with the corresponding solution for the inner or outer fluid region. Matching at the inner boundary yields the vacuum solution whereas matching at the outer boundary gives where Rc and R f are functions of time. These two results have the same functional form inside the vacuum, but differ by a constant. This is problematic: since V is defined through the vacuum, these solutions, unlike those for φ, do need to match. This problem is resolved, however, by appealing to the extra degrees of freedom offered to us by the fact that the velocity potential in each fluid region is fixed only up to the addition of an arbitrary function of time. From equation (41) for the time evolution of φ, we see that for the vacuum solutions for V derived at each end of the boundary to match, we must simply adjust one or both of the velocity potentials φ at the two ends of the vacuum by functions of time alone such that their difference is altered by the following Physically, at the level of the Schrödinger equation, this corresponds to a gauge change, which alters the potential by a function of time, by introducing a phase term (the integral of the change in time) into the wavefunction. This phase term is introduced relative to the cosmological solution wavefunction in the relevant region. The solutions for individual regions can now be put together to form a global solution with all boundary conditions satisfied. The solution (26-27) is taken as the global solution (ψ, V ) for r Rc.
In the vacuum region, between Rc and R f , the velocity potential is undefined, and we match V appropriately at the inner boundary so that the gravitational potential (43) holds throughout the gap. We then employ the gauge freedom in φ to adjust the velocity potential in the outer region, r R f , so that the global gravitational potential is continuous. We do this by subtracting (45) from the phase of the wavefunction in (14), which then becomes the solution for φ in the outer region. Equation (41) can then be used to find the corresponding gravitational potential in the outer region, which is correctly matched at the boundaries. The resulting global solution for ψ = αe iφ/ν and V , defined piecewise in α, φ and V for r Rc, Rc < r < R f and r R f respectively, is thus given by It is worth noting the physically interesting result that the 'correction term' in the velocity potential φ in the outer region r R f is proportional to the difference in the conformal time variables in the inner and outer regions. Also of interest is that, for similar physical systems, it is clear that we will not require any particular conditions on the wavefunction at the boundaries. As illustrated by our (somewhat pathological) example, step discontinuities ensure that the wavefunction and its first derivative certainly are not required to be continuous there, although we can imagine situations where they may be; an example of this might be a smoothly varying density of fluid with a boundary represented by an abrupt change in fluid velocity at a given radius. This behaviour is in marked contrast to familiar systems in quantum mechanics, where it is the gravitational (or other) potential which generally is considered to have discontinuities, but the wavefunction must satisfy physical boundary conditions.

Comparison with quantum Schrödinger systems
It is interesting to consider in more detail how our results compare with what may be considered a similar system in standard quantum mechanics. Surprisingly, the problem of which boundary conditions to apply at a moving boundary is a matter of debate in the literature: an example of the confusion is provided by Samura and Ohmukai (2006). The authors consider reflection and transmition of an incident wavefunction at a moving potential step using the standard time-dependent Schrödinger equation.
They claim that the boundary conditions one would expect to apply, namely that the wavefunction and its first derivative are con- tinuous, are not in fact valid. Since the position of the potential step is a function of t, they claim that one cannot use the continuity of the first derivative because in taking this condition, the boundary is fixed in time. Instead, they argue, one should use an alternative boundary condition, that the phase of the wavefunction is continuous at the boundary, and they proceed to use this to derive expressions for the reflected and transmitted waves.
Although this condition yields the correct final results, it is straightforward to show that the authors are incorrect in claiming that the continuity of the first derivative is not a valid boundary condition to apply. One can in fact use the standard boundary conditions to derive the correct results, simply applying a Galilean transformation to the stationary barrier case.
In this approach, we take the usual Schrödinger equations left and right of the boundary, but take the three wavefunctions, ψi, ψr, and ψt as functions of x ′ = x−vt. We examine the case in which, if the step were stationary, ψt would be purely evanescent, so that we can take the limit to an infinite step, which is the case most closely resembling our cosmological Schrödinger system. By applying the usual boundary conditions of continuity of the wavefunction and its first derivative at the boundary, x = vt, one finds the following results for the wavenumbers kr and kt in terms of ki (setting = 1): These yield wavefunctions which satisfy the Schrödinger equations in each region, and they reduce to the correct results in the case where v = 0. A series expansion of the first derivative of the wavefunction in the case that V0 → ∞ shows that there is a finite slope at the boundary, demonstrating that this method works for both the finite and infinite potential step. This ambiguity in the published literature highlights the fact that determining which boundary conditions to apply to Schrödinger wavefunctions, even in the usual quantum mechanical context, is not a trivial matter. It is interesting that in the case of our classical gravitational system, neither the matter wavefunction nor its derivative need be continuous, in contrast to the usual quantum case. Instead, we require that conditions be placed on the gravitational potential which is coupled to the wavefunction.

MULTIPLE FLUIDS
The fluid description of cosmological dark matter breaks down when multistreaming occurs, most notably at shell-crossing. In general, the onset of shell-crossing occurs after the onset of non-linearity; shell-crossing thus represents a natural 'barrier' beyond which fluid-based methods cannot follow the evolution of structure into the fully nonlinear regime. Attempts to describe structure in the quasi-linear regime using a fluid description therefore resort to approximate methods to extend the approach beyond shell-crossing, e.g. the Zeldovich, adhesion, and free-particle approximation methods.
The Schrödinger formalism is also based on a fluid description, and consequently shares these limitations. Nonetheless, as in other fluid-based approaches, one can model more general physical systems by considering them to be made up of multiple fluids. One still cannot accommodate multi-streaming within each individual fluid, but one can model some multi-streaming situations by considering each stream as a separate fluid. In this section, we therefore extend our previous description by moving from a single fluid to the multiple fluid case, inspired by the application of techniques from multiparticle quantum mechanics. We do this by assigning a distinct wavefunction to distinct regions of the CDM fluid, each described by its own wave-equation, which interact collectively via a common coupled Poisson equation. The toy cosmological system illustrated in Fig. 1 can always be successfully described using the Schrödinger formulation described in the previous section because the inner and outer regions of fluid never overlap, i.e. its inner and outer 'shells' never cross. There is thus a single unique fluid velocity for each point in the system, throughout its evolution. This formulation would break down at multistreaming because the Madelung transformation ansatz (5) for the form of the fluid wavefunction assumes a single-valued velocity field. A good toy model for illustrating multi-streaming in cosmology can be constructed, however, by considering a physical system similar to that discussed in Section 4, but with appropriate initial and boundary conditions to ensure that the two fluid regions do overlap at some stage. We now describe how to construct a new, two-fluid description for such a system that accommodates this multi-streaming.

Multiparticle wavefunctions and multiple fluids
Since we have previously adopted a wavefunction approach to modelling a single fluid, it is natural to consider whether multiparticle methods in quantum mechanics can assist in providing a wavefunction approach to multiple fluids. For simplicity, we will confine our attention to just two fluids.
Several questions immediately arise, such as what is the correct form of the quantum pressure and potential terms that should be used in the 'two-particle' case, and whether 'entangled' states have any role to play in the classical evolution of two coupled fluids, or whether factorisable states are the only ones relevant.
In the single-fluid case, the approach we have been using is strongly related to the de Broglie-Bohm formulation of quantum mechanics (for a full account see, for example, Holland 1989). It is natural, therefore, to turn first to the many-particle version of de Broglie-Bohm theory. We now wish to adapt this approach to obtain an equation that describes the evolution of two fluids, i.e. describes their velocities and particle densities in time and space, using a Schrödinger formalism. Of course, we also need to show how one can obtain individual densities and velocities for each of the fluids from the resulting equation.
In the two-fluid case, the overall wavefunction now has the form ψ = ψ(r1, r2, t), where r1 and r2 are spatial positions for an element within fluid 1 or fluid 2 respectively, and t is the overall common time variable; we have thus introduced a (6 + 1)-dimensional parameter space. As for the single fluid case, we spilt ψ into a magnitude and phase written as We also introduce the gradient operators ∇i (i = 1, 2) for each fluid, where we adopt the convention that bold quantities 'live' in the individual fluid spaces, and non-bold quantities refer to the total space (with 6 'spatial' dimensions). Using Cartesian coordinates ri = (xi, yi, zi) in each fluid space, for example, the corresponding Laplacian operators (for i = 1, 2) are We now assert that the overall two-fluid wave-equation is where V1,2 is the common gravitational potential and the quantum pressure term has the form This expression for P is certainly the natural generalisation to the two-fluid case of the equivalent P -term for a single fluid. It also agrees with the 'many-body quantum potential' Q given in equation (7.1.4) of Holland (1989), although our quantum pressure term is the negative of this, since here we are seeking to remove the term that would otherwise appear in the fluid equations of motion. A useful quantity relating fluid density and velocity in de Broglie-Bohm theory is the probability current. The natural generalisation of the probability current to two fluids, is where ∇ = ∇1 + ∇2. The interpretation of J as a current in six dimensions, corresponding to a flow of 'probability' into and out of regions, just as for the single-fluid case, follows from which is easy to prove by using the two-fluid wave-equation (54). Using the split (52) into amplitude and phase terms, we find and then identifying the individual space velocities we have The conservation equation (57) can now be written The other equation we require is that relating the common gravitational potential V1,2 to the densities of the two fluids sources. The term V1,2 (for which we effectively include a copy for each fluid in the wave-equation) is taken to satisfy the (single-fluid) equation i.e. the potential under which both fluids move is generated by the sum of their densities. The wave-equation (54) and the Poisson equation (62) give what seems a plausible 'two-particle' quantum description of two interacting classical fluids. The most basic test that any many-body equation must satisfy is that it works for a factored state, for which ψ(r1, r2, t) = ψ1(r1, t)ψ2(r2, t).
Inserting this form of ψ into the wave-equation (54), dividing through by ψ1ψ2, and applying the usual separation of variables argument (here to the coordinates r1 and r2), we obtain the two individual fluid equations ψ2.
The wave-equation (54) thus makes sense for factored states. However, generalizing to 'entangled states' appears problematic. This can be seen immediately in equation (62) for the gravitational potential, which requires knowledge of (the sum of) the separate fluid densities ρ1 and ρ2. Given a general two-fluid wavefunction ψ, it is not clear how these would be obtained. Given ψ, we can successfully form the two separate particle velocities, via the definitions (59), which only require that the total phase, φ, be available. However, knowledge of the overall magnitude, α 2 = |ψ| 2 , is not enough to deduce the separate densities. Thus, in the general entangled case, we have no way to generate ρ1 and ρ2 separately from ψ, and furthermore, it appears that we cannot find their sum either (to put into the right hand side of (62)), so that the equation set becomes incomplete as it stands. A possible solution to this problem is to assume continuity in the two spaces separately, which seems a reasonable physical requirement. Supposing that ∂ρ1 ∂t + ∇·(ρ1v1) = 0 and ∂ρ2 ∂t + ∇·(ρ2v2) = 0, (65) and given initial density distributions ρ1 and ρ2, there is enough information available to propagate the two densities individually forward in time, since v1 and v2 are separately available. Thus there is in principle an answer available, although it would not be a useful one numerically. We note further that the intuitive identification α 2 = ρ1ρ2 is not even necessarily valid in the non-factored case, since the 6dimensional continuity equation (60), plus the assumption of continuity for each fluid separately, can be used to show that where the d/dt represents a comoving or streamline derivative in the 7-dimensional space: In the factored case, the right-hand side of (66) will vanish, meaning that the identification of α 2 with ρ1ρ2 is maintained along the (joint-fluid) motion. However, in the non-factored case, there is no guarantee that right-hand side of (66) will indeed vanish, and we have a further problem in interpreting the wavefunction.
To summarize our answers to the questions posed at the beginning of this section, we now have a good two-fluid wave-equation description of the motion of two coupled fluids, for which factored states make sense and return sensible equations in each particle space (with the coupling still present due to the potential, but not explicitly at the coordinate level). The notion of non-factorisable states looks problematic, but given that we are attempting to describe classical, rather than quantum, fluids, this is perhaps unsurprising.

Numerical two-fluid evolution with spherical symmetry
We now apply the two-fluid wave-equation formalism developed above to model explicitly the evolution of a spherically symmetric system consisting of two physically identical, mutually selfgravitating fluids, each of which is described by individual wavefunctions ψ1 and ψ2 respectively. As previously, the wavefunctions are given by where α1(r, t) and φ1(r, t) and α2(r, t) and φ2(r, t) are the square-root densities and velocity potentials for each fluid. The equations that these wavefunctions satisfy are given by (64), which must, in general, be solved numerically.
The simplest way to solve the wave-equations numerically is to take their real and imaginary parts, and to solve for the amplitude and phase of each wavefunction. We therefore solve the following system of five coupled equations for the evolution of the two fluids: It is worth noting that equation (74) for the gravitational potential V does not contain a time derivative and hence, for fixed t, it is, in fact, an ordinary differential equation in r. Moreover, we note that the parameter ν does not appear in the above system of equations. To solve the system of equations (70-74) numerically for the five variables αi(r, t), φi(r, t) (i = 1, 2) and V (r, t), one must first specify the appropriate initial and boundary conditions on the solution. The initial conditions at some time t = 0 (say) are assigned by specifying αi(r, 0) and φi(r, 0) (i = 1, 2). At any given t (including t = 0), knowledge of the αi is sufficient to determine the gravitational potential V by solving (74). In this sphericallysymmetric case, the general solution of (74), at a given time, is where c1 and c2 are arbitrary functions of time. Clearly c2 is merely an overall additive term which fixes the zero of gravitational potential, and c1 expresses whether there is a point mass at the origin. We choose c1 = c2 = 0 to fix V in a simple manner. Furthermore, we will assume that the total fluid density is non-zero at the origin. This means that the dominant behavior of the inner integrand is ∝ r 2 , and thus that the dominant behavior of the outer integrand is ∝ r.
With these assumptions, one can perform the double integration in (75) using a set of two interleaved Simpson's rules, which we have found to work to high accuracy. It remains to consider the boundary conditions on αi and φi at each end of the spatial domain, i.e. at r = 0 and r = rmax. We first focus on the conditions that must apply at the origin. There is a physical requirement at r = 0 that, in the absence of a singularity there, it should be neither a source nor a sink. Since, as mentioned above, we are assuming the total fluid density α 2 1 + α 2 2 > 0 at the origin, we thus require ∂φ1/∂r = 0 and ∂φ2/∂r = 0 (i.e. zero fluid velocities) there. There is also an arbitrary additive constant in each φi to fix and we do this by requiring that φi = 0 at the origin. These two boundary conditions for each φi are compatible with the equations of motion, since ∂φi/∂r = 0 and V = 0 at r = 0 combine to mean that ∂φi/∂t is also zero there, so φi = 0 will be maintained. A further detail that has to be addressed at r = 0 is how to evolve each αi forward in time given that 1/r appears in one of the terms of ∂αi/∂t. Since φi and ∂φi/∂r are both zero at r = 0 we can assume that the dominant behavior of φi is ∝ r 2 . In this case, we can combine the last two terms of ∂αi/∂t to obtain for each fluid If φi is in fact proportional to a higher power of r, then this formula is still correct, since it then returns zero. It is worth noting in the above prescriptions, there is no requirement that the αi have zero gradient at r = 0. Thus a cusp in the αi distributions at the origin is physically allowed, and does not mean that there is a singularity there. We now consider the boundary conditions we must apply at the outer edge r = rmax of the spatial domain. It is, in fact, possible to evolve the equations (70-74) numerically without fixing boundary conditions for αi and φi at the outer end point of the spatial grid, provided we take r = rmax to be sufficiently large that the total fluid density is always zero there. However, any increase in the size of the spatial domain comes at a heavy computational cost for a fixed spatial resolution, so we choose to apply suitable boundary conditions at the outer end point of the spatial grid. To do so, we use the 4 rightmost end points of the grid to calculate third-order spatial derivatives of αi and φi at r = rmax and thereby predict to third-order in r the values of these variables at an 'extra' endpoint, one grid step beyond the desired boundary. Since such a scheme is accurate to third order, it is sufficient for the consistent solution of our system of second-order PDEs.
Having specified the above initial and boundary conditions on the solution, the equations (70-74) can then be integrated numerically. The most straightforward numerical approach is to use a simple first-order time-stepping technique in which the initial distributions αi(r, 0) and φi(r, 0) are 'marched forward': at each time-step t → t + ∆t, we calculate the required spatial partial derivatives on the right-hand sides of (70-73) by taking differences and then add ∆t × ∂αi/∂t and ∆t × ∂φi/∂t to αi and φi (respectively) at each grid point. It is well-known, however, that such a scheme is prone to numerical instability (a fact that we unfortunately verified when we implemented it). We therefore instead employed a fourth-order Runge-Kutta integration technique (see e.g. Press et al. 2007), which calculates spatial derivatives (again by simple differencing) at several points throughout the time interval through which the solution is being advanced. This offers a significant in-crease in accuracy over a standard Euler, or simple finite differencing, technique, which calculates spatial derivatives only at the beginning of the time interval.

Toy physical system
To demonstrate that our approach can indeed be used to describe the evolution of a (spherically-symmetric) two-fluid system after the fluids overlap, we apply it to a simple physical system similar to the spherical overdensity model considered analytically in Section 4. In that case, however, the inner and outer fluid regions did not overlap, the outer fluid region was of infinite spatial extent, and the density profiles possessed sharp step-function discontinuties. We consider here instead a modified physical system, as follows. First, we consider both inner and outer fluid regions to be of finite spatial extent. Second, we specify initial density and velocity distributions for both fluids that are consistent with flat-universe evolution, since it is straightforward to construct such a system in which overlap of fluids is guaranteed to occur. Finally, since our numerical solution scheme relies on obtaining accurate spatial derivatives, it is clearly feasible only for initial density profiles that are smoothed to some degree. Nonetheless, the smoothing is sufficiently slight that we can compare the results we obtain numerically to known analytic results for the evolution of an annulus of a flat universe, as found in Section 4.
The assumed initial αi and φi distributions for each fluid are shown in Fig. 3, together with the resulting initial gravitational potential V derived from (75) using a single application of the Simpson's rule method described above. The initial αi profiles are step functions smoothed with a relatively narrow tanh function. The initial profile φ1 for the inner fluid is that of a flat universe, i.e. quadratic in r, and the velocity potential for the outer fluid was taken to be φ2(r, 0) = 0.01 − 0.1r 2 . We set the outer boundary at rmax = 5, and used 401 grid points to sample the spatial domain. The equations were evolved over 651 time steps, with a total integration time of 3.5. It was checked that the ratio of spatial and time resolutions ∆x/∆t remained greater than or equal to the greatest fluid velocity in the system throughout the entire integration time.
It is worth mentioning how the accuracy of the numerical integration is affected by the choice of the initial conditions. As expected, very sharp edges on the αi profiles lead to inaccuracies in the calculations spatial derivatives via finite differencing when the 'width' of the profile edges is close to the spatial grid-spacing; these inaccuracies can then propagate into the rest of the spatial domain. Also, if either fluid velocity is too fast at any point, such that it exceeds ∆x/∆t, then inaccuracies again propagate through the numerical solutions. This may occur directly by setting the initial velocity too high, but high velocities can also occur later in the evolution if the initial density of either fluid is set very large, resulting in faster gravitational collapse.
Given the initial conditions illustrated in Fig. 3, the resulting αi and φi distributions for each fluid after evolving the system to t = 3.5, together with resulting gravitational potential V , are shown in Fig. 4. We see that, as the two fluids overlap, new features develop in the profiles. In particular, we note a pronounced caustic in fluid 2 at a finite distance from the origin, which signals that fluid is 'piling up' at this radius as a result of shell-crossing. Fluid 1 also possess a corresponding 'hump' at radii larger than the position of the caustic in fluid 2.
Finally, to check our code, we also considered a model in which the initial density distribution of fluid 1 was set to zero everywhere, but the initial density of fluid 2 has the same profile as in t = 0 t = 0 t = 0 r r r Figure 3. The initial square-root densities, α 1 and α 2 (left), initial velocity potentials, φ 1 and φ 2 (middle) and the initial gravitational potential V (right) for the numerical integration of the spherically-symmetric two-fluid system. The solid lines refer to the fluid 1 and the dashed lines to fluid 2. The outer fluid (fluid 2) is given a small inwards initial velocity, as shown. t = 3.5 t = 3.5 t = 3.5 r r r Figure 4. As for Fig. 3, but after evolving the system to a time t = 3.5. Fig. 3 (left panel) and its velocity potential was taken to be that of a spatially-flat universe. In this case, aside from minor differences resulting from the smooth edges of the density distribution of fluid 2, we expect this fluid to evolve as (an annulus of) a spatially-flat universe (with Λ = 0), for which the analytic solution is given by (14) and (15). Our numerical results were found to agree well with this analytic solution. This example also demonstrates that, because of the (slight) smoothing of the initial density profile of fluid 2, the total density at the origin r = 0 was non-zero to machine precision, thereby satisfying our earlier assumption in solving equation (75) for the gravitational potential. Hence our scheme can still be used in many cases to evolve systems which 'appear' to have a hole at the centre.

CONCLUSIONS
We have considered the use of the Schrödinger formalism for fluid dynamics to model the evolution of cosmological dark matter. Our approach differs from previous work in this area in that we consider the full non-linear Schrödinger-Poisson (SP) system, without making any approximations. In particular, we retain the 'quantum pressure' term, which previous authors have discarded to regain a linear wave-equation. By so doing previous authors have had to work in the so-called correspondence limit ν → 0, where ν is a parameter that changes the spatial and velocity resolution of the solutions.
We have shown that cosmological solutions to the full SP system can be readily obtained for closed and spatially-flat universes. Moreover, these solutions can be combined to obtain a global solution to the SP system for a spherical overdensity. This necessitated determining the appropriate boundary conditions on the solutions for each homogeneous region of the model, which differ from those usually assumed in quantum mechanics. We also highlight that the issue of what conditions to apply to wavefunctions at a moving boundary remains a matter of debate in the quantum mechanics literature.
We consider how the Schrödinger formalism for classical fluids can be extended to multiple fluids and derive an appropriate wave-equation for the two-fluid case, which could be easily generalised to more fluids. In particular, we find that for consistency of the equations the multi-fluid wavefunction must be factorisable into wavefunctions for each fluid separately; 'entangled' states appear to be problematic. Nonetheless, our approach provides a new and consistent method for dealing with multiple fluids. We illustrate our method for describing multi-streaming fluids by applying it to the evolution of a modified spherical overdensity toy model, in which the two fluids eventually overlap. We integrate the SP equations numerically to obtain physically sensible results, which illustrate the formation of cusps.
We consider our investigations to show that the full non-linear SP system has the potential to be a useful tool for theoretical investigations of cosmological fluid evolution. The non-linearity of the equations can be troublesome for analytical investigations (although not always), but poses few problems for following the fluid evolution numerically. In any case, we believe these difficulties are offset by the advantages it offers in ensuring that the fluid density is always positive. In future work, we plan to investigate linear perturbation methods within the Schrödinger formalism to determine whether such an approach offers any advantages over traditional first-order perturbation methods applied to the standard fluid equations in following cosmological structure formation in the linear regime.