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 ~10⁻²²eV and the corresponding de-Broglie wave length is ~1kpc. 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.

1. 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 (Bahcall et al., 1995), the rotation curves of galaxies (Einasto et al., 1974), the Bullet Cluster (Clowe et al., 2006), the cosmic microwave background (CMB) (Ade et al., 2016) and the large-scale structure of the universe (Tegmark et al., 2004). However, despite its success on large scales, the CDM paradigm faces three problems on small scales, known as the “small-scale crisis” (Weinberg et al., 2015): (i) the missing satellite problem, (ii) the cusp-core problem, and (iii) the too-big-to-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 (Pontzen and Governato, 2014), and the other is to introduce alternative dark matter models like warm dark matter (WDM) (Colin et al., 2000), decaying dark matter (DDM) (Cheng et al., 2015), self-interacting dark matter (SIDM) (Spergel and Steinhardt, 2000) and fuzzy dark matter (FDM) (Hu et al., 2000).

In the FDM model, the dark matter particles are made of ultra-light bosons in Bose-Einstein condensate (BEC) state (Marsh, 2016). As an alternative to CDM, it suppresses small-scale structures while keeps the success of CDM on large scales (Du et al., 2017; Mocz et al., 2017; Zhang et al., 2018a). The FDM model is phenomenologically different from the CDM model due to its effective "quantum pressure" (QP) which originates from the uncertainty principle (Hu et al., 2000). 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. There are quite a few theoretical studies of such models (Khlopov et al., 1985; Sahni and Wang, 2000; Chavanis, 2012; Mishra et al., 2017). 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 Lee (2017).

The predictions of FDM with mass ~10⁻²²eV are consistent with observations of the large-scale structure (Hlozek et al., 2015), high-z galaxies, CMB optical depth (Bozek et al., 2015), and the density profiles of dwarf spheroidal galaxies (Schive et al., 2014a). The tightest constraints come from the comparison of the recent Lyman-alpha forest observations with FDM hydrodynamic simulations (Armengaud et al., 2017; Iršič et al., 2017; Kobayashi et al., 2017). These works claimed that FDM model with particle mass less than 10⁻²¹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 (Zhang et al., 2018a).

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 (Schive et al., 2014a; Mocz and Succi, 2015; Schwabe et al., 2016; Veltmaat and Niemeyer, 2016; Mocz et al., 2017; Edwards et al., 2018; Nori and Baldi, 2018; Zhang et al., 2018b). 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 section 2, the simulation treatments and code comparison in section 3, the current status and challenges of FDM simulation in section 4, and finally we discuss about possible smoking-gun signatures for the FDM model.

2. 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

$\begin{array}{ll}S={\displaystyle \int }\frac{d^{4}x}{\hslash c^2}\sqrt{-g}\left\{\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -\frac{1}{2}\frac{m^2{c}^2}{{\hslash }^2}{\varphi }^2\right\},& (1)\end{array}$

where we follow the convention in Hui et al. (2017). The related de Broglie wavelength of particles with rest mass m is

$\begin{array}{ll}\frac{\lambda }{2\pi }=\frac{\hslash }{mv}=1.92\text{ kpc }\left(\frac{10^{-22}\text{eV}}{m}\right)\left(\frac{10\text{ km\hspace{0.17em}s}{-}^1}{v}\right).& (2)\end{array}$

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,

$\begin{array}{ll}\begin{array}{c}\text{i}\hslash \frac{\text{d}\Psi }{\text{d}t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\Psi +mV\Psi ,\end{array}& (3)\end{array}$

where Ψ is the plane wave description of the scalar field ϕ,

$\begin{array}{ll}\begin{array}{c}\varphi =\sqrt{\frac{{\hslash }^{3}c}{2m}}(\Psi {\text{e}}^{-\text{i}m{c}^{2}t/\hslash }+{\Psi }^{*}{\text{e}}^{\text{i}m{c}^{2}t/\hslash }),\end{array}& (4)\end{array}$

and V is gravitational potential,

$\begin{array}{ll}{\nabla }^{2}V=4\pi Gm|\text{Ψ}{|}^2.& (5)\end{array}$

The wave function Ψ can be written as

$\begin{array}{ll}\Psi =\sqrt{\frac{\rho }{m}}\text{exp}\left(\frac{\text{i}S}{\hslash }\right)& (6)\end{array}$

in terms of the number density of FDM particles ρ/m, while we can define the gradient of S to be the momentum,

$\begin{array}{ll}\nabla S=mv.& (7)\end{array}$

After transforming the wave function, the Schrödinger-Poisson equations can be written in an equivalent fluid dynamics form with the continuity equation,

$\begin{array}{ll}\frac{\text{d}\rho }{\text{d}t}+\nabla ·\left(\rho v\right)=0,& (8)\end{array}$

and the Euler equation,

$\begin{array}{ll}\frac{\text{d}v}{\text{d}t}+\left(v·\nabla \right)v=-\nabla \left(Q+V\right),& (9)\end{array}$

where the quantum pressure Q is defined as

$\begin{array}{ll}Q=-\frac{{\hslash }^2}{2m^2}\frac{{\nabla }^2\sqrt{\rho }}{\sqrt{\rho }}=-\frac{{\hslash }^2}{2m^2}\left(\frac{{\nabla }^2\rho }{2\rho }-\frac{|\nabla \rho {|}^2}{4{\rho }^2}\right).& (10)\end{array}$

Equations (8) and (9) are known as the Madelung equations (Spiegel, 1980; Uhlemann et al., 2014; Marsh, 2015). In cosmological simulations, we also need to consider the expansion of the universe. Equation (3) should be rewritten as:

$\begin{array}{ll}\text{i}\hslash (\frac{\text{d}\Psi }{\text{d}t}+\frac{3}{2}H\Psi )=-\frac{{\hslash }^2}{2m}{\nabla }^2\Psi +mV\Psi ,& (11)\end{array}$

where H = ȧ/a is the Hubble parameter and a is the scale factor of the universe. The Madelung equations change accordingly

$\begin{array}{ll}\frac{\text{d}\rho }{\text{d}t}+3H\rho +\frac{1}{a}\nabla ·\left(\rho v\right)=0,& (12)\end{array}$

$\begin{array}{ll}\frac{\text{d}v}{\text{d}t}+Hv+\frac{1}{a}(v·\nabla )v=-\frac{1}{a}\nabla V-\frac{1}{a^3}\nabla Q.& (13)\end{array}$

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⁻², 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⁻²² eV.

From the linear perturbation of the Equations (12) and (13), the density contrast evolves according to

$\begin{array}{ll}\ddot{\delta }+2H\dot{\delta }+\left(\frac{{\hslash }^2{k}^4}{4m^2{a}^4}-\frac{4\pi G\bar{\rho }}{a^3}\right)\delta =0.& (14)\end{array}$

A solution is given by a plane wave with wave number

$\begin{array}{ll}{k}_J\left(a\right)={\left(\frac{16\pi G\bar{\rho }{a}^3{m}^2}{{\hslash }^2}\right)}^{1/4}{a}^{1/4}.& (15)\end{array}$

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

$\begin{array}{ll}\begin{array}{l}{D}_{+}\left(k,a\right)=\left[\left(3-x^2\right)cosx+3xsinx\right]/x^2\\ D_{-}\left(k,a\right)=\left[\left(3-x^2\right)cosx-3xsinx\right]/x^2\end{array},x\left(k,a\right)=\sqrt{6}{k}^2/k_J^2\left(a\right).& (16)\end{array}$

For k ≪ k_J(a), the two modes return to the CDM solutions D₊∝a and $D_{-}\propto 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_{+}\propto k^{-4}$. But the Jeans wave number $k_J\left(a\right)\propto a^{1/4}$ is growing over time, and so the small-scale structures will eventually start growing: the smaller the scale (the larger the wave number), the later this mode started growing.

3. 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 (Hlozek et al., 2015), or given by the empirical transfer function (Hu et al., 2000),

$\begin{array}{ll}{P}_{FDM}\left(k\right)=T_F^2\left(k\right)P_{CDM}\left(k\right),T_F\left(k\right)\approx \frac{cos{x}^3}{1+x^8},& (17)\end{array}$

where $x=1.61{\left(\frac{m}{10^{-22}\text{eV}}\right)}^{1/18}\frac{k}{k_{Jeq}},k_{Jeq}=9{\left(\frac{m}{10^{-22}\text{eV}}\right)}^{1/2}\text{Mp}{\text{c}}^{-1}$. It has been shown that using these two methods makes little difference, and the empirical transfer function is a good approximation (Armengaud et al., 2017).

3. Perturb the distribution of particles according to the matter power spectrum;

4. 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.

TABLE 1

Table 1. Summary of the Lagriangian based simulation codes for FDM model (Madelung solvers).

The first three methods in Table 1 are based on the traditional SPH method. The essence of the SPH method is to first assign 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 equation and the equation of state. The kernel smoothing of the field is simplly:

$\begin{array}{ll}{O}_i=\sum m_j\frac{O_j}{{\rho }_j}W\left(\frac{r_{ij}}{h}\right),& (18)\end{array}$

where W is a spherical function with finite support, r_ij is the distance between two particles, and h is the smoothing parameter, notice that it is not the dimensionless Hubble parameter. 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 use three different method to calculate derivatives.

Veltmaat and Niemeyer (2016) used the particle-in-cell method:

1. Assign the physical quantities of each simulation particle onto an auxiliary cubic grid;

2. Calculate the derivatives of the physical fields with the finite difference method;

3. Interpolate the derivatives of physical fields back to the positions of the simulation particles.

The additional force coming from quantum pressure is given by

$\begin{array}{ll}-\nabla Q_i=\frac{{\hslash }^2}{2m^2}{\left(\Delta x\right)}^3{\displaystyle \sum _{j,k,l}\frac{{\left({\nabla }^2\sqrt{\rho }\right)}_{j,k,l}}{\sqrt{{\rho }_{j,k,l}}}}{m}_i\nabla W\left(\frac{\left|r_i-x_{j,k,l}\right|}{h}\right),& (19)\end{array}$

where the seven-point stencil is used to calculate the Laplacian: ${\left({\nabla }^2\sqrt{\rho }\right)}_{j,k,l}=\sqrt{{\rho }_{j+1,k,l}}+\sqrt{{\rho }_{j-1,k,l}}+\sqrt{{\rho }_{j,k+1,l}}+\sqrt{{\rho }_{j,k-1,l}}+\sqrt{{\rho }_{j,k,l+1}}+\sqrt{{\rho }_{j,k,l-1}}-6\sqrt{{\rho }_{j,k,l}}$.

Mocz and Succi (2015) and Nori and Baldi (2018) used similar formulae:

$\begin{array}{l}\nabla O_i=\sum m_j\frac{O_j-O_i}{{\rho }_j}\frac{{\Theta }_j}{{\Theta }_i}\nabla W\left(\frac{r_{ij}}{h}\right),\\ {\nabla }^2{O}_i=\sum m_j\frac{O_j-O_i}{{\rho }_j}\frac{{\Theta }_j}{{\Theta }_i}{\nabla }^{2}W\left(\frac{r_{ij}}{h}\right)-\frac{2}{{\Theta }_i}\nabla O_i·\nabla {\Theta }_i,& (20)\end{array}$

with different choices of the auxiliary function: Θ_i = 1 for Mocz and Succi (2015) and ${\Theta }_i=\sqrt{{\rho }_i}$ for Nori and Baldi (2018). The force contributed by the quantum pressure is given by

$\begin{array}{ll}-\nabla Q_i=\frac{{\hslash }^2}{2m^2}{\displaystyle \sum _j}\frac{m_j}{f_j{\rho }_j}\left(\frac{{\nabla }^2{\rho }_j}{2{\rho }_j}-\frac{|\nabla {\rho }_j{|}^2}{4{\rho }_j^2}\right)\nabla W\left(\frac{r_{ij}}{h}\right),& (21)\end{array}$

where $f_j=1+\frac{h_j}{3{\rho }_j}\sum _k{m}_k\frac{\partial W\left(r_{jk}/h_j\right)}{\partial h_j}$ is a correcting factor when variable smoothing length is used (Springel and Hernquist, 2002).

All the three methods above involve the estimation of density and its derivatives on the grids or at the positions of the particles, and the force (19) or (21) has the form of many body interaction. Therefore, their computational costs are relatively high. Zhang et al. (2018b) improves the SPH method by reducing the quantum pressure to a 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 Table 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:

• The Schrödinger-Poisson equations and Madelung equations are not strictly equivalent. As proved in Wallstrom (1994), a quantization condition ∮_Lv · dl = 2πj (j ∈ ℤ and L is any closed loop.) is necessary to recover the Schrödinger-Poisson equations from the Madelung equations, which is not checked and possibly not obeyed in 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 Silverman (1978), 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 (Merritt, 1996; Dehnen, 2001), 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 (Schive et al., 2014b). Although the core-like structures appear in the Lagrangian based simulations, they are not 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

$\begin{array}{ll}\Psi \left(x,t+\Delta t\right)=T\text{\hspace{0.17em}exp\hspace{0.17em}}\left[-\frac{\text{i}\Delta t}{\hslash }{\displaystyle \int \text{d}}{t}^{\prime }\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+mV\left(x,t^{\prime }\right)\right)\right]\Psi \left(x,t\right)& (22)\end{array}$

where T is the time-ordering symbol. For a sufficiently small time step, it can be approximated as

$\begin{array}{l}\Psi \left(x,t+\Delta t\right)=\text{exp\hspace{0.17em}}(\frac{\text{i}\hslash \Delta t}{2m}{\nabla }^2-\frac{\text{i}m\Delta t}{2\hslash }V\left(x,t+\Delta t\right)\\ -\frac{\text{i}m\Delta t}{2\hslash }V\left(x,t\right))\Psi \left(x,t\right),& (23)\end{array}$

which can be further splitted into three operations according to the Baker–Campbell–Hausdorff formula:

$\begin{array}{l}\Psi \left(x,t+\Delta t\right)=\text{exp\hspace{0.17em}}\left(-\frac{\text{i}m\Delta t}{2\hslash }V\left(x,t+\Delta t\right)\right)\text{exp\hspace{0.17em}}\left(\frac{\text{i}\hslash \Delta t}{2m}{\nabla }^2\right)\\ \text{ exp }\left(-\frac{\text{i}m\Delta t}{2\hslash }V\left(x,t\right)\right)\Psi \left(x,t\right).& (24)\end{array}$

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:

$\begin{array}{ll}D\left(\frac{\Delta t}{2}\right)\Psi \left(x,t\right)=\text{IFFT}\left\{-\frac{\text{i}\hslash }{m}\frac{\Delta t}{2}{k}^2\text{FFT}\left[\Psi \left(x,t\right)\right]\right\}.& (25)\end{array}$

The main differences between the Eulerian based and Lagrangian based methods are that the original Schrödinger-Poisson equations are solved in the former, not the transformed Madelung equations as in the latter, and the Eulerian method can be used to reliably estimate second order derivatives of the fields. From Table 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 10Mpc/h.

TABLE 2

Table 2. Summary of the Eulerain based simulation codes for FDM model (Schrödinger-Poisson solvers).

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 gun features of the FDM model, and it is very important to understand their properties. On the other hand, cosmological scale simulations are needed whenever large scale survey data are used to constrain the FDM model or the properties of galaxy cluster are studied.

To simplify the generation of the granular structures, a self-consistent method was introduced in Lin et al. (2018). 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}_{\odot }$) 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}_{\odot }$), 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 De Martino et al. (2017), Khmelnitsky and Rubakov (2014), and Porayko and Postnov (2014) with Pulsar Timing Array. By studying the modulation of the arriving time of pulses from many pulsars, the Pulsar Timing Array can be used to directly detect the granularity of dark matter distribution in the Milky Way galaxy. The Parkes Pulsar Timing Array collaboration has obtained the first constraints of FDM as m > 10⁻²³eV (Porayko et al., 2018). This method can not only set constraints on the rest mass of FDM particles, but also confidently claim the existence of FDM. 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 other smoking gun feature of the FDM model is the existence of solitonic cores in the dark matter halos. Different simulation groups reported that they have found solitonic cores (Schive et al., 2014a,b; Schwabe et al., 2016; Mocz et al., 2017). The solitonic core solution can also be obtained analytically (Chavanis, 2011; Chavanis and Delfini, 2011). The numerical simulations (Schive et al., 2014a,b) provide the empirical core-halo mass relation,

$\begin{array}{l}{x}_c\approx 160{(\frac{M_h}{{10}^{12}{M}_{\odot }})}^{-\frac{1}{3}}{(\frac{m}{{10}^{-22}\text{eV}})}^{-1}\text{pc},\\ \rho (x)\approx \frac{190{(\frac{m}{{10}^{-22}\text{eV}})}^{-2}{(\frac{x_c}{100\text{pc}})}^{-4}}{{(1+0.091{(\frac{x}{x_c})}^2)}^8}{\text{M}}_{\odot }{\text{pc}}^{-3},& (26)\end{array}$

where x_c is the solitonic core radius, m is the FDM particle mass and M_h is the halo mass. The radius of the solitonic core is defined as the radius where the mass density drops by a factor of 2 from its value at the origin. This relation can be used to compare with observations. It is believed that the dark matter halos made of FDM should host a solitonic core in the center and follow usual NFW profile in the outer region. However, how the density profile transfers from soliton to the outer NFW profile is not consistent among the simulations (Bar et al., 2018). It was claimed that the central solitonic core profile with outer NFW profile provide better fits than CDM predictions for the dwarf galaxies (especially for Fornax) (Schive et al., 2014a). However, it faces challenge from the analysis of rotation curves of many other galaxies (Bar et al., 2018).

FDM simulations with the largest box size (50Mpc/h) were performed in Zhang et al. (2018a). 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 Zhang et al. (2018a) and Nori and Baldi (2018), 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 can 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 (~ 500Mpc/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 De Martino et al. (2018) that the dynamical evidence of the existence of the solitonic core in the center of Milky Way was found. It was also claimed in Schive et al. (2014b) that the existence of the solitonic core can solve the cusp-core problem. Both of these two studies favor a FDM with m ~ 10⁻²²eV. However, the recent constraints from Lyman-alpha forest implied that FDM with m < 10⁻²¹eV was ruled out (Armengaud et al., 2017; Iršič et al., 2017). Such a tension is a problem for FDM models. Zhang et al. (2018a) suggested that considering of the important role of quantum pressure and systematic uncertainties may relieve this tension, but Nori et al. (2018) with further analysis concluded that the quantum pressure can not affect the lyman-alpha forest significantly. Hui et al. (2017) suggested that other astrophysical processes like patchy reionization can relieve the tension. Therefore, FDM with m ~ 10⁻²²eV is still not conclusively ruled out by the observations of Lyman-alpha forest, while 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⁻²¹eV and m > 1.5 × 10⁻²²eV (Amorisco and Loeb, 2018; Lidz and Hui, 2018), respectively. A recent study of rotation curves of near-by galaxies also claimed that FDM with m < 10⁻²¹eV is not favored (Bar et al., 2018). These independent constraints using different methods are also not supporting the value m ~ 10⁻²²eV needed to solve the small-scale crisis. If more arguments and modifications to the current FDM model are made, we can relieve the tensions but lose the beauty and simplicity of the FDM model.

We notice that with m > 10⁻²²eV, FDM model may still be able to solve the small-scale crisis due to the runaway tidal disruption. The runaway tidal disruption 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 (Du et al., 2018). This effect has not been studied by other codes yet, whether it is a physical effect or numerical illusion is unknown. A systematic study however, considering such mechanism needs simulations with sufficiently large box size (> 100Mpc/h) and high resolution ($<10^9{M}_{\odot }/h$) at the same time. This is very difficult to reach now.

4. Summary

In this paper, we have reviewed the basic idea of Fuzzy Dark Matter (FDM) model and the current status of simulations for this model. As a mini review, this work provides a short summary for the readers to follow the state-of-the-art research of FDM.

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 systematic uncertainties 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 ~ 500Mpc/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 observable 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 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. The FDM model is a beautiful model with no more free parameters than the CDM model together with the WIMP assumption. If we can determine the mass range of FDM particles, it will significantly improve our understanding of dark matter and the universe.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

JZ acknowledges the support from China Postdoctoral Science Foundation 2018M632097.

References

Ade, P., Aghanim, N., Arnaud, M., Ashdown, M., Aumont, J., Baccigalupi, C., et al. (2016). Planck 2015 results-xiii. cosmological parameters. Astron. Astrophys. 594:A13. doi: 10.1051/0004-6361/201525830

CrossRef Full Text | Google Scholar

Amorisco, N. C., and Loeb, A. (2018). First constraints on Fuzzy Dark Matter from the dynamics of stellar streams in the Milky Way. arXiv:1808.00464.

Google Scholar

Armengaud, E., Palanque-Delabrouille, N., Yèche, C., Marsh, D. J., and Baur, J. (2017). Constraining the mass of light bosonic dark matter using sdss lyman-α forest. arXiv 1703.09126. doi: 10.1093/mnras/stx1870

CrossRef Full Text | Google Scholar

Bahcall, N. A., Lubin, L. M., and Dorman, V. (1995). Where is the dark matter? Astrophys. J. Lett. 447:L81. doi: 10.1086/309577

CrossRef Full Text | Google Scholar

Bar, N., Blas, D., Blum, K., and Sibiryakov, S. (2018). Galactic rotation curves vs. ultra-light dark matter: implications of the soliton – Host Halo relation. Phys. Rev. D 98:083027. doi: 10.1103/PhysRevD.98.083027

CrossRef Full Text | Google Scholar

Bozek, B., Marsh, D. J., Silk, J., and Wyse, R. F. (2015). Galaxy uv-luminosity function and reionization constraints on axion dark matter. Month. Notices R. Astron. Soc. 450, 209–222. doi: 10.1093/mnras/stv624

CrossRef Full Text | Google Scholar

Chavanis, P.-H. (2011). Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. I. Analytical results. Phys. Rev. D 84:043531. doi: 10.1103/PhysRevD.84.043531

CrossRef Full Text | Google Scholar

Chavanis, P.-H., and Delfini, L. (2011). Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. ii. numerical results. Phys. Rev. D 84:043532. doi: 10.1103/PhysRevD.84.043532

CrossRef Full Text | Google Scholar

Chavanis, P. H. (2012). Growth of perturbations in an expanding universe with Bose-Einstein condensate dark matter. Astron. Astrophys. 537:A127. doi: 10.1051/0004-6361/201116905

CrossRef Full Text | Google Scholar

Cheng, D., Chu, M.-C., and Tang, J. (2015). Cosmological structure formation in Decaying Dark Matter models. J. Cosmol. Astropart. Phys. 7:9. doi: 10.1088/1475-7516/2015/07/009

CrossRef Full Text | Google Scholar

Clowe, D., Bradac, M., Gonzalez, A., Markevitch, M., Randall, S., Jones, C., et al. (2006). A direct empirical proof of the existence of dark matter. Astrophys. J. Lett. 648, L109–L113. doi: 10.1086/508162

CrossRef Full Text | Google Scholar

Colin, P., Avila-Reese, V., and Valenzuela, O. (2000). Substructure and halo density profiles in a warm dark matter cosmology. Astrophys. J. 542:622. doi: 10.1086/317057

CrossRef Full Text | Google Scholar

De Martino, I., Broadhurst, T., Tye, S.-H. H., Chiueh, T., and Schive, H.-Y. (2018). Dynamical evidence of a solitonic core of 10⁹M_⊙ in the Milky Way. arXiv:1807.08153.

Google Scholar

De Martino, I., Broadhurst, T., Tye, S.-H. H., Chiueh, T., Schive, H.-Y., and Lazkoz, R. (2017). Recognising axionic dark matter by compton and de-Broglie scale modulation of pulsar timing. Phys. Rev. Lett. 119:221103. doi: 10.1103/PhysRevLett.119.221103

CrossRef Full Text | Google Scholar

Dehnen, W. (2001). Towards optimal softening in three-dimensional N-body codes - I. Minimizing the force error. Month. Notices R. Astron. Soc. 324, 273–291. doi: 10.1046/j.1365-8711.2001.04237.x

CrossRef Full Text | Google Scholar

Du, X., Behrens, C., and Niemeyer, J. C. (2017). Substructure of fuzzy dark matter halos. Month. Notices R. Astron. Soc. 465, 941–951. doi: 10.1093/mnras/stw2724

CrossRef Full Text | Google Scholar

Du, X., Schwabe, B., Niemeyer, J. C., and Bürger, D. (2018). Tidal disruption of fuzzy dark matter subhalo cores. Phys. Rev. D 97:063507. doi: 10.1103/PhysRevD.97.063507

CrossRef Full Text | Google Scholar

Edwards, F., Kendall, E., Hotchkiss, S., and Easther, R. (2018). PyUltraLight: a pseudo-spectral solver for ultralight dark matter dynamics. J. Cosmol. Astropart. Phys. 10:027 doi: 10.1088/1475-7516/2018/10/027

CrossRef Full Text | Google Scholar

Einasto, J., Kaasik, A., and Saar, E. (1974). Dynamical evidence on massive coronas of galaxies. Nature 250:309. doi: 10.1038/250309a0

CrossRef Full Text | Google Scholar

Hlozek, R., Grin, D., Marsh, D. J., and Ferreira, P. G. (2015). A search for ultralight axions using precision cosmological data. Phys. Rev. D 91:103512. doi: 10.1103/PhysRevD.91.103512

CrossRef Full Text | Google Scholar

Hu, W., Barkana, R., and Gruzinov, A. (2000). Cold and fuzzy dark matter. Phys. Rev. Lett. 85, 1158–1161. doi: 10.1103/PhysRevLett.85.1158

PubMed Abstract | CrossRef Full Text

Hui, L., Ostriker, J. P., Tremaine, S., and Witten, E. (2017). Ultralight scalars as cosmological dark matter. Phys. Rev. D 95:043541. doi: 10.1103/PhysRevD.95.043541

CrossRef Full Text | Google Scholar

Iršič, V., Viel, M., Haehnelt, M. G., Bolton, J. S., and Becker, G. D. (2017). First constraints on fuzzy dark matter from lyman-α forest data and hydrodynamical simulations.

PubMed Abstract | Google Scholar

Khlopov, M. I., Malomed, B. A., and Zeldovich, I. B. (1985). Gravitational instability of scalar fields and formation of primordial black holes. Month. Notices RAS 215, 575–589. doi: 10.1093/mnras/215.4.575

CrossRef Full Text | Google Scholar

Khmelnitsky, A., and Rubakov, V. (2014). Pulsar timing signal from ultralight scalar dark matter. J. Cosmol. Astropart. Phys. 2:019. doi: 10.1088/1475-7516/2014/02/019

CrossRef Full Text | Google Scholar

Kobayashi, T., Murgia, R., De Simone, A., Iršič, V., and Viel, M. (2017). Lyman-alpha constraints on ultralight scalar dark matter: implications for the early and late universe. Phys. Rev. D 96:123514. doi: 10.1103/PhysRevD.96.123514

CrossRef Full Text | Google Scholar

Lee, J.-W. (2017). “Brief history of ultra-light scalar dark matter models” in EPJ Web of Conferences. doi: 10.1051/epjconf/201816806005

CrossRef Full Text | Google Scholar

Lidz, A., and Hui, L. (2018). Implications of a prereionization 21-cm absorption signal for fuzzy dark matter. Phys. Rev. D 98:023011. doi: 10.1103/PhysRevD.98.023011

CrossRef Full Text | Google Scholar

Lin, S.-C., Schive, H.-Y., Wong, S.-K., and Chiueh, T. (2018). Self-consistent construction of virialized wave dark matter halos. Phys. Rev. D 97:103523. doi: 10.1103/PhysRevD.97.103523

CrossRef Full Text | Google Scholar

Marsh, D. J. (2016). Axion cosmology. Phys. Rep. 643, 1–79. doi: 10.1016/j.physrep.2016.06.005

CrossRef Full Text | Google Scholar

Marsh, D. J. E. (2015). Nonlinear hydrodynamics of axion dark matter: relative velocity effects and quantum forces. Phys. Rev. D 91:123520. doi: 10.1103/PhysRevD.91.123520

CrossRef Full Text | Google Scholar

Merritt, D. (1996). Optimal smoothing for N-body codes. Astron. J. 111:2462. doi: 10.1086/117980

CrossRef Full Text | Google Scholar

Mishra, S. S., Sahni, V., and Shtanov, Y. (2017). Sourcing dark matter and dark energy from α-attractors. J. Cosmol. Astropart. Phys. 6:045. doi: 10.1088/1475-7516/2017/06/045

CrossRef Full Text | Google Scholar

Mocz, P., and Succi, S. (2015). Numerical solution of the nonlinear schrödinger equation using smoothed-particle hydrodynamics. Phys. Rev. E 91:053304. doi: 10.1103/PhysRevE.91.053304

PubMed Abstract | CrossRef Full Text | Google Scholar

Mocz, P., Vogelsberger, M., Robles, V., Zavala, J., Boylan-Kolchin, M., and Hernquist, L. (2017). Galaxy formation with BECDM: I. Turbulence and relaxation of idealised haloes. Month. Notices R. Astron. Soc. 471, 4559–4570. doi: 10.1093/mnras/stx1887

CrossRef Full Text

Nori, M., and Baldi, M. (2018). AX-GADGET: a new code for cosmological simulations of Fuzzy Dark Matter and Axion models. arXiv: 1801.08144.

Google Scholar

Nori, M., Murgia, R., Iršič, V., Baldi, M., and Viel, M. (2018). Lyman-α forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies. Month. Notices RAS 482, 3227–3243. doi: 10.1093/mnras/sty2888

CrossRef Full Text | Google Scholar

Pontzen, A., and Governato, F. (2014). Cold dark matter heats up. Nature 506, 171–178. doi: 10.1038/nature12953

PubMed Abstract | CrossRef Full Text | Google Scholar

Porayko, N. K., and Postnov, K. A. (2014). Constraints on ultralight scalar dark matter from pulsar timing. Phys. Rev. D 90:062008. doi: 10.1103/PhysRevD.90.062008

CrossRef Full Text | Google Scholar

Porayko, N. K., Zhu, X., Levin, Y., Hui, L., Hobbs, G., Grudskaya, A., et al. (2018). Parkes Pulsar Timing Array constraints on ultralight scalar-field dark matter. Phys. Rev. D 98:102002. doi: 10.1103/PhysRevD.98.102002

CrossRef Full Text | Google Scholar

Sahni, V., and Wang, L. (2000). New cosmological model of quintessence and dark matter. Phys. Rev. D 62:103517. doi: 10.1103/PhysRevD.62.103517

CrossRef Full Text | Google Scholar

Schive, H.-Y., Chiueh, T., and Broadhurst, T. (2014a). Cosmic structure as the quantum interference of a coherent dark wave. Nat. Phys. 10, 496–499. doi: 10.1038/nphys2996

CrossRef Full Text | Google Scholar

Schive, H. Y., Liao, M. H., Woo, T. P., Wong, S. K., Chiueh, T., Broadhurst, T., et al. (2014b). Understanding the core-halo relation of quantum wave dark matter from 3d simulations. Phys. Rev. Lett. 113:261302. doi: 10.1103/PhysRevLett.113.261302

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwabe, B., Niemeyer, J. C., and Engels, J. F. (2016). Simulations of solitonic core mergers in ultralight axion dark matter cosmologies. Phys. Rev. D 94:043513. doi: 10.1103/PhysRevD.94.043513

CrossRef Full Text | Google Scholar

Silverman, B. W. (1978). Choosing the window width when estimating a density. Biometrika 65, 1–11. doi: 10.1093/biomet/65.1.1

CrossRef Full Text | Google Scholar

Spergel, D. N., and Steinhardt, P. J. (2000). Observational evidence for self-interacting cold dark matter. Phys. Rev. Lett. 84, 3760–3763. doi: 10.1103/PhysRevLett.84.3760

PubMed Abstract | CrossRef Full Text | Google Scholar

Spiegel, E. (1980). Fluid dynamical form of the linear and nonlinear schrödinger equations. Physica D 1, 236–240. doi: 10.1016/0167-2789(80)90015-9

CrossRef Full Text | Google Scholar

Springel, V., and Hernquist, L. (2002). Cosmological smoothed particle hydrodynamics simulations: the entropy equation. Month. Notices RAS 333, 649–664. doi: 10.1046/j.1365-8711.2002.05445.x

CrossRef Full Text | Google Scholar

Tegmark, M., Blanton, M. R., Strauss, M. A., Hoyle, F., Schlegel, D., Scoccimarro, R., et al. (2004). The three-dimensional power spectrum of galaxies from the sloan digital sky survey. Astrophys. J. 606:702. doi: 10.1086/382125

CrossRef Full Text | Google Scholar

Uhlemann, C., Kopp, M., and Haugg, T. (2014). Schrödinger method as N-body double and UV completion of dust. Phys. Rev. D 90:023517. doi: 10.1103/PhysRevD.90.023517

CrossRef Full Text | Google Scholar

Veltmaat, J., and Niemeyer, J. C. (2016). Cosmological particle-in-cell simulations with ultra-light axion dark matter. Phys. Rev. D 94:123523. doi: 10.1103/PhysRevD.94.123523

CrossRef Full Text | Google Scholar

Wallstrom, T. C. (1994). Inequivalence between the Schrodinger equation and the Madelung hydrodynamic equations. Phys. Rev. A 49, 1613–1617. doi: 10.1103/PhysRevA.49.1613

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberg, D. H., Bullock, J. S., Governato, F., de Naray, R. K., and Peter, A. H. (2015). Cold dark matter: controversies on small scales. Proc. Natl. Acad. Sci. U.S.A. 112, 12249–12255. doi: 10.1073/pnas.1308716112

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Kuo, J.-L., Liu, H., Sming Tsai, Y.-L., Cheung, K., and Chu, M.-C. (2018a). The importance of quantum pressure of fuzzy dark matter on Lyα forest. Astrophys. J. 863:73. doi: 10.3847/1538-4357/aacf3f

CrossRef Full Text | Google Scholar

Zhang, J., Sming Tsai, Y.-L., Kuo, J.-L., Cheung, K., and Chu, M.-C. (2018b). Ultralight axion dark matter and its impact on dark halo structure in N-body simulations. Astrophys. J. 853:51. doi: 10.3847/1538-4357/aaa485

CrossRef Full Text | Google Scholar