Warm Dark matter constraints from the joint analysis of CMB, Ly α, and global 21 cm data Open Access

ABSTRACT

1 INTRODUCTION

The concordance Lambda cold dark matter (ΛCDM) model is extremely successful in explaining the Universe’s large-scale structure, e.g. extremely accurate prediction of the cosmic microwave observations (Peebles 1982) and large-scale distribution of the galaxies (Blumenthal et al. 1984). Interestingly, the same cosmological model fails to match with some of the galactic and subgalactic scale observations such as – (i) dearth of low-mass galaxies in low-mass haloes (Peebles 2001), (ii) ΛCDM model’s prediction of cuspy core in the dark matter (DM) haloes opposes the observationally preferred constant density cores (Moore et al. 1999; Subramanian, Cen & Ostriker 2000), (iii) Too big to fail problem for the field galaxies (Oman et al. 2016). The root cause behind all of these problems is the abundance of the small-scale structure due to the very cold nature (with mass ∼100 GeV) of the constituent DM particle in the ΛCDM model.

Recently, a number of hydrodynamical simulations (Kravtsov, Gnedin & Klypin 2004; Governato et al. 2010; Trujillo-Gomez et al. 2015; Sawala et al. 2016; Garrison-Kimmel et al. 2019; Applebaum et al. 2021; Engler et al. 2021; Giri & Schneider 2021) have been trying to solve this issue by considering the baryonic feedback in the form of AGN, or stellar feedback to inhibit the overproduction of the small-scale structure. Nevertheless, incorporating the baryonic feedback self-consistently with the DM-only simulation is extremely non-trivial, and so far, the success is limited (Okamoto, Gao & Theuns 2008; Boylan-Kolchin, Bullock & Kaplinghat 2011; Governato et al. 2012; Teyssier et al. 2013). An alternate solution to these crises that has been proposed is to assume that the DM is ‘non-cold’ (Hu, Barkana & Gruzinov 2000; Bœhm, Fayet & Schaeffer 2001; Wang et al. 2014). The generic feature of such DM candidates is that the small-scale fluctuations in the matter distribution are suppressed relative to the standard CDM. There are several examples, for instance, warm dark matter (WDM)-like sterile neutrinos (Dodelson & Widrow 1994; Laine & Shaposhnikov 2008; Lovell et al. 2016), ultra-light scalars or axions also known as fuzzy DM (Hu, Barkana & Gruzinov 2000; Marsh & Silk 2014; Du, Behrens & Niemeyer 2017; Giri & Schneider 2022), self-interacting DM (Spergel & Steinhardt 2000; Vogelsberger et al. 2014), and others (Bœhm, Fayet & Schaeffer 2001; Dvorkin, Blum & Kamionkowski 2014; Wang et al. 2014).

Among the above, one of the most extensively studied candidates is the WDM with particle masses |$m_x \sim \mathcal {O}$| keV (see e.g. Blumenthal et al. 1984; Bode, Ostriker & Turok 2001; de Vega, Salucci & Sanchez 2012; Lovell et al. 2012). These particles are essentially thermal relics; hence, the small-scale suppression is entirely determined by particle mass (Lovell et al. 2012). Since these models have been widely studied in the literature, there exist straightforward methods to compute the abundance of DM haloes in addition to the modifications in the DM power spectrum (Bode, Ostriker & Turok 2001; Viel et al. 2005; Lovell et al. 2014; Schneider et al. 2014; Lovell 2020), both of which are crucial for our work.

Given the lower value of the WDM particle mass, these models erase the small-scale substructure, delay the structure formation and therefore solve the small-scale problems arising in the ΛCDM model.

As the WDM models delay the structure formation, it consequently delays the formation of the first stars. Therefore any observation related to the formation of the first stars could be used to constrain the mass of the WDM particles. As reionization is believed to start from the first generation of stars, a number of studies (Barkana, Haiman & Ostriker 2001; Somerville, Bullock & Livio 2003; Yoshida et al. 2003; Yue & Chen 2012; Pacucci, Mesinger & Haiman 2013; Dayal, Mesinger & Pacucci 2015; Dayal et al. 2017; Lopez-Honorez et al. 2017; Rudakovskyi et al. 2021; Schneider, Schaeffer & Giri 2023) used reionization related observations to put lower limits on the mass of the WDM in the range of 1.3–5 keV. Another observation related to the first generation of stars is the global 21 cm signal coming from cosmic dawn. A growing body of studies (Safarzadeh, Scannapieco & Babul 2018; Boyarsky et al. 2019; Chatterjee et al. 2019; Leo et al. 2020; Rudakovskyi, Savchenko & Tsizh 2020; Hibbard et al. 2022) exploited this signal to put lower limits on the WDM particle mass in the range of 3–6.6 keV. Furthermore, Viel et al. (2013); Iršič et al. (2017) and more recently Murgia, Iršič & Viel (2018) used Lyman Alpha (Ly α) forest power spectrum measurement using MIKE/HIRES spectrograph coming from high-resolution quasar spectra at redshifts z ∼2–5 to constrain the WDM mass in the range of 2.2−3.6 keV.¹ Very recently, the high-redshift observations coming from the JWST have been used to rule out WDM models with m_X < 1.5–2.0 keV (Maio & Viel 2023; Dayal & Giri 2023). Other than these high-redshift observations, Kennedy et al. (2014) used the count of dwarf galaxies to rule out WDM mass m_X < 2.3 keV, and finally, the most stringent constraint on the WDM mass comes from Nadler et al. (2021) ruling out m_X < 9.7 keV from a combined analysis of strong gravitational lenses and the Milky Way satellite galaxy population.

However, the aforementioned works have two limitations: (i) All of these works use either the reionization-related observations, the global 21 cm signal, or other high-redshift observations, but none of them combines all the data to put constraints on the WDM particles. (ii) Some of the above-mentioned works, where hydrodynamical or semi-numerical simulation is used, did not employ Markov chain Monte Carlo (MCMC)-based methods to quantify degeneracies (if any) between the mass of WDM particle and other cosmological and/or astrophysical parameters. For example, if we change any cosmological parameter that can delay the timing of structure formation of the Universe, that can, in principle, imitate the effect of lowering the mass of a WDM particle (remember that the lighter the mass of a WDM particle, the harder it is to start the structure formation).

To overcome both these issues, one has to first combine cosmic microwave background (CMB), reionization-related observations and a hypothetical data set of 21 cm signal and then vary all the free parameters (mass of WDM particles along with all the other cosmological and astrophysical parameters) simultaneously to put constraints on the mass of WDM particles and quantify, if any, degeneracy between different free parameters. In Chatterjee, Choudhury & Mitra (2021) (referred to as CCM21 hereafter), we have introduced an advanced MCMC-based parameter estimation package called cosmoreionmc which has all the above-mentioned features and therefore provides an ideal opportunity to carry out this investigation.

The rest of the paper is organized as follows. We describe in Section 2 the effect of incorporating the WDM in our galaxy formation model, theoretical modelling for reionization and 21 cm signal. Section 3 describes the findings of this work, and finally, Section 4 summarizes the work.

2 THEORETICAL MODELLING

2.1 Warm dark matter

It is well known that the effect of introducing the mass of WDM in our reionization and global 21 cm signal modelling will be manifested in the DM power spectrum, the halo mass function and consequently on any quantity that depends on either or both of them.

Following Bode, Ostriker & Turok (2001), the DM power spectrum of the WDM can be expressed as

where P_CDM is the usual CDM power spectrum and T_WDM is the transfer function given by Viel et al. (2005)

where μ = 1.2 and α is given by (Viel et al. 2005)

Following Lovell (2020), we write the halo mass function of the WDM as

with β = 2.3 γ = 0.8, δ = −1.0 N_CDM and N_WDM are the number of CDM and WDM haloes, respectively. The half-mode mass M_hm is given by

where |$\bar {\rho }_{m}$| is the background matter density, and the half-mode scale is given by

Following (Sheth & Tormen 1999), the halo mass function for the CDM is given by

where

with A = 0.3222, q = 0.707, p = 0.3, and ν is defined as

where D(z) is the well-known growth function.

2.2 Modelling reionization

For this work, we use the reionization model (CF model hereafter) developed in Choudhury & Ferrara (2005, 2006); Mitra, Choudhury & Ferrara (2012). As the detailed discussion of this model is beyond the scope of this paper, we briefly summarize here the main characteristics of this model.

• In this model, the overdensity of the intergalactic medium (IGM) is described using a lognormal distribution in the low-density regions and as a power-law distribution in the high-density region following the treatment presented in Miralda-Escudé (2003). We write the probability density function (PDF) of the overdensity Δ as

  $$\begin{eqnarray} P(\Delta) &=& \frac{A}{\sigma _b \Delta \sqrt{2 \pi }} \exp \left[-\frac{(\ln \Delta -\mu)^2}{2 \sigma _b^2}\right] \mathrm{if} \Delta \lt \Delta _V, \\ &=& \frac{A}{\sigma _b \Delta _V \sqrt{2 \pi }} \exp \left[-\frac{(\ln \Delta _V-\mu)^2}{2 \sigma _b^2}\right] \left(\frac{\Delta }{\Delta _V} \right)^{\beta } \mathrm{if} \Delta \gt \Delta _V, \\ \end{eqnarray}$$

  (10)

  where the parameters A, μ, and Δ_V are determined by demanding continuity of the derivative of P(Δ) at the transition overdensity Δ_V, and by normalizing the volume and mass to unity. We choose β = −2.5, appropriate for high redshifts. The quantity σ_b is the rms linear mass fluctuations in baryons and is related to the WDM power spectrum as

  $$\begin{eqnarray} \sigma _b^2 = \int _0^{\infty } \mathrm{d}k~k^2~ \frac{P_\mathrm{WDM}(k)}{\left(1 + x_J^2 k^2\right)^2}, \end{eqnarray}$$

  (11)

  where x_J is the Jeans length, which depends on the IGM temperature. Note that the density PDF is sensitive to the value of m_X through σ_b. The most important feature of this model is its ability to calculate the ionization and thermal state of the IGM in the neutral and ionized regions for different species (i.e. hydrogen and helium) separately, simultaneously and self-consistently. Moreover, once all the low-density regions of the IGM are ionized, this model assumes the Universe to be completely ionized.
• The original CF reionization model (Choudhury & Ferrara 2005, 2006; Mitra, Choudhury & Ferrara 2012) assumes the source of reionization to be quasars, PopII and PopIII stars. Therefore, the total photon production rate at a redshift z is given by

  $$\begin{eqnarray} \dot{n}_{\rm ph} (z) = \dot{n}_{\rm ph, \rm stellar}(z)+ \dot{n}_{\rm ph, \rm QSO}(z) . \end{eqnarray}$$

  (12)

  While the quasar contribution can be calculated easily by computing their ionizing emissivities from the observed quasar luminosity function (LF) at z < 7.5 (Kulkarni, Worseck & Hennawi 2019), the calculation for stellar contribution becomes non-trivial if we consider both PopII and PopIII. However, CCM21 shows that the contribution of PopIII stars is negligible as long as we use CMB and quasar absorption-related observations to constrain different parameters. Also, as discussed later, while simulating the global 21 cm signal, we take into account the contribution only from PopII stars. Therefore, in this work, we take stellar contributions only from PopII stars and completely ignore the contributions from PopIII stars. The number of ionizing photon from stellar sources are hence computed using

  $$\begin{eqnarray} \dot{n}_{\rm ph, \rm stellar}(z)= \rho _b \mathbf {\epsilon } \frac{{\rm d}f_{\mathrm{coll}}}{{\rm d}t} \int ^{\infty }_{\nu _H} \left(\frac{{\rm d}N_{\nu }}{{\rm d}M} \right) {\rm d}\nu, \end{eqnarray}$$

  (13)

  where ν_H is the threshold frequency for hydrogen photoionization, ρ_b is the mean comoving density of baryons in the IGM, and ϵ = f_* × f_esc, where f_* and f_esc, respectively, denotes the star formation efficiency and the escape fraction of the ionizing photons. The quantity dN_ν/dM denoting the number of photons emitted per frequency range per unit mass of the star, depends on the stellar spectra and IMF of the stars (Choudhury & Ferrara 2005). Using a standard Salpeter IMF in the mass range |$1\!-\!100 \, {\rm M}_{\odot }$| with a metallicity of |$0.05 \, {\rm M}_{\odot }$|⁠, dN_ν/dM has been computed from the stellar synthesis models of Bruzual & Charlot (2003). We consider ϵ as a free parameter in our model and later constrain it using MCMC (discussed in Sections 3 and 4).
• Two of the observables that the CF model can predict and will be later used in our MCMC analysis are (i) the redshift distribution of Lyman-limit system (dN_LL/dz) and (ii) the hydrogen photoionization rate (Γ_PI). To calculate both the observable, the CF model first calculates the mean free path of the photons using

  $$\begin{eqnarray} \lambda _{\mathrm{mfp}}=\frac{\lambda _{0}}{[1-F_{v}(\Delta _{i})]^{2/3}} \end{eqnarray}$$

  (14)

  where λ₀ is a free parameter of the reionization model, and

  $$\begin{eqnarray} F_{V}(\Delta _i) = \int _0^{\Delta _i} \mathrm{d}\Delta ~ P(\Delta) \end{eqnarray}$$

  (15)

  is the volume fraction of the ionized region as a function of the overdensity Δ_i. It is clear that the mean free path is sensitive to m_X through P(Δ). The dependence of λ_mfp on the mass of the WDM particle is described in detail in the Appendix. Once, we compute λ_mfp, it is straightforward to calculate dN_LL/dz (Choudhury & Ferrara 2005, CCM21). Similarly, the Γ_PI can be calculated using

  $$\begin{eqnarray} \Gamma _{\mathrm{PI}}(z)=(1+z)^{3} \int ^{\infty }_{\nu _H} {\rm d}\nu ~ \lambda _{\mathrm{mfp}}(\nu ,z)~\dot{n}_{\mathrm{ph}}(z)~\sigma _{H}(\nu), \end{eqnarray}$$

  (16)

  where σ_H(ν) is the hydrogen photoionization cross-section and |$\dot{n}_{\mathrm{ph}}(z)$| is the photon production rate as described in equation (12).

Although the introduction of the WDM models does not change the mathematical framework used for the CF reionization model, it affects any quantity of the model that depends on the halo mass function (hmf). For example, the collapse fraction of the DM halo (appeared in equation 13), which depends on the hmf, will change and the modified form will be

where |$\bar {\rho }_{m}$| is the mean comoving density of DM, M_min(z) is the minimum mass for star-forming haloes which is determined by different cooling processes (such as atomic cooling, molecular cooling) and feedback processes (radiative feedback, mechanical feedback, chemical feedback, Lyman Warner feedback, etc.). In the reionization model considered here, we consider only the atomic cooling and on top of that, radiative feedback is incorporated using a Jeans mass prescription described in detail in Choudhury & Ferrara (2005).

2.3 Global 21 cm modelling

The sky averaged 21 cm global differential brightness temperature can be written as (Furlanetto, Oh & Briggs 2006a; Chatterjee et al. 2019)

where T_γ is the background radiation temperature, T_S is the neutral hydrogen spin temperature and x_H i is the neutral hydrogen fraction in the IGM. Under the assumption that the optical depth of the Ly α is very high in the epoch of cosmic dawn and the redshift range we are interested in (discussed later in Section 3.2), the spin temperature T_S is computed using

where T_K is the kinetic temperature of the IGM, x_α is the Ly α is the coupling coefficient.²

Although the kinetic temperature of the IGM computed in the CF reionization model is decided mainly by two processes, namely, the adiabatic cooling and the photoheating from UV photons, the moment we consider a 21 cm signal coming from cosmic dawn, we have to include the X-ray heating term in the temperature evolution equation. However, as we go towards lower redshift during EoR, the X-ray heating can be ignored, and the UV heating becomes dominant once the reionization starts; therefore, we turn off the X-ray heating in the reionization epoch.

The X-ray heating can be computed using (Mineo, Gilfanov & Sunyaev 2012)

where f_Xh,* = f_X × f_h × f_*. f_X is an unknown normalization parameter of our model. It takes into account any discrepancy between the properties of the locally observed galaxy and yet-to-observe high-redshift galaxy. f_h is the parameter corresponding to the fraction of the total X-ray photons that heat the IGM.

To calculate x_α, we first determine the background Ly α flux using

where f_{α, *} = f_α × f_*. f_α is an unknown efficiency parameter such that any uncertainties in the properties of the high-redshift galaxies can be absorbed in this. Further, the effect of any radiative cascading, generating any additional Ly α photons, will also be absorbed in this factor.³ To determine the upper limit z_max of the integral, we assume that all the continuum ionizing photons would be absorbed in the IGM and will not play any part in determining the Ly α radiation. z_max is calculated using (Chatterjee et al. 2020)

where ν_α is the Ly α frequency.⁴ The quantity |$\dot{n}_{\rm \nu ^{\prime }}(z^{\prime })$| is given by

with |$\frac{{\rm d}N_{\nu ^{\prime }}}{{\rm d}M}$| denoting number of photon per unit stellar mass at frequency ν′. Once, we determine the background Ly-α flux, the coupling coefficient can be computed using

where |$\rm S_{\alpha }$| accounts for the detailed atomic physics involved in the scattering process, and we take S_α = 1 (Furlanetto, Oh & Briggs 2006a).

2.3.1 Effect of WDM on 21 cm signal

As discussed earlier, any quantity that depends on the mass of the WDM particle will change due to the change in the mass of the WDM particle and, therefore, the global 21 cm signal as a whole will depend on the value of the WDM mass. In Fig. 1, we have shown the effect of changing WDM mass on the global 21 cm signal. It is evident from Fig. 1 that as we keep decreasing the mass of DM particles, the absorption trough of the global signal continues to shift towards lower redshifts. This is expected as the smaller value of DM particles will lead to a delay in structure formation and therefore cause the absorption trough of the global signal to occur in lower redshifts.

Figure 1.

The global 21 cm differential brightness temperature for the CDM, 7 keV, 3 keV, and 1.5 keV WDM models. As is evident, if we lower the mass of the DM particles, the absorption trough in the signal shifts to lower redshifts.

Open in new tabDownload slide

2.4 The CMB anisotropies

While describing the calculation of the CMB anisotropies in CCM21, we mentioned that we modify the publicly available python-wrapped CAMB Lewis (2013)⁵ to incorporate the reionization history implied by the CF reionization model rather than using the default redshift symmetric tanh model in the camb code. For this work also, we use the same modified camb code to generate the CMB anisotropy data.

2.5 The cosmoreionmc package

As the mass of the DM particle is considered as a free parameter for this work, we have modified our previously developed MCMC-based parameter estimation based package cosmoreionmc according to the methods described in Sections 2.1–2.2. Moreover, the version of the cosmoreionmc used here is more flexible compared to the original version as the inverse of the mass of WDM particles i.e. |$m^{-1}_{X}$| is treated as a free parameter.

Next, we will describe the results obtained from this study.

3 RESULT

Here, we present the results of our analysis on the parameter constraints obtained using cosmoreionmc.

3.1 Constraining mass of the WDM particles using CMB and Quasar absorption data

Using cosmoreionmc, we first obtain constraints on the mass of WDM particles using CMB and reionization-related observations while simultaneously varying the cosmological and astrophysical parameters. The free parameters for this analysis (referred to as CMB + Quasar hereafter) are

where the first five parameters are the usual cosmological parameters, ϵ, λ₀ are the free parameters of our reionization model, and |$m^{-1}_{X}$| is the inverse of the mass of the DM particles in the unit of |$\rm keV^{-1}$|⁠.

In this analysis, we include the reionization-related observations from quasar absorption spectra and the Planck 2018 observations. Data sets related to reionization used in this analysis are (i) photoionization rate Γ_PI data obtained from the combined analysis of quasar absorption spectra and hydrodynamical simulations (Calverley et al. 2011; Becker & Bolton 2013; D’Aloisio et al. 2018; Becker et al. 2021). (ii) The redshift distribution of Lyman-limit system dN_LL/dz (Prochaska, O’Meara & Worseck 2010; Songaila & Cowie 2010; Ribaudo, Lehner & Howk 2011; Fumagalli et al. 2013; O’Meara et al. 2013; Crighton et al. 2019), (iii) Measurement of the upper limit on the neutral hydrogen fractions coming from the dark fractions in quasar spectra (Jin et al. 2023) have been used as priors while calculating the likelihood. On top of that, with the recent studies of the large-scale fluctuations of the effective Ly α optical depth from high redshift quasar spectra (Becker et al. 2015; Eilers et al. 2017; Bosman et al. 2018; Eilers, Davies & Hennawi 2018; Choudhury, Paranjape & Bosman 2020), we put a prior that reionization has to be completed (Q_H ii = 1) at z ≥ 5.3.

The total Likelihood function for this analysis is given by

where

Here, |$\zeta ^{\mathrm{obs}}_{\alpha }$| represents the set of N_obs observational data related to photoionization rates and the distribution of the Lyman-Limit system whereas |$\zeta ^{\mathrm{th}}_{\alpha }$| represents the values from the theoretical model. The σ_α denotes the observational error bars. |$\mathcal {L}_{\mathrm{Pl}}$| is the log-likelihood function corresponding to the Planck 2020 observations (Planck Collaboration 2020).

We assume a broad flat prior for all the eight free parameters. For |$m^{-1}_{X}$| (in unit of |$\rm keV^{-1}$|⁠), we take the flat prior range to be [0.0, 1.0] which allows us to explore m_X in the range [∞, 1.0] keV. In order to explore the parameter space with MCMC chains, we use 32 walkers taking 10⁶ steps. Before producing any result, the convergence of the MCMC chains is ensured using a detailed autocorrelation analysis as described in CCM21 and Foreman-Mackey et al. (2013).

The one-dimensional (1D) marginalized posterior distribution of the |$\left(\frac{m_{X}}{keV}\right)^{-1}$| is shown (in red) in Fig. 2. This figure shows that this analysis rules out WDM particles with m_X < 2.8 keV at 95 per cent confidence level. The result is in close agreement with the constraints coming from Baur et al. (2016) (they constrain m_X < 2.96 keV at 95 per cent confidence level). However, it is slightly weaker than that of Viel et al. (2013) (they rejected WDM particles with m_X < 3.3 eV) and is comparable with the constraints coming from Murgia, Iršič & Viel (2018). Based on the combined observations of medium-resolution spectra of the XQ-100 sample observed with the X-shooter spectrograph (z ∼ 4) and high-resolution spectra of the z ∼ 5 QSOs obtained with the HIRES/MIKE spectrographs, Murgia, Iršič & Viel (2018) rejected the WDM particles with m_X < 2.2−4.1 keV (at 95 per cent confidence level). The exact constraint on WDM mass depends on their assumption regarding the temperature evolution of the IGM. As discussed in Murgia, Iršič & Viel (2018), one has to keep in mind that these constraints on the WDM mass also depend on the choices of the priors on the IGM thermal history and that the different priors can significantly alter these mass limits.

Figure 2.

1D marginalized posterior distribution of the constraints on |$m^{-1}_{X}$| from different scenarios. The red curve represents the case when CMB and Quasar observations are used. The magenta curve depicts the scenario when the added 21 cm signal is simulated from the CDM model, whereas the green curve represents the case with the hypothetical 21 cm signal computed from the 7 keV WDM model.

Open in new tabDownload slide

For a detailed understanding of the constraints on different parameters and their correlation, the posterior distribution of different parameters is shown in Fig. 3. The most important point to note from this figure is that the inverse of WDM mass |$m^{-1}_{X}$| has a correlation with ϵ (and anticorrelation with λ₀). It is clear from |$m^{-1}_{X}-\epsilon$| subplot in Fig. 3 that a larger value of |$m^{-1}_{X}$| (lower value of m_X) requires a larger value of ϵ. This is only to be expected because a larger value of |$m^{-1}_{X}$| implies a delayed structure formation, and the only way to compensate for this is to have a higher ϵ to enhance the reionization process and therefore match with the observations. It is also important to note that this correlation is more prominent for a higher value of |$m^{-1}_{X}$| (smaller value of m_X). This is due to the fact that a smaller value of |$m^{-1}_{X}$| is practically indistinguishable from CDM. Apart from this, ϵ and λ₀ shows strong anticorrelation, this is because λ₀ and ϵ comes as a product at the time of calculating Γ_PI. So, to keep their product unchanged (necessary to match with the observations), if one parameter increases, the other parameter has to decrease and vice versa. As λ₀ and ϵ are strongly correlated and we have already seen that |$m^{-1}_{X}$| and ϵ are correlated, it is only to be expected that λ₀ and |$m^{-1}_{X}$| will be anticorrelated.

Figure 3.

The marginalized posterior distribution of 8 free parameters obtained for the CMB + quasar case. Two-dimensional plots in the figure show the joint probability distribution (confidence contours at 68 per cent and 95 per cent) of any two parameters. It is also clear from the two-dimensional plots between |$m^{-1}_{X}$| and other free parameters that there exists no correlation between them.

Open in new tabDownload slide

3.2 Constraining WDM with CMB, quasar, and a hypothetical global 21 cm signal

Next, we focus on seeing the potential of the global 21 cm signal to put tighter constraints on the mass of WDM particles when added along with the CMB and quasar absorption data. To this aim, we first generate a mock 21 cm signal and then add it with the CMB and quasar data. While constructing the hypothetical signal, we follow the procedure outlined in CCM21. To briefly summarize, we generate the mock signal in the frequency range of 55–235 MHz with frequency channels of width 0.5 MHz. In each frequency channel, we add a Gaussian noise of zero mean and standard deviation σ_i = 10mK to the theoretical signal. Note that this assumed noise is lower compared to what was found in the EDGES or SARAS experiment but is certainly achievable with a longer integration time. As pointed out in CCM21, with a noise similar to these experiments, constraints on the cosmological parameters remain similar to that of the Planck limit. Also, the frequency coverage of our hypothetical signal, especially on the higher frequency end, is wider than that of the EDGES or SARAS-3 experiment. The wider frequency coverage of our hypothetical signal is essential so that the signal be present during the EoR epoch (i.e. the 21 cm signal becomes sensitive to the reionization history), which will provide a more degeneracy-breaking potential to this signal while constraining different parameters. The input parameters (common to both reionization and 21 cm signal) used while generating the mock 21 cm signal are the best-fitting values of different cosmological and reionization model parameters from the CMB + quasar analysis. As is evident, the input values of other free parameters related to 21 cm estimation, which did not appear in the CMB + quasar analysis, are f_Xh,* and f_α,*. To be consistent with the values estimated from the low-redshift observations (Furlanetto, Oh & Briggs 2006b), we take both f_X and f_α to be equal to 1.0 and f_h = 0.2. We take f_* to be 0.01 consistent with our earlier works (see e.g. Mitra, Choudhury & Ferrara 2015; Mitra, Choudhury & Ratra 2018; Chatterjee, Choudhury & Mitra 2021). Note that with f_X = 1.0, f_h = 0.2, and f_* = 0.01, the input parameter f_Xh,* becomes 0.002 and f_α,* becomes 0.001. Since the actual value of m_X is not known, we explore two scenarios to simulate the future data and for making the forecasts: one where the DM is the usual CDM (referred to as CMB + quasar + 21cm_CDM hereafter) and another where m_X = 7 keV (referred to as CMB + quasar + 21cm_{7 keV WDM} hereafter).

The ten free parameters for the joint analysis, including the 21 cm signal, are

In the presence of the mock 21 cm observations, the log-likelihood becomes

where |$\mathcal {L}_{21}$| is the loglikelihood corresponding to mock observational data. Of course, the likelihood corresponding to 21 cm signal, |$\mathcal {L}_{21}$|⁠, will depend on whether the mock data is generated with |$m_{X} \xrightarrow []{} \infty$| or m_X = 7 keV as discussed below.

3.2.1 CMB+Quasar + 21cm_CDM

In this case, the likelihood term |$\mathcal {L}_{21}$| in equation (28) will become

where |$\delta T_b^{\mathrm{mock, CDM}}(\nu _i)$| is the mock brightness temperature data generated using the CDM model.

Once the MCMC run fulfils the convergence criteria, the 1D marginalized distribution of the quantity |$\left(\frac{m_{X}}{keV}\right)^{-1}$| is shown (in magenta) in the Fig. 2. It is evident from this figure that the inclusion of 21 cm data forecasts the constraints to be m_X > 7.7 keV (95  per cent confidence level), which is much more stringent than that derived from the CMB + quasar case. It is slightly weaker than the constraints coming from Nadler et al. (2021) (they constrain m_X < 9.7 keV) and is even stronger than the result obtained in Murgia, Iršič & Viel (2018) (as mentioned earlier, their most stringent constraint on the WDM particles comes out to be m_X < 4.1 keV). Meanwhile, Fig. 4 shows the constraints and the posterior distribution of all the free parameters used in this analysis. Unlike CMB + quasar case, here, |$m^{-1}_{X}$| does not show any correlation with any of the free parameters. This is because, in this case, the allowed range of values of m_X is very high (i.e. m_X > 7.7 keV), and these high m_X WDM models are practically indistinguishable from CDM. Note that in CMB + quasar case, only the low value of m_X shows the correlation/anticorrelation.

Figure 4.

The marginalized posterior distribution of 10 free parameters obtained for the CMB + quasar + 21cm_CDM case. Two-dimensional plots in the figure show the joint probability distribution (confidence contours at 68 per cent and 95 per cent) of any two parameters.

Open in new tabDownload slide

We also note from the 1D posterior distribution of f_Xh,* and f_α,* (bottom row of Fig. 4) that their best-fittng values are 0.002 and 0.01, respectively. This is expected as those were the input values of these two parameters at the time of creating the hypothetical signal. This result also shows that our MCMC analysis with cosmoreionmc is excellent at recovering the ‘true’ parameters of the mock signal.

3.2.2 CMB + quasar + 21cm_{7 keV WDM}

For mock data generated with m_X = 7 keV, the likelihood term |$\mathcal {L}_{21}$| in equation (28) will become

where |$\delta T_b^{\mathrm{mock, 7keV}}(\nu _i)$| is the mock brightness temperature data generated using 7keV WDM model.

After the completion of the MCMC run, the 1D posterior distribution of |$(m_X/\rm {keV})^{-1}$| is shown in green in Fig. 2. As shown from the posterior distribution, the 95 per cent confidence level of |$(m_X/\rm {keV})^{-1}$| comes out to be [0.1, 0.2] implying that the future 21 cm data should allow detection of the WDM particles if m_X ∼ 7 keV.

The free parameters’ posterior distribution is shown in Fig. 5. It is evident that |$m^{-1}_{X}$| is strongly correlated with both f_Xh,* and f_α,*. It is because with smaller and smaller value of m_X (higher and higher in |$m^{-1}_{X}$|⁠), structure formation gets delayed making the appearance of absorption trough at later and later redshifts, the only way to keep both the redshift and the depth of the absorption trough unchanged is to increase the value of f_Xh,* and f_α,*. It is also clear from the plot that f_Xh,* and f_α,* are strongly correlated with each other. It is due to the fact that the increase in f_α,* will try to make the absorption trough deeper, and the only way to compensate for that is to increase the value of f_Xh,*.

Figure 5.

The marginalized posterior distribution of 10 free parameters obtained for the CMB + quasar + 21cm_{7 keV WDM} case. Two-dimensional plots in the figure show the joint probability distribution (confidence contours at 68  per cent and 95  per cent) of any two parameters.

Open in new tabDownload slide

In Fig. 6, we demonstrate the comparison between the mock signal and the recovered signal coming from the MCMC chain after its convergence. It is clear that the best-fitting signal recovered from the MCMC run is in excellent agreement with the mock data in both cases.

Figure 6.

Comparison between the mock 21 cm data and the signal recovered from the MCMC run. The left panel shows the case when the mock signal is produced from the CDM model, and the right panel represents the scenario with the mock data produced from the 7 keV WDM model. In both the panels, the black, magenta and cyan curves represent the mock data, the best-fitting model and models corresponding to 1000 random samples from the MCMC chain, respectively.

Open in new tabDownload slide

4 CONCLUSION AND DISCUSSION

In this work, with the help of our previously developed MCMC-based parameter estimation pipeline cosmoreionmc, we explored three different scenarios to constrain the mass of WDM particles and also shed light on the long-overlooked issue of the degeneracy between cosmological and astrophysical parameters while constraining the mass of the WDM particles. First, we demonstrate that when CMB and quasar absorption-related observations are used along with the CMB angular power spectrum observations, the WDM particles with m_X < 2.8 keV can readily be ruled out. In the next step, we add a hypothetical 21 cm global signal along with the CMB and quasar absorption-related observations to check if adding a 21 cm signal can put more stringent constraints on the mass of WDM compared to the already existing constraints. To this end, we generate two mock 21 cm signals, one with m_X = 7 keV and the other assuming the usual CDM model (i.e. m_X = ∞). For the first scenario, the forecasts give (m_X/keV)⁻¹ in the range [0.1, 0.2] (95 per cent confidence level) implying that the future 21 cm data should allow detection of the WDM particles if m_X ∼ 7 keV. In the second case, the inclusion of 21 cm data forecasts the constraints to be m_X > 7.7 keV (95 per cent confidence level), much stronger than the present ones.

Finally, we will discuss some of the caveats of our analysis presented here. First, we take all the free astrophysical efficiency parameters, e.g. the escape fraction, the X-ray heating efficiency, and the Ly α flux efficiency, to be constants i.e. they do not change with redshift or halo mass. But it is entirely possible that they are not constants as we have assumed here, and therefore the constraints on m_X could be different from the analysis presented here. However, in spite of the simple assumptions, our work highlights the importance of the global 21 cm experiments in constraining the WDM particle mass.

In future work, we are planning to include more observational data e.g. UVLF data from the JWST observations (Bouwens et al. 2021, 2023; Harikane et al. 2022; Naidu et al. 2022; Harikane et al. 2023) and then revisit the constraints on different parameters. We are also considering using a more accurate reionization model to eliminate some of the simplified assumptions used in the code. For example, we ignore molecular cooling completely in our model, despite the fact that molecular cooling for DM halos is an important mechanism that can change the PopIII star formation rate inside a halo. In addition to that, we are also planning to include redshift/halo mass dependency in the efficiency parameters used in this analysis.

ACKNOWLEDGEMENTS

AC and TRC acknowledge the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700. AC would also like to thank Aditya Chowdhury for his constant inspiration to write this paper.

DATA AVAILABILITY

The observational data used here are taken from the literature and the code underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: DEPENDENCE OF MEAN FREE PATH OF PHOTONS ON THE WDM PARTICLE MASS

We compare in Fig. A1 the dependence of the mean free path (λ_mfp) on the mass of the WDM particle. For this exercise, we fix the hydrogen photoionization rate Γ_PI to a value of |$0.3 \times 10^{-12} \, {\rm sec^{-1}}$| and the IGM temperature to 10⁴ K for all the DM models. We study the redshift evolution of λ_mfp only in the redshift range where the Universe is fully ionized in all the DM models. Fixing these quantities ensures that any variation in (λ_mfp) arising from differences in the ionization and thermal histories for the different DM models is absent, the mean free path depends only on the density distribution of the IGM. It is evident from the figure that the λ_mfp obtained in our analysis increases with the decreasing mass of the WDM particle. Recently, Cain et al. (2023) have presented similar analysis using hydrodynamic simulations. We would like to point out that our result regarding the dependence of λ_mfp on the WDM particle mass, i.e. increase of λ_mfp with decrease in WDM particle mass, from our extremely efficient semi-analytic model indeed matches the overall trend obtained from their detailed hydrodynamical simulation.

Figure A1.

Redshift evolution of the mean free path of the ionizing photon (λ_mfp) for CDM and WDM models with different particle mass. The green, orange, and blue curves, respectively, denote CDM, 7keV and 3keV WDM. To make sure that the λ_mfp depends only on the density of the IGM, we kept the temperature of the IGM fixed at 10⁴ K and |$\Gamma _{\rm PI} = 0.3 \times 10^{-12} \, {\rm sec^{-1}}$| in a completely ionized Universe. As is obvious, the mean free path λ_mfp decreases with increasing m_X.

Open in new tabDownload slide