The Importance of Quantum Pressure of Fuzzy Dark Matter on Lyα Forest

Dark matter is one of the intriguing mysteries of modern cosmology. Currently, the leading paradigm of dark matter is the cold dark matter (CDM), supported by the majority of the observations like 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). Despite its success on large scales, the CDM paradigm faces three problems on small scales, dubbed the "small scales crisis" (Weinberg et al. 2015): (i) the missing satellite problem, (ii) the cusp-core problem, and (iii) the too-big-to-fail problem. The essence of these problems is that CDM predicts an excess amount of dark matter on small scales, hence the key to addressing them is to smooth out the small-scale structures by astrophysical processes (Pontzen & Governato 2014) or invoking alternative dark matter models such as warm dark matter (WDM) (Colin et al. 2000), self-interacting dark matter (Spergel & Steinhardt 2000), and fuzzy dark matter (FDM; Hu et al. 2000).

The FDM paradigm, in which the dark matter is made of ultra-light bosons in Bose–Einstein condensate state, is an ideal alternative of CDM, as it suppresses small-scale structures while it inherits the success of CDM on large scales (Marsh 2016; Du et al. 2017; Mocz et al. 2017). For the detailed history and implementation of FDM, one can see  (Lee 2018) and the reference therein. The suppression effect, arising from the effective "quantum pressure" (QP) of FDM, is directly connected to the mass of the ultra-light axion. The predictions of FDM with mass $\sim {10}^{-22}\,\mathrm{eV}$ are consistent with observations of the CMB and large-scale structures (Hlozek et al. 2015); high-z galaxies and CMB optical depth (Bozek et al. 2015); and the density profiles of dwarf spheroidal galaxies (Schive et al. 2014). However, recent results claimed that FDM with mass below ${10}^{-21}\,\mathrm{eV}$ has already been ruled out at 95% confidence level by comparing the results of hydrodynamic simulations to the Lyα forest data (Armengaud et al. 2017; Iršič et al. 2017; Kobayashi et al. 2017).

Lyα forest, a series of absorption lines in the Lyα emission spectrum from distant galaxies and quasars by neutral hydrogen (H i) gas at different redshift, provides the information about the the spatial distribution of H i gas at high redshift. In the leading theory of Lyα forest—the gravitational confinement model—H i clumps are confined by the gravity provided by dark matter halos (Ikeuchi 1986; Rees 1986). Thus, the flux power spectrum of Lyα forest is a biased representation of the underling DM density field power spectrum (Weinberg et al. 2003; Viel et al. 2004; Arinyo-i Prats et al. 2015). Recently, two collaborations, the Baryon Oscillation Spectroscopic Survey (BOSS; Palanque-Delabrouille et al. 2013) and XQ-100 (López et al. 2016), have announced their analyses between redshift z = 2–5 for the flux power spectrum, thus providing tools to constraint dark matter models on an unprecedented high level of precision.

To robustly exclude such a mass range, two tasks have to be carried out before the experimental data analysis. The first is numerical simulations of the FDM system on large scales, which have been performed in a number of approaches: directly solving the Schrödinger–Poission system (Schive et al. 2014), smoothed particle hydrodynamic (SPH) simulation (Mocz & Succi 2015), and N-body simulation with particle-mesh (PM) method (Veltmaat & Niemeyer 2016). However, all of these simulations were restricted in simulation scale and suffered from the singularity problem at zero-density points in the calculation of QP. In this paper, we adopt an independent simulation scheme developed in Zhang et al. (2018), which provides DM-only simulations with SPH and particle–particle (PP) interactions to account for QP, is feasible for implementation in cosmological scale simulations, and avoids the singularity problem.

The second is the uncertainties in hydrodynamic simulations. In the N-body simulation, the uncertainty in the matter power spectrum is ${ \mathcal O }(10 \% )$ originating from the use of different initial conditions, the adoption of different orders of Lagrangian perturbation theory (L'Huillier et al. 2014), finite-box effect (Bagla & Prasad 2006) and the usage of different N-body simulation codes (Heitmann et al. 2010; Kim et al. 2013; Sembolini et al. 2016). These uncertainties can be precisely estimated and controlled because their properties are well studied. In contrast, the uncertainties from hydrodynamic simulations can be of orders of magnitude (Oshea et al. 2005; Vazza et al. 2011; Kim et al. 2013; Sembolini et al. 2016). Different treatments introduced in hydrodynamic simulations in handling the gaseous part induce huge discrepancy in gas density, gas temperature, and galaxy structure. Moreover, different hydrodynamic simulation codes contribute uncertainties to the 1D flux power spectrum at the level of 5%. As stated in Hui et al. (2017), several additional astrophysical processes may also alter the 1D flux power spectrum, such as patchy reionization and smoothing.

We estimate the constraint on the mass range of FDM from Lyα forest data with DM-only cosmological simulations in four different settings: (1) standard CDM simulation, (2) simulation using FDM-modified initial condition and CDM dynamics, (3) simulation using FDM-modified initial condition and dynamics with FDM mass equals to $2.5\times {10}^{-22}\,\mathrm{eV}$, (4) $2.5\times {10}^{-23}\,\mathrm{eV}$. By utilizing the 1D flux power spectrum of the CDM hydrodynamic simulation in Iršič et al. (2017), one can obtain the 1D flux power spectrum for the corresponding condition using the 3D power spectrum ratio from our simulations.

Our results are shown in the aspects of density field, halo mass function, 3D power spectrum, and 1D flux power spectrum with Lyα forest data. We show that the effect of QP in structure formation is non-trivial and introduce important suppression in the growth of structures. Inspired by the linear theory summarized in Mo et al. (2010), we have calculated the uncertainty range of the 1D flux power spectrum based on the unknown hydrogen gas temperature properties. The suppression of the 1D flux power spectrum from higher gas temperature behaves like the suppression caused by FDM.

This paper is arranged as follows. In Section 2, we discuss the methodology of the simulation, simulation set-up and initial conditions. We also review the linear theory of 1D flux power spectrum calculation. In Section 3, we present our results in the aspects of density field, halo mass function, 3D power spectrum, and 1D flux power spectrum. We compare our 1D flux power spectrum to the Lyα forest data. In Section 4, we discuss the gas temperature uncertainty. Finally, we summarize our outcomes in Section 5.

2.1. QP as PP Interaction

Standard N-body simulations have problems in calculating the QP because of the discretization of the density field by the delta function. From the expression of QP, obviously, what we need in the calculation is a smooth density field. There are already some reliable smoothing particle methods used in the N-body simulation of the FDM system with limited box sizes ($\lt 10\,{h}^{-1}\,\mathrm{Mpc}$), (e.g., Mocz & Succi 2015; Veltmaat & Niemeyer 2016). For the sake of Lyα forest, it is important to have a cosmological scale simulation so that the study of structure formation and matter power spectrum on large scales is possible.

We proposed a novel N-body simulation scheme for FDM in a previous work (Zhang et al. 2018), in which the delta function is replaced by a smooth Gaussian kernel function to solve the singularity problem and avoid the simulation crash. This scheme can be used to perform simulations at scales no smaller than $50\,{h}^{-1}\,\mathrm{Mpc}$ as we shall prove later and give a coarse-grained description on small scales. However, the original scheme requires some modifications for cosmological simulations. In this section, we will demonstrate how to embed QP in the cosmological simulations.

For cosmological simulations, one has to consider the transformation from physical to comoving coordinates. Under the transformation described in Appendix A, there is an additional pre-factor ${a}^{-2}$ for the original QP defined in Zhang et al. (2018), where a is the cosmological scale factor. The QP for a cosmological simulation becomes

$$\begin{eqnarray}&&Q=-\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{\chi }^{2}{a}^{2}}\displaystyle \frac{{{\boldsymbol{\nabla }}}^{2}\sqrt{\rho }}{\sqrt{\rho }},\end{eqnarray} \tag{ 1 }$$

where ℏ, ${m}_{\chi }$, and ρ are the reduced Planck constant, FDM particle mass, and the mass density of FDM, respectively. The corresponding acceleration can be written as

$$\begin{eqnarray}\ddot{{\boldsymbol{r}}} & = & \displaystyle \frac{4M{{\hslash }}^{2}}{{M}_{0}{m}_{\chi }^{2}{\lambda }^{4}{a}^{2}}\displaystyle \sum _{j}{{ \mathcal B }}_{j}\exp \left[-\displaystyle \frac{2| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{j}{| }^{2}}{{\lambda }^{2}}\right]\\ & & \times \,\left(1-\displaystyle \frac{2| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{j}{| }^{2}}{{\lambda }^{2}}\right)({{\boldsymbol{r}}}_{j}-{\boldsymbol{r}}),\end{eqnarray} \tag{ 2 }$$

where M, M₀, λ, and ${{ \mathcal B }}_{j}$ are the mass of the simulation particle, a normalization factor accounting for the volume ${\rm{\Delta }}{V}_{j}$ occupied by simulation particles, the de Broglie wavelength of FDM particles, and the correction factor for high-density regions, respectively. For more detailed explanations of ${{ \mathcal B }}_{j}$, please refer to Appendix A of Zhang et al. (2018).

To demonstrate the effect of QP, we consider a two-particle system separated by a distance of order ${ \mathcal O }(\mathrm{kpc})$, and the acceleration caused by QP will be ${ \mathcal O }({{\hslash }}^{2}\,{m}_{\chi }^{-2}\,{\lambda }^{-3})\sim { \mathcal O }({10}^{-10}\,{\rm{m}}\,{{\rm{s}}}^{-2})$. In Figure 1, the QP effect is presented in the plane of (r, $\ddot{r}$). The acceleration from QP, gravity, and their sum are shown by the black dashed line, blue solid line, and red solid line, respectively. Apparently, the effect of QP is attractive if the separation between the two particles is shorter than $\lambda /\sqrt{2}$, otherwise repulsive. To understand this phenomenon, we refer back to the definition of QP in Equation (1): Q is proportional to the curvature of the density, which can be negative, positive, or zero, physically corresponding to negative, positive, and zero forces, respectively. The reason of the difference between Figure 1 here and the Figure 1 of Zhang et al. (2018) is that the mass of the simulation particle is changed from ${10}^{6}{M}_{\odot }$ to $\sim {10}^{8}{M}_{\odot }$. The effect of gravity is directly related to the mass of the simulation particles while the QP is not affected due to the fact that the normalization factor M₀ in Equation (2) is chosen to be the same as the mass of the simulation particles. In other words, the strength of QP is a quantity independent of the simulation particle mass, M. Therefore, the acceleration of gravity is considerably larger in a cosmological simulation, due to the larger mass of the simulation particles. Thus, the structure formation in large scales is still dominated by gravity. However, the QP plays a non-trivial role in the structure formation, particularly in highly nonlinear regions such as dark matter halos.

Zoom In Zoom Out Reset image size

2.2. Simulation Settings

We use the code Gadget2 (Springel 2005), which is a TreePM hybrid N-body code, to perform our simulations. To describe QP, we discretize the interaction term and modify Gadget2⁶ to calculate QP in the same way as in Zhang et al. (2018). (The PM method is helpful for cosmological simulations with periodic boundary conditions.) Because QP behaves like a short-range force, we adopt the original PM code to compute the long-range force and modify the Tree code which takes care of the short-range force to include the calculation of QP. In addition, it is not necessary to set softening length for QP since QP is finite in the $({{\boldsymbol{r}}}_{j}-{\boldsymbol{r}})\simeq 0$ region.

We start our simulation from the redshift z = 99. The related cosmological parameters are DM energy density ${{\rm{\Omega }}}_{m}=0.3$, cosmological constant ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$, baryon energy density ${{\rm{\Omega }}}_{b}=0.04$, dimensionless Hubble parameter h = 0.7, scalar spectral index n_s = 0.96, and the power spectrum normalization factor σ₈ = 0.8.

We generated the initial CDM power spectrum following  Eisenstein & Hu (1998). The suppression of the FDM power spectrum relative to the CDM one on small scales can be characterized by a transfer function ${ \mathcal T }(k,z)$ (Hu et al. 2000)

$$\begin{eqnarray}&&{P}_{F}(k,z)={\left[\displaystyle \frac{{{ \mathcal T }}_{{\mathtt{FDM}}}(k,z)}{{{ \mathcal T }}_{{\mathtt{CDM}}}(k,z)}\right]}^{2}{P}_{C}(k,z)={{ \mathcal T }}^{2}(k,z){P}_{C}(k,z),\end{eqnarray} \tag{ 3 }$$

where k is the wavenumber and ${P}_{F}(k,z)$ and ${P}_{C}(k,z)$ are the 3D power spectra of FDM and CDM, respectively. One has to bear in mind that the physical difference between ${{ \mathcal T }}_{{\mathtt{CDM}}}(k,z)$ and ${{ \mathcal T }}_{{\mathtt{FDM}}}(k,z)$ ⁷ is characterized by the transfer function, ${ \mathcal T }(k,z)$, which transforms the power spectrum from CDM to FDM and so is deduced as the ratio of FDM transfer function ${{ \mathcal T }}_{{\mathtt{FDM}}}(k,z)$ to CDM transfer function ${{ \mathcal T }}_{{\mathtt{CDM}}}(k,z)$.

Following the arguments in Schive et al. (2016), we can well approximate ${ \mathcal T }(k,z)$ by the redshift-independent expression (Hu et al. 2000)

$$\begin{eqnarray}{ \mathcal T }(k) & = & \displaystyle \frac{\cos {x}^{3}}{1+{x}^{8}},\,\,\mathrm{where}\,x=1.61\\ & & \times \,{\left(\displaystyle \frac{{m}_{\chi }}{{10}^{-22}\mathrm{eV}}\right)}^{1/18}\times \displaystyle \frac{k}{{k}_{J}}.\end{eqnarray} \tag{ 4 }$$

The parameter ${k}_{J}=9{({m}_{\chi }/{10}^{-22}\mathrm{eV})}^{1/2}\,{\mathrm{Mpc}}^{-1}$ is the critical scale of Jeans wavenumber at matter-radiation equality. In Figure 2, we present the square of the transfer function ${{ \mathcal T }}^{2}(k)$ with three different FDM masses to demonstrate the suppression of FDM power spectrum in small scales relative to CDM power spectrum. The blue, black, and red solid lines correspond to the masses of the FDM being $2.5\times {10}^{-21}\,\mathrm{eV}$, $2.5\times {10}^{-22}\,\mathrm{eV}$, $2.5\times {10}^{-23}\,\mathrm{eV}$, respectively. The vertical green dashed line represents the Nyquist limit which is the resolution limit of our simulations and the corresponding wavenumber k_Ny is computed to be

$$\begin{eqnarray*}&&{k}_{{Ny}}=\pi {\left(\displaystyle \frac{{N}_{0}}{{V}_{0}}\right)}^{1/3},\end{eqnarray*}$$

where V₀ is the volume of the simulation box, and N₀ is the total number of simulation particles, which are ${(50{h}^{-1}\mathrm{Mpc})}^{3}$ and 512³ respectively in our simulations. The origin of this limit is due to the fact that each simulation particle has an approximate average volume ${({V}_{0}/{N}_{0})}^{1/3}$ and we cannot know what happens inside the simulation particle. In other words, the result is not reliable for $k\gt {k}_{{Ny}}$. One can see that there is a sharp break of FDM power spectrum at $k\sim {k}_{J}$, and severe oscillations occur with suppression for $k\gt {k}_{J}$ for all three different FDM masses. For a smaller mass, the suppression is tremendous compared to a larger mass which behaves almost the same as the CDM case at high k.

Zoom In Zoom Out Reset image size

We modify the code 2LPTic (Crocce et al. 2006) for generating the initial conditions of the CDM power spectrum for cosmological simulations to incorporate the transfer function so as to generate initial conditions with the FDM power spectrum.

2.3. 1D Flux Power Spectrum

The calculation of the 1D flux power spectrum in the linear regime is well summarized in Mo et al. (2010),

$$\begin{eqnarray}&&{P}_{b}=\displaystyle \frac{{P}_{\mathrm{DM}}}{{(1+{k}^{2}/{k}_{J}^{2})}^{2}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{P}_{F}({k}_{z},z)={\int }_{{k}_{z}}^{\infty }\displaystyle \frac{{kdk}}{2\pi }{P}_{b}(k,z)W(k,{k}_{z}).\end{eqnarray} \tag{ 6 }$$

First, the linear dark matter power spectrum P_DM was calculated by CAMB (Lewis & Bridle 2002), with the same parameters used in the simulations. Then, we calculate the baryon power spectrum from Equation (5), in which k_J is the Jeans wavenumber related to the Jeans length ${\lambda }_{J}$ by

$$\begin{eqnarray}&&{k}_{J}=\displaystyle \frac{2\pi }{{\lambda }_{J}}.\end{eqnarray} \tag{ 7 }$$

Finally, we integrate Equation (6) to get the 1D flux power spectrum, in which we need the bias function

$$\begin{eqnarray}W & = & A\exp \left(-\displaystyle \frac{{k}_{z}^{2}{b}_{0}^{2}}{2{H}^{2}}\right)\left[1+\displaystyle \frac{{{\rm{\Omega }}}_{m}^{0.6}}{2+0.7(1-\gamma )}\right.\\ & & \times \,{\left.\displaystyle \frac{{k}_{z}^{2}}{{k}^{2}}-\displaystyle \frac{\gamma -1}{4[2+0.7(1-\gamma )]}\displaystyle \frac{{k}_{z}^{2}{b}_{0}^{2}}{{H}^{2}}\right]}^{2}.\end{eqnarray} \tag{ 8 }$$

Here,

$$\begin{eqnarray*}&&{b}_{0}(z)=\sqrt{\displaystyle \frac{2{k}_{B}{T}_{0}(z)}{{m}_{H}}}\end{eqnarray*}$$

is a parameter related to the velocity dispersion of the H i gas at redshift z, and ${T}_{0}(z)$ is the average temperature of H i gas at redshift z. $\gamma (z)$ is the polytropic index in the equation of state of H i gas.⁸ A is a normalizing factor, whose effect is the same as the parameter σ₈ in P_DM.

For further calculation, we use a power-law parametrization of T(z) and $\gamma (z)$:

$$\begin{eqnarray}&&{T}_{0}(z)={T}_{0}^{A}{\left[\displaystyle \frac{1+z}{5.5}\right]}^{{T}_{0}^{S}},\quad \gamma (z)={\gamma }^{A}{\left[\displaystyle \frac{1+z}{5.5}\right]}^{{\gamma }^{S}},\end{eqnarray} \tag{ 9 }$$

where the 1σ uncertainty ranges are

$$\begin{eqnarray}{T}_{0}^{A} & = & {9.2}_{-0.1}^{+1.2}{10}^{3}\,{\rm{K}},\quad {T}_{0}^{S}=-{2.5}_{-0.5}^{+0.45},\\ {\gamma }^{A} & = & {1.64}_{-0.26}^{+0.01},\quad {\gamma }^{S}=-{0.15}_{-0.61}^{+1.25}.\end{eqnarray} \tag{ 10 }$$

These are copied from the Table  II in Viel et al. (2013). As for σ₈ and ${\lambda }_{J}$, we treat them as redshift-independent constants with the 1σ uncertainty ranges as in Ade et al. (2016) and Rorai et al. (2017):

$$\begin{eqnarray}&&{\sigma }_{8}\in [0.78,0.88],\quad {\lambda }_{J}=100\pm 80\,\mathrm{kpc}.\end{eqnarray} \tag{ 11 }$$

The formulae discussed above are only applicable to linear perturbation, and the insertion of nonlinear matter power spectrum in Equation (6) is problematic. However, for the sake of simplicity, the 1D flux power spectrum is assumed to be related to the DM density power spectrum with the following equation:

$$\begin{eqnarray}&&{P}_{F}({k}_{z},z)={\int }_{{k}_{z}}^{\infty }\displaystyle \frac{{kdk}}{2\pi }{P}_{\mathrm{DM}}(k,z)W(k),\end{eqnarray} \tag{ 12 }$$

which is inspired by the linear theory. Here, W represents the bias introduced by the thermal properties of intergalatic medium (IGM), but the distortions in the direction of line of sight is neglected. Given the form of Equation (12), the 1D flux power spectrum is the 3D matter power spectrum convolved with the bias W along the line of sight. On one hand, we assume that the dynamics of dark matter does not affect the gaseous component, and that the bias W is the same in all four different simulations. On the other hand, we reweigh the matter power spectrum, P_DM, with the 3D matter power spectrum ratio of the simulations FDM/FIC/F23 to the simulation CDM. That is, we first take the derivative of the 1D flux power spectrum P_F, from the mock 1D flux power spectrum in hydrodynamic simulation shown in Iršič et al. (2017), with respect to k_z,

$$\begin{eqnarray}&&-\displaystyle \frac{2\pi }{{k}_{z}}\displaystyle \frac{{{dP}}_{F}}{{{dk}}_{z}}={P}_{\mathrm{DM}}W,\end{eqnarray} \tag{ 13 }$$

and multiply the right-hand side of this equation by the ratio ${P}_{i}/{P}_{\mathrm{CDM}}$, where i = FDM, FIC, and F23. Then, we integrate this expression to get the modified 1D flux power spectrum

$$\begin{eqnarray}&&{P}_{F,i}({k}_{z},z)={\int }_{{k}_{z}}^{\infty }\displaystyle \frac{{kdk}}{2\pi }\displaystyle \frac{{P}_{i}(k,z)}{{P}_{\mathrm{CDM}}(k,z)}{P}_{\mathrm{DM}}(k,z)W(k).\end{eqnarray} \tag{ 14 }$$

By comparing the modified ${P}_{F,i}$ and the original P_F, we can tell the effect of FDM model on Lyman alpha forest.

We have performed four different kinds of simulations for comparison as listed in Table 1. To avoid words cluttering in the following presentation, we use abbreviations for these four simulations. In this section, we first present the density field of the simulation FDM, as well as the density field difference between the simulations FDM and FIC. We then compare the 3D power spectra and 1D flux power spectra of all four simulations to investigate the nonlinear effect of the QP. After showing the 1D flux power spectrum, a numerical chi-square test of Lyα forest data from BOSS and XQ-100 is performed. Finally, we discuss the constraint on the FDM mass.

3.1. Density Field and Halo Mass Function

In the left panel of Figure 3, we depict the density field taken from a slice of the simulation cube. The slice is $0.5\,{h}^{-1}\,\mathrm{Mpc}$ thick, and the density field is calculated from the particle distribution by the triangular shaped cloud (TSC) scheme, but further smoothed by a Gaussian filter (variance σ = 0.15 h⁻¹ Mpc) for better illustration. With the color scale, one can clearly see that the large scales structures include voids, filaments, and knots. The density contrast ${\delta }_{q}$ is defined as

$$\begin{eqnarray*}&&{\delta }_{q}=\displaystyle \frac{{\rho }_{q}}{\bar{{\rho }_{q}}}-1,\end{eqnarray*}$$

where the subscript q denotes that it is from the simulation FDM. The color scale of Figure 3 represents the logarithm of ${\delta }_{q}+1$, which is the ratio of local density, ${\rho }_{q}$, to the average density, $\bar{{\rho }_{q}}$. In the right panel of Figure 3, we compare the density field of the simulation FDM and the simulation FIC. The color scale shows the difference between the logarithm of ${\rho }_{q}/{\bar{\rho }}_{q}$ and the logarithm of ${\rho }_{n}/{\bar{\rho }}_{n}$.

Zoom In Zoom Out Reset image size

We found that the difference between the density fields of the two simulations FDM and FIC is small, but the effect of QP can produce granular structures close to the high-density regions. Note that these granular structures are indeed due to QP, but not the difference in the initial conditions by comparing the density fields of the simulations FIC and CDM.

To illustrate more details in small scale, we select a 5 h⁻¹ Mpc sub-box projected in y direction and zoom in for better resolution. In Figure 4, the left panel is from the simulation FDM and the right panel is from the simulation FIC. The density field is calculated from the particle distribution by the nearest grid point scheme and further smoothed by bilinear algorithm. We note that the cluster in the simulation FDM looks much fuzzier than the one in the simulation FIC. Qualitatively, the formation of low-mass halos is further suppressed when taking QP into consideration.

Zoom In Zoom Out Reset image size

We identified the halos using the package ${\mathtt{AHF}}$ (Knollmann & Knebe 2009), which can build up the hierarchical structure for the halos and sub-halos in the snapshots of our simulations. The halos are identified if the average density of the halo is over 200 times the critical density of the universe. In Figure 5, the halo mass function is presented with the FDM mass ${m}_{\chi }=2.5\times {10}^{-22}\,\mathrm{eV}$. The colors represent different redshifts, and the solid, dotted–dashed and dashed lines represent the halo mass function of the simulations CDM, FIC, and FDM, respectively. The break of the halo mass functions at $M=2.3\times {10}^{9}\,{h}^{-1}\,{M}_{\odot }$ (the green dashed line) is due to our limited resolution. By definition, we cannot identify halos whose mass is below $2.3\times {10}^{9}\,h{}^{-1}{M}_{\odot }$. From halo masses $5\times {10}^{11}\,h{}^{-1}{M}_{\odot }$ to $2\times {10}^{13}\,h{}^{-1}{M}_{\odot }$, there is no recognizable difference among the simulations FDM, FIC, and CDM. However, for halo mass below $5\times {10}^{11}\,h{}^{-1}{M}_{\odot }$, the difference becomes noticeable.

Zoom In Zoom Out Reset image size

To quantitatively see the difference, one can refer to the right panel of Figure 5, which manifestly demonstrates the suppression caused by the QP and the modified initial condition. We can see that QP (FDM simulation) introduces 20% more suppression on the number density of $M\lt 2\,\times {10}^{11}\,{h}^{-1}{M}_{\odot }$ halos than that of FIC simulation, with modified initial condition only. There is no identifiable effect in the simulations FDM and FIC for the halo mass function with $M\gt 5\times {10}^{11}\,{h}^{-1}{M}_{\odot }$.

It is clear that the difference in modified initial conditions and QP start to have significant influence on the formation of halos. Furthermore, it is worth mentioning that by using the code established in Zhang et al. (2018), we are able to explore halo mass smaller than $2\times {10}^{13}\,h{}^{-1}{M}_{\odot }$ in cosmological simulations.

3.2. Impacts on Lyα Forest

Here, we investigate the difference in the 1D flux power spectrum among our four simulations with the method mentioned in Section 2.3, which enables us to quantitatively study the effect of the modified initial condition and dynamics introduced by FDM. Furthermore, our result of 1D flux power spectrum will be used to compare with the BOSS/XQ-100 data.

In Figure 6, we compare the 3D power spectra at different redshifts among different simulations: CDM, FIC, FDM, and F23. In the left panel of Figure 6, the effect of FDM initial condition is almost negligible at low redshifts, but plays an important role at high redshifts. However, the effect of QP is non-negligible, even at low redshifts. During the nonlinear gravitational evolution, even without the QP, the power spectrum of the simulation FIC grows differently from the simulation CDM. However, at z = 0 the power spectra of the simulations, FIC and CDM almost overlap, and one can barely see the difference between them. The effect of the FDM initial condition turns out to be a tiny suppression on power spectrum at low redshifts.

Zoom In Zoom Out Reset image size

On the other hand, the power spectrum of the simulation FDM is clearly different from the simulations FIC and CDM at low redshifts because of the non-trivial QP effect. For the power spectrum of the simulation F23, the suppression on small scales is even more significant than all other simulations due to its larger QP.

In the right panel of Figure 6, we show the ratio R(k) of the 3D power spectra from the simulations FDM, FIC, and F23 to that of simulation CDM. One can see that the QP suppresses the power spectrum by 2%–5% relative to the simulation FIC at $k\lt 10\,h\,{\mathrm{Mpc}}^{-1}$ for three different redshifts: z = 3.0, 3.6, 4.2. However, for $k\lt 1\,h\,{\mathrm{Mpc}}^{-1}$, the 3D power spectra from the simulations FDM, FIC. and F23 are identical. Therefore, we can conclude that, QP has no effect on large scale, which is well expected. For $k\gt 10\,h\,{\mathrm{Mpc}}^{-1}$, as it is approaching our resolution limit, the results should be taken carefully and critically. The effect of QP in the simulation F23 is clearly more significant than the simulation FDM.

The 1D flux power spectra are shown in the left panel of Figure 7 for comparison with the data of BOSS and XQ-100 at three different redshifts z = 3.0, 3.6, 4.2. The solid, dotted–dashed, dashed and dotted lines represent 1D flux power spectra from the simulations CDM, FIC, FDM, and F23, respectively. The dots and error bars are the data from BOSS in darker colors and XQ-100 in lighter colors. At the small-scale region, ${10}^{-2}\,{\mathrm{km}}^{-1}\,{\rm{s}}\lt k\lt {10}^{-1}\,{\mathrm{km}}^{-1}\,{\rm{s}}$, which corresponds to $1\,h\,{\mathrm{Mpc}}^{-1}\lt k\lt 10\,h\,{\mathrm{Mpc}}^{-1}$ in the discussion of the matter power spectrum, the 1D flux power spectra for the simulations FDM and F23 are relatively more suppressed than those of CDM and FIC. Although the simulation FIC differs from CDM in the initial condition, the difference in their spectra is still small.

Zoom In Zoom Out Reset image size

The right panel of Figure 7 shows the ratio R_F(k) of the 1D flux power spectra of the simulations FIC, FDM, and F23 to that of the simulation CDM. The degree of suppression is up to 10% at $k\simeq {10}^{-1}\,{\mathrm{km}}^{-1}\,{\rm{s}}$ for the simulation FDM. Additionally, the 1D flux power spectrum of the simulation F23 is suppressed even more than the simulation FDM. On the contrary, the degree of suppression is smaller than 4% at $k\simeq {10}^{-1}\,{\mathrm{km}}^{-1}\,{\rm{s}}$ for the simulation FIC. For $k\lt {10}^{-2}\,{\mathrm{km}}^{-1}\,{\rm{s}}$, the difference in the 1D flux power spectra for the simulations FDM, FIC, and CDM is less than 2%. In fact, the QP contribution cannot be neglected even if the QP is sub-dominant as shown in Figure 1.

Finally, it is clear that a full investigation of hydrodynamic simulations with the QP is necessary to robustly constrain the mass region of FDM.

It has been reported in Armengaud et al. (2017), Iršič et al. (2017), and Kobayashi et al. (2017) that the FDM mass can be excluded up to ${10}^{-21}\,\mathrm{eV}$ at 2σ significance. However, the systematic uncertainties arising from both simulation (N-body and hydrodynamics) and gas properties (such as the temperature of gas) were not discussed adequately in these works.

Note that the well-known discrepancies in simulation will introduce uncertainties on the observables, e.g., the code-to-code inconsistency in hydrodynamic simulations introduces at least ∼5% uncertainty in 1D flux power spectrum as shown in Regan et al. (2007) and Bird et al. (2013). A brief discussion of the uncertainties in simulation is given in Appendix B for interested readers.

In the previous section, we consider the effects of FDM. However, the thermal properties of IGM also impact the 1D flux power spectrum, as shown in Equations (6) and (12). As these two effects are degenerate in the linear power spectrum, one should consider them simultaneously before reaching any conclusion. In this subsection, we will discuss the possible uncertainty caused by the gas temperature and compare it with that cause by FDM mass to see their degeneracy in the power spectrum.

First, as an illustration of the potential effects of gas, we consider the prediction of the linear theory and its uncertainty. There are several parameters in Equation (6), but for simplicity we only consider the uncertainty of FDM particle mass, m_χ, and average gas temperature, T₀. The error of the 1D flux power spectrum can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{P}_{F}}{{P}_{F}}=\displaystyle \frac{\partial \mathrm{ln}{P}_{F}}{\partial {m}_{\chi }}{\rm{\Delta }}{m}_{\chi }+\displaystyle \frac{\partial \mathrm{ln}{P}_{F}}{\partial {T}_{0}}{\rm{\Delta }}{T}_{0}.\end{eqnarray} \tag{ 15 }$$

We apply this formula to the linear 1D power spectrum at redshift z = 3.0, and set $\gamma =1.64{(1+3/5.5)}^{-0.15}=1.72$, ${\lambda }_{J}=100\,\mathrm{kpc}$ in the calculation of partial derivatives. For the most optimistic estimation, we set ${\rm{\Delta }}{P}_{F}=0$, so that we can express ${\rm{\Delta }}{m}_{\chi }$ in terms of ${\rm{\Delta }}{T}_{0}$. In a most conservative precision, ${\rm{\Delta }}{m}_{\chi }/{m}_{\chi }\lt 100 \%$, one can obtain a maximum temperature uncertainty, ${\rm{\Delta }}{T}_{0,\max }$, and compare it with the current uncertainty of T₀. From Equation (10), the best-fit and 1σ uncertainty of T₀ at redshift z = 3 is about

$$\begin{eqnarray*}&&{T}_{0}={20}_{-3}^{+7}\times {10}^{3}\,{\rm{K}}.\end{eqnarray*}$$

For FDM mass ${m}_{\chi }=2.5\times {10}^{-23}\,\mathrm{eV}$, ${\rm{\Delta }}{T}_{0,\max }$ is larger than ${10}^{4}\,{\rm{K}}$. Therefore, we can actually exclude this kind of FDM model robustly with Lyα observation. However, for FDM mass in the range from ${10}^{-22}\,\mathrm{eV}$ to ${10}^{-21}\,\mathrm{eV}$, we need much higher precision of the temperature to draw the conclusion (see Table 2).

Table 2.  For Each FDM Mass m_χ, the Gas Temperature Uncertainty Should be Less than ΔT_0,max for any Confirmation or Exclusion to be Valid

mχ	${\rm{\Delta }}{T}_{0,\max }$
2.5 × 10−22 eV	1.5 × 103 K
5 × 10−22 eV	2 × 102 K
10−21 eV	1 × 102 K

Download table as:  ASCIITypeset image

Second, we use a similar method as in Section 3.2, with one difference in the RHS of Equation (13); we multiply the ratio of bias function W with different gas temperature rather than the ratio of matter power spectra,

$$\begin{eqnarray}{P}_{F,\mathrm{thermal}}({k}_{z},z) & = & {\int }_{{k}_{z}}^{\infty }\displaystyle \frac{{kdk}}{2\pi }{P}_{\mathrm{DM}}(k,z)W(k)\\ & & \times \,\displaystyle \frac{W({k}_{z},k,{T}_{0})}{W({k}_{z},k,{T}_{0,\mathrm{bestfit}})}.\end{eqnarray} \tag{ 16 }$$

Here, we used the bias function in the linear theory Equation (8). Even though at the redshift range z = 3 to 4, the structure in the universe is nonlinear, we still use the ratio predicted by linear theory because our purpose here is to demonstrate the uncertainty of gas temperature instead of studying the thermal properties of gas which requires full hydrodynamics simulation.

As in Section 3.2, we will calculate the ratio ${R}_{F}={P}_{F,\mathrm{thermal}}/{P}_{F}$ again and compare with with the suppression induced by FDM. The best-fit value ${T}_{0,\mathrm{bestfit}}(z)$ is calculated from Equation (9), with ${T}_{0,\mathrm{bestfit}}^{A}=9.2\times {10}^{3}\,{\rm{K}}$ and ${T}_{0,\mathrm{bestfit}}^{S}\,=-2.5$. The upper limit in the red shaded region in Figure 8 is given by ${T}_{0}^{A}=9.1\times {10}^{3}\,{\rm{K}}$ and ${T}_{0}^{S}=-2.05$, while the lower limit corresponds to ${T}_{0}^{A}=10.4\times {10}^{3}\,{\rm{K}}$ and ${T}_{0}^{S}=-3$, which are chosen according to the 1σ uncertainty Equation (10). From Figure 8, it is obvious that the current constraint on the temperature of gas is not good enough to exclude FDM mass ${m}_{\chi }=2.5\times {10}^{-22}\,\mathrm{eV}$, let alone the fact that we have not even included the uncertainties of other thermal parameters of the gas. A detailed study of the thermal properties of the gas is needed to settle the issue.

Zoom In Zoom Out Reset image size

Finally, as pointed out in Hui et al. (2017), several additional astrophysical processes may also alter the 1D flux power spectrum, such as the effects of fluctuations in the ionizing background (Hui et al. 2017), patchy reionization (Cen et al. 2009), modification of the thermal history (Garzilli et al. 2017), and galactic outflows (Viel et al. 2013). These additional factors are yet absent in many of the present cosmological simulations. To sum up, the uncertainties from both numerical and physical factors prevent us from constraining the mass of FDM precisely, unless models and analyses with percent level precisions are available. A hydrodynamic simulation handling all of the uncertainties is needed in the future to set the correct constraint on the FDM mass.

In this paper, we have extended our FDM smoothed particle hydrodynamics methodology to cosmological N-body simulation. Unlike previous works in the literature, we have implemented not only the FDM initial condition, but also the QP effect to our cosmological N-body simulations. The correct transformation of QP from physical to comoving coordinates has been derived in this work. With this new technique, we have performed four different simulations: CDM, FIC, FDM, and F23. We have shown the difference of over-density between the FIC and FDM simulations. We have found that some granular structures located at higher-density regions can be produced by QP. Remarkably, we are able to probe halos with mass smaller than $2\times {10}^{13}\,h{}^{-1}{M}_{\odot }$ in cosmological simulations based on our methodology.

The matter power spectra from our four simulations tell us that the impact from QP is non-trivial, as shown by the difference between FIC and FDM. Compared with CDM in the region with wavenumber $k\lt 10\,{h}^{-1}$Mpc, the power spectrum suppression due to the effect of initial condition is less than 1% at redshift z = 4.2, but the QP effect can cause <5% suppression in the same region. Hence, the impact from QP on the power spectrum is more significant than that from initial conditions at low redshifts. Moreover, the QP effect also depends on the redshift. At high redshifts, $z\sim 10$ the effect of modified initial condition is more important than QP, but vice versa at low redshifts $z\sim 0$. Considering the DM mass around ${10}^{-23}\,\mathrm{eV}$, the matter power spectrum of the simulation F23 shows a large deviation from that of CDM in the wavenumber $k\gtrsim 2\,{h}^{-1}\,\mathrm{Mpc}$ region.

Using the results of these four different simulations, we then further studied the flux power spectrum of Lyα forest. We obtain the 1D flux power spectrum by integrating the 3D matter power spectrum taken from our simulation result. There is still suppression on the flux power spectra of FDM and F23 on small scales compared to that of CDM. However, the difference between FIC and CDM is small, which indicates that the suppression due to the transfer function of FDM is consumed by the nonlinear evolution. To summarize, by comparing the flux power spectra of different simulations, we demonstrate that the QP causes non-trivial effect on small scales, which could be important in the study of Lyα forest.

Finally, we discussed the uncertainties in simulations and models. As an example, we did a rough estimation and found that the average gas temperature T₀ and FDM mass m_χ are degenerate because a smaller m_χ or a higher T₀ both suppress the 1D flux power spectrum. We found that the current constraint on the average temperature of the gas is not accurate enough to exclude FDM mass ${m}_{\chi }\gt {10}^{-22}\,\mathrm{eV}$, which is summarized in Table 2 and Figure 8.

We conclude that the QP plays an important role in structure formation and affects the prediction for Lyα forest significantly. A further comprehensive hydrodynamic simulation including the QP, and a precise constraint on the gas temperature of Lyα forest are needed to solidly set a lower bound on the FDM particle mass.

We acknowledge Eric Armengaud and Vid Iršič for their discussion and Lachlan Lancaster for his useful suggestions. We also acknowledge Lam Hui and Tom Broadhust for their useful comments. This work is partially supported by a CUHK Discretionary Fund and by the MoST of Taiwan under the grant no. 105-2112-M-007-028-MY3.

We show the full transformation of the QP from physical to comoving coordinates. The Lagrangian can be expressed as

$$\begin{eqnarray}&&L=\displaystyle \frac{1}{2}M{{\boldsymbol{v}}}^{2}-\displaystyle \frac{{{\hslash }}^{2}}{2{{m}_{\chi }}^{2}}\displaystyle \frac{M}{\bar{\rho }}{({\rm{\nabla }}\sqrt{\rho })}^{2}-M\phi ({\boldsymbol{r}}),\end{eqnarray} \tag{ 17 }$$

where M, v, m_χ, ρ, $\bar{\rho }$, and ϕ are mass and velocity of simulation particles, mass of FDM, density in physical coordinate, average density in physical coordinate, and gravitational potential, respectively. The first, second, and the third terms in Equation (17) are the kinetic energy, the potential energy arising from the QP, and the gravitational potential energy, respectively. We follow the basic transformations

$$\begin{eqnarray}{\boldsymbol{r}} & = & a{\boldsymbol{x}},\quad {\boldsymbol{v}}=\dot{{\boldsymbol{r}}}=a\dot{{\boldsymbol{x}}}+\dot{a}{\boldsymbol{x}},\quad \rho ={\rho }_{x}/{a}^{3},\\ {{\rm{\nabla }}}_{x} & = & a{{\rm{\nabla }}}_{r},\quad {{\rm{\Delta }}}_{x}={a}^{2}{{\rm{\nabla }}}_{r}\cdot {{\rm{\nabla }}}_{r}={a}^{2}{{\rm{\Delta }}}_{r},\end{eqnarray} \tag{ 18 }$$

which incorporate the scale factor a to account for the expansion of the universe. Here, r denotes the physical coordinates and x denotes the comoving coordinates.

After including the transformations, the Lagrangian Equation (17) becomes

$$\begin{eqnarray}&&L=\displaystyle \frac{1}{2}M{(a\dot{{\boldsymbol{x}}}+\dot{a}{\boldsymbol{x}})}^{2}-\displaystyle \frac{1}{{a}^{2}}{K}_{\rho x}-M\phi ({\boldsymbol{x}}),\end{eqnarray} \tag{ 19 }$$

where we define

$$\begin{eqnarray}&&{K}_{\rho x}=\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{\chi }^{2}}\displaystyle \frac{M}{\bar{{\rho }_{x}}}{({{\rm{\nabla }}}_{x}\sqrt{{\rho }_{x}})}^{2}.\end{eqnarray} \tag{ 20 }$$

The transformation of the Poisson equation can be written as

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{r}\phi =4\pi G(\rho -{\rho }_{{\rm{\Lambda }}})\to {{\rm{\Delta }}}_{x}\phi =\displaystyle \frac{4\pi G}{a}({\rho }_{x}-{\rho }_{{\rm{\Lambda }}x}),\end{eqnarray} \tag{ 21 }$$

where ${\rho }_{{\rm{\Lambda }}}$ is the density of the cosmological constant. To simplify the equation, we perform a canonical transformation on the Lagrangian

$$\begin{eqnarray}L & = & L-\displaystyle \frac{{dF}({\boldsymbol{x}},t)}{{dt}},\\ F & = & \displaystyle \frac{1}{2}{Ma}\dot{a}{{\boldsymbol{x}}}^{2},\end{eqnarray} \tag{ 22 }$$

which does not change the equation of motion. Now, the Lagrangian can be written as

$$\begin{eqnarray}&&L=\displaystyle \frac{1}{2}{{Ma}}^{2}{\dot{{\boldsymbol{x}}}}^{2}-\displaystyle \frac{1}{{a}^{2}}{K}_{\rho x}-M{\rm{\Phi }},\end{eqnarray} \tag{ 23 }$$

where we define ${\rm{\Phi }}=\phi +(1/2)a\ddot{a}{{\boldsymbol{x}}}^{2}$ for simplicity. Consequently, the Poisson equation is converted to

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{x}{\rm{\Phi }}={{\rm{\Delta }}}_{x}\left(\phi +\displaystyle \frac{1}{2}a\ddot{a}{{\boldsymbol{x}}}^{2}\right)=\displaystyle \frac{4\pi G}{a}({\rho }_{x}-{\rho }_{{\rm{\Lambda }}x})+3a\ddot{a}.\end{eqnarray} \tag{ 24 }$$

From the second equation of the Friedmann equations,

$$\begin{eqnarray}&&\displaystyle \frac{\ddot{a}}{a}=-\displaystyle \frac{4\pi G}{3}(\bar{\rho }-{\rho }_{{\rm{\Lambda }}})=-\displaystyle \frac{4\pi G}{3{a}^{3}}(\bar{{\rho }_{x}}-{\rho }_{{\rm{\Lambda }}x}),\end{eqnarray} \tag{ 25 }$$

one can obtain

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{x}{\rm{\Phi }}=\displaystyle \frac{4\pi G}{a}({\rho }_{x}-\bar{{\rho }_{x}}).\end{eqnarray} \tag{ 26 }$$

Additionally, we define Ψ = aΦ and acquire ${{\rm{\Delta }}}_{x}{\rm{\Psi }}\,=4\pi G({\rho }_{x}-\bar{{\rho }_{x}})$. Now, the Lagrangian can be written as

$$\begin{eqnarray}&&L=\displaystyle \frac{1}{2}{{Ma}}^{2}{\dot{{\boldsymbol{x}}}}^{2}-\displaystyle \frac{1}{{a}^{2}}{K}_{\rho x}-\displaystyle \frac{1}{a}M{\rm{\Psi }}.\end{eqnarray} \tag{ 27 }$$

Hence, we can define the canonical momentum as

$$\begin{eqnarray}&&{\boldsymbol{p}}=\displaystyle \frac{\partial L}{\partial \dot{{\boldsymbol{x}}}}={{Ma}}^{2}\dot{{\boldsymbol{x}}},\end{eqnarray} \tag{ 28 }$$

and write down the Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2{{Ma}}^{2}}{{\boldsymbol{p}}}^{2}+\displaystyle \frac{1}{{a}^{2}}{K}_{\rho x}+\displaystyle \frac{1}{a}M{\rm{\Psi }}.\end{eqnarray} \tag{ 29 }$$

As a consequence, we obtain the equations of motion

$$\begin{eqnarray}\dot{{\boldsymbol{x}}} & = & \displaystyle \frac{\partial H}{\partial {\boldsymbol{p}}}=\displaystyle \frac{{\boldsymbol{p}}}{{{Ma}}^{2}},\\ \dot{{\boldsymbol{p}}} & = & -\displaystyle \frac{\partial H}{\partial {\boldsymbol{x}}}=-\displaystyle \frac{1}{{a}^{2}}{{\rm{\nabla }}}_{x}{K}_{\rho x}-\displaystyle \frac{M}{a}{{\rm{\nabla }}}_{x}{\rm{\Psi }}.\end{eqnarray} \tag{ 30 }$$

In the equations of motion, the terms for gravity and QP have prefactors $1/a$ and $1/{a}^{2}$, respectively, hence we need to treat them separately in the comoving coordinates. In the physical coordinates, we can deal with these two terms simultaneously since a = 1.

The systematic uncertainties in N-body simulation are well studied. The matter power spectrum is affected by the methodology of generating initial condition for simulations, such as the choice between the first-order Lagrangian perturbation theory and the second-order Lagrangian perturbation theory can introduces ∼6% difference in the matter power spectrum, as reported in L'Huillier et al. (2014). The halo mass function in the large mass end is sensitive at ∼7% level to the choice of the finite simulation box size, which can be corrected owing to their clear nature (Bagla & Prasad 2006).

The uncertainty resulting from different N-body simulation codes is numerically hard to estimate, and the origin of the discrepancies among different codes is also barely known. Such kinds of issues are comprehensively discussed in Heitmann et al. (2010). Roughly speaking, the discrepancy between different codes is within 10% in the matter power spectrum and halo mass function. However,  Kim et al. (2013) and Sembolini et al. (2016) reported a 20% discrepancy in the center and 10% in the outskirts of the halo by comparing the halo density profile from several different codes.

Unlike the ${ \mathcal O }(10 \% )$ errors in the N-body simulations, the uncertainties involved in hydrodynamic simulations are much larger. Usually, a considerable uncertainty can be introduced due to the treatments of the gaseous component. As demonstrated in Oshea et al. (2005), Vazza et al. (2011), Kim et al. (2013), and Sembolini et al. (2016), a large inconsistency of the galaxy structure can be caused by using different hydrodynamic codes. Quantitatively, the errors of the gas density and temperature in the centers of the galaxies are about one to two orders of magnitude, estimated by comparing different hydrodynamic codes.

The 1D flux power spectrum of Lyα forest can also contain some level of uncertainties, e.g., see the code comparison in Regan et al. (2007) and Bird et al. (2013). For hydrodynamic simulations using ENZO (Bryan et al. 2014) and using Gadget2, there is also a ∼5% difference in the 1D flux power spectrum and ∼10% difference in the probability distribution of density and temperature. Moreover, for simulations using AREPO (Springel 2010) and Gadget2, there are also 5% difference in their 1D flux power spectra.