Evolution of Cosmic Voids in the Schrodinger-Poisson Formalism

We investigate the evolution of cosmic voids in the Schrodinger Poisson formalism, finding wave mechanical solutions for the dynamics in a standard cosmological background with appropriate boundary conditions. We compare the results in this model to those obtained using the Zel'dovich approximation. We discuss the advantages of studying voids in general and the advantages of the Schrodinger Poisson description over other approaches. In particular, emphasizing the utility of the free particle approximation. We also discuss a dimensionless number, similar to the Reynolds number, for this system which allows our void solutions to be scaled to systems of different physical dimensions.


INTRODUCTION
Over the past few decades, galaxy redshift surveys, exemplified by those resulting from the Sloan Digital Sky Survey described by York et al. (2000), Eisenstein et al. (2011) and Blanton et al. (2017), have revealed that the present-day Universe displays a rich hierarchical pattern of structure -the "Cosmic Web" -that encompasses a vast range of length scales, from individual galaxies to clusters and filaments surrounding enormous voids. On the other hand, observations of the cosmic microwave background, most recently from Planck Collaboration (2020), suggest that the early Universe was almost homogeneous, with only slight temperature fluctuations seen in the cosmic microwave background radiation. Models of structure formation link these observations through the effect of gravity, relying on the fact that small initially over-dense regions accrete additional matter as the universe expands (a mechanism known as gravitational instability). The growth of density perturbations via gravitational instability is well understood in the linear regime, but the non-linear regime is much more complicated and generally not amenable to analytic solution. Numerical N-body simulations such as those discussed in Springel et al. (2005) have led the way towards a fuller understanding of the structure formation process, but although such calculations have been priceless in establishing quantitative predictions of the large-scale structure expected to arise in a particular cosmology, it remains important to develop as full an understanding as possible using analytical models. This is not only because simulation is not the same as comprehension but also because N -body methods are computational expensive which makes it difficult to explore a large parameter space using them alone.
To complement the "brute force" approach of purely computational techniques, a variety of analytic techniques has been developed for describing various aspects of the evolution of density perturbations to various levels of approximation; for a good review of many of these approaches,see Bernardeau et al. (2002). None of these capture all the details of the structure formation process, but each can provide useful information in particular circumstances. * Aoibhinn. Gallagher@mu.ie In this paper we study a particular formulation of the gravitational instability problem -the Schrödinger-Poisson (SP) approach (explained below) -which is based on wave mechanics rather than fluid mechanics. This approach was first proposed in the cosmological context by Widrow & Kaiser (1993), and has subsequently generated a great deal of interest; see e.g. Coles (2002), Szapudi & Kaiser (2003), Coles & Spencer (2003 Short & Coles (2006b) and Johnston et al. (2010). It turns out that this approach, while initially developed as an approximate method for handling the evolution of cold dark matter, is also applicable to a scenarios in which the behaviour of dark matter is inherently wave-like rather than particle-like, such as if the dark matter is a very light bosonic particle, perhaps in the form of a condensate, or some other form of "fuzzy" dark matter; see for example Brook & Coles (2009), Schwabe et al. (2020, Dentler et al. (2021), Hui (2021) and references therein.
The behaviour of cosmic voids provides an interesting test case for the SP system that complements the case of spherical collapse studied by Johnston et al. (2010). While gravitational collapse tends to amplify asymmetry, Icke (1984) proved a "Bubble Theorem" according to which an initially asymmetric void becomes more spherical as the Universe evolves; see also numerical computations by Bertschinger (1985). The study of a spherically symmetric void is therefore arguably more realistic than the collapse of a spherical over-density. Moreover, because matter tends to become compressed in sheets or filaments as the Cosmic Web develops, voids actually dominate space in terms of volume fraction. This is true for isolated voids, which is what we have done in this paper. As well as these physical aspects of void evolution there is an important statistical property shown by White (1979) that voids contain information about the entire hierarchy of n-point correlation functions at all orders. For these reasons and others there has been a great deal of recent interest in the behaviour of cosmic voids, see e.g. Sheth & van de Weygaert (2004), Bos et al. (2012), Achitouv et al. (2015), Demchenko et al. (2016) and Pisani et al. (2020).
Pulling together these chains of thought, in this paper we apply the Schrödinger-Poisson system to the case of isolated spherical voids. Although of course in a realistic context void expansion would be hindered by neighbouring structures, it is important to obtain as full an understanding as possible of the behaviour of the simplest case within the SP system. The layout of the paper is as follows: in Section 2, we outline the background cosmology needed to build the results of this paper; in Section 3, we outline the basics of the SP system; in Sections 4 and 5, we look directly at the void and apply what was covered in previous sections to solve for the dynamics of the void. In Section 6 we discuss scaling properties of the solutions in terms of an effective Reynolds number. In Section 7 we present our conclusions and suggest future applications of our results and the approach that leads to them.

BACKGROUND COSMOLOGY
The current standard model of cosmology assumes the Cosmological Principle, according to which the universe is assumed to be homogeneous and isotropic, at least on large scales. Space-times consistent with this are described by the Robertson-Walker metric: (1) where κ is the spatial curvature: κ = 0 represents flat space; κ = +1 represents constant positively curved space (closed universe); and κ = −1 represents constant negatively curved space (open universe). The time coordinate t is cosmological proper time and a(t) is the cosmic scale factor. The time evolution of the cosmic scale factor, a(t), is determined via Einstein's gravitational field equations through the Friedman equation, and the density-pressure relation, where the dots denote derivatives with respect to cosmological proper time t, thus describing the global expansion or contraction of the universe. These models can be further parameterised by the Hubble parameter H =ȧ/a and the density parameter Ω = 8πGρ/3H 2 . As usual the present epoch is defined by t = t 0 , when H = H 0 and Ω = Ω 0 . Throughout this paper we assume Newtonian gravity, since the scale of the perturbations are much smaller than the effective cosmological horizon d H = c/H. We also do all calculations in co-moving coordinates, where x is the co-moving spatial coordinate (x is fixed in the frame of the Hubble expansion), r is the spatial coordinate in our reference frame, and a(t) is the cosmic scale factor.
2.1. The Zel'dovich Approximation A simple and elegant approximation to the non-linear stage of gravitational evolution is the Zel'dovich approximation, introduced by Zel'Dovich (1970). The initial density distribution is considered homogeneous and collisionless. The initial (Lagrangian) coordinates are q, then the Eulerian coordinates at time t are given by where a(t) is the cosmic scale factor, b(t) describes the evolution of a perturbation in the linear regime, and is defined byb In a flat matter-dominated universe b ∝ t 2/3 . s(q) here is the displacement term, described by the initial velocity potential φ 0 by, ∇φ 0 (q) = s(q).
Mass conservation requires ρ(r, t)dr = ρ 0 dq, so the density field as a function of Lagrangian coordinates is where ρ = (a 0 /a) 3 ρ 0 . Here ∂r ∂q is the deformation tensor, which accounts for the gravitational evolution of the density field. Due to the nature of s(q), the deformation tensor is a real symmetric matrix with eigenvectors that define a set of three principal (orthogonal) axes. Rewriting Equation (9) in terms of the eigenvalues −α(q), −β(q) and −γ(q), we get (10) If we use comoving coordinates (Equation (5)) in 1 dimension, we get the much simpler expression for ρ The Zel'dovich approximation is an excellent tool for studying structure formation and is highly accurate up to the point of shell crossing Shandarin & Zel'dovich (1989), at which point it breaks down. To combat this we need a model that goes beyond shell crossing, this is one reason why we look to the SP model. Other methods for studying multi-stream flow can be seen in Buchert & Domínguez (1998) and Gough & Uhlemann (2022).

SCHRÖDINGER POISSON REGIME
We start with the Newtonian equations for a selfgravitating perfect fluid. The Euler Equation the Continuity Equation and the Poisson Equation Here v, p and ρ are the (peculiar) fluid velocity, pressure and density, respectively, and V here is the gravitational potential. We simplify our treatment by assuming a cold (pressureless), irrotational fluid moving under the influence of a general potential V . This means that Equation (12) now becomes the much simpler Bernoulli equation where v = ∇φ, where φ is the velocity potential. We now make the Madelung Transformation from Equations (15) and (13). Equation (16) is Schrödinger Equation, with potential V , ν acting as with the addition of non-linear term P . Although ψ is governed by the same equation as the evolution of a single-particle wave function, that is not how it should be interpreted. In particular, |ψ| 2 represents a physical density not a probability density, and its evolution is completely unitary -there is nothing like the wave-function collapse that occurs in standard quantum mechanics in this system.
It's important to note a crucial advantage of this description, namely that because ρ = |ψ| 2 the condition that ρ ≥ 0 is automatically enforced if one applies, e.g., perturbation theory to ψ. This is not the case for approaches based on standard Eulerian perturbation theory applied to δ = (ρ − ρ 0 )/ρ 0 which predict ρ < 0 when δ < −1 and are therefore very unsuitable for describing voids. Note also that because the wavefunction describes a delocalized particle there are no singularities analogous to the caustics that form in the Zel'dovich approximation.
In the cosmological context we take V to be the gravitational potential determined via the Poisson Equation and we get a system of coupled Partial Differential Equations (PDEs). We set ν to be an adjustable parameter so we have some control over the frequency of any oscillatory solutions that could arise from this system. ν has dimensions such that φ/ν is dimensionless, φ is a velocity potential so has dimensions L 2 T −1 , these are also the dimensions of kinematic viscosity, so ν can be thought of as a viscosity parameter for the model; we return to this later.
The original intentions of this model were to find a fluid interpretation of a quantum system, Widrow & Kaiser (1993) first proposed this be used to simulate cold dark matter (CDM). Widrow & Kaiser (1993) start with just the gravitational potential and not yet the fluid approach. It has been used predominantly in a cosmological context since. The quantum nature of this model is also an advantage as it gives us a new perspective on the behaviour of large scale structure. This is especially useful when one seeks to look beyond shell crossing, as this model can handle multiple streams.

COSMOLOGICAL SOLUTIONS
Analytic solutions to the SP system can be obtained straightforwardly for the simplest case of a homogeneous and isotropic fluid. For this P ≡ 0, since |ψ| 2 = α 2 = ρ is defined to be homogeneous. We can start our study of a void by using solutions for an open and flat universe for the interior and exterior of the void in order to make a snapshot of a simple void configuration as follows.
4.1. Spatially flat universe For the following we assume a flat universe (with κ = 0). We begin here by taking the well known cosmological result for the evolution of the density field under a globally defined cosmic time, t, given by Under spherical symmetry the Laplacian operator is simply where r here is now the radial variable r = |x|. This density evolution obtained previously can be substituted into Equation (14) to obtain an equation for the gravitational potential V . This allows the calculation of φ by substitution into Equation (16), resulting in the wavefunction ψ, as shown in Johnston et al. (2010):  (20) where λ ≡ 3 2 Λc 2 /3. We therefore also find φ = Λc 2 12 coth(λt)r 2 .
Also following Johnston et al. (2010), a series expansion of ρ = α 2 and φ/r 2 about Λ = 0 provides the well known results of ρ ∝ t 2 and V = r 2 /9t 2 for this particular case.
4.2. Open universe Working in a similar way to Section 4.1, we find solutions for a universe with negative curvature, κ = −1, an open universe, to be where A ≡ 4πGρ 0 /3, C = γ 1/2 /4A, µ = 3M/4π and γ is a constant. For a simpler derivation we have converted to conformal time by using adη = dt where a is the scale factor of the open universe.
Having established this simple solution we can now proceed to study the dynamics of a void, taking into account the fact that the density surrounding the void will not remain homogeneous as the interior expands.

MODELLING AN UNDERDENSITY
In this Section we look at simple cases of a void growing under the gravity from its boundary. We will use the solutions found in Section 4 to form an analytic model of a void, in analogous fashion to the spherical collapse model studied by Johnston et al. (2010). We will also discuss the results of the numerical solutions to the SP system of equations for a simple one dimensional void model, both 'free-particle' and including the gravitational potential in the evolution. These solutions will be compared to the Zel'dovich approximation for the same initial conditions, which is exact up to the moment of shell crossing, Shandarin & Zel'dovich (1989). We see that shell crossing happens rather quickly in the case of a void though the formation of multi-stream regions forming at the edge of the void. This means that the Zel'dovich approximation applies for a shorter portion of the evolution in the void case than in the collapse case.

Analytic solutions to the void problem
Using the results from Section 4 and following the same method as described by Johnston et al. (2010), we can construct a simple model of a void. We construct the model of the void with an inner region and an outer fluid separated by a shell. The inner fluid is modeled by the open universe solutions and the outer fluid is modeled by the flat universe solutions. Both fluids are homogeneous but the fluid in the inner region is less dense than the outer fluid so expands more quickly, pushing into the surrounding medium and evacuating the void. This creates a small region on the cusp where the density is higher than the outer region. An illustration of how this would look is shown in Figure 1. This figure shows the inside of the void in white, this is the area that gets emptied out. The yellow area is the shell of the void, this is where the denser region is created at the edge of the void. The pink area is the the rest of the universe which extends to infinity.
Using the cosmological solutions shown in Section 4, and following the same procedure as Johnston et al. (2010), we compute an expression for the expansion of the radii of both sections of the void, in conformal time η. We have R o for the open-universe section of the void and R f for flat universe section of the void. and where B is a free parameter corresponding to the initial conformal time η 0 when the configuration is set up, at which point R o = R f . The evolution is illustrated by Figure 2 which shows the evolution of each of these radii  with respect to conformal time with an arbitrary scaling on the vertical axis just to show that the inner region of the void expands much more quickly than the outer region, and an overlapping region develops after some time as the inner void expands into a region that was previously in the space exterior to the void.

The 'free particle' void
In this Section we will present numerical solutions of the void evolution in one dimension using the SP system. We start by studying the approximate case of a void where the gravitational potential is calculated at each time-step, but not included in the Schrödinger equation for the evolution of the void. In other words, the gravitational potential is updated as the fluid moves but the effect of this evolution on the fluid motion is ignored. This has the effect that the fluid moves as a collection of 'free particles', a shortcut that has been shown to be remarkably accurate for some applications; see Short & Coles (2006a). The Schrödinger -Poisson system in this approach simplifies to These equations thus decouple, making the computation much simpler. We begin with a compensated void, shown in Figure 3. A compensated void, is such that the mean density of the density field ρ is set to be 1. To get this shape we use a smoothed tanh function. We need to set up an initial velocity field in which the fluid particles are moving away from the centre of the void, quickest at the centre and more slowly towards the edges of the void. This behaviour can be seen in Figure 4a. To get φ, as shown in Figure 4b, we solve ∇φ = v. This solution for φ is only determined up to a constant, which we fix by setting the potential at infinity to be zero (the edge of the box). These initial conditions were chosen to match the phenomenology from the previous solution in Section 4.
We performed these calculations using periodic boundary conditions. These boundary conditions work quite well with these initial conditions as we have constructed them in such away that the boundary is far enough away from any structure so do not generate any significant cross-boundary interaction as the boundaries are effectively infinitely far away. The Schrödinger-Poisson solution is seen in Figure 5; this solution is analogous to the Zel'dovich approximation if we work in comoving coordinates x, defined by Equation (5), and time coordinate b, as defined by Equation (7), so these two solutions are comparable.
The Zel'dovich approximation is calculated using Equation (11), for ρ 0 in Figure 3 and s(q) = ∇φ 0 (q), where φ 0 is that in Figure 4b. Shell-crossing occurs shortly after the second plot in Figure 5, where multistreaming begins and the fluid density field is no longer well defined under Zel'dovich, as it cannot account for multi-streaming. Examining Figure 5 we see the evolution of the void at three different time-steps in the Schrödinger-Poisson and two different time-steps in the Zel'dovich approximation. At early time-step, we see the void begins to expand slightly and a small amount of matter accumulates at the boundary of the void, as seen in Section 5.1. Just before shell-crossing we see that the two models deviate, with Zel'dovich void boundaries not expanding as rapidly, and the SP void is further expanded. In the final time-step, we see that the inter-  ference effect caused by multi-streaming is in full effect as the void continues to expand. This interference effect occurs after shell crossing, and is also seen in gravitational collapse described by SP dynamics in, e.g., Coles & Spencer (2003) Johnston et al. (2010) and Uhlemann et al. (2019). We think the best way to see what is happening with this model is to watch a movie of this, where the density field changes with time, see ancillary files, which contains this movie.

Void Expansion Speed
It can be seen in Figure 5 that the edges of the void move outwards in comoving coordinates as time goes on. Since the peaks in ρ are sharply defined at least until the strong interference develops after shell-crossing, their positions can be used to define an expansion speed that can be extracted easily from the numerical calculations. Doing this, we found the group velocity of this peak to be generally constant in time (as defined by b), with some jumps once the interference pattern starts to develop. This can be seen in Figure 6, which shows the position of the highest peak at the cusp of the void changes with respect to time. The resulting plot is roughly linear in the coordinates used, which indicates a free expansion  of the void region discussed by e.g. Bertschinger (1985) and Coles & Barrow (1990)). At very late times one expects a void embedded in a baryonic fluid to enter an adiabatic phase described by the Sedov solution for blast waves (Sedov (1946)), but the physical behaviour of the quantum fluid described here would probably be quite different and in any case we only look at the initial phase of the expansion here.

Including the potential term
In this section we look at including a static potential and a time-varying potential term in the Schrödinger-Poisson system of equations. For the static potential, the potential is calculated by Equation (29) from the initial conditions, and then used as a time-independent constant in For the time-varying potential, both Equations (28) and (29) are solved simultaneously as a coupled system of PDEs for a time-dependent ψ and V . Since the evolution of the gravitational potential is quite small, we did not expect either of these extensions to have a drastic effect on the results, but we include them here for completeness. The effect of the potential on the system is only seen at late time-steps and only has an effect on the interference pattern seen after shell crossing. These effects can be seen by comparing Figures 5 and 7. The initial conditions show in Figures 3 and 4 were used for these calculations also, so they can be compared. Comparing Figures 5 and 7, we don't see a large difference between the images, although it is worth noting that for the same time-step in Figure 7 we see the interference pattern coming into effect sooner than in Figure 5. This is expected, as we expect the inclusion of the evolution of the potential to speed up the gravitational effects of the system. It is also worth noting that the interference pattern is almost identical between the static and time-varying potential, 7. This shows us that the potential does not change much over this time-scale which Fig. 7.-Comparing the solution when we include a static potential (dashed orange) and a time dependent potential (solid blue). These plots are the same time-steps as those given in Figure 5 so that we can compare in turn validates the use of the static potential or even the free-particle approximation. We do not see a notable difference in the behaviour of the outward peak velocity, when the potential is included.

VISCOSITY IN A SCHRÖDINGER FLUID
It is interesting to remark that the parameter ν involved in the SP formalism has the same dimensions as a viscosity. This was discussed for example in Short & Coles (2006b) in connection with the adhesion model, an extension of the Zel'dovich approximation; see Gurbatov et al. (1989).
In order to discuss kinematic viscosity in the SP system in more detail one should compare it not with the Euler Equation (12) but with a Navier-Stokes Equation in which case the kinematic viscosity η arises naturally as a consequence of the quantum potential term discussed above in Section 3, though its behaviour is however somewhat different from that of a classical viscous fluid, specifically because the viscosity varies spatially and can even be negative. We will defer a fuller discussion of the origin and role of viscosity in the SP system to future work but, continuing with the simple fluid analogy of the Schrödinger-Poisson model, one can suggest defining something analogous to a Reynolds number for this model as In normal fluid mechanics (e.g. Landau & Lifshitz (1959)), the Reynolds number is a dimensionless quantity formed from u (a velocity), l (a reasonable length-scale for the system) and ν is the kinematic viscosity of the fluid. We see in Section 3 that ν is essentially a kinematic viscosity. Although ν is well defined, it is not clear how one would define the other two quantities. There are many velocity quantities present in this model and it is not exactly clear which one we should choose. There is a momentum associated with the wave-function, but this is not uniform in space and it is not clear how one could convert this to a single number. There is also the velocity of the outward peak of the void, v peak , a sort of expansion velocity. Although this is not well-defined, it is a single number with dimensions of velocity that captures the behaviour of the void. There is also the issue of the length-scale to choose for this Reynolds number. It is obvious to choose the size of the void, however, since the void is expanding, this length is not constant. For this we have decided to choose the initial size of the void for our length. This is still somewhat arbitrary as the initial time and therefore size of the void is arbitrarily chosen.
To test the reliability of this pseudo-Reynolds number, we looked at the relationships between the quantities used to define it. Since u is the calculated value v peak , we started by fixing l, and comparing various ν values with the v peak produced by this. As seen in Figure 8a, it is clear that ν and v peak are proportional to one another. Then we fixed ν and calculated v peak for various l values and as seen in Figure 8b, it is clear that l and v peak have   an inversely proportional relationship. Due to the nature of how v peak is calculated it is not possible to fix u and look at the relationship between ν and l, but it suffices for a quick illustration.
The next obvious step is to compute a Reynolds number for the system described above. The initial size of the void can be seen in Figure 3. However, it is not very clear in Figure 3 where the edge of the void is,or what the edge of a void would mean, but we have chosen it, in this case, to be where the density begins to 'drop off'. This was defined to be at ±1, therefore making this length l = 2.0. The frequency parameter, or viscosity, was set to be ν = 0.05. As seen above the peak velocity was calculated and found to be v peak = 1.55. This would give us a Reynolds number R = 62.0. This is equivalent to a void 30 Mpc wide, with expansion velocity of 27 km s −1 and a viscosity of 4 × 10 10 km 2 s −1 . The actual value for the effective viscosity is not well known, so we used the value stated and discussed in Gibson (2006) for illustration.

CONCLUSIONS AND DISCUSSION
In this paper we have explored the dynamical evolution of cosmic voids through the lens of the Schrödinger-Poisson approach. Following Widrow & Kaiser (1993), we argue that whether or not dark matter does behave quantum-mechanically, this is an advantageous approximate approach that avoids the formation of caustics in the Zel'dovich approximation outlined in Section 2.1. In an era dominated by massive N -body simulations there is still space for analytical approaches such as these because they are much less expensive from a computational perspective and allow us to explore parameter space much more quickly.
Combined with other works using the SP model to study gravitational collapse such as Johnston et al. (2010) and Uhlemann et al. (2019), demonstrates that this has much greater potential for the evolution of voids than the the Zel'dovich approximation; see Figure 5. The efficacy of this approach is further enhanced by the fact that this model has scaling properties so that computations for voids of different sizes and initial expansion rates can be computed cheaply.
Obviously isolated voids such as those we have studied here are idealized so their main function is as a testing ground to enable us to develop insight into the behaviour of the SP system. In future work we plan to explore the evolution of more realistic 3D structures. A particularly interesting potential application of the approximate dynamics deployed here is the problem of cosmic reconstruction, for example from redshift-space to real space, or evolving structure back in time from the present epoch to some earlier state, discussed for example by Brenier et al. (2003).
The above comments apply to our use of the SP dynamics as approximate model for the behaviour of collisionless cold dark matter. If it turns out that dark matter is not cold but instead is for example ultra-light axionic matter then ν takes on physical relevance as it is inversely proportional to the mass of the light particle. It is important to stress however that the SP approach has a much wider range of potential use than this specific scenario.

ACKNOWLEDGMENTS
PC is grateful to William Hunter for discussions during a preliminary investigation of this topic. We thank Cora Uhlemann and Alex Gough for their constructive comments and interesting conversation regarding the Schrödinger-Poisson formalsim. We also thank the referees for their very helpful review which aided the presentation of the results in this paper.