# Dark matter-baryon scattering effects on temperature perturbations and
implications for cosmic dawn

Kathleen Short  José Luis Bernal  Kimberly K. Boddy  Vera Gluscevic  and Licia
Verde

###### Abstract

The nature of dark matter remains unknown, but upcoming measurements probing
the high-redshift Universe may provide invaluable insight. In the presence of
dark matter-baryon scattering, the suppression in the matter power spectrum
and the colder mean gas temperature are expected to modify the evolution of
cosmic dawn and reionization. However, the contributions from such
interactions to the baryon and dark matter temperature perturbations have been
neglected thus far. In this work, we derive these contributions, evolve the
cosmological perturbations until the end of the dark ages and show that they
may have a significant impact in the beginning of cosmic dawn. In particular,
we find that the amplitude of the temperature power spectrum at large scales
can change by up to 1–2 orders of magnitude and that the matter power spectrum
is further suppressed with respect to $\Lambda$CDM by $5$–$10\%$ at $k\sim
200\,{\rm Mpc^{-1}}$ compared to the computation ignoring these contributions
for scattering cross sections at current CMB limits. As a case example, we
also compute the HI power spectrum from the dark ages, finding significant
differences due to the changes in the temperature and ionization fraction
power spectra. We argue that these new contributions must be included in
studies of this dark matter model relying on cosmic dawn and reionization
observables.

##  1 Introduction

Cosmological observables provide very tight constraints on the abundance of
matter in the Universe and in particular on the abundance of the dark matter
(DM) [1, 2, 3, 4, 5]. Understanding the microscopic nature of DM by exploring
phenomenology beyond the known gravitational interactions is one of the
primary goals in cosmology and astroparticle physics. Regarding particle DM,
different strategies include direct detection [6] and collider [7]
experiments, but cosmological and astrophysical searches often provide
sensitivities at different masses ranges [8]. However, despite decades of
searches yielding progressively stronger bounds on the nature of the DM, the
solution remains elusive.

Many well-motivated particle candidates for DM, including weakly-interacting
massive particle (WIMP) models [9, 10], predict elastic scattering between the
DM and baryons which has cosmological consequences. While direct detection
experiments are highly sensitive to DM-nucleon interactions for DM masses
larger than a few GeV, cosmological probes cover complementary regions of
parameter space and are particularly sensitive to scattering of sub-GeV DM
with baryons [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29]. In a cosmological context, the interaction between DM and baryons
leads to transfer of heat and momentum between the two fluids which can affect
a range of physical scales and leave imprints on cosmological observables
throughout cosmic history. A key signature of the interactions is a
suppression of power in cosmological perturbations on small physical scales.

In cosmology, DM-baryon scattering is often characterised by a power law
dependence of the momentum-transfer cross section on the relative particle
velocity, $\sigma(v)=\sigma_{0}v^{n}$, with a power law index $n$ and an
unknown coefficient $\sigma_{0}$ that quantifies the strength of the
interaction at hand. If the dependence of the scattering cross section on the
relative velocity is steep enough ($n<-2$), its effect is larger later on in
cosmic history, thus impacting late-time observables more than the cosmic
microwave background (CMB) anisotropy. Moreover, the impact of the
interactions on the baryon temperature and matter clustering at small scales
is expected to affect the formation of the first stars.

Upcoming cosmological probes, in particular line-intensity mapping, will
provide new avenues to investigate the nature of DM. While post-reionization
measurements are promising in this regard (see e.g., [30, 31]), the HI
intensity mapping signal from cosmic dawn is expected to be especially
sensitive to the DM microphysics. In general, any model with clustering
suppression at small scales delays the collapse of the first halos and
therefore the formation of the first stars, affecting the time evolution of
cosmic dawn and reionization [32, 33, 34, 35, 36, 37, 38, 39, 40].
Furthermore, since the HI signal is highly sensitive to the gas temperature
through the HI spin temperature, any modification to the gas temperature
(either due to energy injection or heat transfer) is expected to greatly
modify the HI signal [41, 42, 43, 44, 45, 46, 47, 48, 49]. This has motivated
studies of the signatures of DM-baryon scattering in the HI intensity mapping
from the dark ages and cosmic dawn [50, 51, 46, 18, 44, 45, 52, 53, 54],
especially in attempts to explain the anomalous EDGES measurement [55].
Experiments such as EDGES [55], SARAS [56], and LEDA [57], which measure the
sky-averaged brightness temperature of the HI line at cosmic dawn, and HERA
[58] and SKA [59], which will provide measurements of its power spectrum, are
thus expected to return invaluable insight on the nature of the dark matter.

However, the direct contribution from DM-baryon scattering to the gas
temperature perturbations and the actual dark matter temperature perturbations
have been thus far ignored. Temperature perturbations are usually negligible
for CMB observations [60] but are critical for signals like the HI intensity
mapping, and they may affect the clustering of the first stars. Furthermore,
temperature perturbations appear in the sound speed terms of the DM and baryon
fluid evolution equations, modifying the clustering at scales below the
corresponding Jeans scales. Therefore, we deem the inclusion of these
contributions and variables in the evolution equations critical for a detailed
assessment of the effects of DM-baryon scattering at low redshift.

In this work, we compute for the first time the effects of DM-baryon elastic
scattering in the baryon temperature and ionization fraction perturbations, as
well as considering the DM temperature perturbations. We implement this
derivation in a modified version of the linear Boltzmann solver
class111https://github.com/lesgourg/class_public [61], developed for DM-baryon
scattering cosmologies in ref. [16], and self-consistently evolve the
perturbation equations up to the beginning of cosmic dawn (which we consider
to be at $z\sim 30$).222The specific start of cosmic dawn depends on the
cosmological and astrophysical model. However, this does not affect our main
conclusions. We show numerical results for the $n=-4$ case with cross section
values set to the 95% confidence level (C.L.) upper limits from CMB analyses
when all the dark matter interacts with baryons [16, 17]. As expected, we find
significant changes in the gas temperature and ionization fraction
perturbations due to the additional terms sourced by DM-baryon scattering.
Furthermore, there is an increased suppression of the baryon and total matter
power spectra relative to $\Lambda$CDM when compared to the prediction for DM-
baryon scattering models neglecting the newly-derived contributions to the
temperature perturbations. We argue that these changes are large enough to
have a sizeable impact in the beginning of cosmic dawn and must be taken into
account when setting the initial conditions for cosmic dawn. Once the first
stars form, the heating of the gas will be dominated by these sources and a
dedicated study, beyond the scope of this work, is required to follow the
evolution of the gas and dark matter temperature perturbations coherently with
the effects of the scattering in the matter clustering. As an example of the
impact of the new terms in an observable, we compute the HI power spectrum
from the dark ages.

This paper is organised as follows. In section 2, after reviewing the modified
Boltzmann equations for DM-baryon scattering cosmologies, we present the main
result of this work: we derive the evolution of the temperature and ionization
fraction perturbations sourced by the DM-baryon interactions. In section 3 we
present our numerical results, quantifying the effects on the DM and baryon
temperature and ionization fraction perturbations and the consequent effect on
the baryon and total matter power spectra. Finally, in section 4 we discuss
the phenomenological impact of the new temperature perturbation contributions
on cosmological observables and as a case example, we compute the effects on
the HI line-intensity mapping signal during the dark ages. We discuss and
conclude in section 5. Hereinafter we focus on the power spectrum at $z=30$,
but relevant results at $z=50$ can be found in appendix A, while the effects
of varying the scattering cross section and the fraction of interacting DM are
discussed in appendix B.

##  2 Evolution of perturbations with DM-baryon scattering

In this section, we derive the contributions to the perturbed DM and baryon
temperatures and ionization fraction arising from the DM-baryon collision
term, and incorporate these contributions into the modified Boltzmann
equations for DM-baryon interacting cosmologies. We consider DM-baryon
scattering of the form $\sigma(v)=\sigma_{0}v^{n}$, where $v$ is the relative
velocity between scattering particles.

Perturbations in the ionization fraction during recombination and the dark
ages affect the evolution of the baryon temperature and density perturbations,
which has a direct impact on the HI signal and structure formation. Here we
follow the nomenclature introduced by ref. [60] and refer to this effect as
“perturbed recombination”.

###  2.1 Modified Boltzmann equations

We review the modified Boltzmann equations accounting for scattering between
the DM and baryons, following the treatment presented in ref. [16]. Note that
here we work in the Newtonian gauge, but the equations in the synchronous
gauge are identical up to terms involving the metric perturbations. The
evolution equations for the DM and baryon density perturbations,
$\delta_{\chi}$ and $\delta_{b}$, and velocity divergences, $\theta_{\chi}$
and $\theta_{b}$, respectively, are given by

$\displaystyle\begin{split}\dot{\delta_{\chi}}={}&-\theta_{\chi}+3\dot{\phi}\\\
\dot{\delta_{b}}={}&-\theta_{b}+3\dot{\phi}\\\
\dot{\theta_{\chi}}={}&-\frac{\dot{a}}{a}\theta_{\chi}+c^{2}_{\chi}k^{2}\left(\delta_{\chi}+\delta_{T_{\chi}}\right)+k^{2}\psi+\widetilde{R}_{\chi}\left(\theta_{b}-\theta_{\chi}\right)\\\
\dot{\theta_{b}}={}&-\frac{\dot{a}}{a}\theta_{b}+c^{2}_{b}k^{2}\left(\delta_{b}+\delta_{T_{b}}\right)+k^{2}\psi+R_{\gamma}\left(\theta_{\gamma}-\theta_{b}\right)+\frac{\rho_{\chi}}{\rho_{b}}\widetilde{R}_{\chi}\left(\theta_{\chi}-\theta_{b}\right),\end{split}$
(2.1)

where $k$ is the wave number of a given Fourier mode; $a$ is the scale factor;
$\phi$ and $\psi$ are the Newtonian scalar metric potentials; $c_{\chi}$ and
$c_{b}$ are the DM and baryon sound speeds, respectively; $\rho_{\chi}$ and
$\rho_{b}$ are their energy densities; $R_{\gamma}$ is the momentum-transfer
rate for baryon-photon coupling due to Compton scattering; and
$\widetilde{R}_{\chi}$ is the (modified) momentum-transfer rate arising from
the DM-baryon interactions. The overdot denotes a derivative with respect to
conformal time. We have added the temperature perturbations
$\delta_{T_{\chi}}$ and $\delta_{T_{b}}$ ($\delta_{T_{\chi,b}}\equiv\delta
T_{\chi,b}/\bar{T}_{\chi,b}$, where $\bar{T}_{\chi,b}$ denotes the mean
temperatures) of the DM and baryon gas to the sound speed terms in the
momentum equations, as is done for the baryons in refs. [62, 63]. Previously
the DM temperature perturbations $\delta_{T_{\chi}}$ have been neglected
entirely; in this work, as will become clear in the following, we include it
for the first time. As we will see below, accounting for the DM temperature
perturbations will not only change the gas temperature power spectrum but also
the growth of matter overdensities at small scales, due to this additional
contribution to the DM sound speed.

At early times ($z>10^{4}$), the DM-baryon relative bulk velocity is small
compared to the relative thermal velocity dispersion, and the momentum-
transfer coefficient $R_{\chi}$ is approximately independent of the velocity
divergences such that

$R_{\chi}=a\rho_{b}\frac{Y_{H}\sigma_{o}}{m_{\chi}+m_{b}}\mathcal{N}_{n}\bar{v}_{\rm{th}}^{(n+1)},$
(2.2)

where $\mathcal{N}_{n}\equiv 2^{(5+n)/2}\Gamma(3+n/2)/(3\sqrt{2})$,
$\bar{v}_{\rm{th}}^{2}=T_{\chi}/m_{\chi}+T_{b}/m_{b}$ is the relative thermal
velocity dispersion squared, $T_{\chi}\,(m_{\chi})$ and $T_{b}\,(m_{b})$ are
the DM and baryon temperatures (masses), respectively, $Y_{H}$ is the mass
fraction of hydrogen, and we neglect helium scattering. This approximation
breaks down at lower redshifts for $n\leq-2$ as the relative bulk velocity
becomes non-negligible, resulting in non-linear Boltzmann equations and mode-
coupling. In order to restore linearity to the evolution equations, we follow
the prescription of ref. [16] in which a modified momentum-transfer rate is
introduced to include the effects of mode mixing. Defining

$\displaystyle\begin{split}V^{2}_{\rm{flow}}\equiv{}&\int_{0}^{k}\,\frac{dk^{\prime}}{k^{\prime}}\Delta^{2}_{\zeta}\left[\frac{\theta_{b}(k^{\prime},z)-\theta_{\chi}(k^{\prime},z)}{k^{\prime}}\right]^{2}\\\\[5.0pt]
V^{2}_{\rm{RMS}}\equiv{}&\int_{k}^{\infty}\,\frac{dk^{\prime}}{k^{\prime}}\Delta^{2}_{\zeta}\left[\frac{\theta_{b}(k^{\prime},z)-\theta_{\chi}(k^{\prime},z)}{k^{\prime}}\right]^{2},\end{split}$
(2.3)

where $\Delta^{2}_{\zeta}$ is the primordial curvature perturbation variance
per $\ln k$, the modified momentum-transfer rate is given by

$\widetilde{R}_{\chi}=R_{\chi}\left[1+\frac{V^{2}_{\rm{RMS}}(k,z)/3}{\bar{v}_{\rm{th}}^{2}}\right]^{(n+1)/2}\times{}_{1}F_{1}\left(-\frac{n+1}{2},\frac{5}{2},-\frac{V^{2}_{\rm{flow}}(k,z)}{2[\bar{v}^{2}_{\rm{th}}+V^{2}_{\rm{RMS}}(k,z)/3]}\right)\,,$
(2.4)

where ${}_{1}F_{1}$ is the confluent hypergeometric function of the first
kind. In the limit where $\langle V^{2}_{\chi b}\rangle\equiv
V^{2}_{\rm{flow}}+V^{2}_{\rm{RMS}}\ll\bar{v}^{2}_{\rm{th}}$,333$\langle
V^{2}_{\chi b}\rangle$ is the full variance of the relative bulk velocity,
integrated over all $k$. we recover $\widetilde{R}_{\chi}=R_{\chi}$, as
desired.

Before recombination, supersonic relative velocities between the DM and
baryons are generated, leading to a root-mean-square velocity of $\sim 30$ km
s-1 at the time of recombination. This supersonic relative velocity permits
the baryons to freely stream out of the DM potential wells, suppressing
density perturbations at scales below the characteristic advection scale over
a Hubble time at the time of decoupling ($k\sim 40\,\rm Mpc^{-1}$) [63, 64,
65, 66]. The relative velocity is coherent over scales smaller than $k\sim
0.3\,\rm Mpc^{-1}$, which allows the use of moving-background perturbation
theory to account for a formally non-linear term, by perturbing around a local
background value of the relative bulk velocity within coherent patches.
Following ref. [64], the evolution of small-scale perturbations after
recombination in the presence of a local relative bulk velocity $\bm{v}_{\chi
b}^{\rm bg}$ can be written to first order as444Note that, contrary to ref.
[64], we use derivatives with respect to conformal time and use the standard
definition of $\theta$ as the velocity divergence with respect to comoving
space.

$\displaystyle\begin{split}\dot{\delta}_{\chi}={}&i({\bf\it
v}^{\,\rm{bg}}_{\chi b}\cdot{\bf\it k})\,\delta_{\chi}-\theta_{\chi}\,\\\
\dot{\theta}_{\chi}={}&i({\bf\it v}^{\,\rm{bg}}_{\chi b}\cdot{\bf\it
k})\theta_{\chi}-\frac{\dot{a}}{a}\theta_{\chi}+k^{2}\phi+c^{2}_{\chi}k^{2}(\delta_{\chi}+\delta_{T_{\chi}})+\widetilde{R}_{\chi}(\theta_{b}-\theta_{\chi})\,\\\
\dot{\delta}_{b}={}&-\theta_{b}\,\\\
\dot{\theta}_{b}={}&-\frac{\dot{a}}{a}\theta_{b}+k^{2}\phi+c^{2}_{b}k^{2}(\delta_{b}+\delta_{T_{b}})+\frac{\rho_{\chi}}{\rho_{b}}\widetilde{R}_{\chi}(\theta_{\chi}-\theta_{b})\,\\\
\textrm{with}\quad\frac{k^{2}}{a^{2}}\phi={}&-\frac{3}{2}\frac{H^{2}_{0}}{a^{3}}(\Omega^{0}_{b}\delta_{b}+\Omega^{0}_{\chi}\delta_{\chi})\,.\end{split}$
(2.5)

In this regime, we are considering scales much smaller than the horizon and
times after recombination, and so we neglect the relativistic contributions to
the evolution equations and the baryon-photon coupling term. We assume that
${\bf\it v}^{\,\rm{bg}}_{\chi b}$ follows an isotropic Gaussian distribution
with variance given by $\langle V^{2}_{\chi
b}\rangle=V^{2}_{\rm{flow}}+V^{2}_{\rm{RMS}}$ and place ourselves in the local
baryon rest frame.

The small-scale patches with coherent relative velocities are correlated as a
relic of the correlation of the relative velocities after recombination.
Therefore, the small-scale fluctuations modulated by the relative velocities
also affect the large-scale fluctuations through uncorrelated, additive
quadratic terms [63, 65]. The quadratic terms at small scales are averaged
over regions of the size of the coherence scale of the relative velocity and
contribute to large-scale perturbations, following the correlation function of
the local relative velocities.

###  2.2 Temperature and ionization fraction perturbations

To derive the evolution equations for the DM and baryon temperature
perturbations, we begin from the first law of thermodynamics

$\frac{3}{2}dT_{\chi,b}-T_{\chi,b}d\log\rho_{\chi,b}=dQ_{\chi,b},$ (2.6)

where $dQ_{\chi,b}$ is the heating rate per particle. In the absence of DM-
baryon interactions, the only relevant source of heating is the standard
Compton collision heating rate for baryons scattering with CMB photons, given
by

$\dot{Q}_{C}=3\frac{\mu_{b}}{m_{e}}R_{\gamma}(T_{\gamma}-T_{b})=\frac{4a\sigma_{T}\rho_{\gamma}x_{e}}{m_{e}(1+x_{\rm
He}+x_{e})}(T_{\gamma}-T_{b})$ (2.7)

where $\mu_{b}=m_{H}(1+4x_{\rm He})/(1+x_{\rm He}+x_{e})$ is the mean
molecular weight of the baryons, $m_{e}$ is the electron mass, $T_{\gamma}$ is
the photon temperature, $x_{e}=n_{e}/n_{\rm H}$ is the free electron fraction,
$x_{\rm He}=n_{\rm He}/n_{\rm H}$ is the constant helium to hydrogen number
density ratio, $\rho_{\gamma}=a_{r}T_{\gamma}^{4}$ is the photon energy
density with radiation constant $a_{r}$, and $\sigma_{T}$ is the Thompson
cross section.

In our case, there is an additional contribution sourced by DM-baryon
interactions $\dot{Q}^{\rm iDM}$, for both the baryons and DM. Scattering of
the baryons with the colder DM leads to heat exchange between the two species,
cooling down the baryons and heating the DM. Moreover, the non-negligible
relative bulk velocity induces a drag force between the two fluids which
generates an additional temperature-independent heating term that dominates
over the baryon cooling at high redshift for masses $m_{\chi}\gtrsim 1$ GeV
[51]. The heat exchange rates due to DM-baryon interactions, accounting for a
non-negligible relative bulk velocity, were generalised for all $n>-5$ in ref.
[16]; in what follows, we substitute $V^{2}_{\chi b}$ with its averaged value
$\langle V^{2}_{\chi b}\rangle$ to maintain linearity, according to the
prescription set out in that work. Let us define

$\displaystyle\begin{split}\Gamma_{\chi}\equiv{}&2\frac{\bar{\rho}_{b}}{\rho_{b}}R^{\prime}_{\chi}\left[{}_{1}F_{1}\left(-\frac{n+3}{2},\frac{3}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{\rm{th}}}\right)-\frac{\langle
V^{2}_{\chi
b}\rangle}{3\bar{v}^{2}_{\rm{th}}}{}_{1}F_{1}\left(-\frac{n+1}{2},\frac{5}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{\rm{th}}}\right)\right]\,,\\\\[10.0pt]
\Gamma_{b}\equiv{}&\frac{2\mu_{b}}{m_{\chi}}\frac{\bar{\rho}_{\chi}}{\rho_{b}}R^{\prime}_{\chi}\left[{}_{1}F_{1}\left(-\frac{n+3}{2},\frac{3}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{\rm{th}}}\right)-\frac{\langle
V^{2}_{\chi
b}\rangle}{3\bar{v}^{2}_{\rm{th}}}{}_{1}F_{1}\left(-\frac{n+1}{2},\frac{5}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{\rm{th}}}\right)\right]\,,\\\\[10.0pt]
\mathcal{H}_{\chi}\equiv{}&2\frac{\bar{\rho}_{b}}{\rho_{b}}R^{\prime}_{\chi}\frac{m_{b}}{3}\langle
V^{2}_{\chi
b}\rangle{}_{1}F_{1}\left(-\frac{n+1}{2},\frac{5}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{th}}\right)\,,\\\\[10.0pt]
\mathcal{H}_{b}\equiv{}&\frac{2\mu_{b}}{m_{\chi}}\frac{\bar{\rho}_{\chi}}{\rho_{b}}R^{\prime}_{\chi}\frac{m_{\chi}}{3}\langle
V^{2}_{\chi
b}\rangle{}_{1}F_{1}\left(-\frac{n+1}{2},\frac{5}{2},-\frac{\langle
V^{2}_{\chi b}\rangle}{2\bar{v}^{2}_{th}}\right),\end{split}$ (2.8)

where the heat-transfer rate coefficient is given by
$R^{\prime}_{\chi}=R_{\chi}m_{\chi}/(m_{\chi}+m_{b})$. The heating rates due
to DM-baryon interactions may then be written as

$\displaystyle\begin{split}\dot{Q}^{\rm
iDM}_{\chi}={}\frac{3}{2}\Gamma_{\chi}\frac{\rho_{b}}{\bar{\rho}_{b}}(T_{b}-T_{\chi})+\mathcal{H}_{\chi}\frac{\rho_{b}}{\bar{\rho}_{b}}\,,\qquad\dot{Q}^{\rm
iDM}_{b}={}\frac{3}{2}\Gamma_{b}\frac{\rho_{\chi}}{\bar{\rho}_{\chi}}(T_{\chi}-T_{b})+\mathcal{H}_{b}\frac{\rho_{\chi}}{\bar{\rho}_{\chi}}\,.\end{split}$
(2.9)

Combining eqs. (2.6), (2.7) and (2.9), we obtain the full evolution equations
for the DM and baryon temperatures

$\displaystyle\begin{split}\dot{T_{\chi}}-\frac{2}{3}\frac{\dot{\rho}_{\chi}}{\rho_{\chi}}T_{\chi}={}&\Gamma_{\chi}\frac{\rho_{b}}{\bar{\rho}_{b}}(T_{b}-T_{\chi})+\mathcal{H}_{\chi}\frac{\rho_{b}}{\bar{\rho}_{b}}\,,\\\\[10.0pt]
\dot{T_{b}}-\frac{2}{3}\frac{\dot{\rho}_{b}}{\rho_{b}}T_{b}={}&\Gamma_{b}\frac{\rho_{\chi}}{\bar{\rho}_{\chi}}(T_{\chi}-T_{b})+\mathcal{H}_{b}\frac{\rho_{\chi}}{\bar{\rho}_{\chi}}+\frac{2\mu_{b}}{m_{e}}R_{\gamma}\left(T_{\gamma}-T_{b}\right).\end{split}$
(2.10)

Expanding eq. (2.10) and linearising using
$\dot{\rho}_{i}/\rho_{i}=-3\dot{a}/a+\dot{\delta}_{i}$ for species $i$, we
find the homogeneous solution gives the evolution of the background
temperatures

$\displaystyle\begin{split}\dot{\bar{T}}_{\chi}+2\frac{\dot{a}}{a}\bar{T}_{\chi}={}&\Gamma_{\chi}(\bar{T}_{b}-\bar{T}_{\chi})+\mathcal{H}_{\chi}\,,\\\\[10.0pt]
\dot{\bar{T}}_{b}+2\frac{\dot{a}}{a}\bar{T}_{b}={}&\Gamma_{b}(\bar{T}_{\chi}-\bar{T}_{b})+\mathcal{H}_{b}+\frac{2\mu_{b}}{m_{e}}R_{\gamma}\left(\bar{T}_{\gamma}-\bar{T}_{b}\right)\,,\end{split}$
(2.11)

which matches the results found in previous work [51, 16]. In figure 1 we show
the evolution of the background temperatures as function of redshift for
different DM masses and $n=-4$ when setting $\sigma_{0}$ to current 95% C.L.
upper limits from CMB power spectrum analyses as reported in table 1.

Figure 1:  Evolution of the mean free-electron fraction over the $\Lambda$CDM
prediction (left) and evolution of the mean baryon and DM temperatures (right)
as a function of redshift for several DM masses for $n=-4$. The $\Lambda$CDM
mean baryon temperature is shown in black, and the CMB temperature is shown in
dotted grey. For each DM mass, the coefficient of the momentum-transfer cross
section $\sigma_{0}$ can be found in table 1, consistent with $95\%$ C.L.
upper limits derived in refs. [16, 17].

At first order, the perturbations
$\delta_{T_{i}}=(T_{i}-\bar{T}_{i})/\bar{T}_{i}$ to the DM and gas
temperatures then evolve according to

$\displaystyle\dot{\delta}_{T_{\chi}}-\frac{2}{3}\dot{\delta}_{\chi}={}$
$\displaystyle\Gamma_{\chi}\left[\left(\frac{\bar{T}_{b}-\bar{T}_{\chi}}{\bar{T}_{\chi}}\right)\delta_{b}+\frac{\bar{T}_{b}}{\bar{T}_{\chi}}(\delta_{T_{b}}-\delta_{T_{\chi}})\right]+\mathcal{H}_{\chi}\left(\frac{\delta_{b}-\delta_{T_{\chi}}}{\bar{T}_{\chi}}\right)$
(2.12)
$\displaystyle\begin{split}\dot{\delta}_{T_{b}}-\frac{2}{3}\dot{\delta}_{b}={}&\Gamma_{b}\left[\left(\frac{\bar{T}_{\chi}-\bar{T}_{b}}{\bar{T}_{b}}\right)\delta_{\chi}+\frac{\bar{T}_{\chi}}{\bar{T}_{b}}(\delta_{T_{\chi}}-\delta_{T_{b}})\right]+\mathcal{H}_{b}\left(\frac{\delta_{\chi}-\delta_{T_{b}}}{\bar{T}_{b}}\right)\,\\\\[10.0pt]
&+\frac{2\mu_{b}}{m_{e}}R_{\gamma}\left[\left(\frac{\bar{T}_{\gamma}}{\bar{T}_{b}}-1\right)\left\\{\delta_{x_{e}}-\frac{\bar{x}_{e}\delta_{x_{e}}}{1+x_{\rm{He}}+\bar{x}_{e}}+\delta_{\gamma}\right\\}-\frac{\bar{T}_{\gamma}}{\bar{T}_{b}}(\delta_{T_{b}}-\delta_{T_{\gamma}})\right]\,,\end{split}$
(2.13)

where $\delta_{T_{\gamma}}$ is the photon temperature fluctuation,
$\delta_{\gamma}=4\delta_{T_{\gamma}}$ is the photon density fluctuation,
$\delta_{x_{e}}\equiv\delta x_{e}/\bar{x}_{e}$ is free electron fraction
fluctuation, and helium fraction perturbations are neglected. In both
equations the first two terms correspond to the contributions from the DM-
baryon scattering, while the contributions from Compton heating are encoded in
the second line of the baryon temperature perturbation equation. A detailed
derivation of the fluctuations sourced by the standard Compton heating term
can be found in e.g., refs. [60, 62, 67]. Eqs. (2.12) and (2.13) are the main
theoretical result of this work, and we study the effects of including the
temperature perturbations in the following sections.

At small scales and after recombination, where the supersonic relative
velocities between DM and baryons are coherent and the growth of perturbations
is affected, we solve the system using the moving-background perturbation
theory as described in eq. (2.5). Temperature (and ionization fraction, see
below) perturbations are affected by the coherent bulk velocities indirectly
through the DM and baryon perturbations, hence we need to take this into
account. Furthermore, we neglect the photon temperature perturbations in this
regime as they are negligible.

####  2.2.1 Ionization fraction fluctuations

The evolution of the ionization fraction perturbations is modified by DM-
baryon interactions indirectly through the dependence on the gas temperature
fluctuations. After recombination, the free-electron fraction obeys [60]

$\dot{x}_{e}\approx-a\alpha_{\rm{B}}n_{H}x^{2}_{e},$ (2.14)

where $\alpha_{B}$ is the case-B recombination coefficient used in recfast
[68],

$\alpha_{B}=F\frac{a(T_{b}/10^{4}\,\mathrm{K})^{b}}{1+c(T_{b}/10^{4}\,\mathrm{K})^{d}}\,\mathrm{m}^{3}\mathrm{s}^{-1}\,,$
(2.15)

with $a=4.309\times 10^{-19}$, $b=-0.6166$, $c=0.6703$, $d=0.5300$ and
$F=1.125$ is a fudge factor to reproduce the result of a multi-level atom
calculation [69, 70]. The linear ionization fraction perturbations then evolve
as

$\dot{\delta}_{x_{e}}=\frac{\dot{\bar{x}}_{e}}{\bar{x}_{e}}(\delta_{x_{e}}+\delta_{b}+\delta_{\alpha_{\rm{B}}}),$
(2.16)

where $\delta_{n_{H}}\approx\delta_{b}$ up to very small corrections, and

$\delta_{\alpha_{\rm{B}}}=\frac{\partial\ln\mathcal{\alpha}_{\rm{B}}}{\partial\ln
T_{b}}\delta_{T_{b}}\,.$ (2.17)

####  2.2.2 Large-scale enhancement of temperature and ionization
perturbations

At small scales, the relative velocities only affect the temperature and
ionization perturbations indirectly through their effect on the matter. Over
large scales, the small-scale relative velocities have no dynamical effect on
the growth of overdensities, since the non-linear terms are full divergences
that integrate to zero. However, the quadratic cooling of temperatures at
small scales is non-adiabatic and therefore the temperature and ionization
fraction need to be evolved as a coupled system up to second order, which
affects the large-scale temperature and ionization perturbations differently
[71].

Perturbing the temperature evolution equations (2.10) to second order – e.g.,
$T_{i}=\bar{T}_{i}(1+\delta_{T_{i}}^{\rm I}+\delta_{T_{i}}^{\rm II})$ – we
find

$\displaystyle\dot{\delta}^{\rm{II}}_{T_{\chi}}={}$
$\displaystyle\frac{2}{3}\dot{\delta}_{\chi}(\delta^{\rm{I}}_{T_{\chi}}-\delta_{\chi})+\Gamma_{\chi}\left[\frac{\bar{T}_{b}}{\bar{T}_{\chi}}\left(\delta^{\rm{II}}_{T_{b}}-\delta^{\rm{II}}_{T_{\chi}}+\delta^{\rm{I}}_{T_{b}}\delta_{b}\right)-\delta^{\rm{I}}_{T_{\chi}}\delta_{b}\right]+\mathcal{H}_{\chi}\frac{\delta^{\rm{II}}_{T_{\chi}}}{\bar{T}_{\chi}}$
(2.18)
$\displaystyle\begin{split}\dot{\delta}^{\rm{II}}_{T_{b}}={}&\frac{2}{3}\dot{\delta}_{b}(\delta^{\rm{I}}_{T_{b}}-\delta_{b})+\Gamma_{b}\left[\frac{\bar{T}_{\chi}}{\bar{T}_{b}}\left(\delta^{\rm{II}}_{T_{\chi}}-\delta^{\rm{II}}_{T_{b}}+\delta^{\rm{I}}_{T_{\chi}}\delta_{\chi}\right)-\delta^{\rm{I}}_{T_{b}}\delta_{\chi}\right]+\mathcal{H}_{b}\frac{\delta^{\rm{II}}_{T_{b}}}{\bar{T}_{b}}\\\\[10.0pt]
&\qquad\qquad\quad\ \
+\frac{2\mu_{b}}{m_{e}}R_{\gamma}\left[\frac{\bar{T}_{\gamma}-\bar{T}_{b}}{\bar{T}_{b}}\delta^{\rm{II}}_{x_{e}}-\delta^{\rm{I}}_{x_{e}}\delta^{\rm{I}}_{T_{b}}-\frac{\bar{T}_{\gamma}}{\bar{T}_{b}}\delta^{\rm{II}}_{T_{b}}\right],\end{split}$
(2.19)

where we have neglected photon temperature fluctuations and assumed
$1+x_{\rm{He}}+x_{e}\approx 1+x_{\rm{He}}$ to be locally homogeneous, valid at
small (sub-horizon) scales. The second-order perturbations to the ionization
fraction $\delta^{\rm{II}}_{x_{e}}$ evolve according to [71]

$\displaystyle\begin{split}\dot{\delta}^{\rm{II}}_{x_{e}}=\frac{\dot{x}_{e}}{\bar{x}_{e}}\left(\delta^{\rm{II}}_{x_{e}}+\frac{d\log\mathcal{A}_{B}}{d\log
T_{b}}\delta^{\rm{II}}_{T_{b}}+\delta^{\rm{II}}_{\dot{x}_{e}}\right),\end{split}$
(2.20)

where $\delta^{\rm{II}}_{\dot{x}_{e}}$ represents the part of
$\delta\dot{x}_{e}/\dot{\bar{x}}_{e}$ quadratic in the perturbations. Each of
the source terms in the above equations are previously smoothed over the small
volumes with coherent relative velocities.

##  3 Numerical results

We solve the perturbation equations in different regimes, accounting for all
new terms in the temperature perturbations sourced by DM-baryon scattering,
and also for the coherent relative velocities at small scales and the
consequent quadratic contributions at large scales, as required in each
regime. Therefore, we use a different approach in each case.

We have implemented the new contributions to the perturbed recombination
arising from DM-baryon interactions derived in section 2.2 (see eqs. (2.12),
(2.13) and (2.16)) into a modified version of the linear Boltzmann solver
class that has been developed for DM-baryon scattering models in ref. [16]. We
use this version of class to evolve the linear perturbation equations at all
scales before recombination and only large scales ($k\lesssim 10$ Mpc-1) after
that. We separately solve the evolution equations at small scales post-
recombination as function of $k$ and $\bm{v}_{b\chi}\cdot\bm{k}$ (the local
relative velocity in each patch) and average the results using the probability
distribution function of $\bm{v}_{b\chi}\cdot\bm{k}$. For this we take the
initial conditions from the modified version of class at $z=1010$ and evolve
them until $z=20$.555At $z\lesssim 40\text{--}30$ we expect the first stars to
form, triggering cosmic dawn, which will then dominate the temperature
perturbations; the effects of DM-baryon scattering on the temperature
perturbations during cosmic dawn are beyond the scope of this study and left
for future work. With these results, we follow ref. [71] to compute the
quadratic corrections at large scales, using the smoothed perturbations at
small scales and the solution of the system in eqs. (2.18), (2.19), and
(2.20). For $z\lesssim 50$, non-linear effects become important at the several
percent level at high-$k$. Due to this, we restrict our $k_{\rm max}=500\,\rm
Mpc^{-1}$ when computing the large-scale modulation of small-scale
perturbations at these redshifts, which should result in a conservative
estimate.

For all results presented in this work we consider DM-baryon interactions with
$n=-4$ velocity dependence and adopt the best-fit Planck 2018 [1] values for
the $\Lambda$CDM cosmological parameters. We show results for selected masses
in the range 1 MeV–10 GeV, setting $\sigma_{0}$ for each DM mass to the
corresponding 95% C.L. upper limit found in ref. [16] where available. The
mass-scaling relationship $\sigma_{0}\propto(m_{\chi}+m_{H})$ of refs. [17,
18] is used to obtain our benchmark 10 GeV limit. For reference these values
are reported in table 1. As anticipated in the introduction, the perturbed
recombination – even with the new contributions derived in the previous
section – has no significant effect on the CMB power spectra.

DM mass $m_{\chi}$ | $\sigma_{0}\,[\rm{cm}^{2}]$
---|---
1 MeV | $1.7\times 10^{-41}$
0.1 GeV | $1.9\times 10^{-41}$
1 GeV | $3.5\times 10^{-41}$
10 GeV | $2.0\times 10^{-40}$
Table 1: Benchmark values of the momentum-transfer cross section coefficient
$\sigma_{0}$ for each DM mass used in all of our numerical results. Consistent
with the 95% C.L. upper limits found in refs. [16, 17, 18].

###  3.1 Temperature and ionization fraction perturbations

We first focus on the perturbations of the baryon temperature and ionization
fraction. We compare the results considering DM-baryon scattering with respect
to that of $\Lambda$CDM, both when including and neglecting the DM-baryon
interactions contribution to the perturbed recombination (in the latter case,
this amounts to setting $\Gamma_{b}=\mathcal{H}_{b}=0$ in eq. (2.13) and
ignoring $\delta_{T_{\chi}}$). We refer to these two cases as $\delta^{\rm
iDM}$ and $\delta^{\rm std}$ and represent them in figures with solid (or
dark) and dashed (or light), respectively, unless otherwise stated.

Figure 2:  The evolution of the gas temperature and ionization fraction
fluctuations at $k=1\,\mathrm{Mpc}^{-1}$ as a fraction of the respective
fluctuation in $\Lambda$CDM, including (dark) and neglecting (light) the new
temperature perturbations sourced by DM-baryon interactions. In the bottom
panel of each figure, we show the relative contribution of each of the three
terms in eq. (2.13) (labelled by its coefficient) to the temperature
fluctuations, scaled to the range $[-1,1]$. When neglecting the fluctuations
from DM-baryon scattering, there is only a contribution from the term
proportional to $R_{\gamma}$. Note the change of scale in the $y$-axis of each
panel.

The evolution of the baryon temperature and ionization fraction perturbations
as a function of redshift for $k=1$ Mpc-1 is shown in figure 2. In the bottom
panels, we explicitly show the relative contributions of each term in eq.
(2.13) to the evolution of the perturbations, with each term labelled by its
rate coefficient respectively: $R_{\gamma}$, $\Gamma_{b}$, or
$\mathcal{H}_{b}$. For $m_{\chi}\gtrsim 1$ GeV, the drag heating term
($\mathcal{H}_{b}$) dominates at early times, increasing temperature
fluctuations and therefore indirectly enhancing the standard Compton heating
term through its dependence on $\delta_{T_{b}}$, analogous to the effect on
the mean temperatures. These two terms have opposite impact in the evolution
of $\delta_{T_{b}}$, but their balance results in enhanced temperature
perturbations for all redshifts considered. On the other hand, the ionization
fraction fluctuations are suppressed relative to $\Lambda$CDM, as
recombination is less efficient in the overdensities due to the higher gas
temperature, i.e. there are more free electrons. A similar effect – albeit
with smaller amplitude – is seen for $m_{\chi}=1$ GeV, with the exception that
the contribution from the cooling term ($\Gamma_{b}$) becomes important at
low-$z$, and $\delta_{T_{b}}$ is suppressed for $z\lesssim 60$. For
$m_{\chi}<1$ GeV, while the $\mathcal{H}_{b}$ contribution can enhance
temperature perturbations at $z\sim 10^{3}$, the cooling term dominates at
late times; this results first in suppressed temperature fluctuations for
$70\lesssim z\lesssim 200$, and later enhanced negative fluctuations
$z\lesssim 70$ compared to $\Lambda$CDM. For $m_{\chi}=0.1$ GeV, the DM and
baryon temperatures reach equilibrium at $z\sim 50$, and so the dominant
contribution to the cooling term becomes negligible. The remaining terms add
up to increase the temperature perturbations, and hence the gas temperature
grows again in over-dense regions. For $m_{\chi}=1$ MeV, the DM and baryon
temperatures reach equilibrium later, at $z\lesssim 30$. Finally, for these
two cases, the lower temperature in overdensities results in an enhancement of
the (negative) ionization fraction fluctuations, due to more efficient
recombination. In all cases, the interdependence of the three contributions to
the temperature perturbations results in a richer redshift evolution of the
perturbed temperatures with respect to neglecting the terms sourced by DM-
baryon scattering.

Figure 3:  Ratio of the linear power spectra (with respect to $\Lambda$CDM) of
gas temperature fluctuations at $z=30$ for several DM masses, including
(solid) and neglecting (dashed) fluctuations sourced by DM-baryon
interactions, shown until $k=10\,\rm Mpc^{-1}$. At smaller scales, other
effects become important, which are considered in section 3.2.

In figure 3, we demonstrate how the baryon temperature power spectrum changes
due to the DM-baryon interactions and in particular, the importance of the new
terms derived in the previous section. We plot the ratio of the baryon
temperature power spectrum for different DM masses with respect to the
$\Lambda$CDM prediction at $z=30$. The contributions from the DM-baryon
scattering to the baryon temperature power spectrum modify the amplitude by a
factor of a few up to several tens (even hundreds for $m_{\chi}=1$ MeV at
$z=50$, see figure A.1), in a roughly scale-independent manner. Except in the
case of 1 GeV, for which temperature fluctuations become suppressed relative
to $\Lambda$CDM for $30\lesssim z\lesssim 60$, the power spectrum is always
larger than the $\Lambda$CDM prediction at the end of the dark ages ($z\sim
30$–$50$). Moreover, we see that including the DM-baryon scattering
contributions to the perturbed recombination can enhance or suppress the
baryon acoustic oscillation (BAO) feature in the temperature power spectrum.
This is a result of the changes, with respect to considering Compton
scattering only, in the terms that dominate the baryon temperature
perturbations and their dependence on the density perturbations.

These results depend on the strength of the interacting cross section, which
we explore for a DM mass of 1 GeV in appendix B. Interestingly, the change in
amplitude of the gas temperature power spectrum does not necessarily decrease
with the scattering cross section; reducing the cross section alters the
balance of competing terms on the right-hand side of eq. (2.13), which can
result in more or less enhancement or suppression of the temperature
perturbations depending on the redshift.

###  3.2 Small-scale effects

Figure 4:  Ratio (with respect to $\Lambda$CDM) of the small-scale power
spectra of baryon temperature (top), baryon density (middle) and matter
density (bottom) fluctuations at $z=30$ for scales $k=40$–$500\,\rm Mpc^{-1}$.
We compare the result when including (solid) and neglecting (dashed) the novel
contributions to perturbed recombination from DM-baryon interactions for each
DM mass, finding additional suppression of small-scale baryon and matter
density perturbations in the former case.

In addition to the almost scale-independent enhancement of the temperature
fluctuations, there are several effects that alter the power spectra only at
scales $k\gtrsim 10$ Mpc-1. First, there is the averaged effect due to the
relative bulk velocities; this effect is qualitatively similar to $\Lambda$CDM
cosmologies, featuring a suppression of fluctuations around the baryonic Jeans
scale. However, DM-baryon interactions lead to a further suppression of power
and a series of damped oscillations at very small scales. We find that
averaging over relative bulk velocities can modify the shape of this
oscillatory feature.

Second, the changes to the temperature fluctuations sourced by DM-baryon
interactions are imprinted on the small-scale density fluctuations through the
sound speed terms. The dark matter fluid sound speed term now has an
additional contribution due to the presence of $\delta_{T_{\chi}}$, which acts
to suppress the DM density perturbations at small scales for all cases since
$\delta_{T_{\chi}}$ is always positive. Similarly, $\delta_{T_{b}}$ has new
contributions, which modify the baryon fluid sound speed term: the cooler
(hotter) baryon temperature in overdensities for sub-GeV masses
($m_{\chi}\gtrsim 1$ GeV) tends to enhance (suppress) the baryon density power
spectrum. However, due to the coupling of the two fluids, the suppressing
effect on the DM due to $\delta_{T_{\chi}}$ dominates in all cases; we find
that the baryon power spectrum is increasingly suppressed due to the
contributions to perturbed recombination from DM-baryon scattering for all
masses considered. The effect is smaller for sub-GeV masses because of the
opposite trends in the effects of $\delta_{T_{b}}$ and $\delta_{T_{\chi}}$
described above. Of course, the new contributions to the sound speed terms
result in a larger suppression of the small-scale total matter power spectrum,
too. Consequently, the temperature power spectrum is also affected at small
scales through its dependence on the density perturbations.

We illustrate this effect in figure 4, showing the ratio of the small-scale
baryon temperature, baryon density and total matter density power spectra with
respect to $\Lambda$CDM at $z=30$, comparing the predictions for DM-baryon
scattering cosmologies with (solid) and without (dashed) including the new
contributions to perturbed recombination. For $k\lesssim 100$ Mpc-1 we find a
similar scale-independent enhancement of the baryon temperature power spectrum
as demonstrated in figure 3. However, this power spectrum drops at higher $k$
due to the DM-baryon interactions. This suppression is faster when all the
terms in the temperature perturbations are included, since the contributions
to the sound speed terms increase. This can be seen comparing the panels
corresponding to the baryon temperature and the baryon density power spectra.
The relative suppression of baryon density perturbations, including all
contributions to perturbed recombination, is greater for DM masses $\gtrsim 1$
GeV due to the combined effect of increased contributions to both the DM and
baryon sound speed terms (e.g., for a 10 GeV DM, there is a further $\sim
10\%$ suppression of the baryon power spectrum at $k\sim 200$ Mpc-1, compared
to a little over $\sim 5\%$ for sub-GeV masses at this scale.)

On the other hand, the DM temperature perturbations $\delta T_{\chi}$ are
larger for lighter DM masses; this results in a greater suppression of the DM
density perturbations, and therefore the total matter power spectrum, at small
scales. For sub-GeV masses, we find that the total matter power spectrum is
suppressed by $\sim 10$–$20\%$ over scales $k\sim 200$–$300$ Mpc-1 relative to
the case when neglecting the temperature perturbations sourced by DM-baryon
scattering. For heavier DM, the small-scale suppression is less steep; over
the same scales, the matter power suppression ranges from $\sim 5$–$15\%$ for
$m_{\chi}=1$ GeV and from $\sim\rm{few}$–$10\%$ for $m_{\chi}=10$ GeV,
compared to neglecting the new contributions.

Overall, depending on the DM mass, we find that including the terms sourced by
DM-baryon scattering in the perturbed recombination leads to an additional
$\sim 5$–$10\%$ suppression of the total matter power spectrum relative to
$\Lambda$CDM for scales $k\gtrsim 200$ Mpc-1. These results will also depend
on the strength of the interaction cross section $\sigma_{0}$. We discuss the
effect varying $\sigma_{0}$ for a 1 GeV mass in appendix B. As one might
expect, the additional small-scale suppression that results from our full
treatment of the perturbed recombination scales inversely with the scattering
cross section, and the effect is less pronounced as the interaction strength
weakens. However, there is still an extra suppression at the few percent level
even for a value of $\sigma_{0}$ one order of magnitude below current limits.

###  3.3 Large-scale enhancement of baryon temperature power spectrum

Finally, we study the large-scale quadratic corrections to the baryon
temperature power spectrum. In the case of $\Lambda$CDM, this correction
amounts to a $\sim 10\%$ enhancement of power at $k\lesssim
0.01\,\rm{Mpc}^{-1}$ for $z=30$ [71]. Scattering between DM and baryons can
enhance the quadratic contribution to the large-scale temperature fluctuations
(both in absolute terms and relative to the linear power spectrum), as we
illustrate in figure 5. For example, for a 1 GeV DM mass, the large-scale
temperature fluctuations are enhanced by $\sim 40\%$ at $z=30$ for $k\sim
0.01\,\rm{Mpc}^{-1}$ and increase at larger scales (becoming comparable to the
linear contribution at $k\sim 0.001\,\rm{Mpc}^{-1}$). Neglecting the
temperature fluctuations sourced by DM-baryon interactions not only
significantly mis-estimates the linear power spectrum, but the second order
contribution at large scales as well.

Figure 5:  Quadratic correction to large-scale gas temperature fluctuations at
$z=30$, including the DM-baryon scattering contribution to perturbed
recombination (iDM) and neglecting them (std.). We show the $\Lambda$CDM
prediction in black for reference. Note that for $m_{\chi}\leq 1$ GeV, the
temperature fluctuations $\delta_{T_{b}}$ can change sign from negative to
positive at certain values of $z$; this explains the “dip” feature in the
monopole evident for the 1 MeV case. For each DM mass we set the momentum-
transfer cross section to the values found in table 1.

##  4 Consequences of temperature perturbations in the presence of DM-baryon
scattering

In the previous sections, we have derived the new contributions to perturbed
recombination arising from DM-baryon scattering and their direct effects on
the evolution of cosmological perturbations, with special attention to the
temperature and density power spectra. While the newly-derived contributions
do not have an impact on CMB power spectra, the modified baryon temperature
and growth of structure at small scales may have sizable signatures in other
observables that are sensitive to the abundance of light collapsed structures
and the baryonic gas.

It is expected that the first stars form at $z\sim 30$–$40$ in the densest
volumes of the Universe where gas collapses into dark matter halos and cools.
As discussed above, accounting for the temperature perturbations arising from
the DM-baryon scattering results in a $\sim 10\%$ suppression of the baryon
and matter density power spectra with respect to the (already suppressed)
previous predictions. This suppression occurs at the scales of interest for
the formation of the dark matter halos that could host the first stars.
Therefore, the DM-baryon interactions may delay cosmic dawn even more than
originally predicted by first estimates that neglected the contributions to
the temperature perturbations, which should have a very strong impact in the
HI line-intensity mapping observables during cosmic dawn [32, 33, 34, 35, 36,
37, 38, 39, 40]. Moreover, the modified temperature perturbations may induce
additional spatial correlations in the clustering of the first stars. At the
very least, our results show that the inclusion of all relevant sources of
temperature perturbations is very likely to modify the initial conditions at
the onset of cosmic dawn. Nonetheless, it will arguably still play an
important role once sources of radiation switch on, modifying the predicted
evolution of the cosmic dawn and the epoch of reionization. Given the
complexity of such processes and the expected highly non-trivial interplay
between the effects of the DM-baryon scattering (both on densities and
temperature) and astrophysical processes, we leave a dedicated study of this
to future work. In any case, we anticipate that these contributions will
increase the sensitivity of experiments targeting both the high-redshift HI
power spectrum – such as HERA [58] and SKA [59] – and the sky-averaged
measurements – e.g. EDGES [55], SARAS [56] and LEDA [57] – to DM-baryon
scattering models.

As an example, we study the effect on the HI line-intensity mapping from the
dark ages, before the formation of the first stars, which can be computed from
the evolution equations in section 2. We first briefly review the computation
of the HI line-intensity mapping observables from the dark ages, drawing
heavily from refs. [71, 72]; we refer the interested reader to these works for
further details. Then, we show the predicted power spectra. We note that for
$z\lesssim 50$, non-linearities become important in the HI fluctuations and
can affect the small-scale power spectrum at the many percent level [72].
However, we emphasise that the aim of this work is to demonstrate the
importance of including DM-baryon interactions in the temperature
perturbations (not just the background) when evaluating the HI signal for
these models. For this purpose it is sufficient to limit ourselves to a linear
treatment, but non-linear effects should be included for a high-precision
prediction of the signal and comparison with data.

###  4.1 HI line intensity and fluctuations

The spin temperature $T_{s}$ of the baryon gas is defined through the ratio of
the populations of neutral hydrogen in triplet to singlet state, and can be
computed, in the limit $T_{\star}\ll T_{b},T_{\gamma}$, as

$T_{s}=T_{\gamma}+(T_{b}-T_{\gamma})\frac{C_{10}}{C_{10}+A_{10}\frac{T_{b}}{T_{\star}}},$
(4.1)

where $C_{01}$ and $C_{10}$ are the upward and downward collision transition
rates, respectively [72], $A_{10}\approx 2.85\times 10^{-15}\,\rm s^{-1}$ is
the spontaneous Einstein decay rate of the HI transition, and
$T_{\star}\approx 0.068\rm K$ is the energy difference between the two spin
states.

The HI brightness temperature $T_{\rm HI}$, defined as the difference between
the brightness of the HI radiation field and the background CMB radiation
field, is equivalent to the line intensity. The observed redshifted brightness
temperature today of the HI radiation emitted at redshift $z$ is given by

$T_{\rm
HI}=\frac{(T_{s}-T_{\gamma})(1-e^{-\tau})}{1+z}\approx\frac{(T_{s}-T_{\gamma})}{1+z}\tau,$
(4.2)

where we have used the optically thin limit, valid during all of the dark
ages, and $\tau$ is the optical depth for the HI radiation.

We compute the HI brightness temperature fluctuations to second order in the
hydrogen density, temperature perturbations, and ionization fraction,
neglecting fluctuations in the photon temperature. Accounting for redshift-
space distortions, the observed brightness temperature fluctuations to second
order are

$\begin{split}\delta T_{\rm HI}^{\rm obs}={}&-\overline{T}_{\rm
HI}\delta_{v}+\mathcal{T}_{H}\delta_{b}+\mathcal{T}_{T}\delta_{T_{b}}+\mathcal{T}_{x}\delta_{x_{e}}\\\
&\quad+\mathcal{T}_{HH}\delta_{b}^{2}+\mathcal{T}_{TT}\delta_{T_{b}}^{2}+\mathcal{T}_{xx}\delta_{x_{e}}^{2}+\mathcal{T}_{HT}\delta_{b}\delta_{T_{b}}+\mathcal{T}_{Hx}\delta_{x_{e}}\delta_{b}+\mathcal{T}_{Tx}\delta_{T_{b}}\delta_{x_{e}},\end{split}$
(4.3)

where $\overline{T}_{\rm HI}$ is the average brightness temperature,
$\delta_{v}\equiv\partial_{\parallel}v_{\parallel}/H$ is a dimensionless small
quantity (where $\partial_{\parallel}v_{\parallel}$ is the line-of-sight
gradient in proper space of the component of the peculiar velocity along the
line of sight), and we assume $\delta_{H}\approx\delta_{b}$. The coefficients
$\mathcal{T}_{i}$ are functions of redshift only (for explicit expressions see
ref. [73]).

The monopole source can be defined as

$\delta_{s}\equiv\frac{\mathcal{T}_{H}\delta_{b}+\mathcal{T}_{T}\delta_{T_{b}}^{\rm{I}}}{\overline{T}_{\rm
HI}}$ (4.4)

and the total contribution of quadratic terms (recalling that
$\delta_{T_{b}}=\delta^{\rm{I}}_{T_{b}}+\delta^{\rm{II}}_{T_{b}}$ contains
quadratic terms itself) is defined as

$\delta T_{\rm
HI}^{\rm{II}}\equiv\mathcal{T}_{HH}\Delta\delta_{b}^{2}+\mathcal{T}_{TT}\Delta\delta_{T_{b}}^{2}+\mathcal{T}_{HT}\Delta\delta_{b}\delta_{T_{b}}+\mathcal{T}_{T}\delta_{T_{b}}^{\rm{II}}\,,$
(4.5)

where we have neglected terms involving the ionization fraction perturbations
because they are much smaller, $\Delta\delta^{2}$ is the fluctuation around
the mean of $\delta^{2}$, and relativistic terms are neglected. Therefore, the
perturbations of the observed brightness temperature are

$\delta T_{\rm HI}^{\rm{obs}}=\overline{T}_{\rm
HI}(\delta_{s}-\delta_{v})+\delta T_{\rm HI}^{\rm{II}}=\overline{T}_{\rm
HI}(\delta_{s}-\frac{\mu^{2}\theta_{b}}{H})+\delta T_{\rm HI}^{\rm{II}}\,,$
(4.6)

where $\mu=\bm{\hat{n}}\cdot\bm{\hat{k}}$. At small scales we only consider
the first term, while at large scales the quadratic corrections are
uncorrelated with the linear terms.

Figure 6:  Ratio (with respect to $\Lambda$CDM) of the HI power spectrum as a
function of redshift at $k=1\,\rm{Mpc}^{-1}$. We show the result when
including the additional temperature fluctuations due to DM-baryon
interactions (blue, iDM) compared to neglecting them (red, std.) For each DM
mass, the value of the cross section $\sigma_{0}$ is found in table 1. Note
that for sub-GeV masses, there is effectively no HI absorption line $z\lesssim
50$, indicated by the shaded region.

In figure 6, we illustrate the effect of DM-baryon scattering on the HI power
spectrum as a function of redshift for each DM mass. We show the ratio of the
power spectra sourced by the monopole at $k=1$ Mpc-1 with respect to
$\Lambda$CDM, comparing the result when including (blue) and neglecting (red)
the new contributions to the gas temperature fluctuations from DM-baryon
scattering. The primary effect of including the additional gas temperature
fluctuations is to modify the overall amplitude of the HI power spectrum over
a wide range of scales and redshifts. This is not surprising, since this was
also the main effect at large scales for the baryon gas temperature power
spectrum which is one of the leading contributions to $\delta T_{\rm HI}$ at
linear order.

For $m_{\chi}<1$ GeV we find a large enhancement at high redshifts, peaking at
$z\sim 70$, compared to neglecting the effects of DM-baryon interactions in
the perturbed recombination. Note that for sub-GeV masses, the cooling of the
baryon gas temperature drives the mean brightness temperature
$\overline{T}_{\rm HI}\rightarrow 0$ much faster, and there is effectively no
HI absorption line below $z\sim 50$. For a 1 GeV DM mass, the amplitude of the
HI power spectrum is enhanced for $z\gtrsim 40$ due to the enhanced baryon
temperature fluctuations, except close to $z\sim 80$ where there is a near
cancellation of monopole terms
($\mathcal{T}_{H}\delta_{b}+\mathcal{T}_{T}\delta^{\rm{I}}_{T}\sim 0$); below
$z\sim 40$, the HI power spectrum is suppressed due to the suppression of the
baryon temperature fluctuations. Finally, we see that for a 10 GeV DM the HI
power spectrum is enhanced $z\gtrsim 50$; below this redshift, even though the
baryon temperature fluctuations are enhanced, this results in a closer
cancellation of the monopole terms and the amplitude of the HI power spectrum
is actually suppressed (for the benchmark $\sigma_{0}$).

In the left panels of figure 7 we show the ratio of the large-scale HI power
spectrum relative to $\Lambda$CDM at $z=30$ and 50 up to $k=10\,\rm Mpc^{-1}$.
At $z=30$, the amplitude of HI power spectrum for $m_{\chi}\geq$ 1 GeV is
suppressed by $\sim 15\%$ compared to the case when considering standard
Compton scattering only in the temperature perturbations. Moreover, we see
that the additional gas temperature fluctuations from DM-baryon interactions
can enhance or modify the BAO feature in the power spectrum depending on which
term in eq. (4.3) dominates. Observing this feature could help to distinguish
a DM-baryon scattering signal from other exotic physics which might modify the
amplitude of the HI power spectrum. The right panels of figure 7 show the
ratio of the HI power spectrum relative to $\Lambda$CDM at $z=30$ and 50 for
scales $k\geq 40\,\rm Mpc^{-1}$. The effects on smaller scales are a
combination of the results at larger scales (i.e., maintaining the amplitude
change shown in figure 6) and the suppression of temperatures and densities at
scales smaller than the Jeans scale. However, the magnitude of the monopole
term grows significantly in some cases due to the averaging over relative bulk
velocities and the non-trivial interplay of the terms in eq. (4.3), contrary
to the density and temperature perturbations.

Figure 7:  Ratio (with respect to $\Lambda$CDM) of the power spectrum of HI
fluctuations at $z=30$ (top) and $z=50$ (bottom). On the left, we show the
large-scale power spectrum both including (solid) and neglecting (dashed) the
DM-baryon scattering contributions to perturbed recombination until $k=10\,\rm
Mpc^{-1}$. On the right, we show the same for scales $k=40$–$500\,\rm
Mpc^{-1}$.

Finally, in figure 8, we show the power spectrum of the quadratic terms
$\delta T^{\rm II}_{\rm HI}$ compared to the monopole fluctuations for both
$\Lambda$CDM and DM-baryon scattering cases at $z=50$. In $\Lambda$CDM, the
relative contribution of quadratic terms is maximal around $z\sim 30$,
contributing greater than $\sim 10\%$ of the monopole fluctuations for
$k\lesssim 0.01$ Mpc-1. At this redshift, the quadratic contribution is
unimportant for sub-GeV masses as the monopole term is negligible, and for
masses 1–10 GeV the relative contribution of quadratic terms is similar to or
less than for $\Lambda$CDM. However, when there are interactions between DM
and baryons, we find considerable modifications to the absolute and relative
amplitudes of the quadratic contributions. In particular, the relative
contribution of the quadratic corrections can grow significantly at higher
redshifts. The different redshift dependence of the relative contribution of
these corrections to the HI and baryon temperature power spectra is due to the
fact that the HI signal depends non-trivially on other quantities too (see eq.
(4.5)).

For both 0.1 and 10 GeV masses, the relative contribution peaks around $z\sim
50$, becoming comparable to if not larger than the signal sourced by the
monopole for scales $k\lesssim 0.01$ Mpc-1. A similar behaviour is found for
$m_{\chi}=1$ MeV but at slightly higher redshift (not shown). In all cases,
neglecting the contributions from DM-baryon interactions in the temperature
fluctuations significantly mis-estimates the amplitude of all contributions to
the HI power spectrum.

Figure 8:  Power spectrum of HI fluctuations at $z=50$ comparing the standard
monopole term (solid) and the quadratic contribution (dashed). We compare our
full calculation (iDM), to the result when ignoring the fluctuations sourced
by DM-baryon scattering (std.), and to $\Lambda$CDM (black). The interaction
cross section $\sigma_{0}$ for each DM mass is set to the values found in
table 1.

##  5 Conclusions

We have considered the cosmological effects of dark matter that features
velocity-dependent elastic interaction with protons, where the momentum-
transfer cross section is parameterised as a power law of the relative
particle velocity, $\sigma(v)=\sigma_{0}v^{n}$. In particular, we focused on
understanding the effects of the perturbed baryon and dark matter
temperatures, deriving for the first time the contributions arising directly
from the DM-baryon collision term and studying their effects on the evolution
of cosmological perturbations. We presented numerical results for Coulomb-like
scattering with $n=-4$. Such interactions have been previously invoked to
explain the anomalous HI signal reported by the EDGES collaboration, but are
also well-motivated from a theoretical standpoint within millicharged models
for dark matter.

While the effect of DM-baryon scattering on the background DM and baryon
temperatures has been well-studied in the literature, this is the first time
the impact at the level of the perturbations has been considered. We derived
the modified perturbed recombination equations for the temperatures and the
ionization fraction and evolved them self-consistently with the linear
Boltzmann equations appropriate for DM-baryon scattering cosmologies until the
end of the dark ages ($z\sim 30$). By doing so, we also computed how the
modified DM and baryon temperature perturbations affect the growth of
structure by properly accounting for their contribution to the sound speed
terms in the momentum evolution of each fluid respectively.

We found that the baryon temperature (and ionization fraction) perturbations
are enhanced by up to 1–2 orders of magnitude when including the novel DM-
baryon scattering contributions to perturbed recombination, for scattering
cross sections within current CMB bounds. These changes are significant over a
wide redshift throughout the dark ages, and largely scale-independent for
$k\lesssim 100\,\rm Mpc^{-1}$. Moreover, the DM-baryon interactions result in
non-negligible DM temperature perturbations as well. These additional
contributions to the DM and baryon temperature perturbations also affect the
growth of density perturbations at the Jeans scale through the sound speed
terms in the Boltzmann equations; we found a further $\sim 5\%$–$10\%$
suppression of the baryon and matter power spectra at $k\sim 200\,{\rm
Mpc^{-1}}$ relative to $\Lambda$CDM at $z=30$, compared to the prediction
without including these effects (depending on the DM mass and scattering cross
section).

The effects of DM-baryon scattering on the temperature perturbations, and the
resulting suppression of small-scale structures, have important implications
for the modelling and interpretation of cosmological observables sensitive to
the gas temperature and the abundance of light collapsed objects. The
evolution of cosmic dawn and reionization is highly sensitive to the early-
time star formation history, which in turn depends upon the fraction of gas
collapsed into halos where it can sufficiently cool and form stars. Therefore,
our results will be critical to accurately predict the initial conditions of
cosmic dawn and we anticipate this to have a significant impact on the
reionization history, delaying the milestone epochs of cosmic dawn, X-ray
heating and reionization. This will provide future observations with
additional sensitivity to DM-baryon scattering models. One example of such an
observable is HI intensity mapping. To illustrate this point, we computed the
HI power spectrum from the dark ages including the newly derived contributions
to the temperature perturbations from DM-baryon scattering. We found that the
amplitude of the HI power spectrum is modified by tens of percent at $z=30$,
and up to factor of a few at higher redshifts. The effect of the interactions
is also expected to alter the global HI signal [74].

We have considered the end of the cosmic dark ages to be at $z\sim 30$.
However, if cosmic dawn begins earlier, our main conclusions still hold; the
novel effects of DM-baryon scattering on the temperature perturbations and the
matter power spectra are significant at higher redshift, which is shown
explicitly for $z=50$ in Appendix A. We emphasise that the (previously
neglected) temperature perturbations sourced by DM-baryon scattering derived
in this work are required to properly model the initial conditions for cosmic
dawn. Once the first stars form, their radiation will become the dominant
heating source over the Compton or DM-baryon scattering contributions, and a
dedicated study of the interplay of these different effects is necessary to
predict the evolution of the cosmic dawn signal. Our results strongly motivate
further investigation into how these effects would impact the cosmic dawn
signal beyond $z=30$, which are likely to further improve the constraining
power of forthcoming probes of the high-redshift Universe to DM-baryon
interactions.

## Acknowledgments

KS acknowledges support from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
713673. KS is also supported by an INPhINIT fellowship from "la Caixa"
Foundation (ID 100010434), grant code LCF/BQ/DI17/11620047. JLB is supported
by the Allan C. and Dorothy H. Davis Fellowship. KB acknowledges support from
the National Science Foundation (NSF) under Grant No. PHY-2112884. VG
acknowledges support from the NSF under Grant No. PHY-2013951. LV acknowledges
support from the European Union’s Horizon 2020 research and innovation
programme ERC (BePreSysE, grant agreement 725327). Funding for this work was
partially provided by project PGC2018-098866-B-I00 MCIN/AEI/
10.13039/501100011033 y FEDER “Una manera de hacer Europa”, and the “Center of
Excellence Maria de Maeztu 2020-2023” award to the ICCUB (CEX2019-000918-M
funded by MCIN/AEI/ 10.13039/501100011033).

##  Appendix A Results at $z=50$

In the main text we show results at $z=30$ which we consider to mark the end
of the dark ages and the beginning of cosmic dawn. However, as demonstrated in
figure 2, there are significant changes to the temperature and ionization
fraction perturbations during much of the redshift range of the dark ages. We
highlight here that there are important effects beyond $z>30$ which, given the
astrophysical uncertainties surrounding early stellar formation, may be
relevant for modelling the initial conditions of cosmic dawn at higher
redshift. Figures A.1 and A.2 are analogous to figures 3 and 4, but shown for
$z=50$.

Figure A.1:  Ratio of the linear power spectra (with respect to $\Lambda$CDM)
of the gas temperature fluctuations at $z=50$ for several DM masses, including
(solid) and neglecting (dashed) fluctuations sourced by DM-baryon
interactions, shown until $k=10\,\rm Mpc^{-1}$. Figure A.2:  Ratio (with
respect to $\Lambda$CDM) of small-scale power spectra of baryon temperature
(top), baryon density (middle) and matter density (bottom) fluctuations at
$z=50$ for scales $k=40$–$500\,\rm Mpc^{-1}$. We compare the result when
including (solid) and neglecting (dashed) additional contributions to
perturbed recombination from DM-baryon interactions for each DM mass.

##  Appendix B Varying scattering cross section and the fraction of
interacting DM

Figure B.1:  Effect of varying the scattering cross section for a DM mass of 1
GeV on the large-scale baryon temperature power spectrum (left) and the small-
scale total matter power spectrum (right) at $z=30$. We compare the result
when including (solid) and neglecting (dashed) the additional contributions to
perturbed recombination resulting from the DM-baryon interactions for several
values of $\sigma_{0}$ below the 95% C.L. upper limit given in table 1.

Throughout the main text we show results for the case where all of the dark
matter interacts with baryons and adopt a benchmark value of $\sigma_{0}$ for
each DM mass, consistent with the 95% C.L. upper limits derived in refs. [16,
17, 18]. In this appendix we vary these parameters, considering different
values of $\sigma_{0}$ and small fractions of strongly interacting DM.

In figure B.1 we show how results change when reducing the strength of the
interaction cross section by a numerical factor below the current upper limit
for $m_{\chi}=1$ GeV. We find that as we reduce $\sigma_{0}$, there is an
intermediate range where the change in the amplitude of the temperature power
spectrum at $z=30$, with respect to neglecting temperature perturbations
sourced by DM-baryon scattering, actually increases significantly before
decreasing again. This is due to how the evolution of the temperature
perturbations $\delta_{T_{b}}$ depends on the different contributions in eq.
(2.13); varying the cross section alters the balance of these competing terms
and can result in either more or less enhancement or suppression of the
temperature perturbations, depending on the redshift. For example, for an
interaction cross section set at 80% of our benchmark $\sigma_{0}$, we find
that the temperature power spectrum is suppressed by almost an order of
magnitude at $z=30$. This also has an impact in the large-scale HI power
spectrum; the amplitude of the HI signal becomes suppressed by $\sim 25\%$ at
$z=30$ relative to neglecting the DM-baryon scattering effects on the
temperature perturbations (compared to $\sim 15\%$ for our benchmark
$\sigma_{0}$).

On the other hand, the amount of extra suppression in the total matter power
spectrum due to the DM-baryon interaction contribution to the perturbations
scales inversely with $\sigma_{0}$. For a cross section one order of magnitude
weaker than current limits, there is an additional suppression of a $\sim$few
% in the total matter power spectrum at $k=200\,\rm Mpc^{-1}$.

Lastly, we also consider how our results change when only a fraction
$f_{\chi}$ of the total DM density is interacting; in this regime larger
values of the cross section are allowed for small interacting fractions. In
figure B.2 we show results for various fractions with $\sigma_{0}$ set to the
corresponding 95% C.L. upper limit derived in ref. [16] for a mass of 1 MeV
(reported in table 2).

For $f_{\chi}\gtrsim 2\%$ the limit on $\sigma_{0}$ roughly scales with
$f_{\chi}$. In this regime, we find that the amount of extra suppression in
the matter power spectrum declines with smaller $f_{\chi}$. On the other hand,
the amplification of the temperature perturbations can be substantially
increased; e.g., for $f_{\chi}=0.1$, the large-scale temperature power
spectrum is enhanced by almost an order of magnitude at $z=30$ when fully
accounting for all contributions to perturbed recombination. The difference in
scale of the impact on the baryon temperature and the matter power spectra can
be explained as follows. As shown in figure 2, each contribution to the gas
temperature depends indirectly on the rest and the evolution with redshift is
much richer than when only Compton heating is considered. Therefore, even if
only a small part of the DM interacts, the change can be very large. However,
the effect on the total matter density decreases significantly because the
temperature of the (interacting) dark matter now only affects a small fraction
of the total dark matter.

The sensitivity of CMB constraints to DM-baryon interactions degrades
significantly for very small interacting fractions (below $f_{\chi}\sim 0.4\%$
[16]), and large values of $\sigma_{0}$ are permitted for which the
interacting sub-component of DM is tightly coupled to the baryons. In this
regime, we find there can still be considerable modifications to the baryon
temperature perturbations (for $f_{\chi}=0.3\%$, suppressed by almost a factor
of $\sim 100$ at $z=30$), while the BAO feature is markedly enhanced.
Moreover, for such small tightly coupled fractions, the additional suppression
in the small-scale matter power spectrum can increase again.

These results motivate further investigation and modelling of the effects of
these contributions to the full cosmic dawn signal, especially in the context
of millicharged DM models invoked to explain the EDGES measurement.

Figure B.2:  Effect of varying the interacting fraction of DM $f_{\chi}$ for a DM mass of 1 MeV, with $\sigma_{0}$ set to the respective 95% C.L. upper limits derived in ref. [16] (see table 2). We show the large-scale baryon temperature power spectrum (left) and the small-scale total matter power spectrum (right) at $z=30$, comparing the result when including (solid) and neglecting (dashed) the additional contributions to perturbed recombination resulting from the DM-baryon interactions. $f_{\chi}$ | $\sigma_{0}\,[\rm{cm}^{2}]$
---|---
1 | $1.7\times 10^{-41}$
0.1 | $1.9\times 10^{-40}$
0.05 | $6.0\times 10^{-40}$
0.01 | $5.5\times 10^{-39}$
0.003 | $5.0\times 10^{-35}$
Table 2: Benchmark values of the momentum-transfer cross section coefficient
$\sigma_{0}$ for various interacting DM fractions $f_{\chi}$ for $m_{\chi}=1$
MeV, consistent with 95% C.L. upper limits derived in ref. [16].

## References

  * [1] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [2] eBOSS Collaboration, S. Alam et al., “Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory,” Phys. Rev. D 103 no. 8, (2021) 083533, arXiv:2007.08991 [astro-ph.CO].
  * [3] DES Collaboration, T. M. C. Abbott et al., “Dark Energy Survey Year 3 Results: Cosmological Constraints from Galaxy Clustering and Weak Lensing,” arXiv:2105.13549 [astro-ph.CO].
  * [4] C. Heymans et al., “KiDS-1000 Cosmology: Multi-probe weak gravitational lensing and spectroscopic galaxy clustering constraints,” Astron. Astrophys. 646 (2021) A140, arXiv:2007.15632 [astro-ph.CO].
  * [5] D. Brout et al., “The Pantheon+ Analysis: Cosmological Constraints,” arXiv:2202.04077 [astro-ph.CO].
  * [6] J. Billard et al., “Direct Detection of Dark Matter – APPEC Committee Report,” arXiv:2104.07634 [hep-ex].
  * [7] A. Boveia and C. Doglioni, “Dark Matter Searches at Colliders,” Ann. Rev. Nucl. Part. Sci. 68 (2018) 429–459, arXiv:1810.12238 [hep-ex].
  * [8] K. K. Boddy et al., “Astrophysical and Cosmological Probes of Dark Matter,” in 2022 Snowmass Summer Study.  3, 2022.  arXiv:2203.06380 [hep-ph].
  * [9] G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: evidence, candidates and constraints,” Physics Reports 405 no. 5-6, (Jan, 2005) 279–390. http://dx.doi.org/10.1016/j.physrep.2004.08.031.
  * [10] L. Roszkowski, E. M. Sessolo, and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81 no. 6, (2018) 066201, arXiv:1707.06277 [hep-ph].
  * [11] X.-l. Chen, S. Hannestad, and R. J. Scherrer, “Cosmic microwave background and large scale structure limits on the interaction between dark matter and baryons,” Phys. Rev. D 65 (2002) 123515, arXiv:astro-ph/0202496.
  * [12] K. Sigurdson, M. Doran, A. Kurylov, R. R. Caldwell, and M. Kamionkowski, “Dark-matter electric and magnetic dipole moments,” Phys. Rev. D 70 (2004) 083501, arXiv:astro-ph/0406355. [Erratum: Phys.Rev.D 73, 089903 (2006)].
  * [13] C. Dvorkin, K. Blum, and M. Kamionkowski, “Constraining Dark Matter-Baryon Scattering with Linear Cosmology,” Phys. Rev. D 89 no. 2, (2014) 023519, arXiv:1311.2937 [astro-ph.CO].
  * [14] V. Gluscevic and K. K. Boddy, “Constraints on Scattering of keV–TeV Dark Matter with Protons in the Early Universe,” Phys. Rev. Lett. 121 no. 8, (2018) 081301, arXiv:1712.07133 [astro-ph.CO].
  * [15] K. K. Boddy and V. Gluscevic, “First Cosmological Constraint on the Effective Theory of Dark Matter-Proton Interactions,” Phys. Rev. D 98 no. 8, (2018) 083510, arXiv:1801.08609 [astro-ph.CO].
  * [16] K. K. Boddy, V. Gluscevic, V. Poulin, E. D. Kovetz, M. Kamionkowski, and R. Barkana, “Critical assessment of CMB limits on dark matter-baryon scattering: New treatment of the relative bulk velocity,” Phys. Rev. D98 no. 12, (2018) 123506, arXiv:1808.00001 [astro-ph.CO].
  * [17] W. L. Xu, C. Dvorkin, and A. Chael, “Probing sub-GeV Dark Matter-Baryon Scattering with Cosmological Observables,” Phys. Rev. D 97 no. 10, (2018) 103530, arXiv:1802.06788 [astro-ph.CO].
  * [18] T. R. Slatyer and C.-L. Wu, “Early-Universe constraints on dark matter-baryon scattering and their implications for a global 21 cm signal,” Phys. Rev. D 98 no. 2, (2018) 023013, arXiv:1803.09734 [astro-ph.CO].
  * [19] Z. Li, V. Gluscevic, K. K. Boddy, and M. S. Madhavacheril, “Disentangling Dark Physics with Cosmic Microwave Background Experiments,” Phys. Rev. D 98 no. 12, (2018) 123524, arXiv:1806.10165 [astro-ph.CO].
  * [20] E. O. Nadler, V. Gluscevic, K. K. Boddy, and R. H. Wechsler, “Constraints on Dark Matter Microphysics from the Milky Way Satellite Population,” Astrophys. J. Lett. 878 no. 2, (2019) 32, arXiv:1904.10000 [astro-ph.CO]. [Erratum: Astrophys.J.Lett. 897, L46 (2020), Erratum: Astrophys.J. 897, L46 (2020)].
  * [21] DES Collaboration, E. O. Nadler et al., “Milky Way Satellite Census. III. Constraints on Dark Matter Properties from Observations of Milky Way Satellite Galaxies,” Phys. Rev. Lett. 126 (2021) 091101, arXiv:2008.00022 [astro-ph.CO].
  * [22] D. V. Nguyen, D. Sarnaaik, K. K. Boddy, E. O. Nadler, and V. Gluscevic, “Observational constraints on dark matter scattering with electrons,” Phys. Rev. D 104 no. 10, (2021) 103521, arXiv:2107.12380 [astro-ph.CO].
  * [23] M. A. Buen-Abad, R. Essig, D. McKeen, and Y.-M. Zhong, “Cosmological constraints on dark matter interactions with ordinary matter,” Phys. Rept. 961 (2022) 1–35, arXiv:2107.12377 [astro-ph.CO].
  * [24] K. Maamari, V. Gluscevic, K. K. Boddy, E. O. Nadler, and R. H. Wechsler, “Bounds on velocity-dependent dark matter-proton scattering from Milky Way satellite abundance,” Astrophys. J. Lett. 907 no. 2, (2021) L46, arXiv:2010.02936 [astro-ph.CO].
  * [25] K. Bechtol et al., “Dark Matter Science in the Era of LSST,” arXiv:1903.04425 [astro-ph.CO].
  * [26] J. Ooba, H. Tashiro, and K. Kadota, “Cosmological constraints on the velocity-dependent baryon-dark matter coupling,” JCAP 09 (2019) 020, arXiv:1902.00826 [astro-ph.CO].
  * [27] Y. Ali-Haïmoud, J. Chluba, and M. Kamionkowski, “Constraints on Dark Matter Interactions with Standard Model Particles from Cosmic Microwave Background Spectral Distortions,” Phys. Rev. Lett. 115 no. 7, (2015) 071304, arXiv:1506.04745 [astro-ph.CO].
  * [28] Y. Ali-Haïmoud, “Testing dark matter interactions with CMB spectral distortions,” Phys. Rev. D 103 no. 4, (2021) 043541, arXiv:2101.04070 [astro-ph.CO].
  * [29] K. K. Rogers, C. Dvorkin, and H. V. Peiris, “New limits on light dark matter \- proton cross section from the cosmic large-scale structure,” arXiv:2111.10386 [astro-ph.CO].
  * [30] J. B. Bauer, D. J. E. Marsh, R. Hložek, H. Padmanabhan, and A. Laguë, “Intensity Mapping as a Probe of Axion Dark Matter,” Mon. Not. Roy. Astron. Soc. 500 no. 3, (2020) 3162–3177, arXiv:2003.09655 [astro-ph.CO].
  * [31] J. L. Bernal, A. Caputo, and M. Kamionkowski, “Strategies to Detect Dark-Matter Decays with Line-Intensity Mapping,” Phys. Rev. D 103 no. 6, (2021) 063523, arXiv:2012.00771 [astro-ph.CO].
  * [32] S. Das, R. Mondal, V. Rentala, and S. Suresh, “On dark matter - dark radiation interaction and cosmic reionization,” JCAP 08 (2018) 045, arXiv:1712.03976 [astro-ph.CO].
  * [33] M. Escudero, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz, and P. Villanueva-Domingo, “A fresh look into the interacting dark matter scenario,” JCAP 06 (2018) 007, arXiv:1803.08427 [astro-ph.CO].
  * [34] M. Safarzadeh, E. Scannapieco, and A. Babul, “A limit on the warm dark matter particle mass from the redshifted 21 cm absorption line,” Astrophys. J. Lett. 859 no. 2, (2018) L18, arXiv:1803.08039 [astro-ph.CO].
  * [35] L. Lopez-Honorez, O. Mena, and P. Villanueva-Domingo, “Dark matter microphysics and 21 cm observations,” Phys. Rev. D 99 no. 2, (2019) 023522, arXiv:1811.02716 [astro-ph.CO].
  * [36] J. B. Muñoz, C. Dvorkin, and F.-Y. Cyr-Racine, “Probing the Small-Scale Matter Power Spectrum with Large-Scale 21-cm Data,” Phys. Rev. D 101 no. 6, (2020) 063526, arXiv:1911.11144 [astro-ph.CO].
  * [37] J. B. Muñoz, S. Bohr, F.-Y. Cyr-Racine, J. Zavala, and M. Vogelsberger, “ETHOS - an effective theory of structure formation: Impact of dark acoustic oscillations on cosmic dawn,” Phys. Rev. D 103 no. 4, (2021) 043512, arXiv:2011.05333 [astro-ph.CO].
  * [38] D. Jones, S. Palatnick, R. Chen, A. Beane, and A. Lidz, “Fuzzy Dark Matter and the 21 cm Power Spectrum,” Astrophys. J. 913 no. 1, (2021) 7, arXiv:2101.07177 [astro-ph.CO].
  * [39] S. K. Giri and A. Schneider, “Imprints of fermionic and bosonic mixed dark matter on the 21-cm signal at cosmic dawn,” arXiv:2201.02210 [astro-ph.CO].
  * [40] D. Sarkar, J. Flitter, and E. D. Kovetz, “Exploring delaying and heating effects on the 21-cm signature of fuzzy dark matter,” arXiv:2201.03355 [astro-ph.CO].
  * [41] C. Evoli, A. Mesinger, and A. Ferrara, “Unveiling the nature of dark matter with high redshift 21 cm line experiments,” JCAP 11 (2014) 024, arXiv:1408.1109 [astro-ph.HE].
  * [42] A. Fialkov, R. Barkana, and A. Cohen, “Constraining Baryon–Dark Matter Scattering with the Cosmic Dawn 21-cm Signal,” Phys. Rev. Lett. 121 (2018) 011101, arXiv:1802.10577 [astro-ph.CO].
  * [43] A. Lidz and L. Hui, “Implications of a prereionization 21-cm absorption signal for fuzzy dark matter,” Phys. Rev. D 98 no. 2, (2018) 023011, arXiv:1805.01253 [astro-ph.CO].
  * [44] J. B. Muñoz and A. Loeb, “A small amount of mini-charged dark matter could cool the baryons in the early Universe,” Nature 557 no. 7707, (2018) 684, arXiv:1802.10094 [astro-ph.CO].
  * [45] J. B. Muñoz, C. Dvorkin, and A. Loeb, “21-cm Fluctuations from Charged Dark Matter,” Phys. Rev. Lett. 121 no. 12, (2018) 121301, arXiv:1804.01092 [astro-ph.CO].
  * [46] E. D. Kovetz, V. Poulin, V. Gluscevic, K. K. Boddy, R. Barkana, and M. Kamionkowski, “Tighter limits on dark matter explanations of the anomalous EDGES 21 cm signal,” Phys. Rev. D 98 no. 10, (2018) 103529, arXiv:1807.11482 [astro-ph.CO].
  * [47] O. Mena, S. Palomares-Ruiz, P. Villanueva-Domingo, and S. J. Witte, “Constraining the primordial black hole abundance with 21-cm cosmology,” Phys. Rev. D 100 no. 4, (2019) 043540, arXiv:1906.07735 [astro-ph.CO].
  * [48] K. Short, J. L. Bernal, A. Raccanelli, L. Verde, and J. Chluba, “Enlightening the dark ages with dark matter,” JCAP 07 (2020) 020, arXiv:1912.07409 [astro-ph.CO].
  * [49] J. L. Bernal, A. Raccanelli, L. Verde, and J. Silk, “Signatures of primordial black holes as seeds of supermassive black holes,” JCAP 05 (2018) 017, arXiv:1712.01311 [astro-ph.CO]. [Erratum: JCAP 01, E01 (2020)].
  * [50] H. Tashiro, K. Kadota, and J. Silk, “Effects of dark matter-baryon scattering on redshifted 21 cm signals,” Phys. Rev. D90 no. 8, (2014) 083522, arXiv:1408.2571 [astro-ph.CO].
  * [51] J. B. Muoz, E. D. Kovetz, and Y. Ali-Hamoud, “Heating of Baryons due to Scattering with Dark Matter During the Dark Ages,” Phys. Rev. D92 no. 8, (2015) 083528, arXiv:1509.00029 [astro-ph.CO].
  * [52] R. Barkana, “Possible interaction between baryons and dark-matter particles revealed by the first stars,” 1803.06698. https://arxiv.org/abs/1803.06698.
  * [53] R. Barkana, N. J. Outmezguine, D. Redigolo, and T. Volansky, “Strong constraints on light dark matter interpretation of the EDGES signal,” Phys. Rev. D98 no. 10, (2018) 103005, arXiv:1803.03091 [hep-ph].
  * [54] A. Berlin, D. Hooper, G. Krnjaic, and S. D. McDermott, “Severely Constraining Dark Matter Interpretations of the 21-cm Anomaly,” Phys. Rev. Lett. 121 no. 1, (2018) 011102, arXiv:1803.02804 [hep-ph].
  * [55] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen, and N. Mahesh, “An absorption profile centred at 78 megahertz in the sky-averaged spectrum,” Nature 555 no. 7694, (2018) 67–70, arXiv:1810.05912 [astro-ph.CO].
  * [56] S. Singh, J. Nambissan T., R. Subrahmanyan, N. Udaya Shankar, B. S. Girish, A. Raghunathan, R. Somashekar, K. S. Srivani, and M. Sathyanarayana Rao, “On the detection of a cosmic dawn signal in the radio background,” Nature Astronomy (Feb., 2022) , arXiv:2112.06778 [astro-ph.CO].
  * [57] D. C. Price, L. J. Greenhill, A. Fialkov, G. Bernardi, H. Garsden, B. R. Barsdell, J. Kocz, M. M. Anderson, S. A. Bourke, J. Craig, M. R. Dexter, J. Dowell, M. W. Eastwood, T. Eftekhari, S. W. Ellingson, G. Hallinan, J. M. Hartman, R. Kimberk, T. J. W. Lazio, S. Leiker, D. MacMahon, R. Monroe, F. Schinzel, G. B. Taylor, E. Tong, D. Werthimer, and D. P. Woody, “Design and characterization of the Large-aperture Experiment to Detect the Dark Age (LEDA) radiometer systems,” Mon. Not. Roy. Astron. Soc. 478 no. 3, (Aug., 2018) 4193–4213, arXiv:1709.09313 [astro-ph.IM].
  * [58] D. R. DeBoer et al., “Hydrogen Epoch of Reionization Array (HERA),” Publ. Astron. Soc. Pac. 129 no. 974, (2017) 045001, arXiv:1606.07473 [astro-ph.IM].
  * [59] SKA Collaboration, D. J. Bacon et al., “Cosmology with Phase 1 of the Square Kilometre Array: Red Book 2018: Technical specifications and performance forecasts,” Publ. Astron. Soc. Austral. 37 (2020) e007, arXiv:1811.02743 [astro-ph.CO].
  * [60] A. Lewis, “Linear effects of perturbed recombination,” Phys. Rev. D 76 (2007) 063001, arXiv:0707.2727 [astro-ph].
  * [61] D. Blas, J. Lesgourgues, and T. Tram, “The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes,” JCAP 07 (2011) 034, arXiv:1104.2933 [astro-ph.CO].
  * [62] S. Naoz and R. Barkana, “Growth of linear perturbations before the era of the first galaxies,” Mon. Not. Roy. Astron. Soc. 362 (2005) 1047–1053, arXiv:astro-ph/0503196.
  * [63] D. Tseliakhovich, R. Barkana, and C. Hirata, “Suppression and Spatial Variation of Early Galaxies and Minihalos,” Mon. Not. Roy. Astron. Soc. 418 (2011) 906, arXiv:1012.2574 [astro-ph.CO].
  * [64] D. Tseliakhovich and C. Hirata, “Relative velocity of dark matter and baryonic fluids and the formation of the first structures,” Phys. Rev. D 82 (2010) 083520, arXiv:1005.2416 [astro-ph.CO].
  * [65] N. Dalal, U.-L. Pen, and U. Seljak, “Large-scale BAO signatures of the smallest galaxies,” JCAP 11 (2010) 007, arXiv:1009.4704 [astro-ph.CO].
  * [66] A. Fialkov, R. Barkana, D. Tseliakhovich, and C. M. Hirata, “Impact of the Relative Motion between the Dark Matter and Baryons on the First Stars,” Mon. Not. Roy. Astron. Soc. 424 (2012) 1335–1345, arXiv:1110.2111 [astro-ph.CO].
  * [67] L. Senatore, S. Tassev, and M. Zaldarriaga, “Cosmological Perturbations at Second Order and Recombination Perturbed,” JCAP 08 (2009) 031, arXiv:0812.3652 [astro-ph].
  * [68] S. Seager, D. D. Sasselov, and D. Scott, “A new calculation of the recombination epoch,” Astrophys. J. Lett. 523 (1999) L1–L5, arXiv:astro-ph/9909275.
  * [69] J. A. Rubino-Martin, J. Chluba, W. A. Fendt, and B. D. Wandelt, “Estimating the impact of recombination uncertainties on the cosmological parameter constraints from cosmic microwave background experiments,” Mon. Not. Roy. Astron. Soc. 403 (2010) 439, arXiv:0910.4383 [astro-ph.CO].
  * [70] N. Lee and Y. Ali-Haïmoud, “HYREC-2: a highly accurate sub-millisecond recombination code,” Phys. Rev. D 102 no. 8, (2020) 083517, arXiv:2007.14114 [astro-ph.CO].
  * [71] Y. Ali-Haïmoud, P. D. Meerburg, and S. Yuan, “New light on 21 cm intensity fluctuations from the dark ages,” Phys. Rev. D 89 no. 8, (2014) 083506, arXiv:1312.4948 [astro-ph.CO].
  * [72] A. Lewis and A. Challinor, “The 21cm angular-power spectrum from the dark ages,” Phys. Rev. D 76 (2007) 083005, arXiv:astro-ph/0702600.
  * [73] A. Pillepich, C. Porciani, and S. Matarrese, “The bispectrum of redshifted 21-cm fluctuations from the dark ages,” Astrophys. J. 662 (2007) 1–14, arXiv:astro-ph/0611126 [astro-ph].
  * [74] T. Driskell et al., “Bounds on dark matter-baryon scattering from the global 21 cm signal and the suppressed matter power spectrum,” In prep. (2022) .