Clues to the nature of dark matter from first galaxies

TL;DR

This study uses high-resolution galaxy formation simulations to explore how a 3 keV warm dark matter model affects early universe galaxy properties, reionization, and observable signatures compared to cold dark matter.

Contribution

It provides the first detailed simulation-based analysis of warm dark matter effects on high-redshift galaxy formation and reionization, highlighting key differences from cold dark matter.

Findings

WDM suppresses low-mass halo formation.

WDM results in fewer stars at fixed halo mass.

WDM leads to higher dark halo fractions below 10^9 M_sun.

Abstract

We use thirty-eight high-resolution simulations of galaxy formation between redshift 10 and 5 to study the impact of a 3 keV warm dark matter (WDM) candidate on the high-redshift Universe. We focus our attention on the stellar mass function and the global star formation rate and consider the consequences for reionization, namely the neutral hydrogen fraction evolution and the electron scattering optical depth. We find that three different effects contribute to differentiate warm and cold dark matter (CDM) predictions: WDM suppresses the number of haloes with mass less than few $10^9$ M$_{\odot}$; at a fixed halo mass, WDM produces fewer stars than CDM; and finally at halo masses below $10^9$ M$_{\odot}$, WDM has a larger fraction of dark haloes than CDM post-reionization. These three effects combine to produce a lower stellar mass function in WDM for galaxies with stellar masses at and below $\sim 10^7$ M$_{\odot}$. For $z > 7$, the global star formation density is lower by a factor of two in the WDM scenario, and for a fixed escape fraction, the fraction of neutral hydrogen is higher by 0.3 at $z \sim 6$. This latter quantity can be partially reconciled with CDM and observations only by increasing the escape fraction from 23 per cent to 34 per cent. Overall, our study shows that galaxy formation simulations at high redshift are a key tool to differentiate between dark matter candidates given a model for baryonic physics.

Tables

Table 1. Table 1: Size and resolution of the 6 N 𝑁 N -body simulation boxes.

Box	Size	N	mdm	${ϵ}_{\mathrm{dm}}$
	Mpc		M⊙	kpc
1	$5$	${250}^3$	$1.00\times {10}^6$	$0.50$
2	$15$	${250}^3$	$2.70\times {10}^7$	$1.50$
3	$20$	${500}^3$	$8.00\times {10}^6$	$1.00$
4	$40$	${600}^3$	$3.70\times {10}^7$	$1.67$
5	$100$	${500}^3$	$1.00\times {10}^9$	$5.00$
6	$200$	${500}^3$	$8.00\times {10}^9$	$10.0$

Table 2. Table 2: Resolution levels used in our zoom-in simulations.

Level	$m_{\mathrm{DM}}$	$m_{\mathrm{gas}}$	${ϵ}_{\mathrm{DM}}$	${ϵ}_{\mathrm{gas}}$
	M⊙	M⊙	pc	pc
1	$1.95\times {10}^3$	$3.56\times {10}^2$	$63$	$27$
2	$1.56\times {10}^4$	$2.85\times {10}^3$	$125$	$53$
3	$1.25\times {10}^5$	$2.28\times {10}^4$	$250$	$107$
4	$4.21\times {10}^5$	$7.70\times {10}^4$	$375$	$160$
5	$5.78\times {10}^5$	$1.06\times {10}^5$	$417$	$178$

Table 3. Table 3: Best-fitting parameter values for global SFR density, given by equation ( 6 ).

Cosmology	$\kappa$	$\lambda$	$\mu$	$z_0$
CDM	0.00131	0.254	6.03	0.0937
WDM	0.676	0.377	6.16	0.381

Table 4. Table 4: Best-fitting parameter values for halo mass functions.

$z$	Cosmology	$A$	$B$	$C$	$M_0$	$\alpha$
5	CDM	52.93	3.641	15.31	11.05	0.1315
	WDM	48.99	3.863	16.06	8.574	0.2371
6	CDM	53.23	3.851	17.12	9.568	0.1462
	WDM	53.60	4.331	14.72	9.092	0.2569
7	CDM	47.86	3.819	14.86	8.572	0.1862
	WDM	62.80	5.0117	17.48	9.201	0.2411
8	CDM	54.53	4.338	16.04	8.975	0.1837
	WDM	60.06	4.987	14.13	9.495	0.2628
9	CDM	51.62	4.274	15.18	8.557	0.1963
	WDM	58.73	4.960	15.83	8.878	0.2616
10	CDM	54.11	4.652	11.09	9.837	0.2199
	WDM	65.11	5.537	21.02	8.175	0.2643

Table 5. Table 5: The median parameter values, including the 1 σ 𝜎 \sigma uncertainties, derived from the MCMC sample for the M ⋆ subscript 𝑀 ⋆ M_{\star} - M 200 subscript 𝑀 200 M_{200} relation in log-log space.

$z$	Cosmology	slope	intercept
5	CDM	${1.536}_{-0.051}^{+0.049}$	$-{8.24}_{-0.47}^{+0.50}$
	WDM	${1.683}_{-0.062}^{+0.058}$	$-{9.73}_{-0.58}^{+0.61}$
6	CDM	${1.490}_{-0.052}^{+0.050}$	$-{7.67}_{-0.47}^{+0.50}$
	WDM	${1.621}_{-0.070}^{+0.065}$	$-{9.01}_{-0.63}^{+0.67}$
7	CDM	${1.448}_{-0.060}^{+0.056}$	$-{7.25}_{-0.51}^{+0.56}$
	WDM	${1.483}_{-0.082}^{+0.078}$	$-{7.61}_{-0.74}^{+0.78}$
8	CDM	${1.358}_{-0.069}^{+0.064}$	$-{6.37}_{-0.58}^{+0.62}$
	WDM	${1.444}_{-0.100}^{+0.094}$	$-{7.23}_{-0.87}^{+0.93}$
9	CDM	${1.354}_{-0.083}^{+0.078}$	$-{6.25}_{-0.96}^{+0.73}$
	WDM	${1.350}_{-0.13}^{+0.13}$	$-{6.28}_{-1.16}^{+1.21}$
10	CDM	${1.43}_{-0.12}^{+0.11}$	$-{6.96}_{-0.99}^{+1.06}$
	WDM	${1.41}_{-0.20}^{+0.19}$	$-{6.89}_{-1.71}^{+1.81}$

Table 6. Table 6: The median parameter values, including the 1 σ 𝜎 \sigma uncertainties, resulting from the MCMC procedure for f ⋆ subscript 𝑓 ⋆ f_{\star} , which follows the functional form of equation ( 5 ).

$z$	Cosmology	$\beta$	$M_1$
5	CDM	${1.45}_{-0.32}^{+0.63}$	$-{8.73}_{-0.17}^{+0.16}$
	WDM	${1.64}_{-0.47}^{+1.04}$	$-{8.96}_{-0.17}^{+0.19}$
6	CDM	${1.82}_{-0.42}^{+0.88}$	$-{8.46}_{-0.12}^{+0.15}$
	WDM	${1.82}_{-0.50}^{+1.05}$	$-{8.62}_{-0.16}^{+0.19}$
7	CDM	${2.02}_{-0.48}^{+0.96}$	$-{8.35}_{-0.12}^{+0.14}$
	WDM	${2.08}_{-0.62}^{+1.22}$	$-{8.59}_{-0.16}^{+0.22}$
8	CDM	${2.43}_{-0.63}^{+1.15}$	$-{8.15}_{-0.13}^{+0.13}$
	WDM	${2.50}_{-0.76}^{+1.31}$	$-{8.20}_{-0.22}^{+0.27}$
9	CDM	${2.36}_{-0.58}^{+1.11}$	$-{7.99}_{-0.14}^{+0.10}$
	WDM	${2.74}_{-0.88}^{+1.30}$	$-{7.99}_{-0.30}^{+0.30}$
10	CDM	${2.55}_{-0.69}^{+1.19}$	$-{8.03}_{-0.15}^{+0.12}$
	WDM	${2.78}_{-0.94}^{+1.29}$	$-{7.97}_{-0.38}^{+0.37}$

Table 7. Table 7: The best-fitting parameters derived from the MCMC procedure for the SFR- M ⋆ subscript 𝑀 ⋆ M_{\star} relation, which is linear in log-log space.

$z$	slope	intercept
5	0.999	$-8.947$
6	1.019	$-8.493$
7	1.070	$-8.802$
8	1.030	$-8.408$
9	0.9930	$-8.053$
10	0.9994	$-8.026$

Equations

$$P_{\mathrm{WDM}}(k)=[T_{\mathrm{WDM}}(k){]}^2{P}_{\mathrm{CDM}}(k),$$

$$T_{\mathrm{WDM}}(k)={\left[1+(\alpha k{)}^{2\nu }\right]}^{-5/\nu }.$$

$$\alpha =0.049{\left(\frac{m_{\mathrm{x}}}{1\mathrm{keV}}\right)}^{-1.11}{\left(\frac{{\Omega }_{\mathrm{x}}}{0.25}\right)}^{0.11}{\left(\frac{h}{0.7}\right)}^{1.22}{h}^{-1}\mathrm{Mpc},$$

$$\log \left(\frac{dN}{d\log M}\right)=A-B\log M-C{\mathrm{e}}^{{\left(M_0-\log M\right)}^{\alpha }},$$

$$f_{\star }=\frac{1+\tanh \left[\beta (\log M-M_1)\right]}{2},$$

$$\log \left(\frac{{\rho }_{\mathrm{SFR}}(z)}{{\mathrm{M}}_{\odot }\mathrm{Mp}{\mathrm{c}}^{-3}\mathrm{y}{\mathrm{r}}^{-1}}\right)=\kappa -\lambda z-\mu e^{z_0-z},$$

$${\dot{Q}}_{\text{Hii}}=\frac{{\dot{n}}_{\mathrm{ion}}}{⟨n_{\mathrm{H}}⟩}-\frac{Q_{\text{Hii}}}{t_{\mathrm{rec}}}.$$

$$t_{\mathrm{rec}}$$

$$\approx 0.97\mathrm{Gyr}{\left(\frac{C_{\text{Hii}}}{3}\right)}^{-1}{\left(\frac{T}{2\times 10^4\mathrm{K}}\right)}^{0.7}{\left(\frac{1+z}{7}\right)}^{-3},$$

$${\dot{n}}_{\mathrm{ion}}=f_{\mathrm{esc}}{\xi }_{\mathrm{ion}}{\rho }_{\mathrm{SFR}},$$

$$\log \left(\frac{{\xi }_{\mathrm{ion}}}{\mathrm{photons}{\mathrm{s}}^{-1}{\mathrm{M}}_{\odot }\mathrm{yr}}\right)=53.14,$$

$$\tau (z)=c⟨n_{\mathrm{H}}⟩{\sigma }_{\mathrm{T}}{\int }_0^z{f}_{\mathrm{e}}{Q}_{\text{Hii}}(z^{\prime })H^{-1}(z^{\prime })(1+z^{\prime }{)}^{2}d{z}^{\prime },$$

Full text

1 Introduction

In the next several years, we expect a large amount of data on the high-redshift ($z>6$) Universe to become available to the scientific community. Facilities like Atacama Large Millimetre/sub-millimetre Array (ALMA) and James Webb Space Telescope (JWST) will open a completly new observational window on the first billion years of the life of our Cosmos. These data will help us understand the early phases of galaxy formation, but they might also guide us towards a better understanding of the dark side of our Universe, as for example the nature of dark matter.

The current leading model for dark matter is based on a Cold candidate (CDM), initially motivated by the possible existence of Weakly Interactive Massive Particles (WIMPs) predicted by some extensions of the Standard Particle Models (Bertone et al., 2005). Such candidates should be in the mass range of several GeV (Bergström, 2012) and have been extensively searched for in underground labs, but with very little success to date (see Roszkowski et al., 2018, for a recent review). The lack of direct detection of a possible CDM particles raises the question if the dark matter particle might have a much lower mass than initially thought, and possibly lie in the keV range, being what is usually dubbed Warm Dark Matter (WDM).

WDM was initially put forward as a possible solution of the so called ‘small-scale crisis’ of CDM (e.g Moore, 1994; Klypin et al., 1999; Moore et al., 1999; Colín et al., 2000). More recently, there has been a mounting evidence that the solution of these problems most likely lies in the baryonic sector. Recent simulations and semi-analytical models have been able to reconcile a CDM universe with observations both on the number of galactic satellites (Bullock et al., 2000; Macciò et al., 2010; Sawala et al., 2016; Buck et al., 2018) and on the dark matter distribution at small scales by simply providing a more accurate treatment of galaxy formation and baryonic physics within a CDM model (Governato et al., 2010; Zolotov et al., 2012; Di Cintio et al., 2014; Oñorbe et al., 2015; Frings et al., 2017).

Moreover, current limits from the Lyman-$\alpha$ (Ly$~{}\alpha$) forest on the matter power spectrum (Iršič et al., 2017; Yèche et al., 2017) have set very stringent limits on the mass of a possible (thermal) WDM candidate, which is now constrained to $m_{\nu}>3.0$ keV. With such a mass, a WDM candidate will be practically indistinguishable from a cold one regarding the dark matter distribution on small scales (Colín et al., 2008; Macciò et al., 2012; Shao et al., 2013; Herpich et al., 2014) or the number of satellites (Macciò & Fontanot, 2010; Lovell et al., 2016). The similarity of currently allowed WDM models to CDM on small scales calls for a new venue on where to look for more clues on the actual mass of the dark matter candidate, and this venue can lie in the high-redshift Universe.

In WDM, due to the lack of power on small scales, structure formation is delayed with respect to a cold scenario. At present time, the effect of such an initial delay has been erased by the nonlinearity of galaxy formation, but it can still present itself in the very early stages of the Universe, when the first galaxies were formed. This signature can manifest itself in the formation of first stars (e.g. Gao & Theuns, 2007), or first galaxies (e.g. Dayal et al., 2015). Recently, Corasaniti et al. (2017) have used measurements of the galaxy luminosity function at $z=$ 6, 7, and 8 to derive constraints on WDM. They have combined very high resolution $N$-body (pure gravity) simulations with an empirical approach based on halo abundance matching to derive predictions for the luminosity function in different dark matter models. They have obtained a lower limit on the WDM thermal relic particle mass of 1.5 keV at a 2$\sigma$ level. This study relies on dark matter only simulations that might miss some of the nonlinear effects between the nature of dark matter and galaxy formation, which is a possible limitation.

The first galaxies are also thought to be the main source of the ionizing photons responsible for cosmic reionization, which may then provide another constraint to the nature of dark matter. Many recent studies have studied the connection between WDM and reionization (e.g. Lapi & Danese, 2015; Tan et al., 2016; Lopez-Honorez et al., 2017; Das et al., 2018). In particular, Carucci & Corasaniti (2019) extend the work of Corasaniti et al. (2017) to derive the requisite escape fractions of ionizing photons to match recent observations, but given the uncertain nature of this quantity, the results largely independent of underlying dark matter scenario. Two studies employing semi-analytic methods to study WDM in the context of high redshift and reionization also find that with a dark matter scenario consistent with observational constraints the ionizing population is shifted to lower halo masses, but that distinguishing CDM and WDM with current observations would be difficult (Bose et al., 2016; Dayal et al., 2017). Furthermore, Villanueva-Domingo et al. (2018) have used hydrodynamical simulations in a box of 20 $h^{-1}$ Mpc to derive constraints on WDM from galaxy luminosity functions, the ionization history, and the Gunn-Peterson effect. They have concluded that while some effects are indeed present, until the modelling of baryonic effects (star formation and feedback mainly) are constrained better, any conclusions on the nature of dark matter derived from reionization observables remain model-dependent. One possible way to anchor the parameters describing the baryonic physics (besides having better constraints from observations) is to use observations at one redshift (usually $z=0$) as a constraint and then check predictions at a much higher redshift.

In this study, we attempt to connect the high- and low-redshift Universe via hydrodynamical simulations of galaxy formation. We want to extend to the high-redshift Universe ($z>6$) the baryonic physics used in the NIHAO project (Wang et al., 2015), which has proven very successful in reproducing key observations in the low-redshift Universe ($0<z<2$), as for example the stellar mass - halo mass relation (Wang et al., 2015), the local velocity function (Macciò et al., 2016), and galaxies scaling relations (Dutton et al., 2017).

We use 19 high-resolution, zoom-in simulations in each cosmology (for a total of 38) to accurately sample the stellar mass function in the redshift range $5<z<10$, and we study the impact of a 3 keV dark matter candidate on several current and future observables, like the stellar mass function, the global star formation rate, the hydrogen neutral fraction, and the Cosmic Microwave Background (CMB) optical depth.

The paper is organized as follows. In Section 2, we present the set up of our galaxy formation simulations, in Section 3, we present our results and forecasts for current and future observations. Finally in Section 4, we discuss the implications or our results and present our conclusions.

2 Simulations

All simulations herein were ran in a flat $\rm\Lambda CDM$/$\rm\Lambda WDM$ cosmology, with parameters from the Planck Collaboration et al. (Planck Collaboration et al., 2014): Hubble Parameter $H_{0}=67.1$ km s*-1* Mpc*-1*, matter density $\Omega_{\rm m}=0.3175$, radiation density $\Omega_{\rm r}=0.00008$, dark energy density $\Omega_{\Lambda}=1-\Omega_{\rm m}-\Omega_{\rm r}=0.6824$, power spectrum normalization $\sigma_{\rm s}=0.8344$, and power spectrum slope $n=0.9624$.

2.1 Initial Conditions

Initial conditions for both the $N$-body and hydrodynamical, zoom-in simulations of galaxy formation were generated using a modified version of the grafic2 package (Bertschinger, 2001; Penzo et al., 2014). For WDM, the process is identical to CDM except for a modification to the power spectrum. The same random seed was used to generate all pairs of CDM/WDM 3 keV initial conditions.

WDM power spectra were computed with the use of a relative transfer function, as in Bode et al. (2001):

[TABLE]

where

[TABLE]

Here, $\nu=1.12$ (Viel et al., 2005) and $\alpha$ is the length scale of the break in the WDM power spectrum:

[TABLE]

where $m_{\rm x}$ is the mass of the WDM particle and $\Omega_{\rm x}$ the density.

2.2 $N$-body Simulations

First, we run a series of $N$-body simulations that only include dark matter, detailed in Table 1. We used six boxes that range from 5 to 200 Mpc length on a side with varying particle number. Softening was set to 1/40 of the intra-particle distance on the initial conditions grid. From this set of simulations, we are then able to sample a wide range of halo masses at $z=5-10$.

2.3 Hydrodynamical Simulations

From the dark matter-only simulations, 19 haloes were chosen to be re-simulated at higher resolution with baryons included. Table 2 shows the five resolution levels used in our simulations. These levels were chosen to have approximately one million particles inside the virial radius of the galaxy across the whole range of halo masses.

All simulations were performed with the SPH code gasoline (Wadsley et al., 2004, 2017). The code was setup in the framework of the NIHAO project (Wang et al., 2015), including metal cooling, chemical enrichment, star formation, and feedback from massive stars and supernovae (SN).

Stars are formed from gas cooler than $T$ = 15,000 K and denser than $n_{\rm th}=10{\rm~{}cm}^{-3}$. The star formation efficiency used was $c_{\star}$ = 0.1. Cooling via hydrogen, helium, and various metal-lines is included as in Shen et al. (2010), including photoionization and heating from ultraviolet (UV) background in Haardt & Madau (2012).

SN feedback is implemented using the blastwave approach as described in Stinson et al. (2013), which relies on delaying the cooling of nearby particles to a SN event. We also include what we dubbed Early Stellar Feedback: we inject 13 per cent of the UV luminosity of the stars as thermal energy before any SN events take place without disabling the cooling (see Stinson et al., 2013, for more details).

NIHAO galaxies have been shown to be consistent with a wide range of low-redshift galaxy properties, such as the stellar mass - halo mass relation and star formation rates (Wang et al., 2015); stellar disk kinematics (Obreja et al., 2016); cold and hot gas content (Stinson et al., 2015; Wang et al., 2016; Gutcke et al., 2016); and resolve the too-big-to-fail problem (Dutton et al., 2016). We are therefore confident that the simulations presented herein can be used as plausible tools to study the effect of a WDM cosmology on the high-redshift Universe.

Haloes in all simulations were identified using the MPI+OpenMP hybrid Amiga Halo Finder 111http://popia.ft.uam.es/AMIGA (AHF; Knollmann & Knebe, 2011). The virial masses of the haloes are defined as the mass within a sphere containing $\Delta=200$ times the cosmic critical matter density. The virial (total) mass is denoted as $M_{200}$, and $M_{\star}$ indicates the total stellar mass within the virial radius.

Beyond the 19 target galaxies, we include any galaxy within our zoom simulations that has greater than 10,000 particles. We also only consider central haloes, i.e. not satellites, and require no pollution from larger low-resolution dark matter particles, i.e. any particle larger than the lowest level in Table 2. The result is of order 100 galaxies at $z=5$, and tens of galaxies at $z=10$. Simulation analysis was done using the Python package pynbody (Pontzen et al., 2013). All parameter fits were performed using Markov Chain Monte Carlo (MCMC) via the Python package emcee (Foreman-Mackey et al., 2013).

3 Results

In the following, we will present the results of our $N$-body and galaxy simulations in CDM and WDM cosmologies at $z=5-10$. In all plots, CDM results are shown in blue with circles, while 3 keV WDM results are shown in red with squares. Section 3.1 shows the halo mass functions obtained from the $N$-body simulations. In Section 3.2, we find the stellar mass - halo mass relation obtained from our simulations. Section 3.3 shows our results for the fraction of galaxies that form stars as a function of halo mass. Section 3.4 shows the stellar mass functions obtained by convolving all prior results. Finally in Section 3.6, we compare our results for the reionization history and CMB optical depth to observational constraints.

3.1 Halo Mass Function

The $N$-body boxes were used to construct and fit halo mass functions. We used a five-parameter, modified Schechter function (Schechter, 1976):

[TABLE]

where $A$, $B$, $C$, $\alpha$, and $M_{0}$ are fitting parameters that are given in Appendix A and $M$ represents $M_{200}/{\rm M}_{\odot}$. The original Schechter function struggled to fit the WDM results, so we adopted this modified version for both cosmologies. The uncertainty for each halo mass bin is considered Poissonian.

The halo mass function represents the abundance of haloes at a given mass. As shown in Fig. 1, the halo mass functions exhibit a suppression of the number density of low-mass haloes in a WDM cosmology, a direct consequence of the lack of power on small scales in the WDM power spectrum. The lines represent the MCMC parameter fit, the values of which are given in Appendix A. The relations for the two cosmologies follow the same trend at $z=5-10$. Importantly, the difference between CDM and WDM is larger at higher redshift. At $z=0$, the halo mass functions would be identical for the two scenarios.

3.2 Stellar Mass - Halo Mass Relation

Fig. 2 shows the stellar mass - halo mass relation obtained from the set of zoom-in simulations run in CDM and WDM cosmologies.222To be included in Fig. 2 and the fit used in Section 3.6, we require that the haloes have at least 10 stellar particles. Based on this relation, CDM and WDM galaxies cannot be distinguished from each other. While our simulations closely follow the abundance matching relation(s) at lower redshift (see Wang et al., 2015), our predicted relation at these high redshifts consistently lies above the extrapolation from Moster et al. (2013).333The original NIHAO simulations of Wang et al. (2015) use the stellar mass within 20 percent of the virial radius, which would bring our results even closer to the relation. At high redshift, this cut is too conservative. This offset (minor in our simulations) suggests that high-redshift extrapolations of relations calibrated at $z<4$ should be used with caution. Many studies at high redshift exhibit similar offsets, see e.g. Rosdahl et al. (2018), but also see Ma et al. (2018). This figure also indicates that a larger fraction of WDM galaxies remain dark at high redshift compared to CDM, as seen by the smaller number of red points relative to blue. This effect is further explored in the next section.

The solid lines represent a linear fit to the logarithmic data, which is done independently for each cosmology and redshift. The median (1$\sigma$) parameters derived from the MCMC sample generate the solid lines (shaded regions), the values of which are located in Appendix A These (small) uncertainties are propagated to the stellar mass function results. As shown in Fig. 2, CDM and WDM behave similarly. We find weak evidence for redshift evolution of the normalization in that lower redshift has a higher value for the intercept of the linear fit, which is consistent with Ceverino et al. (2017). However, Ma et al. (2018) find no redshift evolution.

3.3 Dark Fraction

Star formation becomes more inefficient at low masses, to the point where no gas is able to collapse to the center and halo and form stars due to the UV background (Gnedin, 2000). Fig. 2 seems to suggest a larger fraction of ‘dark’ galaxies in WDM w.r.t. CDM, given that there are fewer red squares than blue circles. To quantify this effect, we want to compute the fraction of objects able to form stars as a function of their virial mass. A halo is considered dark if it contains no stellar particles.

Results are presented in Fig. 3 for $z=6$ and $z=10$. The fitting curve is a hyperbolic tangent two-parameter function

[TABLE]

where $\beta$ and $M_{1}$ are parameters and $M$ once again represents $M_{200}/{\rm M}_{\odot}$. The points represent the histogram of the $f_{\star}=1-N_{\rm dark}/N_{\rm bin}$, where $N_{\rm dark}$ are the number of dark haloes and $N_{\rm bin}$ are the total number of haloes in the bin. When fitting this curve, the uncertainties on $f_{\star}$ are derived from the number of haloes in each bin. The resultant median (1$\sigma$) uncertainties of the fitted parameters are generated from the MCMC sample and are represented by the solid line (shaded region), all values are found in Appendix A. We find weak evidence that WDM exhibits a steeper dark fraction fraction transition, which is consistent with Macciò et al. (2019).

In these simulations, reionization occurs $z\gtrsim 9$, and the shift of minimum halo mass for star formation is evident. At earlier times (higher redshift), the $M_{200}$ at which stars form is lower due to a lower UV background. At $z>8$, WDM has a lower minimum mass for star formation than CDM, and vice versa for lower redshift. Another interesting result is that only approaching $M_{200}>5\times 10^{10}{\rm M}_{\odot}$ guarantees that stars will form in galaxies at this redshift. The main caveat for these observations is that the uncertainties for $f_{\star}$ are large given our small halo numbers, and distinguishing with confidence between the two cosmologies is beyond our statistical reach.

3.4 Stellar Mass Function

The next step is to convolve the results from previous sections in order to obtain the stellar mass function. We started from the fitted stellar mass - halo mass relation, including the uncertainty in the parameters, for our set of galaxy simulations (Fig. 2) to convert the halo masses in Fig. 1 into stellar masses. We then multiplied by the fraction, including the uncertainty of the fit, of galaxies that actually formed stars (Fig. 3) to obtain the differential number density of galaxies as a function of their stellar mass, which is shown in Fig. 4. The solid line represents the best-fitting (and median) parameters and the shaded region is the propagated uncertainties, where the $M_{\star}$-$M_{200}$ relation dominates at high stellar masses (few high-mass haloes) and $f_{\star}$ dominates at low stellar masses (suppression of star formation becomes important). The $M_{\star}$-$M_{200}$ uncertainty is reflected the stellar massCeverino2017 uncertainty ($x$ axis), and $f_{\star}$ drives the uncertainty in the number density ($y$ axis).

As expected, the stellar mass functions substantially differ only at low masses, while the small differences at high masses are due to cosmic variance and low number statistics. WDM predicts fewer low-mass galaxies than CDM, with a larger difference at high redshift. For our choice of warm candidate mass (3 keV), the relative difference between CDM and WDM is only a factor of a few stellar masses around $10^{6}$ ${\rm M}_{\odot}$ (see also Villanueva-Domingo et al., 2018), below what is currently observable (Bouwens et al., 2015; Bouwens et al., 2017; Ceverino et al., 2017). On the other hand, our results are in agreement with Corasaniti et al. (2017), who found a lower limit of 2 keV from analysis of the luminosity function. Future observations and facilities might improve this limit (e.g. JWST), but then a very careful understanding of the baryonic physics implemented in the simulations will be needed, as discussed in Villanueva-Domingo et al. (2018).

3.5 Star Formation Rates

The star formation rate (SFR) of a galaxy is not only of interest on its own, the reionization history of the Universe is crucially dependent on this quantity. First, we compute the SFR-$M_{\star}$ relation for our galaxies as a function of redshift, where the SFR includes star formation occurring within the past 100 Myr at the redshift in question. As shown in Fig. 5, CDM and WDM follow the same relation at $z=7$, and the same is true for lower and higher redshifts. With this in mind, we fit this relation with a single power law for both models at each redshift (shown by the dashed line). The usual MCMC procedure was performed with each SFR value weighted by the number of stellar particles that formed in the past 100 Myr. The best-fitting slope and intercept values444As a check, the $1\sigma$ uncertainties were calculated as detailed in previous relation fittings, but the resultant values were very small and would be entirely subdominant to the uncertainties in $M_{\star}$-$M_{200}$ and $f_{\star}$. As such, the uncertainty in this relation is not presented nor propagated to the stellar mass functions. are found in Appendix A and are similar to those found in Ma et al. (2018). Contrary to Ma et al. (2018), we find evidence of redshift evolution in the intercept or normalization, though no evolution in the slope. Therefore, at lower redshift, a lower SFR is expected for the same stellar mass.

We then combine our stellar mass functions with the fits of SFR versus $M_{\star}$ to obtain an integrated global SFR density, $\rho_{\rm SFR}$, as a function of redshift (Fig. 6). Stellar mass functions were integrated between $10^{3}-10^{12}~{}{\rm M}_{\odot}$, and the results are insensitive to order-of-magnitude changes in these limits. The values of $\rho_{\rm SFR}$ obtained for $z=5-10$ were fitted following:

[TABLE]

where $\kappa$, $\lambda$, $\mu$, and $z_{0}$ are parameters, the best-fitting values derived from the usual MCMC procedure for cosmologies are found in Table 7. Note that the errorbars represent the propagated uncertainties from the previous relations and decrease at lower redshift as we have more statistical power. At $z=10$, WDM is a factor of two lower than the CDM case, and this factor decreases as redshift decreases. By $z=7$, the two cosmologies are very similar and essentially indistinguishable at even lower redshift. Both scenarios lie within the range of values between the uncorrected and dust-corrected results from Bouwens et al. (2015) as indicated by the shaded region in the figure. Note that we do not correct for dust or consider any effects beyond our intrinsic star formation.

3.6 Reionization History

In this section, we use the predicted galaxy populations in the CDM and WDM models to explore the effects of the WDM cosmology on cosmic reionization, using a commonly used, simple ODE model. We emphasize that this model is a rough estimate of the reionization history, but is a sufficient first step given the uncertainties in reionization modelling. In particular, we treat $f_{\rm esc}$, the fraction of ionizing photons that escape from galaxies into the IGM, as a free parameter that is tuned to match current observational measurements.

We follow the approach described in Kuhlen & Faucher-Giguère (2012), starting with this differential equation for the ionized fraction:

[TABLE]

Here, $\langle n_{\rm H}\rangle=X_{\rm p}\Omega_{\rm b}\rho_{\rm c}$ is the comoving hydrogen density, in terms of the hydrogen mass fraction $X_{\rm p}=0.75$, the baryon density $\Omega_{\rm b}=0.049$, and the critical density $\rho_{\rm c}$. The IGM recombination time is:

[TABLE]

where $\alpha_{\rm B}(T)$ is the case-B recombination coefficient at an IGM temperature $T=2\times 10^{4}$ K, with the value of $1.6\times 10^{-13}$ cm3 s*-1* (Storey & Hummer, 1995), and $Y_{\rm p}=0.25$ being the helium mass fraction. $C_{\text{H\,{ii}}}$ is the effective clumping factor of ionized gas in the IGM; we adopt a constant value of 3 in this work. Of course, a lower value would decrease the resultant $f_{\rm esc}$ and vice versa. For example, $C_{\text{H\,{ii}}}$ = 1 requires an $f_{\rm esc}$ that is approximately 10 per cent lower.

Finally, $\dot{n}_{\rm ion}$ is the global production rate of ionizing photons:

[TABLE]

where $f_{\rm esc}$ is the effective fraction of photons that escape from galaxies to reionize the IGM, and $\xi_{\rm ion}$ is the ionizing photon production efficiency for a typical stellar population per unit time per unit SFR:

[TABLE]

where the right-hand value is chosen to be consistent with recent works, e.g. Lovell et al. (2018). In fact, some observational work points to a higher value (e.g. Topping & Shull, 2015). The exact value is unimportant as $\xi_{\rm ion}$ is completely degenerate with $f_{\rm esc}$, and a higher $\xi_{\rm ion}$ merely indicates a lower $f_{\rm esc}$ and vice versa.

The differential equation was evolved backwards in time, starting with an ionized fraction of $0.1$ at $z=10$, a value consistent with results from simulations, e.g. Dixon et al. (2016), and observations of the CMB optical depth measurements (Planck Collaboration et al., 2016). This simplification was necessary due to our small sample of sufficiently massive haloes at $z>10$, limiting the accuracy of predicted stellar mass functions (see Fig. 1). The $f_{\rm esc}$ was treated as a free parameter and tuned separately for the CDM and WDM scenarios with the goal of matching observations and reaching $Q_{\text{H\,{ii}}}(z=6)=1$. We found that in CDM $f_{\rm esc}=0.23$ matched the observed data well, while WDM can be forced to match CDM with an increased $f_{\rm esc}=0.34.$ Fig. 7 summarizes our ionization history results. The initial ionization fraction assumption does not significantly impact our results. For example, doubling the initial ionized fraction decreases the resultant $f_{\rm esc}$ at the percent level.

The physical value of $f_{\rm esc}$, including the potential redshift and halo mass dependence, is a hotly debated topic. Measurements of nearby galaxies have yielded very small values of a few per cent (e.g. Rutkowski et al., 2017; Steidel et al., 2018). Faisst (2016) have observed a somewhat larger $f_{\rm esc}$ and evidence of increasing values with redshift. Simulations have not reached a consensus, but $f_{\rm esc}$ appears to vary significantly from halo to halo (e.g. Paardekooper et al., 2015; Trebitsch et al., 2017; Rosdahl et al., 2018). Although our values are generally larger, the high redshift Universe remains uncertain, and $f_{\rm esc}$ and $\xi_{\rm ion}$ are degenerate, meaning a larger $\xi_{\rm ion}$ could bring our results more in line with observations.

Finally, we used the resultant ionized fractions to compute the CMB optical depth:

[TABLE]

where $c$ is the speed of light; $\sigma_{\rm T}$ is the Thompson cross section; $\rm H(z)$ is the Hubble parameter. $f_{\rm e}=1\ +\ \eta\ Y_{\rm p}/4X_{\rm p}$ is the number of free electrons per hydrogen nucleus. We consider helium to be singly ionized $(\eta=1)$ at the same rate as hydrogen for $z>3$ and doubly ionized $(\eta=2)$ at $z\leq 3$, which is consistent with recent observations (Worseck et al., 2018). Fig. 8 shows our optical depth results.

The WDM 3 keV model struggles to reach the $1\sigma$ confidence level of the Planck measurement for optical depth with the CDM-inspired $f_{\rm esc}=0.23$. This result is a direct consequence of the delay in reionization shown in Fig. 7. The low optical depth even for our CDM results can be partially explained by the lack of simulation data at $z>10$. Furthermore, reionization modelling that includes inhomogeneity and more physical processes will generically produce slower reionization histories and therefore a larger $\tau$.

4 Conclusion

With upcoming facilities such as ALMA and JWST, we will have access to unprecedented data for the high-redshift Universe and the formation and properties of the first galaxies. The data are sensitive to the early phase of structure formation and hence might help shed light on the nature of dark matter. In this paper, we present a set a large set of high-resolution cosmological hydrodynamical simulations of galaxy formation performed in two different cosmological models, Cold Dark Matter (CDM) and Warm Dark Matter (WDM). For the WDM simulations, we have chosen a particle mass of 3 keV, in agreement with current constraints (Iršič et al., 2017).

Our simulations cover the redshift range $5<z<10$ and span three orders of magnitude in halo mass, from $10^{8}$ to $10^{11}$ ${\rm M}_{\odot}$ (at $z=6$) and more than four in stellar mass, for a total of 38 individual zoom-in simulations among the two cosmologies. In addition to the main zoom-in galaxy, we include all haloes with more than 10,000 particles, resulting in over 100 galaxies at $z=5$. Simulations in the two different cosmologies have been performed with the same implementation of the baryonic physics as in the NIHAO project (Wang et al., 2015), which has proven to be extremely successful in reproducing several observed properties of galaxies at low and intermediate redshift (Wang et al., 2015; Dutton et al., 2017) and for a very extended mass range (e.g. Macciò et al., 2017; Buck et al., 2018).

WDM simulations have a systematically lower stellar mass function at low masses. This deficit is due to a combination of three effects: a suppression of the number of haloes with mass less than $\approx 10^{9}~{}{\rm M}_{\odot}$, a slightly lower star formation efficiency at a fixed halo mass, and a larger fraction of dark haloes for masses below $\approx 10^{9}~{}{\rm M}_{\odot}$. At all redshifts, the two mass functions start to depart at a stellar mass of $10^{7}~{}{\rm M}_{\odot}$. The difference peaks at about an order of magnitude at $M_{\star}\approx 10^{5.5}~{}{\rm M}_{\odot}$ and is more pronounced at higher redshift. Unfortunately, for a realistic candidate of 3 keV, these differences are on the mass scales not yet accessible by current facilities, such as Hubble Space Telescope.

The delayed structure formation in WDM models also manifests itself in the global star formation density rate (star formation rate per unit volume) as a function of redshift, which is about a factor of a few lower in WDM w.r.t. CDM at high redshift. By $z=6$, this quantity is very similar in these two models.

Finally, the different star formation histories in the two cosmologies also leave an imprint in the reionization history of the Universe. For a fixed photon escape fraction ($f_{\rm esc}=0.23$), the fraction of neutral hydrogen is higher by  0.3 in WDM for a redshift range between $5<z<7$; a larger escape fraction ($f_{\rm esc}=0.34$) is needed in order to reconcile WDM with current data and have a completely ionized Universe by $z\approx 6$. The same effect is visible in the CMB optical depth, where WDM is inconsistent with Planck constraints for $f_{\rm esc}=0.23$, and a larger value than even 0.34 is needed to match the data.

Our results suggest that high redshift is one of the best places to look for the signature of a possible warm candidate on structure formation, in agreement with recent studies from Corasaniti et al. (2017). Even for a quite warm-ish mass of 3 keV, the delayed structure formation of WDM leaves a distinct imprint on observable quantities, such as the stellar mass function, and on the evolution of cosmic reionization. On the other hand, as recently shown by Villanueva-Domingo et al. (2018), a firm comparison with current data is at the moment prevented by our lack on knowledge on how to parametrize the ‘baryonic’ side of structure formation, as for example the UV photon escape fraction. Our results, combined with previous works, strengthen the need to combine galaxy formation and cosmology if we want to shed light on the nature of the dark components of our Universe.

Acknowledgments

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at the Leibniz Supercomputing Centre (www.lrz.de) and the High Performance Computing resources at New York University Abu Dhabi.

Appendix A Parameter fitting values

In this appendix, we present the parameters derived from the MCMC fitting procedure for the majority of relations presented in the text. Table 4 shows our best-fittingting parameters for the halo mass function at $z=5-10$ in CDM and 3 keV WDM cosmologies, according to equation (4). Table 5 includes the median and 1$\sigma$ parameter values for the $M_{\star}$-$M_{200}$ relation, which in linear in log-log space. Table 6 gives the same quantities for the fraction of dark galaxies, as described by equation (5).