Cosmological Simulation for Fuzzy Dark Matter Model

Fuzzy Dark Matter (FDM), motivated by string theory, has recently become a hot candidate for dark matter. The rest mass of FDM is believed to be $\sim 10^{-22}$eV and the corresponding de-Broglie wave length is $\sim 1$kpc. Therefore, the quantum effect of FDM plays an important role in structure formation. In order to study the cosmological structure formation in FDM model, several simulation techniques have been introduced. We review the current status and challenges in the cosmological simulation for the FDM model in this paper.


INTRODUCTION
The nature of dark matter is one of the key mysteries of modern cosmology and physics. Dark matter is widely believed to be dominated by cold dark matter (CDM), supported by different observations such as the mass-to-light ratio of clusters of galaxies [1], the rotation curves of galaxies [2], the Bullet Cluster [3], the cosmic microwave background (CMB) [4] and the large-scale structure of the universe [5]. However, despite its success on large scales, the CDM paradigm faces three problems on small scales, known as the "small-scale crisis" [6]: (i) the missing satellite problem, (ii) the cusp-core problem, and (iii) the too-bigto-fail problem. The key point of these problems is that CDM model predicts too much or too compact structures on small scales. Two approaches are under discussion to solve these problems. One is to smooth out the small-scale structure by astrophysical processes [7], and the other is to introduce alternative dark matter models like warm dark matter (WDM) [8], decaying dark matter (DDM) [9], self-interacting dark matter (SIDM) [10] and fuzzy dark matter (FDM) [11].
In the FDM model, the dark matter particles are made of ultra-light bosons in Bose-Einstein condensate (BEC) state [12]. As an alternative to CDM, it suppresses small-scale structures while keeps the success of CDM on large scales [13,14,15]. The FDM model is phenomenologically different from the CDM model due to its effective "quantum pressure" (QP) which originates from the uncertainty principle [16]. Apart from FDM, this model has many other names, such as wave dark matter (ΦDM), ultra-light axion (ULA), scalar field dark matter (SFDM), which is mainly due to historical reasons. These models have slightly different self-interactions and theoretical considerations. However, models in which dark matter has no self-interaction are phenomenologically the same as FDM. The history of FDM and implementation is summarized in Ref. [17].
The predictions of FDM with mass ∼ 10 −22 eV are consistent with observations of the large-scale structure [18], high-z galaxies, CMB optical depth [19], and the density profiles of dwarf spheroidal galaxies [20]. The tightest constraints come from the comparison of the recent Lyman-alpha forest observations with FDM hydrodynamic simulations [21,22,23]. These works claimed that FDM model with particle mass less than 10 −21 eV is ruled out at 95% confidence level. However, it has been pointed out that the quantum pressure plays quite non-trivial role in structure formation, which is neglected in the hydrodynamic simulations for Lyman-alpha forest. The simulation uncertainties are also important issues for making such tight constraints [15].
In order to constrain the parameter space of the FDM model, or to look for smoking-guns for it, simulation is extremely important. There have been eight different codes proposed to perform simulations for the FDM model [20,24,25,26,14,27,28,29]. They can be classified into two major approaches: solving the Schrödinger-Poission equation or the "equivalent" Madelung equations. We reviewed these works and summarized them into a table. The pros and cons of these different simulation methods were clearly stated. We gave some comments on the current status and challenges for FDM simulation.
The paper is organized in the following sections: we review the basic equations necessary for the FDM model in Sec. 2, the simulation treatments and code comparison in Sec. 3, the current status and challenges of FDM simulation in Sec. 4, and finally we discuss about possible smoking-gun signatures for the FDM model.

BASIC EQUATIONS
To study the structures on galactic scales in the low red-shift universe, it is safe to ignore the self-interaction of the scalar field describing the FDM. The action has the following form where we follow the convention in Ref. [30]. The related de Broglie wavelength of particles with rest mass m is λ 2π = mv = 1.92 kpc 10 −22 eV m Using the least action principle and WKB approximation in the non-relativistic limit, one can simplify the governing equations of the scalar field to the Schrödinger-Poisson equations, where Ψ is the plane wave description of the scalar field φ, and V is gravitational potential,

FDM Simulation
The wave function Ψ can be written as Ψ = ρ m exp iS (6) in terms of the number density of FDM particles ρ/m, while we can define the gradient of S to be the momentum, ∇S = mv.
After transforming the wave function, the Schrödinger-Poisson equations can be written in an equivalent fluid dynamics form with the continuity equation, and the Euler equation, dv dt where the quantum pressure Q is defined as Eqs. (8) and (9) are known as the Madelung equations [31,32,33]. In cosmological simulations, we also need to consider the expansion of the universe. Eq. 3 should be rewritten as: where H =ȧ/a is the hubble parameter and a is the scale factor of the universe. The Madelung equations change accordingly dρ dt In short, the difference between the FDM and CDM models lies in the existence of the quantum pressure Q. The derivation of the quantum pressure is well described in text books about BEC. Since Q ∝ m −2 , the quantum pressure in lab BEC systems is negligible. However, this effect is important in the FDM model whose particle mass is around m ∼ 10 −22 eV.
From the linear pertubation of the equations (12) and (13), the density contrast evolves according tö A solution is given by a plane wave with wavenumber If k < k J (a), gravity dominates and the structure will collapse, while modes with k > k J (a) will expand due to the repulsive quantum pressure. So this is the Jeans wavenumber of the FDM model. The growing mode D + (k, a) and the decaying mode D − (k, a) of equation (14) are For k ≪ k J (a), the two modes return to the CDM solutions D + ∝ a and D − ∝ a −2/3 , meaning that the FDM and CDM have the exact same behavior on the large scales. On the other hand, for k ≫ k J (a), the growth of the structure in FDM is suppressed because D + ∝ k −4 . But the Jeans wave number k J (a) ∝ a 1/4 is growing over time, and so the small-scale structures will eventually start growing: the smaller the scale (the larger the wavenumber), the later this mode started growing.

SIMULATION REVIEW
A typical N-body cosmological simulation contains the following steps: 1. Distribute the simulation particles in the simulation box homogeneously and isotropically. The FDM model has the same preparation of these pre-initial conditions as the CDM model. 2. Calculate the matter power spectrum at a relatively high redshift, such as z = 99, according to the prediction of the linear perturbation theory. The modification brought by the FDM model can be either calculated by AxionCAMB [34], or given by the empirical transfer function Ref. [16], where x = 1.61 m Mpc −1 . It has been shown that using these two methods makes little difference, and the empirical transfer function is a good approximation [22]. 3. Perturb the distribution of particles according to the matter power spectrum;

Solve the Euler equation and Poisson equation iteratively (continuity equation is naturally obeyed
using N-body simulation) until the desired redshift, such as z = 0. To incorporate the quantum pressure into the Lagrangian particle tracking simulation scenario, there are four different codes available, summarized in Table 1.
The first three methods in Tab.1 are based on the traditional SPH method. The essence of the SPH method is to first assigning all physical quantities (like density ρ, velocity v and pressure P ) on each simulation particles, then calculate the physical fields by a special interpolation method -kernel smoothing, which in turn give rise to the time evolution of the simulation particles through the Euler where W is a spherical function with finite support. The quantum pressure (10), however, not only depends on the field itself but also its derivatives up to second order. The three different implementations of the SPH methods listed above used three different method to calculate derivatives.
Ref. [26] used the particle-in-cell method: 1. Assign the physical quantities of each simulation particle onto an auxiliary cubic grid; 2. Calculated the derivatives of the physical fields with the finite difference method; 3. Interpolate the derivatives of physical fields back to the position of the simulation particles.
The additional force coming from quantum pressure is given by where the seven-point stencil is used to calculate the Laplacian: Ref. [24] and [28] used similar formulae: with different choices of the auxiliary function: Θ i = 1 for Ref. [24] and Θ i = √ ρ i for Ref. [28]. The force contributed by the quantum pressure is given by where f j = 1 + h j 3ρ j k m k ∂W r jk /h j ∂h j is a correcting factor when variable smoothing length is used [35].
All the three methods above involve the estimation of density and its derivatives on the grids or at the position of the particles, and the the force (20) and (21) have the form of many body interaction. Therefore, their computational costs are relatively high. Ref. [27] improves the SPH method by reducing the quantum pressure to two body particle-particle interaction, hence the additional force can be easily added to the tree algorithm in the TreePM method without the need to resort to the SPH method, and the computational time are greatly reduced.
From Tab.1, we conclude that all these Lagrangian based simulations cannot produce granular structures which are expected to appear as the result of quantum interference. There are two possible explanations that may be viewed as the fundamental flaws of Lagrangian based simulations of the FDM model:

FDM Simulation
• The Schrödinger-Poisson equations and Madelung equations are not strictly equivalent. As proved in Ref. [36], a quantization condition L v · dl = 2πj (j ∈ Z and L is any closed loop.) is necessary to recover the Schrödinger-Poisson equations from the Madelung equations, which is neither checked nor possibly obeyed in the Lagrangian based simulations. • The smoothing kernel method which is indispensable in Lagrangian based simulations cannot accurately estimate the matter density field and its second order derivative simultaneously if merely a single smoothing length is used. As proved in Ref. [37], the relative error of the estimation of the second order derivative of the density field could be as large as 100% when the smoothing length is chosen to minimize the error of density estimation [38,39], to solve the Poisson equation.
It is a consensus that in the center of a virialized FDM halo, there is a solitonic core made of wave function in the ground state with the same phase [40]. Although the core-like structures appear in the Lagrangian based simulations, they are not very trustworthy due to the two reasons listed above.
Apart from Lagrangian based simulations, FDM model can also be studied by Eulerian based simulations summarized in Table 2. The physical fields on the grid also need to be suitably set at the initial moment according to the cosmological linear theory prediction. The time evolution of the wave function is given by where T is the time-ordering symbol. For a sufficiently small time step, it can be approximated as which can be further splitted into three operations according to the Baker-Campbell-Hausdorff formula: This formula has a close resemblance to the kick-drift-kick time evolution in the particle method. The "kick" step is done in real space, which effectively just changes the phase angle at each point. The "drift" step is completed in the Fourier space: The main differences the Eularian based and the Lagrangian based methods are that the original Schrödinger-Poisson equations are solved in the former, not the transformed Madelung equations in the latter, and the Eulerian method can be used to reliably estimate second order derivatives of the fields. From Tab. 2, we find that most of the Eulerian based simulations can produce granular structures and solitonic cores. The sizes of the simulation boxes in these simulations, however, are not large enough to be considered as cosmological scale simulations. The daunting computational costs make it too difficult to perform simulations with box size larger than 10 Mpc/h For current simulation codes, there have to be a trade-off between the fidelity of the simulations and the scale of the simulations. On one hand, both the granular structures and the solitonic cores are smoking To simplify the generation of the granular structures, a self-consistent method was introduced in Ref. [42]. In their simplified model, the halo is composed of "smooth" density distribution along the radial direction and "granular" interference structure along the angular direction, the radial direction density profile is given by a typical guess and the angular direction density distribution is described by spherical harmonics. By fitting to the simulations, they find that the fermionic King model is the best fit energy distribution function and the generated halo is quite similar to that in simulation. With this method, they can generate a halo as massive as Milky Way (∼ 10 12 M ⊙ ) with granular structures. However, it is still very difficult to self-consistently construct a very massive halo, such as a cluster scale halo (∼ 10 14 M ⊙ ), because of the limited computational power and poor generating algorithm. Making use of the information about the granular structures, a very promising smoking gun detection method was introduced in Ref. [43,44,45] with Pulsar Timing Array. By studying the modulation of the arriving time of pulses from many pulsars, people can use Pulsar Timing Array to directly detect the granularity of dark matter distribution in Milky Way galaxy. This method can not only set constraints on the rest mass of FDM particles, but also make it possible to confidently claim the existence of FDM. But only if we understand the granular structures of the Milky Way halo in much detail, we can make correct predictions for the modulation pattern observed in Pulsar Timing Array.
The FDM simulation with the largest box size (50 Mpc/h) was performed in Ref. [15]. The limit of such simulations is their poor resolution, and they are unable to resolve the granular structures. The advantages of the Madelung solvers are their efficiency and mature related data analysis tools. Measured by the two large simulations Ref. [15] and Ref. [28], whose algorithms are different, the matter power spectrum at z = 0 is not only suppressed at small scales by the modification of the initial conditions, but also by the effect of quantum pressure with an additional ∼ 10% suppression at small scales. Such suppression is well expected and confirmed by different simulations. However, none of the Schrödinger-Poisson solvers measured the matter power spectrum due to their small box sizes. Therefore, it is still under investigation about the effect of quantum pressure in the structure formation. It is still inconclusive from simulations how much suppression of the matter power spectrum will be introduced by the quantum pressure. In order to study large scale structure and use observations like weak lensing and red shift distortion to constrain FDM model, a much larger box size (∼ 500 Mpc/h) is necessary. This is still a challenge for all the codes, among which the Madelung solvers are more hopeful to reach such a goal.
Recently, it was claimed in Ref. [46] that the dynamical evidence of the existence of the solitonic core in the center of Milky Way was found. It was also claimed in Ref. [40] that the existence of the solitonic core can solve the cusp-core problem. Both of these two studies favor a FDM with m ∼ 10 −22 eV. However, the recent constraints from Lyman-alpha forest implied that FDM with m < 10 −21 eV was ruled out [22,21]. Such a tension is a problem for FDM models. Ref. [15] suggested that consideration of the important role of quantum pressure and systematics may resolve this tension. Ref. [30] also suggested that other astrophysical processes like patchy reionization can relieve the tension. These works suggest that FDM with m ∼ 10 −22 eV is still possible. More serious studies are clearly needed. Other than the constraints from Lyman-alpha forest, the thickness of the stellar stream and the recent EDGES experiment also set constraints on FDM with m > 5 × 10 −21 eV and m > 1.5 × 10 −22 eV [47,48], respectively. These independent constraints using different methods also seem to be in tension with the requirement to solve the small-scale crisis. If we introduce more arguments and modifications to the current FDM model, we can relieve the tensions but lose the beauty and simplicity of the FDM model.
We notice that a runaway tidal disruption may be very important in the structure formation of FDM halos [49]. It was found in simulations that the solitonic cores in the center of halos can be easily tidal disrupted in a runaway pattern when they rotate around the central massive halo. This means, even with m > 10 −22 eV, FDM model may still be able to solve the small-scale crisis due to the runaway tidal disruption. A systematic study considering such mechanism needs simulations with much larger box size and high resolution at the same time.

SUMMARY
The cosmological simulation is important for understanding the structure formation and looking for smoking-gun signatures for the FDM model. The current simulation codes are not adequate to study the large-scale structure and halo properties under the framework of FDM. The codes designed to solve the SP equations are highly accurate, and many important features of the FDM model such as solitonic cores and granular structures are discovered. But these codes are too computationally heavy to perform simulations with large box size. The codes designed to solve the Madelung equations are less accurate, but the suppression of the matter power spectrum and halo mass function are given by simulations with these codes. These codes are more efficient and compatible with existing data analysis tools, but not accurate enough to resolve granular structures. All of the current FDM cosmological simulations are not large enough in terms of the box size to study large-scale structures systematically. Much efforts are needed to improve these methods.
The FDM model is an interesting alternative to the CDM model. The "small-scale crisis" in CDM might be solved in the FDM model, but more studies are needed to confirm this suggestion. The tensions from different observations on the rest mass of FDM particles may be relieved in many ways, such as considering the systematics in the simulations, invoking astrophysical processes and new mechanisms in the FDM model. In order to understand the structure formation under the framework of FDM model, more and better simulations are important.
With the Madelung solvers, the simulations with box size ∼ 500 Mpc/h can be expected in the near future, which is sufficient to constrain FDM models with observations such as weak lensing and red shift distortion. We need to first make sure that all different codes draw to the same conclusion about these observations. We also need to make sure that all the observables we measured both in observations and in simulations are consistent with the FDM framework. The effect of quantum pressure on the large scale structures can be degenerate with other models such as Warm Dark Matter (WDM), Self-Interacting Dark Matter (SIDM), Decaying Dark Matter (DDM) and so on. Therefore, even if we find conclusive evidence that the small scale structures are suppressed from observations, it is still not conclusive to claim that the FDM model is the correct model of the dark matter. On the other hand, the Schrödinger-Poisson Solvers disclose the possibility of looking for smoking gun signals of the FDM model. Both the existence of a Zhang FDM Simulation solitonic core in the center of a dark matter halo and the granular structure of dark matter halos are unique features of the FDM model, different from all the other models. It is possible to find granular structures for the the FDM model, using Pulsar Timing Array. It is also possible to rule out the FDM model by the next generation observations and more careful data analysis. FDM model is a beautiful model with no more free parameters than CDM model together with WIMP assumption. If we can determine the mass range of FDM particles, it will significantly improve our understanding of dark matter and the universe.

ACKNOWLEDGEMENT
J. Z acknowledges the support from China Postdoctoral Science Foundation 2018M632097.