Construction and evolution of equilibrium configurations of the Schr\"odinger-Poisson system in the Madelung frame

We present the construction of ground state equilibrium configurations of the Schr\"odinger-Poisson (SP) system in the Madelung frame and evolve such configuration using finite volume methods. We compare the behavior of these configurations when evolved within the SP and Madelung frames, in terms of conservation of mass and energy. We also discuss the issues of the equations in the Madelung frame and others inherent to the numerical methods used to solve them.


Introduction
The Schrödinger-Poisson (SP) system of equations has recently called attention because these equations rule the dynamics of the Fuzzy Dark Matter model (FDM), which proposes that dark matter is an ultralight spinless boson (e.g. [1][2][3][4]). In this model, Schrödinger equation plays the role of the Gross-Pitaevskii equation for a Bose gas whereas Poisson equation provides the potential trap that contains the bosons, whose source is the bosonic gas density itself. The solution of the SP system of equations is the essence of the analyses of this dark matter candidate.
Among the most interesting discoveries within this model is that structures resulting from the evolution of dark matter fluctuations, accommodate in density configurations with a core and tail (e.g. [5,6]). The core happens to have the profile of ground state equilibrium configurations of the SP system [7,8], which have been found, in isolated scenarios, to be late-time attractor solutions [9,10]. Now, on the numerical aspect of the simulations based on the numerical solution of the SP system under various astrophysical scenarios, it happens that different numerical methods and approaches are used, see e. g. [11,12] for reviews on the subject. The various analyses use two main frames for the solution of the SP system, some of them consist on the direct solution of the SP equations, whereas others solve the hydrodynamical, Madelung version of the equations [13,14].
Examples solving the SP system include [15], where the first structure formation with ultralight dark matter was experimented with, [5] where the structure formation within the FDM model using high resolution simulations receives an important boost, in [16] the oscillation mode spectrum of cores is studied, in [17] the solitonic behavior of cores is revised, the attractor nature of cores is first presented in [10] for the collapse of non-spherical fluctuations, in [18] the star forming structures in FDM Filaments is presented, in [19] zoom in simulation of a galactic halo formation solves the SP equations together with SPH methods, simulations of core mergers are developed in [20] and in [21] tidal disruption of FDM subhalo cores is studied.
On the other hand equations in the madelung frame have been used to study various other problems, for example, the core-halo mass relations and the possible formation of supermassive black holes [22], in [23] the existence of vorticity is analyzed in core structures, in [24] local properties of cores is studied, including random motion and collisions, in [25] the Jeans mass shrinking is studied within the scalar field dark matter model, the merger of cores was seen in [26], in [27] the equations are solved using the SPH method and in [28] the Schrödinger-Poisson Vlasov-Poisson correspondence is proposed.
Interesting issues emerge in the Madelung frame though, first of all is that Schrödinger equation is cast in an analog fashion to Euler equations, and therefore the methods developed for the evolution of fluids can be applied. On the other hand, an important drawback of the Madelung frame is that unlike Euler evolution equations of a compressible fluid, there is not an equation for the balance of energy, and there is no clear Equation of State (EoS) for this quantum fluid analog. Various proposals starting from the Boltzmann equation indicate an EoS for example in [25] and [28] where the SP system is identified with the Poisson-Vlasov system. In our approach below, no EoS or microscopic approach is considered.
Amid the important advances in FDM simulations, the comparison of numerical solutions in the SP and Madelung frames are needed in order to learn about the pros and cons of each approach, and some studies have started in the structure formation scenario [29]. In this paper we go slightly back and practice a comparison in the simplest of scenarios where the two frames can be compared, namely, the construction and evolution of ground state equilibrium configurations, which, as said before, plays an essential role as galactic cores in the FDM astrophysics. We construct the ground state equilibrium configurations directly on the Madelung frame and evolve them with a code that implements methods used in hydrodynamics in order to compare the dynamics in the two frames.
The paper is organized as follows. In Section 2 we rewrite the SP equations in the Madelung frame and construct the ground state equilibrium solutions. In Section 3 we draw a numerical method suitable for the evolution of these configurations and compare the dynamics of equilibrium configurations in the two frames. Finally in Section 4 we draw some conclusions.

Madelung transform
We start by writing the SP equations in the FDM regime: where Ψ is the wave function or parameter order of the boson gas and V is the potential trap, which is the gravitational potential sourced by the gas mass itself. These equations are already scaled so that constants are absorbed. The Madelung transformation defines the wave function in terms of new variables ρ, v, S, Ψ = √ ρe iS , and transforms Schrödinger equation into [13]: equations for ρ and the phase S. By further defining v = ∇S, the SP system is finally transformed into the constrained evolution system where is known as the quantum potential. In this frame ρ is interpreted as the density of a fluid and v as its velocity field. Equations (3) and (4) correspond to the transformed Schrödinger equation, which looks now as a set of two flux balance laws corresponding to the mass conservation and a momentum density balance, analog to Euler equation for a fluid, except by the quantum potential. Notice that unlike the fluid dynamics equations, there is no equation for the balance of an energy neither an equation of state for the fluid, although in some regimes it is possible to identify the fluid with a polytrope with adiabatic index Γ = 2 [30].

Diagnostics
The quantities to monitor in both frames are the mass of the system M, and its total energy E = K + W, where K and W are the kinetic and potential energies defined in each frame by that hold when integrated in the spatial domain.

Ground state equilibrium configurations in the SP frame
We briefly summarize the well known construction of equilibrium configurations in the SP frame, extensively described in e. g. [7,8]. The first assumption is spherical symmetry, and thus spherical coordinates (r, θ, φ, t) are appropriate, so that the wave function depends on the radial coordinate and time Ψ = Ψ(r, t). The second assumption is that the wave function depends harmonically of time, that is Ψ(r, t) = e iωt ψ(r), which ensures |Ψ| 2 is time independent and consequently the gravitational potential as well. Equations (1)-(2) are reduced to which define an eigenvalue problem for ψ(r) provided suitable boundary conditions for the wave function ψ(0) = ψ c and dψ/dr(0) = 0 at the origin, and isolation, which means ψ(r → ∞) → 0 and dψ/dr(r → ∞) → 0. For the gravitational potential the conditions at the origin are V(0) = V c , dV/dr(0) = 0 and at infinity a monopolar condition is used V(r → ∞) = −M/r where M = |ψ| 2 r 2 dr. Fulfillment of these conditions imply that for each central value of the wave function ψ c there is a unique eigenfrequency ω that satisfies the equations. In this sense, one can construct a one parameter family of solutions labeled by the central value ψ c . Finally, unlike excited state solutions, ground state configurations are characterized by the condition that ψ(r) has no zeroes within the domain of integration, otherwise solutions that satisfy the boundary conditions and has zeroes, are known as excited states (see e.g. [7][8][9]).
This eigenvalue problem is solved on a discrete domain r ∈ [0, r max ] using the shooting method as described in [8], with a tolerance on the fulfillment of boundary conditions at large radius. The solution for ψ c = 1 in Fig. 1.

Ground state equilibrium configuration in the Madelung frame
These are spherically symmetric stationary solutions of the system (3)- (5). For the construction of a stationary configuration we also assume spherical symmetry and use spherical coordinates. The flow is assumed to be stationary, which drops the time derivatives in equations (3)-(4). The mass conservation reduces to an identity whereas the momentum density balance becomes the equation which can be integrated immediately to obtain a constraint between the gravitational and quantum potentials Q(r) Considering the quantum potential is given by (6), this equation is a second-order differential equation for the density ρ: Notice that substitution of ρ = √ ψ in this equation gives the stationary Schrödinger equation (10) with ω = −V 0 . In order to solve this equations we write it down as a first-order system by defining u(r) = ρ (r). In this way Eq. (13) becomes Now, Poisson equation (5) is rewritten also as a first order system using the definition of m (r) = 4πr 2 ρ(r) which leads to the equations The system (14)- (17) is the set of equations for equilibrium configurations to be solved. However notice that their form is not suitable for numerical integration at the origin. Expansion of u (r) near zero gives: is the equation to be solved at the origin instead of (14). The boundary conditions imposed to these equations are m(0) = 0, ρ(0) = ρ c , ρ (0) = u(0) = 0 and V(0) = V c , a value that can be arbitrary, and at infinity lim r→∞ ρ(r) = ρ ∞ , lim r→∞ u(r) = 0. With these conditions the system (14)-(18) becomes an eigenvalue problem for the eigenvalue V 0 .
We numerically integrate the system and approximate the value of V 0 using the Shooting method, on a discrete domain r ∈ [0, r max ], with the condition of minimizing the value of the where ρ ∞ is a finite small value that helps approximate the solution to zero at infinity with finite numerical precision.
In Figure 1 we show a zoom of the density profile ρ obtained for the solution with central density ρ c = 1, and ρ ∞ = 0 within a small tolerance. Superposed we show the solution of the eigenvalue problem (10)-(11) for ψ c = 1. The eigenvalue of the solution in the two frames is

Evolution
A further comparison between the solutions in the SP and Madelung frames is the evolution of these configurations. In fact using one or the other frame motivates the use of different numerical methods. In the SP frame it is common to use Finite Difference methods, with a variety of time integrators, explicit or implicit, because Schrödinger equation is dispersive, which prevents the formation of discontinuities. In the contrary, Eq. (4) is quasilinear, which may lead to the formation of discontinuities and shocks even for smooth initial data, and then other numerical methods are needed.
In the SP frame it is possible to show that when solving (1)-(2) the wave function Ψ oscillates with a frequency that coincides with ω obtained from the solution of the eigenvalue problem of Eqs. (10)- (11). Nevertheless in the Madelung frame there is not an equivalent diagnostics.
A very important aspect is the differentiation of the wave function Ψ, since in general it can be considered for an equilibrium configurations as a smooth function. On the other hand, in the Madelung framework some of the variables are not even continue. We can see this from the definition of S, namely the argument of the function Ψ, which is not defined at the origin, and therefore for all cases where S is not spatially constant, it will lead to initial data with discontinuous velocities that will tend to produce shock waves, which do not appear in the SP framework. A discontinuity of the velocity leads to a discontinuity of the density, and consequently Q and ∇Q are undefined for the right hand side of Eq. (4). Only in cases where the function S is constant can the evolution lead to a solution, even a weak solution.
We then consider that problems in which SP and Madelung frames can be compared, need to be those with a constant velocity field. In what follows we describe a comparison of the evolution of the ground state configuration in each of the frames, accompanied with a description of the respective numerical methods.
To start the evolution we interpolate the equilibrium solutions from Figure 1, constructed in spherical coordinates into a 3D cubic domain described with Cartesian coordinates D = [x min , x max ] 3 , uniformly discretized with N cells along each direction that defines the numerical Once the equilibrium configurations in the SP frame is interpolated into the 3D domain described in Cartesian coordinates, we solve the system (1)-(2) to evolve the configuration in the SP frame using appropriate numerical methods. Likewise, when the equilibrium configuration constructed within the Madelung frame is interpolate into the 3D domain we solve equations (3)-(6).

Evolution in the SP frame
. The wave function evolves according to Schrödinger equation (1), whereas Poisson equation (2) is a constraint. For the evolution two methods are common and implemented here: the Method of Lines (MoL) and the implicit Cranck-Nicholson (CN).
For the MoL, the semi-discrete version of the equations uses Finite Differences with second order accurate spatial derivatives in Schrödinger equation: where Ψ i,j,k is the wave function at point (x i , y j , z k ) ∈ D d at time t n . We integrate this equation from time t n to time t n+1 using a third order Runge-Kutta explicit integrator. The Crank-Nicolson method assumes that the evolution of the wave function from time t n to time t n+1 is constructed as follows where we add the upper label to the wave function indicating the time. We solve for Ψ n+1 i,j,k using the Alternating Direction Implicit (ADI) strategy, that splits the application of the Hamiltonian along each of the three spatial dimensions as follows: where R i,j,k , S i,j,k and T i,j,k are auxiliary numbers that store the values of the wave function after applying the derivative operator along each of the spatial directions. Each of the three first equations defines a tridiagonal system of equations when derivatives are discretized: where α = 1 4 i ∆t ∆x 2 and β = 1 2 i∆t, that we solve using forward and backward substitution. Notice that the potential is evaluated at the intermediate time V n+1/2 , which is important because the potential is time-dependent. We calculate this potential as the average V n+1/2 = 1 2 (V n+1 + V n ), where V n+1 is calculated by solving the system (22) for Ψ n+1 , then solving Poisson equation for the source |Ψ n+1 | 2 to obtain V n+1 , going back the time-step, and integrating in time with the averaged potential.
In both MoL and CN methods, Poisson equation is solved using a Multigrid method that uses a 2-levels V-cycle as in [31]. For the evolution of the equilibrium configuration we impose a boundary condition consistent with the isolation of the system and use a monopolar boundary condition for the gravitational potential V ∂D = −M/r ∂D , where r ∂D is the distance from the origin to a given point of the numerical boundary and M is given by Eq. (7).
For the wave function we implement a sponge that absorbs the modes approaching the faces of the boundary and prevents their reflection. This sponge consists in the addition of an imaginary potential that acts as a sink as described in [8,32]. Briefly, the potential V in (1) is redefined as V → V + V im , where V is the solution of (2) and V im is a spherical function with tanh-profile that goes from zero within a sphere containing the region where the interesting physics happens and minus one near the boundary of the domain, which absorbs the wave function if it approaches the boundary.

Methods for the Madelung frame.
In this frame the evolution of the system for density ρ and velocity field v is ruled by equations (3) and (4), whereas again Poisson equation (5) is a constraint that has to be fulfilled at each time step during the evolution of density and velocity.
We use a Finite Volume discretization to solve the system (3)-(4) following [33]. Essential to the method is the appropriate construction of numerical fluxes and the characteristic structure of the equations, which can be cast in the form where u is a vector of conserved variables, F = [ F x , F y , F z ] a vector of fluxes and S a vector of sources. These elements for the system (3)-(4) read as follows We calculate the numerical fluxes using the HLLE formula [34], according to which the fluxes at the left u L and at the right u R from each intercell boundary, along each of the three Cartesian directions are where i = x, y, z, λ − and λ + are the approximations of the characteristic velocities which are calculated as Notice that with this approach we have not involved the pressure of the fluid, in this case the pressure-like tensor associated to the quantum fluid within the fluxes. Instead, the gradient of Q is assumed to contribute as a source.
Atmosphere. An inherent ingredient of FV methods is the use of an atmosphere, which is defined as a minimum value of the fluid density, set to ρ atm . In fluid dynamics it is useful to avoid the divergence of temperature T = p/ρ, with p the pressure, which leads to the divergence of all other involved state variables. In our case there is no pressure, but there is the quantum potential Q in (6), which diverges or is undefined for ρ = 0, unless ∇ 2 √ ρ compensates the divergence.  Figure 2. Central density |Ψ( 0, t)| 2 for the solution in the SP frame and ρ( 0, t) in the Madelung frame, as functions of time. This plot shows that the first oscillation mode is similar, although not perfectly equal in the two frames. In the SP frame the period coincides with that in [16], whereas the period in the Madelung frame shrinks in time.
Boundary conditions. In order to implement the isolation condition of the system, unlike the sponge in the SP frame, we impose outflow boundary conditions on the fluid variables. In order to solve Poisson equation we use a monopolar boundary condition for the gravitational potential V ∂D = −M/r ∂D , where r ∂D is the distance from the origin to a given point of the numerical boundary and M is given by Eq. (7), and solve using the same Multigrid 2-level V-cycle method used for the SP frame.

Evolution of an equilibrium configuration
In order to compare the essentials of ground state configurations we evolve the equilibrium configurations from Section 2 and track their behavior. Using the numerical methods described above for the evolution, we integrate in time the equations for these configurations centered at the coordinate origin, and in Figure 2 we show the central value of the density, |Ψ( 0, t)| 2 in the SP frame and ρ( 0, t) in the Madelung frame. The result indicates that the configuration remains nearly stationary with an oscillation mode consistent with the dominant spherical mode of the configuration with period T = 21.64 as pointed out in [16] using the SP frame, with both the MoL and CN methods. When using the Madelung frame there is a reduction of the period in time.
In Figure 3 we show the evolution of the mass M and total energy E = K + W as functions of time. We can see that in the Madelung frame the mass is slightly bigger than in the SP frame, which is due to the contribution of the atmosphere. Since this atmosphere is of order of ρ atm ∼ 10 −8 , it contributes with a mass of order M atm ∼ 10 −3 . On the other hand, the mass in the SP frame decreases when using the MoL, which is due to the dissipation of the time integration with the explicit RK3 method, whereas the conservation is better when using the CN method, indicating the evolution is closer to unitary.
Concerning the total energy, the values in the two frames disagree due to the difference in mass, which contributes to the gravitational energy W. For this reason in the Madelung frame the total energy E is smaller than in the SP frame.

Boosted equilibrium configuration
A second comparison is dynamical and corresponds to a boosted equilibrium configuration. For this we apply an initial velocity v 0 x = 1 to the equilibrium configuration placed at the initial position (−10, 0, 0), in order to obtain a configuration traveling to the right along the x-axis.  The linear momentum at initial time is applied to the equilibrium configuration in the SP frame as Ψ( x, 0) = ψ 0 e iv 0 x x , where ψ 0 is the wave function obtained from the solution of the eigenvalue problem (10)- (11). In the Madelung frame the density is ρ( x, 0) = ρ 0 , and the initial speed is set to v( x, 0) = (v 0 x , 0, 0) where ρ 0 is that of the solution of the eigenvalue problem (14)- (18).
Snapshots of the density are shown in Figure 4. Notice that the maximum of the configuration using the hydrodynamical method in the Madelung frame locates at the appropriate position, since the velocity field is one. On the other hand, the configurations evolved in the SP frame retard with respect to the appropriate location of the maximum density. Nevertheless, this is due to the implementation of the momentum, which actually is a phase, and not a true velocity field at initial time.  In Figure 5 we show the evolution of the total energy E = K + W and the mass M in the two frames. Notice that the energy is well preserved in the two frames and the three methods. The mass decreases in a small fraction in the SP frame that we attribute to the sponge, whereas in the Madelung frame the mass remains relatively constant.  Figure 5. Evolution of the total energy E = K + Wand the mass M as a function of time in the two frames for the boosted configuration.

Conclusions
We have described the construction of ground state equilibrium configurations of the Schrödinger-Poisson system directly from the equations in the Madelung frame and verified that the solutions in both frames coincide within numerical accuracy.
Since the equations in the Madelung frame are similar to those of a fluid, we implement a simple scheme involving shock capturing methods used in fluid dynamics to simulate the evolution in time. We compare the evolution in SP and Madelung frames using the evolution equilibrium configuration when at rest and when the configuration is boosted.
We find that important limitations of the FV method arise. One of them is the need of an atmosphere used to avoid the singularity of the quantum potential at rarified regions. We notice that this also happens in simple cases within Quantum Mechanics, for example in excited states for a particle in a box or a harmonic oscillator potential, where the density has zeroes, then using an atmosphere is needed to avoid the singularity 1/ √ ρ of Q. However the use of the atmosphere introduces a non-derivable √ ρ that will spoil the factor ∇ 2 √ ρ of Q in the intersection between the atmosphere and higher values of ρ. It seems that using solely the Madelung frame leads to a dilemma.
This type of problems seems to motivate the use of patches to the FV method, like in [35] where in low density regions the Smoothed Particle Hydrodynamics (SPH) method is used; or the use the SPH method to solve the evolution equations in the Madelung frame for ρ and v in [27]. It would be interesting to investigate how SPH develops in these simple scenarios.
Funding: IAR receives support within the CONACyT graduate scholarship program under the CVU 967478. This research is supported by grants CIC-UMSNH-4.9 and CONACyT Ciencias de Frontera Grant No. Sinergias/304001. The runs were carried out in the Big Mamma cluster of the Laboratorio de Inteligencia Artificial y Supercómputo, IFM-UMSNH.

Data Availability Statement:
The data underlying this article will be shared on reasonable request to the corresponding author.