Impacts of new small-scale N-body simulations on dark matter annihilations constrained from cosmological 21cm line observations

In scenarios of dark matter (DM) for the weakly-interacting massive particle (WIMP) Jungman:1995df , a DM particle with mass of the weak scale should have a thermally-averaged annihilation cross-section of the order of $\langle\sigma v\rangle=2\times 10^{-26}{\rm cm}^{3}/{\rm sec}$ Steigman:2012nb ; Saikawa:2020swg . This value is called the “canonical annihilation cross-section” and is obtained to agree with the observed abundance of DM which was produced from the thermal bath in the early Universe by the thermal freezeout mechanism (e.g, see Ref. Saikawa:2020swg and references therein). Nowadays, a lot of new experimental projects have been proposed to observe signatures of annihilating DM, which heightens the momentum toward cosmologically verifying the existence of the WIMP DM.

On the other hand, it is still an open question how density fluctuations of DM can evolve nonlinearly at high redshifts $z={\cal O}(10)-{\cal O}(10^{2})$ and form structures at small scales down to $k\sim 10^{6}-10^{7}{\rm Mpc}^{-1}$ because nonlinear evolution of the density fluctuations becomes important. Actually two of the authors (KK and RT) of this paper have performed the detailed N-body simulation of cold DM (CDM) in a separate paper Takahashi:2021pse . As a result, they are suggesting that the halo formation is really sizable even at such a high redshift epoch. The results of Ref. Takahashi:2021pse enables us to calculate the boost factor as a function of redshift with a sufficient level of precisions. It is a striking point that this result means that DM could have inevitably annihilated at the high redshifts, and the energy injection by the annihilating DM is expected to be the order of ${\cal O}(10^{-21})\mathrm{eV}/\mathrm{sec}/\mathrm{cm}^{3}$ at $z\sim 20$, which affects the absorption feature of the global 21cm line-spectrum (e.g., see Refs. Poulin:2016anj ; DAmico:2018sxd ; Mena:2019nhm ; Liu:2020wqz ; Bolliet:2020ofj and references therein).

So far, the Experiment to Detect the Global Epoch of Reionization Signature (EDGES) collaboration reported the observational data for the absorption feature of the cosmological global 21cm line-spectrum  Bowman:2018yin at around $z\sim 17$. In a pioneering work DAmico:2018sxd , by using this data the authors obtained upper bounds on the annihilation cross-section of DM in order not to reduce the trough of the absorption feature due to the extra heatings by the annihilations. The boost factors adopted in DAmico:2018sxd were taken from the values from the following two papers with assuming somehow aggressive levels of model-dependent approximations: i) Ref Ripamonti:2006gq in which the Press-Schechter mass function formalism Press:1973iz with the Navarro-Frenk-White (NFW) halo profile Navarro:1996gj was assumed, and ii) Ref. Liu:2016cnk in which a halo model with the Einasto profile Einasto:1965czb was assumed. ¹¹1About other related works to constrain DM from the data of EDGES, see also the following papers and references therein, Refs.Yang:2018gjd ; Jia:2019yhr ; Liu:2020wqz for annihilations, Refs.Clark:2018ghm ; Mitridate:2018iag ; Liu:2020wqz ; Bolliet:2020ofj for decays, and Refs. Tashiro:2014tsa ; Xu:2018efh ; Munoz:2018pzp ; Slatyer:2018aqg for DM-baryon interactions. And see also Refs. Pospelov:2018kdh ; Moroi:2018vci ; Yoshiura:2018zts for productions of additional photons to fit the EDGES data.

On contrary to those previous works, we compute the boost factor as a function of $z$ by adopting the raw numerical data of the detailed N-body simulations at the small scales reported in Takahashi:2021pse . By using this boost factor, we update the upper bounds on annihilating DM with the modes into $W^{+}W^{-}$, $b\bar{b}$, $e^{+}e^{-}$ and $\gamma\gamma$ as conservatively as possible.

This paper is organized as follows. In Section II we review how energy injection in general affects the evolution of the intergalactic medium (IGM). In Section III we review how high energetic injection from DM annihilation is deposited into IGM. In Section IV we describe how we estimate the annihilation boost factor based on our dedicated N-body simulation. In Section V we present our results. We conclude in the final section.

In this section, we review how energy injection in general takes part in evolution of IGM. For illustrative purpose, we here follow the simple description of hydrogen in IGM based on the effective three-level atom modelPeebles:1968ja ; Zeldovich:1969en ; Seager:1999km . In numerical calculation we present in Section V, we adopt the recombination code HyRec³³3https://pages.jh.edu/~yalihai1/hyrec/hyrec.html, which is based on the state-of-the-art effective multi-level atom model (See AliHaimoud:2010dx ; Chluba:2010ca for details). We in this paper focus only on hydrogen ionization and recombination assuming helium is neutral, which should be a good approximation as long as we are interested in the Dark Ages Liu:2016cnk (See also Liu:2019bbm ).

The evolution of ionization fraction, $x_{e}$, is then described by the following equation:

where $T_{m}$ and $T_{\gamma}$ are respectively the temperatures of gas and photon, $n_{\rm H}$ is the number density of hydrogen, $E_{0}\simeq 13.6\,{\rm eV}$ is the ionization energy of hydrogen, $E_{\alpha}=3E_{0}/4$ is the energy of Ly-$\alpha$ photon, $\alpha_{\rm H}$ is the case-B recombination coefficient and $\beta_{\rm H}$ is the corresponding ionization rate. The Peeble’s $C$-factor, which represents the probability that a hydrogen atom initially in the $n=2$ shell reaches the ground state without being photoionized, is given by

where $\Lambda\simeq 8.23$s^-1 is the two-photon decay rate of the hydrogen 2$s$-state, and $H(t)$ is the Hubble expansion rate. The last term in (1) represents the effects of energy injection, which we will describe shortly after.

The evolution of the gas temperature $T_{m}$ is described by the following equation:

where $\Gamma_{C}$ is the coupling rate of $T_{m}$ to $T_{\gamma}$, which is predominated by the Compton scattering,

where $\sigma_{T}$ is the Thomson scattering cross section, $a_{r}$ is the radiation constant, $m_{e}$ is the electron mass, and $f_{\rm He}$ is the number ratio of helium to hydrogen.

The last term in each of Eqs. (1) and (3), which is proportional to the energy injection rate per unit volume per time, $dE_{\rm inj}/(dVdt)$, represents the effect of energy injection. As defined in Slatyer:2015jla ; Slatyer:2015kla the coefficients $f_{\rm ion}(t)$, $f_{\rm exc}(t)$, and $f_{\rm heat}(t)$ (collectively denoted by $\{f_{c}(t)\}$ hereafter) are the fractions of injected energy deposited into the hydrogen ionization, the hydrogen excitation and the heating of gas, respectively, which will be discussed in the next section.

In the case of DM annihilation, the energy injection rate is given as⁴⁴4 Here DM is assumed to be self-conjugate.

where $\bar{\rho}_{\rm DM}$ is the mean energy density of DM. $B(z)=\langle\rho_{\rm DM}^{2}\rangle/{\bar{\rho}_{\rm DM}}^{2}$ is the boost factor due to the inhomogeneity of DM distribution, which will be discussed in Section IV, $\langle\sigma v\rangle$ is the annihilation cross-section averaged over the phase space distribution and $m_{\rm DM}$ is the DM mass.

The deposition fractions $\{f_{c}(t)\}$ depend on particle constituents and their energy spectra from DM annihilation as well as their interaction with IGM. We compute $\{f_{c}(t)\}$ stepwise as follows:

1. 1.

   Once the primary annihilation processes (e.g. ${\rm DM}+\overline{\rm DM}\to{\rm SM}+\overline{\rm SM}$) are specified, the standard model (SM) particle constituents and their energy spectra of the final state can be computed based on Monte Carlo event generators (e.g. PYTHIA⁵⁵5http://home.thep.lu.se/~torbjorn/Pythia.html and HERWIG⁶⁶6https://herwig.hepforge.org) which can simulate cascades of primary annihilation products into stable SM particles. In this paper, for $m_{\rm DM}$ above $5\,{\rm GeV}$ we adopt PYTHIA to compute the energy spectra (For details, we refer to Sjostrand:2007gs ; Bahr:2008pv ; Sjostrand:2014zea and reference therein).

   On the other hand, for $m_{\rm DM}$ below $5\,{\rm GeV}$, where we in this paper restrict ourselves to primary annihilation channels into $e^{+}e^{-}$ and $\gamma\gamma$, we omit final state radiations and adopt monochromatic energy spectra from the primary annihilation processes. Since the fraction of energy carried by final state radiations is small and primary annihilation products efficiently deposit their energy into IGM at low energy, our treatment should be a good approximation.

2. 2.

   Energetic electrons, positrons and photons ejected from DM annihilation subsequently interact with IGM. How those energetic particles lose their energy through interaction with IGM and affect ionization and heating of IGM have been studied by many authors, e.g., Shull:1985 ; Chen:2003gz ; Padmanabhan:2005es ; Ripamonti:2006gq ; Kanzaki:2008qb ; Slatyer:2009yq ; Kanzaki:2009hf ; Evoli:2012zz . Energetic electrons/positrons lose their energy on timescales shorter than the Hubble time. Meanwhile, the timescale of energetic photons above $\simeq 10^{3}\,{\rm eV}$ and below $\simeq 10^{11}\,{\rm eV}$ can be longer than the Hubble time, which requires detailed computation of energy deposition over cosmological timescales. In this paper, we adopt the results of Slatyer:2015kla ⁷⁷7https://faun.rc.fas.harvard.edu/epsilon/, which treats the effects of energy injection at linear level. For full treatments including feedback of modification of IGM evolution in computation of $\{f_{c}(t)\}$ we refer to Liu:2019bbm .
   

Analytically, we can estimate the energy injection rate to be

This order-of-magnitude energy injection rate can affect the absorption feature of the global 21cm line spectrum Poulin:2016anj ; DAmico:2018sxd ; Mena:2019nhm ; Liu:2020wqz ; Bolliet:2020ofj .

DM annihilation is enhanced by inhomogeneity in DM distribution, which can be encapsulated in the boost factor $B(z)$. Even at redshifts as high as $z\gtrsim 15$ which we are focusing on in this paper, DM fluctuations at small scales have grown to be nonlinear. Therefore, to estimate $B(z)$, one needs to trace the nonlinear evolution of DM fluctuations. N-body simulations have been a powerful tool for this purpose.

Denoting $\delta(\mbox{\boldmath{$x$}};z)$ as the DM density contrast at a spatial position $x$ at redshift $z$, the boost factor is defined as $B(z)=1+\langle\delta^{2}(\mbox{\boldmath{$x$}};z)\rangle$. It can be recast using the Fourier transform as Serpico2012 ; Sefusatti2014

with $\Delta^{2}(k;z)\equiv k^{3}P(k;z)/(2\pi^{2})$ where $P(k;z)$ is the power spectrum of DM density fluctuations. The dimensionless power spectrum $\Delta^{2}(k;z)$ is almost flat at $k\gtrsim 10\,\,{\rm Mpc}^{-1}$, where horizon crossing takes place during the radiation dominated era. Therefore, the expression indicates that $B(z)$ is contributed from a wide range of scales. This necessitates that N-body simulations should be performed with a variety of box-sizes, enabling DM fluctuations to be resolved at relevant scales. However, so far there have been few studies performing N-body simulations as such focusing on redshifts of our interest. In Takahashi:2021pse , some of the authors of this paper have addressed this issue by performing dedicated N-body simulations.

As presented in Takahashi:2021pse , a suite of cosmological N-body simulations with a variety of box-sizes (i.e., side lengths of cubic boxes) is performed. The box-sizes range from $1\,{\rm kpc}$ to $10\,{\rm Mpc}$ to cover a wide range of scales, $k\simeq 1$–$10^{7}\,\,{\rm Mpc}^{-1}$. The simulations are comprised of $2560^{3}$ collisionless particles. The initial linear power spectrum is prepared using the transfer function Yamamoto1998 with the free streaming damping of DM particles at $k_{\rm fs}=10^{6}\,\,{\rm Mpc}^{-1}$ Green2004 . Initial conditions of the simulations are set at redshift $z=400$ based on the second-order Lagrangian perturbation theory Crocce2006 ; Nishimichi2009 . We employed the gravity solver GreeM Ishiyama2009 to follow the nonlinear gravitational evolution. To compute the boost factor, $\Delta^{2}(k)$ is constructed by connecting estimated DM power spectra at different wave number bands that depend on box-sizes of simulations. For more details of the simulations and analyses, we refer to Takahashi:2021pse .

Figure 1 shows $B(z)$ computed by using data of the simulations. As references, the figure also shows $B(z)$ computed based on linear perturbation theory as well as one based on the halo model in Evoli:2014pva , which is referred to as the ”Boost 1” model in DAmico:2018sxd .

Here, in the linear theory, $B(z)-1$ simply evolves as $(1+z)^{-2}$.

Compared to other estimations, e.g. the ”Boost 1” model adopted in DAmico:2018sxd as a conservative choice, our $B(z)$ is smaller than it at $z\lesssim 50$. This results in suppressed DM annihilation rate and hence may lead to much more conservative upper bounds on dark matter annihilation cross-section.

In figure 2, we demonstrate the evolution of $x_{e}(z)$ and $T_{m}(z)$ in the presence of DM annihilation with DM mass $m_{\rm DM}=100\,{\rm GeV}$. We here assume that there is no significant heating from astrophysical sources.

The differences between the results of the linear theory and the N-body simulation for the $W^{+}W^{-}$ and $b\bar{b}$ emissions are larger than the ones for the line $e^{+}e^{-}$ and $\gamma\gamma$ emissions. That is because more soft daughter electromagnetic particles such as $e^{+}e^{-}$ or $\gamma\gamma$ are produced through cascade decays of unstable mesons and baryons in cases for the $W^{+}W^{-}$ and $b\bar{b}$ emissions, compared with the cases for the line $e^{+}e^{-}$ and $\gamma\gamma$ emissions. Then, the energy-deposition is more efficient for $W^{+}W^{-}$ and $b\bar{b}$ due to the delayed deposition Liu:2018uzy .

$T_{m}$ is increased compared to cases where energy injection from dark matter annihilation is absent, which should result in modified evolution of the spin temperature, $T_{s}$, associated with the hyperfine splitting of neutral hydrogen ground states. This allows us to constrain DM annihilation cross-section from observations of differential brightness temperature of redshifted 21 cm line emission

before reionization (See e.g. Furlanetto:2006jb ). The time-evolution of $T_{s}$ in general depends on relative couplings of $T_{s}$ with $T_{\gamma}$, $T_{m}$ and the color temperature $T_{c}$, which is the effective temperature associated with background Lyman-$\alpha$ radiation. As IGM is always optically thick for Lyman-$\alpha$ radiation during the cosmological epoch we are interested in, it is reasonable to assume $T_{c}\approx T_{m}$. The fact that EDGES has reported a global absorption signal Bowman:2018yin

indicating $T_{m}<T_{\gamma}$, then the upper bounds on the DM annihilation cross-section should exist. This is because DM annihilation in general suppresses the absorption amplitude of global 21cm signals by ionizing and heating IGM. In Figure 2, the differences between the solid and dotted lines look much smaller than the ones directly-expected from the face values of the difference between the two lines in the boost factors from our N-body simulations and linear perturbation calculations in Figure 1. That is because the delayed deposition occurred as was discussed in Ref. Liu:2018uzy ; Basu:2020qoe .

We in particular obtain conservative upper bounds on DM annihilation cross-section by maximizing the absorption amplitude in the absence of DM annihilation DAmico:2018sxd . This can be realized by assuming no heating of IGM other than DM annihilation. We also assume tight coupling of spin temperature to gas temperature via Lyman-$\alpha$ pumping (Wouthuysen-Field effect Wouthuysen:1952 ; Field:1959 ), which can also maximize the absorption depth. Figure 3 shows upper bounds on DM annihilation cross-section based on this strategy. Our baseline calculation in the absence of DM annihilation, gives $T_{\rm 21cm}\simeq-230\,{\rm mK}$. We put upper bounds on the annihilation cross-section of the DM by requiring $T_{\rm 21cm}\leq-75\,{\rm mK}$, which correspond to the 2$\sigma$ bound given uncertainties of EDGES.

From this figure, we find that the upper bounds on the annihilation cross sections, which are conservatively-obtained in this study for the $b\bar{b}$, $e^{+}e^{-}$ and $\gamma\gamma$ modes, are milder than the ones in the results of DAmico:2018sxd . That is because the boost factor we adopted is smaller than the ones in DAmico:2018sxd . When we assume the canonical value of the annihilation cross-section, $\langle\sigma v\rangle=2\times 10^{-26}{\rm cm}^{3}/{\rm sec}$, we can exclude the masses of DM for each mode to be $m_{\rm DM}<15$ GeV ($b\bar{b}$), $m_{\rm DM}<3$ GeV ($e^{+}e^{-}$), and $m_{\rm DM}<1$ GeV ($\gamma\gamma$) at 95$\%$ C.L.

In this paper we have revisited the possible constraints on annihilation cross-sections of DM from the observations on the cosmological global 21cm line-spectrum reported by EDGES. By adopting the latest data of high-redshift dark-matter halo-formations ($z={\cal O}(10)-{\cal O}(10^{2})$) performed by the detailed N-body simulations at the small scales, we have updated the boost factor of the annihilating DM due to the clumpiness.

With this updated value of the boost factor, we obtained the more conservative upper bounds on the annihilation cross-sections than the ones reported in the previous work. In this study, we can exclude the masses of DM for $m_{\rm DM}<15$ GeV ($m_{\rm DM}<3$ GeV) at 95$\%$ C.L. for the mode into $b\bar{b}$ ($e^{+}e^{-}$) by assuming the canonical value of the annihilation cross-section, $\langle\sigma v\rangle=2\times 10^{-26}{\rm cm}^{3}/{\rm sec}$. These bounds obtained from the global 21cm spectrum are cosmologically robust along with the ones from CMB Slatyer:2015jla and BBN Kawasaki:2015yya because they do not depend on local astrophysical uncertainties.

In the future, we can improve sensitivities on the constraints on the annihilation cross-section by adopting more precise observational data which are expected to be reported by new projects such as HERA Beardsley:2014bea , SKA SKAspec , Omniscope Omniscope or DAPPER Burns:2021ndk .