# Anticipating a New-Physics Signal in Upcoming 21-cm Power Spectrum
Observations

Rennan Barkana<EMAIL_ADDRESS>School of Physics and Astronomy, Tel-Aviv
University, Tel-Aviv, 69978, Israel  Institute for Advanced Study, 1 Einstein
Drive, Princeton, New Jersey 08540, USA  Department of Astronomy and
Astrophysics, University of California, Santa Cruz, CA 95064, USA  Anastasia
Fialkov<EMAIL_ADDRESS>Institute of Astronomy, University of
Cambridge, Madingley Road, Cambridge, CB3 0HA, UK  Kavli Institute for
Cosmology, Madingley Road, Cambridge CB3 0HA, UK  Hongwan Liu
<EMAIL_ADDRESS>Center for Cosmology and Particle Physics, Department
of Physics, New York University, New York, NY 10003, U.S.A.  Department of
Physics, Princeton University, Princeton, New Jersey, 08544, U.S.A.  Nadav
Joseph Outmezguine<EMAIL_ADDRESS>Berkeley Center for Theoretical Physics,
University of California, Berkeley, CA 94720, U.S.A.  Theory Group, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, U.S.A.

###### Abstract

Dark matter-baryon interactions can cool the baryonic fluid, which has been
shown to modify the cosmological 21-cm global signal. We show that in a two-
component dark sector with an interacting millicharged component, dark matter-
baryon scattering can produce a 21-cm power spectrum signal with acoustic
oscillations. The signal can be up to three orders of magnitude larger than
expected in $\Lambda$CDM cosmology, given realistic astrophysical models. This
model provides a new-physics target for near-future experiments such as HERA
or NenuFAR, which can potentially discover or strongly constrain the dark
matter explanation of the putative EDGES anomaly.

In recent years, rapid progress has been made toward turning 21-cm cosmology
into reality, opening a valuable window into the Universe at $6\lesssim
z\lesssim 30$. During this epoch, the intergalactic medium (IGM) reached its
lowest temperature in cosmological history, before heating up due to star
formation, which began in dark matter (DM) halos of sufficient mass. As a
result, 21-cm measurements probe the previously unknown period of cosmic dawn
and the first stars Madau  _et al._ (1997); Loeb and Zaldarriaga (2004);
Tseliakhovich and Hirata (2010); Furlanetto  _et al._ (2006); Barkana (2016),
and are particularly sensitive to new-physics processes that impact the
thermal state of the IGM or the process of star formation during this epoch
Sitwell  _et al._ (2014); Sekiguchi and Tashiro (2014); Nebrin  _et al._
(2019); Muñoz  _et al._ (2020); Jones  _et al._ (2021); Hotinli  _et al._
(2021); Flitter and Kovetz (2022); Evoli  _et al._ (2014); Lopez-Honorez  _et
al._ (2016); D’Amico  _et al._ (2018); Liu and Slatyer (2018); Clark  _et al._
(2018); Cheung  _et al._ (2019); Mitridate and Podo (2018); Barkana (2018);
Tashiro  _et al._ (2014); Muñoz  _et al._ (2015); Muñoz and Loeb (2018);
Barkana  _et al._ (2018); Fialkov  _et al._ (2018); Berlin  _et al._ (2018);
Muñoz  _et al._ (2018); Kovetz  _et al._ (2018); Liu  _et al._ (2019a);
Creque-Sarbinowski  _et al._ (2019); Aboubrahim  _et al._ (2021); Adshead  _et
al._ (2022).

The observable in 21-cm cosmology is the brightness temperature of radiation
with wavelength $21\text{\,}\mathrm{cm}$, $T_{21}$, absorbed or emitted by the
hyperfine states of neutral hydrogen atoms, and observed at a redshifted
wavelength. Measurements of the sky-averaged $T_{21}$ (i.e., the global 21-cm
signal) have generated significant excitement in recent years. The EDGES
collaboration reported a detection of the global signal Bowman  _et al._
(2018), finding a large absorption trough of
$-0.5^{+0.2}_{-0.5}~{}$\text{\,}\mathrm{K}$$ at $z\sim 17$ at 99% confidence.
This result is in significant tension with expectations from $\Lambda$CDM
cosmology Barkana (2018). More recently, however, the SARAS experiment found
that the central value of the EDGES absorption profile is inconsistent with
their measurements at 95% confidence Singh  _et al._ (2021). Near-future
global signal experiments such as PRIZM Philip  _et al._ (2019), SCI-HI Voytek
_et al._ (2014), REACH de Lera Acedo  _et al._ (2022), and MIST Liu  _et al._
(2019b), as well as future results from EDGES and SARAS will help to clarify
the situation soon.

Meanwhile, $T_{21}$ power spectrum measurements have been improving steadily.
Over the last decade, experiments such as GMRT Paciga  _et al._ (2013), MWA
Dillon  _et al._ (2014); Beardsley  _et al._ (2016); Li  _et al._ (2019);
Barry  _et al._ (2019); Trott  _et al._ (2020); Yoshiura  _et al._ (2021),
LOFAR Patil  _et al._ (2017); Mertens  _et al._ (2020), PAPER Kolopanis  _et
al._ (2019), LEDA Garsden  _et al._ (2021) and HERA Abdurashidova  _et al._
(2022a, b) have set increasingly strong upper limits on the power spectrum for
comoving wavenumbers $0.033\text{\,}{\mathrm{Mpc}}^{-1}$, in a broad redshift
range of $6\lesssim z\lesssim 17$; current limits are potentially one order of
magnitude away from optimistic $\Lambda$CDM expectations Reis  _et al._
(2021).

Inspired by the EDGES result, recent effort has been directed toward finding
models that enhance the global signal. This can be accomplished by either
increasing the brightness of the background radiation Feng and Holder (2018);
Bowman  _et al._ (2018); Ewall-Wice  _et al._ (2018); Pospelov  _et al._
(2018); Fialkov and Barkana (2019); Reis  _et al._ (2020), or by reducing the
baryon temperature Evoli  _et al._ (2014); Lopez-Honorez  _et al._ (2016);
D’Amico  _et al._ (2018); Liu and Slatyer (2018); Clark  _et al._ (2018);
Cheung  _et al._ (2019); Mitridate and Podo (2018), typically through DM-
baryon scattering Barkana (2018); Tashiro  _et al._ (2014); Muñoz  _et al._
(2015, 2018); Muñoz and Loeb (2018); Fialkov  _et al._ (2018); Berlin  _et
al._ (2018); Barkana  _et al._ (2018); Kovetz  _et al._ (2018); Liu  _et al._
(2019a); Creque-Sarbinowski  _et al._ (2019); Aboubrahim  _et al._ (2021);
Adshead  _et al._ (2022). In this Letter, we revisit the two-fluid dark
sector—comprising a dominant cold DM (CDM) and a subdominant millicharged DM
(mDM)—first presented in Ref. Liu  _et al._ (2019a). This model generates a
large global signal absorption trough by cooling baryons efficiently, without
introducing significant drag on the baryonic fluid in the early Universe.
Here, we show that the same model also leads to an enhanced power spectrum
that can be several orders of magnitude larger than $\Lambda$CDM expectations.
A strong correlation between baryon temperature and the baryon-DM bulk
relative velocity naturally imprints large acoustic oscillations on the signal
Barkana (2018); Muñoz  _et al._ (2015); Fialkov  _et al._ (2018), with a much
larger amplitude than is possible through the effect on galaxy formation in
standard astrophysical models Tseliakhovich and Hirata (2010); Dalal  _et al._
(2010); McQuinn and O’Leary (2012); Visbal  _et al._ (2012); Muñoz  _et al._
(2018); Muñoz (2019). An oscillation signature is possible in principle within
the simpler, non-interacting millicharged DM model Muñoz and Loeb (2018);
Muñoz  _et al._ (2018), but it is erased by drag at early times throughout the
small parameter space of this model that remains consistent with observational
constraints, particularly the constraints from the cosmic microwave background
(CMB) de Putter  _et al._ (2019); Boddy  _et al._ (2018); Kovetz  _et al._
(2018). The predicted power spectrum in the interacting millicharged DM model
provides a near-future, new-physics target for 21-cm power spectrum
experiments as their sensitivity improves.

21-cm cosmology: As background radiation photons pass through a region of
neutral hydrogen, they interact with the hydrogen hyperfine states.
Consequently, absorption, spontaneous emission and stimulated emission of
21-cm photons change the background radiation intensity at that wavelength.
The photon intensity then redshifts; the temperature contrast between the
transmitted radiation and the background radiation, observed today at a
wavelength $21(1+z)$ $\text{\,}\mathrm{cm}$, is referred to as the 21-cm
brightness temperature and denoted by $T_{21}(z)$. We assume throughout this
Letter that the background radiation temperature is the CMB temperature,
$T_{\gamma}(z)=2.725(1+z)\,$\text{\,}\mathrm{K}$$. $T_{21}$ is determined by
the spin temperature $T_{\text{S}}$ of neutral hydrogen gas, via Madau  _et
al._ (1997); Furlanetto  _et al._ (2006); Barkana (2016)

$\displaystyle
T_{21}(z)\simeq\frac{1}{1+z}\left[T_{\text{S}}(z)-T_{\gamma}(z)\right]\left[1-e^{-\tau(z)}\right]\,,$
(1)

where $\tau(z)$ is the effective optical depth of photons with wavelength 21
cm at redshift $z$.

$T_{\text{S}}$ is determined by the interaction of neutral hydrogen atoms
with: 1) the CMB photons at temperature $T_{\gamma}$; 2) other hydrogen atoms
in the gas with temperature $T_{\text{b}}$; and 3) Lyman-$\alpha$ (Ly$\alpha$)
radiation, which influences $T_{\text{S}}$ through the Wouthuysen-Field (WF)
effect Wouthuysen (1952); Field (1958). The first process drives
$T_{\text{S}}\to T_{\gamma}$, while the other two processes pull
$T_{\text{S}}\to T_{\rm b}$ instead; as a result, $T_{\text{S}}$ takes a value
between $T_{\gamma}$ and $T_{\text{b}}$. The relative strength of these
processes can be parameterized by a single coupling coefficient
$x_{\text{tot,eff}}$, encapsulating the collisional and Ly$\alpha$
interactions Barkana (2016),

$\displaystyle
T_{\text{S}}^{-1}(\bm{x},z)=\frac{T_{\gamma}^{-1}(z)+x_{\text{tot,eff}}(\bm{x},z)T_{\text{b}}^{-1}(\bm{x},z)}{1+x_{\text{tot,eff}}(\bm{x},z)}\,,$
(2)

where we have also introduced the spatial ($\bm{x}$) dependence of $T_{\rm
S}$, resulting from inhomogeneities in $T_{\rm b}$ and $x_{\text{tot,eff}}$.
Further details on the value of $x_{\text{tot,eff}}$ and other details of the
calculation of $T_{21}$ are discussed in the Supplemental Material.

There are two main types of 21-cm observables. The first is the global or sky-
averaged brightness temperature, which we denote by $\left\langle
T_{21}\right\rangle(z)$. After the first stars formed and started emitting
significant Ly$\alpha$ radiation at $z\lesssim 25$, but prior to substantial
X-ray heating of the intergalactic medium (IGM) at $z\lesssim 15$,
$T_{\text{b}}\leq T_{\text{S}}\leq T_{\gamma}$, leading to a global signal
that is in absorption, i.e., $\left\langle T_{21}\right\rangle\leq 0$. A value
of $\left\langle T_{21}\right\rangle(z=17)$ no lower than approximately
$-$150\text{\,}\mathrm{mK}$$ is expected within $\Lambda$CDM cosmology Reis
_et al._ (2021).

The second observable is the power spectrum of spatial fluctuations in the
brightness temperature, $P_{T_{21}}(k,z)$, which is a function of comoving
wavenumber $k$ and redshift $z$. The dimensionless power spectrum
$P_{T_{21}}/\left\langle T_{21}\right\rangle^{2}$ is the Fourier transform of
the two-point correlation function (2PCF) of fluctuations,
$\xi_{T_{21}}(\bm{x})=\langle\delta_{T_{21}}(\bm{x}_{1})\delta_{T_{21}}(\bm{x}_{2})\rangle_{\bm{x}}$,
where $\langle\cdots\rangle_{\bm{x}}$ denotes a spatial average over all pairs
of points $\bm{x}_{1}$ and $\bm{x}_{2}$ such that
$\bm{x}=\bm{x}_{1}-\bm{x}_{2}$, and where $\delta_{T_{21}}\equiv
T_{21}/\langle T_{21}\rangle-1$. For our isotropic Universe, $\xi_{T_{21}}$ is
a function only of $x\equiv|\bm{x}|$ and $P_{T_{21}}$ can be written as

$\frac{P_{T_{21}}(k,z)}{\left\langle
T_{21}\right\rangle^{2}(z)}=4\pi\int_{0}^{\infty}dx\,x^{2}\xi_{T_{21}}(x,z)\frac{\sin(kx)}{kx}\,.$
(3)

$P_{T_{21}}$ is often equivalently expressed in terms of the power per $\log
k$, $\Delta_{21}^{2}(k,z)\equiv k^{3}P_{T_{21}}(k,z)/(2\pi^{2})$.

Two-fluid dark sector: The two-fluid interacting millicharged DM model, first
proposed in Ref. Liu  _et al._ (2019a), is capable of cooling baryons during
the cosmic dark ages sufficiently to produce a global signal $\langle
T_{21}\rangle(z=17)$ consistent with the EDGES result, while avoiding
stringent CMB constraints on momentum transfer between baryons and dark matter
Dvorkin  _et al._ (2014); Xu  _et al._ (2018); Gluscevic and Boddy (2018);
Boddy  _et al._ (2018); Kovetz  _et al._ (2018). To accomplish this, the dark
sector comprises two components: a millicharged component (mDM) with mass
$m_{\text{m}}$ and electric charge $Q$ that makes up $f_{\text{m}}\lesssim
4\times 10^{-3}$ of the total DM energy density, and a cold component (CDM)
which accounts for the remainder. A light mediator between mDM and CDM allows
for energy transfer between these two fluids.

Prior to recombination, the mDM electric charge and mDM-CDM couplings are set
such that the mDM fluid is tightly coupled to the baryons. After
recombination, the baryonic fluid becomes mostly neutral; as a result, the mDM
fluid decouples from the baryons, and instead becomes coupled to the CDM
fluid, which has a temperature below the baryon temperature at all times. The
mDM-baryon interactions now transfer heat from the baryons to the entire CDM-
mDM fluid, cooling $T_{\text{b}}$ well below the $\Lambda$CDM expectation. The
tight mDM-CDM coupling allows for heat flowing from the baryonic fluid to be
shared among all dark-sector particles, greatly enhancing the available heat
capacity for cooling. In this Letter, we focus on
$m_{\text{m}}\gtrsim$1\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, which are
sufficiently massive to avoid CMB $N_{\text{eff}}$ limits on a combination of
millicharged particles and light dark photons Adshead  _et al._ (2022).

As in the $\Lambda$CDM paradigm, baryon acoustic oscillations set up a local
bulk relative velocity between the baryons (together with the mDM) and CDM,
$\bm{v}_{\text{bC}}(\bm{x})$, with a root-mean-square velocity of
$v_{\text{rms}}\equiv\langle
v_{\text{bC}}^{2}(\bm{x})\rangle^{1/2}\approx$29\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
at $z=1010$ Tseliakhovich and Hirata (2010). Unlike the non-interacting
millicharged DM model Muñoz and Loeb (2018), the mDM-CDM coupling in the
interacting model restores the mDM-baryon velocity difference after cosmic
recombination. Since the mDM-baryon interaction responsible for heat transfer
from the baryons weakens rapidly with velocity, patches with large initial
$v_{\text{bC}}$ remain hotter than the rapidly-cooling patches with vanishing
$v_{\text{bC}}$ Barkana (2018); Muñoz  _et al._ (2015). This is an important
point: the relative bulk motion at $z\sim 1010$ results in different initial
conditions for the baryon temperature evolution at each spatial location. This
spatial variation in $T_{\rm b}$ leads to spatial variation in the brightness
temperature $T_{21}$ through Eqs. (1) and (2). In particular, the correlation
with $v_{\text{bC}}$ imprints the acoustic oscillations in the $v_{\text{bC}}$
power spectrum onto the 21-cm power spectrum Tseliakhovich and Hirata (2010);
Barkana (2018).

Following Ref. Liu  _et al._ (2019a), we compute $T_{\text{b}}$ as a function
of dark matter parameters $Q$ and $m_{\text{m}}$, fixing the CDM mass at
$m_{\text{C}}=$100\text{\,}\mathrm{Me\kern-0.50003ptV}$$ and also
$f_{\text{m}}=10^{-4}$. We choose the maximal coupling between mDM and CDM,
permitting both tight coupling between mDM and baryons before recombination,
and a sufficiently small drag on the baryonic fluid to avoid CMB power
spectrum constraints (Boddy  _et al._ (2018); Kovetz  _et al._ (2018), also
see Appendix C of Ref. Liu  _et al._ (2019a)). We integrate the differential
equations governing the properties of the mDM, CDM and baryon fluids starting
from photon decoupling at $z=1010$, for various initial bulk velocities
$v_{\rm bC}$, ultimately obtaining
$T_{\text{b}}(Q,m_{\text{m}};v_{\text{bC}};z)$.

Astrophysics Modeling: To relate
$T_{\text{b}}(Q,m_{\text{m}};v_{\text{bC}},z)$ from our model to a value of
$T_{21}(Q,m_{\text{m}};v_{\text{bC}},z)$ through Eqs. (1) and (2), we need an
astrophysical model for the Ly$\alpha$ and X-ray radiation fields. Ly$\alpha$
photons determine the coupling of $T_{\text{S}}$ to $T_{\text{b}}$ in Eq. (2),
while both Ly$\alpha$ and X-ray photons lead to IGM heating, partially
counteracting cooling by mDM. To determine both of these effects, we rely on a
large-scale, semi-numerical 21-cm code based on Refs. Visbal  _et al._ (2012);
Fialkov  _et al._ (2014a); Reis  _et al._ (2021). In this Letter, we aim to
highlight the discovery potential of the two-fluid dark sector model in the
21-cm power spectrum; we therefore account for the minimal effect of realistic
astrophysical models. This contrasts with the typically-adopted simplistic
approach that derives the maximum possible signal by assuming no astrophysical
heating together with full Ly$\alpha$ coupling (i.e.,
$x_{\text{tot,eff}}\rightarrow\infty$); in practice, this limit is not
possible, since strong Ly$\alpha$ coupling brings along with it significant
heating, thus narrowing the range of possible 21-cm signals Reis  _et al._
(2021). We use an ensemble of 140 realistic astrophysical models from the
semi-numerical simulations, chosen to minimize astrophysical heating and
maximize Ly$\alpha$ coupling. As shown below, we must use an ensemble since
the astrophysical model that gives the maximum absorption of $T_{21}$ depends
on the DM parameters and on redshift. Further details on the simulations and
the astrophysical parameters are discussed in the Supplemental Material.

As we focus on models with subdominant X-ray heating, the dominant process
counteracting dark cooling is Ly$\alpha$ heating. Ly$\alpha$ heating results
from scattering of Ly$\alpha$ photons, either directly from atomic recoil Chen
and Miralda-Escude (2004); Chuzhoy and Shapiro (2006); Furlanetto and
Pritchard (2006), or by mediating heat transfer from the CMB to the baryons
Venumadhav  _et al._ (2018). From each simulation, we obtain the spatially
averaged baryon temperature ${T_{\text{b}}^{A}}(z)$, with $A$ indexing the 140
astrophysical models. To account for astrophysical heating in our calculations
we define the excess heating as $\delta
T_{\text{b}}^{A}(z)\equiv{T_{\text{b}}^{A}}(z)-T_{\text{b}}^{0}(z)$, where
$T_{\text{b}}^{0}(z)$ is the $\Lambda$CDM prediction for the baryon
temperature in the absence of heating. The final baryon temperature is then

$T_{\text{b}}^{A}(Q,m_{\text{m}};v_{\text{bC}};z)=T_{\text{b}}(Q,m_{\text{m}};v_{\text{bC}};z)+\delta
T_{\text{b}}^{A}(z)$ (4)

for each astrophysical model. The second ingredient we obtain from each
simulation is the average effective coupling
$x_{\text{tot,eff}}^{A}(T_{\text{b}};z)$ from the spatially averaged
Ly$\alpha$ radiation field. $T^{A}_{21}(Q,m_{\text{m}};v_{\text{bC}};z)$ is
then determined by substituting $x_{\text{tot,eff}}^{A}(T_{\text{b}};z)$ into
Eq. (2) together with Eq. (4).

Note that we have ignored spatial variations in $\delta T_{\text{b}}^{A}$, as
well as in the Ly$\alpha$ radiation field, using only their mean values, and
thus leaving the dependence of $T_{\text{b}}$ on $v_{\rm bC}$ as the only
source of fluctuations. While a fully self-consistent treatment including
astrophysical heating and dark cooling would account for all spatial
variations concurrently, the simplified prescription presented here is
computationally much more feasible, and is a reasonable approximation within
the parameter range of interest, for which the effect of velocity fluctuations
strongly dominate.

Bulk relative velocity: In the two-fluid interacting dark-sector model, the
fluctuations of $T_{\text{21}}$ are set primarily by the dependence of the
baryon temperature on $v_{\text{bC}}$. This correlates the spatial variations
of $T_{\text{21}}$ (at any redshift), to those of $v_{\rm bC}$ at $z\sim
1010$, assuming that any other sources of such fluctuations are subdominant.

Since the drag between the baryons and CDM is small, $\bm{v}_{\text{bC}}$ at
recombination is as it is in the $\Lambda$CDM paradigm: a Gaussian random
field with a power spectrum $P_{v}(k)$ defined through
$\langle\tilde{v}_{\text{bC}}(\bm{k})\tilde{v}_{\text{bC}}(\bm{k}^{\prime})\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}+\bm{k}^{\prime})P_{v}(k)$,
where
$\hat{k}\tilde{v}_{\text{bC}}(\bm{k})\equiv\tilde{\bm{v}}_{\text{bC}}(\bm{k})$,
the Fourier transform of $\bm{v}_{\text{bC}}(\bm{x})$. $P_{v}(k)$ exhibits the
characteristic acoustic oscillations in $k$ Tseliakhovich and Hirata (2010).

Under our assumption that the spatial fluctuations of $T_{21}$ are dominated
by its dependence on $v_{\text{bC}}$, the global signal is evaluated as

$\displaystyle\langle T_{21}\rangle=\left(\frac{3}{2\pi
v_{\text{rms}}^{2}}\right)^{3/2}\int
d^{3}\bm{v}_{\text{bC}}\,T_{21}(v_{\text{bC}})\exp\left(-\frac{3v_{\text{bC}}^{2}}{2v_{\text{rms}}^{2}}\right)\,.$
(5)

Similarly, the $T_{21}$ 2PCF is given by

$\xi_{T_{21}}(\bm{x})=\int d^{3}\bm{v}_{\text{bC},a}\int
d^{3}\bm{v}_{\text{bC},b}\\\
\times\mathcal{P}(\bm{v}_{\text{bC},a},\bm{v}_{\text{bC},b};\bm{x})\delta_{T_{21}}(v_{\text{bC},a})\delta_{T_{21}}(v_{\text{bC},b})\,,$
(6)

where
$\mathcal{P}\left(\bm{v}_{\text{bC},a},\bm{v}_{\text{bC},b};\bm{x}\right)$ is
the joint PDF of 3D bulk relative velocities at points $a$ and $b$ separated
by vector $\bm{x}$. Since $v_{\text{bC}}$ is a Gaussian random field,
$\mathcal{P}$ is completely specified by $P_{v}(k)$. In the Supplemental
Material, we specify the exact structure of $\mathcal{P}$, and explain the
details of our computation of $\xi_{T_{21}}$. In particular, we improve on
previous results Dalal  _et al._ (2010); Ali-Haïmoud  _et al._ (2014) by
reducing the 6D integral in Eq. (6) to a 3D integral instead of a 4D one,
making the integral much easier to evaluate numerically. The power spectrum is
then calculated through Eq. (3).

Figure 1:  Left: The minimum envelope (i.e., maximum absorption) of the
predicted global signal $\left\langle T_{21}\right\rangle$ for
$m_{\text{m}}=$3\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ across all
astrophysical models, for various values of the charge (colored lines). Right:
The maximum envelope of the power spectrum $\Delta_{21}^{2}$ at
$k=$0.13\text{\,}{\mathrm{Mpc}}^{-1}$$ for the same values of the charge
(colored lines). Existing upper limits over a range of redshifts reported by
LOFAR Patil  _et al._ (2017); Mertens  _et al._ (2020) (green), and HERA
Abdurashidova  _et al._ (2022a, b) (orange) are shown as arrows. In both
panels, the signals for all astrophysical models for $Q=2\times 10^{-2}$ are
shown in gray; the fiducial range of both quantities expected in standard
$\Lambda$CDM cosmology is shaded in purple Reis  _et al._ (2021).

Results: Figure 1 shows the predicted $\left\langle T_{21}\right\rangle$ and
$\Delta^{2}_{21}$ for comoving wavenumber
$k=$0.13\text{\,}{\mathrm{Mpc}}^{-1}$$ for the two-fluid interacting
millicharged DM model with
$m_{\text{m}}=$3\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, for a range of $Q$
values that are viable given current constraints. To indicate the most easily
observable models, we plot lines in color representing the minimum $\langle
T_{21}\rangle$ and maximum $\Delta_{21}^{2}$ envelopes, obtained by varying
over all 140 astrophysical models. To give a sense of the variability in these
models, we show in gray lines the 140 different models considered for the
fixed value $Q=2\times 10^{-2}$. The full fiducial range of values that are
possible in standard $\Lambda$CDM cosmology are shaded in purple.

The global signal shown on the left has the characteristic absorption profile,
corresponding to Ly$\alpha$ emission driving $T_{\text{S}}\to T_{\text{b}}$
before heating of the IGM causes $T_{\text{b}}$ to increase; by adding X-ray
heating it is possible to vary the shape further. As pointed out in Ref. Liu
_et al._ (2019a), the global signal from our model can attain the central
value of the EDGES absorption profile at $z=17$, producing signals as large as
$\langle T_{21}\rangle(z=21)\sim-$1250\text{\,}\mathrm{mK}$$.

The power spectrum shown on the right is predicted to reach
$\Delta_{21}^{2}(z=21)\sim${10}^{5}\text{\,}{\mathrm{mK}}^{2}$$ for $Q=5\times
10^{-3}$, several orders of magnitude larger than the conventional
$\Lambda$CDM expectation Reis  _et al._ (2021) indicated by the purple band.
For $9\leq z\leq 30$, 21-cm power spectrum experiments such as MWA Dillon  _et
al._ (2014); Yoshiura  _et al._ (2021), LOFAR Patil  _et al._ (2017); Mertens
_et al._ (2020), and HERA Abdurashidova  _et al._ (2022a, b) have already
reported upper limits; we limit ourselves to $z\geq 9$ to avoid uncertainties
due to reionization, even though stringent upper limits at $z=7.9$ have
recently been reported by HERA Abdurashidova  _et al._ (2022b).

From comparing the $Q=2\times 10^{-2}$ envelopes of both panels of Fig. 1, it
is evident that models leading to a maximal global signal absorption feature
do not correspond to the largest $\Delta_{21}^{2}$ values. In the left panel
of Fig. 2, we show a scatter plot in the space of
$\Delta_{21}^{2}(k=$0.13\text{\,}{\mathrm{Mpc}}^{-1}$)$ and $\langle
T_{21}\rangle$ at $z=17$, for
$m_{\text{m}}=$3\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ and the experimentally
allowed range of $Q$ Liu  _et al._ (2019a). Each gray line is a fixed
astrophysical model, with varying charge $Q$. The colored dots indicated the
maximal $\Delta_{21}^{2}$ and minimal $\left\langle T_{21}\right\rangle$
across all astrophysical models, for a fixed $Q$. We see that increasing $Q$
first leads to gradually increasing values of $\Delta_{21}^{2}$, before
$\langle T_{21}\rangle$ starts decreasing significantly. This behavior stems
from the $\sigma\propto v_{\text{rel}}^{-4}$ relative velocity dependence of
the Rutherford scattering cross section between mDM and baryons—cooling is
most efficient for regions where $v_{\text{bC}}=0$, while for higher
velocities, less efficient cooling or even heating can take place Muñoz  _et
al._ (2015). This leads to enhanced fluctuations in $T_{\text{b}}$ and hence
$T_{21}$. As $Q$ is increased further, these fluctuations are diminished,
since large regions of the IGM—with most values of $v_{\text{bC}}$—experience
strong and early cooling; this behavior is reflected in the right panel of
Fig. 1, where $\Delta_{21}^{2}$ begins to decrease for $Q\gtrsim 5\times
10^{-3}$.

Figure 2:  Left: Scatter plot of the global signal and the power spectrum at
$k=0.13\,\rm Mpc^{-1}$ at $z=17$, for
$m_{\text{m}}=$3\text{\,}\mathrm{Ge\kern-0.50003ptV}$$. Each gray line is a
fixed astrophysical model with varying charge $Q$, while the color of each
large dot represents a value of $Q$ and indicates the maximum fluctuation and
absorption signal for that value. The EDGES central value of $\langle
T_{21}\rangle=-$500\text{\,}\mathrm{mK}$$ with its error range is shown for
reference (brown, horizontal), although the EDGES profile is ruled out at 95%
confidence by SARAS. The projected sensitivities of NenuFAR Mertens  _et al._
(2021) (green, dotted) and HERA Muñoz  _et al._ (2018) (blue, dotted) with
$10^{3}$ hours of observation are shown. Right: Contours of
$\Delta_{T_{21}}^{2}(k,z)$, maximized over both astrophysics and particle
physics parameters.

In the same figure, we show the forecast sensitivity of NenuFAR Mertens  _et
al._ (2021) and HERA Muñoz  _et al._ (2018) after $10^{3}$ hours of
observation, for $k=$0.13\text{\,}{\mathrm{Mpc}}^{-1}$$. We note that the
Square Kilometre Array (SKA) is projected to do even better, reaching a
sensitivity of around $1~{}$mK2 Koopmans  _et al._ (2015). For many choices of
$(Q,m_{\text{m}})$, we see that 21-cm power spectrum experiments are sensitive
to these models, even though the global signal only has a value of $\langle
T_{21}\rangle\approx-$200\text{\,}\mathrm{mK}$$, only slightly larger than the
$\Lambda$CDM expectation shown in Fig. 1 Reis  _et al._ (2021). At present,
measurements by the SARAS experiment have excluded the central value of the
EDGES absorption profile at the 95% confidence level, but with the exclusion
likely depending on the signal shape and not just on the amplitude.

Finally, in the right panel of Fig. 2, we plot the maximum $\Delta_{21}^{2}$
across all astrophysical and new physics models parameters, for experimentally
relevant values of $k$ and $z$. Across a broad range of $k$ and $z$,
$\Delta_{21}^{2}\gtrsim${10}^{3}\text{\,}{\mathrm{mK}}^{2}$$, reaching values
as large as $\Delta_{21}^{2}\sim${10}^{5}\text{\,}{\mathrm{mK}}^{2}$$,
exceeding the fiducial expectation from $\Lambda$CDM cosmology by several
orders of magnitude. Furthermore, large acoustic oscillations are imprinted in
the $T_{21}$ signal through the correlation between $T_{\text{b}}$ and
$v_{\text{bC}}$. This is currently the only viable new-physics model with
large acoustic oscillations in the 21-cm power spectrum, and provides a new-
physics target for experiments like NenuFAR and HERA.

Conclusion: We study the 21-cm power spectrum of the two-fluid interacting
millicharged DM model first studied in Ref. Liu  _et al._ (2019a), finding
that the model can produce a power spectrum that will be readily detectable in
near-future runs of 21-cm power spectrum experiments such as NenuFAR and HERA,
and eventually the SKA. Large power spectra are possible even in models in
which the predicted global signal is close to the standard $\Lambda$CDM range.
This model is currently the only viable model which produces a large 21-cm
power spectrum with acoustic oscillation features, demonstrating the power of
such experiments in searching for new physics. Our results should provide a
useful new-physics benchmark for upcoming 21-cm power spectrum experimental
results.

Acknowledgments: The authors would like to thank Yacine Ali-Haïmoud, Tomer
Volansky and Omer Katz for useful discussions RB acknowledges the support of
the Israel Science Foundation (grant No. 2359/20), The Ambrose Monell
Foundation and the Institute for Advanced Study as well as the Vera Rubin
Presidential Chair in Astronomy at UCSC and the Packard Foundation. AF was
supported by the Royal Society University Research Fellowship. HL is supported
by NSF grant PHY-1915409, the DOE under Award Number DESC0007968 and the
Simons Foundation. The work of NJO was supported in part by the Zuckerman STEM
Leadership Program and by the National Science Foundation (NSF) under the
grant No. PHY-1915314. This work was performed in part at the Aspen Center for
Physics, which is supported by NSF grant PHY-1607611. This work was also
performed using the Princeton Research Computing resources at Princeton
University which is a consortium of groups including the Princeton Institute
for Computational Science and Engineering and the Princeton University Office
of Information Technology’s Research Computing department. This research made
use of the IPython Perez and Granger (2007), Jupyter Kluyver  _et al._ (2016),
matplotlib Hunter (2007), NumPy Harris  _et al._ (2020), pyfftlog Talman
(1978); Hamilton (2000); Werthmüller (2020), seaborn Waskom  _et al._ (2017),
pandas McKinney (2010), SciPy Virtanen  _et al._ (2020), and tqdm da Costa-
Luis (2019) software packages.

## References

  * Madau  _et al._ (1997) Piero Madau, Avery Meiksin, and Martin J. Rees, “21-cm tomography of the intergalactic medium at high redshift,” Astrophys. J. 475, 429 (1997), arXiv:astro-ph/9608010 .
  * Loeb and Zaldarriaga (2004) Abraham Loeb and Matias Zaldarriaga, “Measuring the small - scale power spectrum of cosmic density fluctuations through 21 cm tomography prior to the epoch of structure formation,” Phys. Rev. Lett. 92, 211301 (2004), arXiv:astro-ph/0312134 .
  * Tseliakhovich and Hirata (2010) Dmitriy Tseliakhovich and Christopher Hirata, “Relative Velocity of Dark Matter and Baryonic Fluids and the Formation of the First Structures,” Phys. Rev. D82, 083520 (2010), arXiv:1005.2416 [astro-ph.CO] .
  * Furlanetto  _et al._ (2006) Steven R. Furlanetto, S. Peng Oh, and Frank H. Briggs, “Cosmology at low frequencies: The 21 cm transition and the high-redshift Universe,” Physics Reports 433, 181–301 (2006), arXiv:astro-ph/0608032 [astro-ph] .
  * Barkana (2016) Rennan Barkana, “The Rise of the First Stars: Supersonic Streaming, Radiative Feedback, and 21-cm Cosmology,” Phys. Rept. 645, 1–59 (2016), arXiv:1605.04357 [astro-ph.CO] .
  * Sitwell  _et al._ (2014) Michael Sitwell, Andrei Mesinger, Yin-Zhe Ma, and Kris Sigurdson, “The Imprint of Warm Dark Matter on the Cosmological 21-cm Signal,” Mon. Not. Roy. Astron. Soc. 438, 2664–2671 (2014), arXiv:1310.0029 [astro-ph.CO] .
  * Sekiguchi and Tashiro (2014) Toyokazu Sekiguchi and Hiroyuki Tashiro, “Constraining warm dark matter with 21 cm line fluctuations due to minihalos,” JCAP 08, 007 (2014), arXiv:1401.5563 [astro-ph.CO] .
  * Nebrin  _et al._ (2019) Olof Nebrin, Raghunath Ghara, and Garrelt Mellema, “Fuzzy Dark Matter at Cosmic Dawn: New 21-cm Constraints,” JCAP 04, 051 (2019), arXiv:1812.09760 [astro-ph.CO] .
  * Muñoz  _et al._ (2020) Julian B. Muñoz, Cora Dvorkin, and Francis-Yan Cyr-Racine, “Probing the Small-Scale Matter Power Spectrum with Large-Scale 21-cm Data,” Phys. Rev. D 101, 063526 (2020), arXiv:1911.11144 [astro-ph.CO] .
  * Jones  _et al._ (2021) Dana Jones, Skyler Palatnick, Richard Chen, Angus Beane, and Adam Lidz, “Fuzzy Dark Matter and the 21 cm Power Spectrum,” Astrophys. J. 913, 7 (2021), arXiv:2101.07177 [astro-ph.CO] .
  * Hotinli  _et al._ (2021) Selim C. Hotinli, David J. E. Marsh, and Marc Kamionkowski, “Probing ultra-light axions with the 21-cm Signal during Cosmic Dawn,” (2021), arXiv:2112.06943 [astro-ph.CO] .
  * Flitter and Kovetz (2022) Jordan Flitter and Ely D. Kovetz, “Closing the window on fuzzy dark matter with the 21cm signal,” (2022), arXiv:2207.05083 [astro-ph.CO] .
  * Evoli  _et al._ (2014) Carmelo Evoli, Andrei Mesinger, and Andrea Ferrara, “Unveiling the nature of dark matter with high redshift 21 cm line experiments,” JCAP 11, 024 (2014), arXiv:1408.1109 [astro-ph.HE] .
  * Lopez-Honorez  _et al._ (2016) Laura Lopez-Honorez, Olga Mena, Ángeles Moliné, Sergio Palomares-Ruiz, and Aaron C. Vincent, “The 21 cm signal and the interplay between dark matter annihilations and astrophysical processes,” JCAP 08, 004 (2016), arXiv:1603.06795 [astro-ph.CO] .
  * D’Amico  _et al._ (2018) Guido D’Amico, Paolo Panci, and Alessandro Strumia, “Bounds on Dark Matter annihilations from 21 cm data,” Phys. Rev. Lett. 121, 011103 (2018), arXiv:1803.03629 [astro-ph.CO] .
  * Liu and Slatyer (2018) Hongwan Liu and Tracy R. Slatyer, “Implications of a 21-cm signal for dark matter annihilation and decay,” Phys. Rev. D 98, 023501 (2018), arXiv:1803.09739 [astro-ph.CO] .
  * Clark  _et al._ (2018) Steven Clark, Bhaskar Dutta, Yu Gao, Yin-Zhe Ma, and Louis E. Strigari, “21 cm limits on decaying dark matter and primordial black holes,” Phys. Rev. D 98, 043006 (2018), arXiv:1803.09390 [astro-ph.HE] .
  * Cheung  _et al._ (2019) Kingman Cheung, Jui-Lin Kuo, Kin-Wang Ng, and Yue-Lin Sming Tsai, “The impact of EDGES 21-cm data on dark matter interactions,” Phys. Lett. B 789, 137–144 (2019), arXiv:1803.09398 [astro-ph.CO] .
  * Mitridate and Podo (2018) Andrea Mitridate and Alessandro Podo, “Bounds on Dark Matter decay from 21 cm line,” JCAP 05, 069 (2018), arXiv:1803.11169 [hep-ph] .
  * Barkana (2018) Rennan Barkana, “Possible Interaction Between Baryons and Dark-Matter Particles Revealed by the First Stars,” Nature 555, 71–74 (2018), arXiv:1803.06698 [astro-ph.CO] .
  * Tashiro  _et al._ (2014) Hiroyuki Tashiro, Kenji Kadota, and Joseph Silk, “Effects of Dark Matter-Baryon Scattering on Redshifted 21 Cm Signals,” Phys. Rev. D 90, 083522 (2014), arXiv:1408.2571 [astro-ph.CO] .
  * Muñoz  _et al._ (2015) Julian B. Muñoz, Ely D. Kovetz, and Yacine Ali-Haïmoud, “Heating of Baryons due to Scattering with Dark Matter During the Dark Ages,” Phys. Rev. D92, 083528 (2015), arXiv:1509.00029 [astro-ph.CO] .
  * Muñoz and Loeb (2018) Julian B. Muñoz and Abraham Loeb, “A small amount of mini-charged dark matter could cool the baryons in the early Universe,” Nature 557, 684 (2018), arXiv:1802.10094 [astro-ph.CO] .
  * Barkana  _et al._ (2018) Rennan Barkana, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky, “Strong Constraints on Light Dark Matter Interpretation of the Edges Signal,” Phys. Rev. D 98, 103005 (2018), arXiv:1803.03091 [hep-ph] .
  * Fialkov  _et al._ (2018) Anastasia Fialkov, Rennan Barkana, and Aviad Cohen, “Constraining Baryon–Dark Matter Scattering with the Cosmic Dawn 21-Cm Signal,” Phys. Rev. Lett. 121, 011101 (2018), arXiv:1802.10577 [astro-ph.CO] .
  * Berlin  _et al._ (2018) Asher Berlin, Dan Hooper, Gordan Krnjaic, and Samuel D. McDermott, “Severely Constraining Dark Matter Interpretations of the 21-Cm Anomaly,” Phys. Rev. Lett. 121, 011102 (2018), arXiv:1803.02804 [hep-ph] .
  * Muñoz  _et al._ (2018) Julian B. Muñoz, Cora Dvorkin, and Abraham Loeb, “21-cm Fluctuations from Charged Dark Matter,” Phys. Rev. Lett. 121, 121301 (2018), arXiv:1804.01092 [astro-ph.CO] .
  * Kovetz  _et al._ (2018) Ely D. Kovetz, Vivian Poulin, Vera Gluscevic, Kimberly K. Boddy, Rennan Barkana, and Marc Kamionkowski, “Tighter limits on dark matter explanations of the anomalous EDGES 21 cm signal,” Phys. Rev. D98, 103529 (2018), arXiv:1807.11482 [astro-ph.CO] .
  * Liu  _et al._ (2019a) Hongwan Liu, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky, “Reviving Millicharged Dark Matter for 21-Cm Cosmology,” Phys. Rev. D100, 123011 (2019a), arXiv:1908.06986 [hep-ph] .
  * Creque-Sarbinowski  _et al._ (2019) Cyril Creque-Sarbinowski, Lingyuan Ji, Ely D. Kovetz, and Marc Kamionkowski, “Direct millicharged dark matter cannot explain the EDGES signal,” Phys. Rev. D 100, 023528 (2019), arXiv:1903.09154 [astro-ph.CO] .
  * Aboubrahim  _et al._ (2021) Amin Aboubrahim, Pran Nath, and Zhu-Yao Wang, “A Cosmologically Consistent Millicharged Dark Matter Solution to the Edges Anomaly of Possible String Theory Origin,” JHEP 12, 148 (2021), arXiv:2108.05819 [hep-ph] .
  * Adshead  _et al._ (2022) Peter Adshead, Pranjal Ralegankar, and Jessie Shelton, “Dark radiation constraints on portal interactions with hidden sectors,” (2022), arXiv:2206.13530 [hep-ph] .
  * Bowman  _et al._ (2018) Judd D. Bowman, Alan E. E. Rogers, Raul A. Monsalve, Thomas J. Mozdzen, and Nivedita Mahesh, “An Absorption Profile Centred at 78 Megahertz in the Sky-Averaged Spectrum,” Nature 555, 67–70 (2018), arXiv:1810.05912 [astro-ph.CO] .
  * Singh  _et al._ (2021) Saurabh Singh, Jishnu Nambissan T., Ravi Subrahmanyan, N. Udaya Shankar, B. S. Girish, A. Raghunathan, R. Somashekar, K. S. Srivani, and Mayuri Sathyanarayana Rao, “On the detection of a cosmic dawn signal in the radio background,” (2021), arXiv:2112.06778 [astro-ph.CO] .
  * Philip  _et al._ (2019) L. Philip, Z. Abdurashidova, H. C. Chiang, N. Ghazi, A. Gumba, H. M. Heilgendorff, J. M. Jáuregui-García, K. Malepe, C. D. Nunhokee, J. Peterson, J. L. Sievers, V. Simes, and R. Spann, “Probing Radio Intensity at High-Z from Marion: 2017 Instrument,” Journal of Astronomical Instrumentation 8, 1950004 (2019), arXiv:1806.09531 [astro-ph.IM] .
  * Voytek  _et al._ (2014) Tabitha C. Voytek, Aravind Natarajan, José Miguel Jáuregui García, Jeffrey B. Peterson, and Omar López-Cruz, “Probing the Dark Ages at $z\sim$ 20: The SCI-HI 21 cm All-sky Spectrum Experiment,” Astrophys. J. Lett. 782, L9 (2014), arXiv:1311.0014 [astro-ph.CO] .
  * de Lera Acedo  _et al._ (2022) E. de Lera Acedo, D. I. L. de Villiers, N. Razavi-Ghods, W. Handley, A. Fialkov, A. Magro, D. Anstey, H. T. J. Bevins, R. Chiello, J. Cumner, A. T. Josaitis, I. L. V. Roque, P. H. Sims, K. H. Scheutwinkel, P. Alexander, G. Bernardi, S. Carey, J. Cavillot, W. Croukamp, J. A. Ely, T. Gessey-Jones, Q. Gueuning, R. Hills, G. Kulkarni, R. Maiolino, P. D. Meerburg, S. Mittal, J. R. Pritchard, E. Puchwein, A. Saxena, E. Shen, O. Smirnov, M. Spinelli, and K. Zarb-Adami, “The REACH radiometer for detecting the 21-cm hydrogen signal from redshift z $\approx$ 7.5–28,” Nature Astronomy (2022), 10.1038/s41550-022-01709-9.
  * Liu  _et al._ (2019b) Adrian Liu, H. Cynthia Chiang, Abigail Crites, Jonathan Sievers, and Renée Hložek, “High-redshift 21cm Cosmology in Canada,”  (2019b), 10.5281/zenodo.3756080, arXiv:1910.03153 [astro-ph.CO] .
  * Paciga  _et al._ (2013) Gregory Paciga _et al._ , “A refined foreground-corrected limit on the HI power spectrum at z=8.6 from the GMRT Epoch of Reionization Experiment,” Mon. Not. Roy. Astron. Soc. 433, 639 (2013), arXiv:1301.5906 [astro-ph.CO] .
  * Dillon  _et al._ (2014) Joshua S. Dillon _et al._ , “Overcoming real-world obstacles in 21 cm power spectrum estimation: A method demonstration and results from early Murchison Widefield Array data,” Phys. Rev. D 89, 023002 (2014), arXiv:1304.4229 [astro-ph.CO] .
  * Beardsley  _et al._ (2016) A. P. Beardsley _et al._ , “First Season MWA EoR Power Spectrum Results at Redshift 7,” Astrophys. J. 833, 102 (2016), arXiv:1608.06281 [astro-ph.IM] .
  * Li  _et al._ (2019) W. Li _et al._ , “First Season MWA Phase II EoR Power Spectrum Results at Redshift 7,” Astrophys. J. 887, 141 (2019), arXiv:1911.10216 [astro-ph.CO] .
  * Barry  _et al._ (2019) N. Barry _et al._ , “Improving the Epoch of Reionization Power Spectrum Results from Murchison Widefield Array Season 1 Observations,”  (2019), 10.3847/1538-4357/ab40a8, arXiv:1909.00561 [astro-ph.IM] .
  * Trott  _et al._ (2020) Cathryn M. Trott _et al._ , “Deep multiredshift limits on Epoch of Reionization 21 cm power spectra from four seasons of Murchison Widefield Array observations,” Mon. Not. Roy. Astron. Soc. 493, 4711–4727 (2020), arXiv:2002.02575 [astro-ph.CO] .
  * Yoshiura  _et al._ (2021) S. Yoshiura _et al._ , “A new MWA limit on the 21 cm power spectrum at redshifts $\sim$13–17,” Mon. Not. Roy. Astron. Soc. 505, 4775–4790 (2021), arXiv:2105.12888 [astro-ph.CO] .
  * Patil  _et al._ (2017) A. H. Patil _et al._ , “Upper limits on the 21-cm Epoch of Reionization power spectrum from one night with LOFAR,” Astrophys. J. 838, 65 (2017), arXiv:1702.08679 [astro-ph.CO] .
  * Mertens  _et al._ (2020) F. G. Mertens _et al._ , “Improved upper limits on the 21-cm signal power spectrum of neutral hydrogen at $\bm{z\approx 9.1}$ from LOFAR,” Mon. Not. Roy. Astron. Soc. 493, 1662–1685 (2020), arXiv:2002.07196 [astro-ph.CO] .
  * Kolopanis  _et al._ (2019) Matthew Kolopanis _et al._ , “A simplified, lossless re-analysis of PAPER-64,”  (2019), 10.3847/1538-4357/ab3e3a, arXiv:1909.02085 [astro-ph.CO] .
  * Garsden  _et al._ (2021) Hugh Garsden, Lincoln Greenhill, Gianni Bernardi, Anastasia Fialkov, Daniel C. Price, Daniel Mitchell, Jayce Dowell, Marta Spinelli, and Frank K. Schinzel, “A 21-cm power spectrum at 48 MHz, using the Owens Valley Long Wavelength Array,” Mon. Not. Roy. Astron. Soc. 506, 5802–5817 (2021), arXiv:2102.09596 [astro-ph.CO] .
  * Abdurashidova  _et al._ (2022a) Zara Abdurashidova _et al._ (HERA), “First Results from HERA Phase I: Upper Limits on the Epoch of Reionization 21 cm Power Spectrum,” Astrophys. J. 925, 221 (2022a), arXiv:2108.02263 [astro-ph.CO] .
  * Abdurashidova  _et al._ (2022b) Zara Abdurashidova _et al._ (HERA), “Improved Constraints on the 21 cm EoR Power Spectrum and the X-Ray Heating of the IGM with HERA Phase I Observations,” (2022b), arXiv:2210.04912 [astro-ph.CO] .
  * Reis  _et al._ (2021) Itamar Reis, Anastasia Fialkov, and Rennan Barkana, “The subtlety of Ly-a photons: changing the expected range of the 21-cm signal,” (2021), arXiv:2101.01777 [astro-ph.CO] .
  * Feng and Holder (2018) Chang Feng and Gilbert Holder, “Enhanced global signal of neutral hydrogen due to excess radiation at cosmic dawn,” Astrophys. J. Lett. 858, L17 (2018), arXiv:1802.07432 [astro-ph.CO] .
  * Ewall-Wice  _et al._ (2018) A. Ewall-Wice, T.-C. Chang, J. Lazio, O. Dore, M. Seiffert, and R.A. Monsalve, “Modeling the Radio Background from the First Black Holes at Cosmic Dawn: Implications for the 21 cm Absorption Amplitude,” Astrophys. J. 868, 63 (2018), arXiv:1803.01815 [astro-ph.CO] .
  * Pospelov  _et al._ (2018) Maxim Pospelov, Josef Pradler, Joshua T. Ruderman, and Alfredo Urbano, “Room for New Physics in the Rayleigh-Jeans Tail of the Cosmic Microwave Background,” Phys. Rev. Lett. 121, 031103 (2018), arXiv:1803.07048 [hep-ph] .
  * Fialkov and Barkana (2019) Anastasia Fialkov and Rennan Barkana, “Signature of Excess Radio Background in the 21-cm Global Signal and Power Spectrum,” Mon. Not. Roy. Astron. Soc. 486, 1763–1773 (2019), arXiv:1902.02438 [astro-ph.CO] .
  * Reis  _et al._ (2020) Itamar Reis, Anastasia Fialkov, and Rennan Barkana, “High-redshift radio galaxies: a potential new source of 21-cm fluctuations,” Mon. Not. Roy. Astron. Soc. 499, 5993–6008 (2020), arXiv:2008.04315 [astro-ph.CO] .
  * Dalal  _et al._ (2010) Neal Dalal, Ue-Li Pen, and Uros Seljak, “Large-Scale Bao Signatures of the Smallest Galaxies,” JCAP 1011, 007 (2010), arXiv:1009.4704 [astro-ph.CO] .
  * McQuinn and O’Leary (2012) Matthew McQuinn and Ryan M. O’Leary, “The impact of the supersonic baryon-dark matter velocity difference on the z~20 21cm background,” Astrophys. J. 760, 3 (2012), arXiv:1204.1345 [astro-ph.CO] .
  * Visbal  _et al._ (2012) Eli Visbal, Rennan Barkana, Anastasia Fialkov, Dmitriy Tseliakhovich, and Christopher Hirata, “The signature of the first stars in atomic hydrogen at redshift 20,” Nature 487, 70 (2012), arXiv:1201.1005 [astro-ph.CO] .
  * Muñoz (2019) Julian B. Muñoz, “Robust Velocity-induced Acoustic Oscillations at Cosmic Dawn,” Phys. Rev. D 100, 063538 (2019), arXiv:1904.07881 [astro-ph.CO] .
  * de Putter  _et al._ (2019) Roland de Putter, Olivier Doré, Jérôme Gleyzes, Daniel Green, and Joel Meyers, “Dark matter interactions, helium, and the cosmic microwave background,” Phys. Rev. Lett. 122, 041301 (2019).
  * Boddy  _et al._ (2018) Kimberly K. Boddy, Vera Gluscevic, Vivian Poulin, Ely D. Kovetz, Marc Kamionkowski, and Rennan Barkana, “Critical Assessment of CMB Limits on Dark Matter-Baryon Scattering: New Treatment of the Relative Bulk Velocity,” Phys. Rev. D98, 123506 (2018), arXiv:1808.00001 [astro-ph.CO] .
  * Wouthuysen (1952) S. A. Wouthuysen, “On the excitation mechanism of the 21-cm (radio-frequency) interstellar hydrogen emission line.” Astron. J. 57, 31–32 (1952).
  * Field (1958) George B. Field, “Excitation of the Hydrogen 21-cm Line,” Proceedings of the IRE 46, 240–250 (1958).
  * Dvorkin  _et al._ (2014) Cora Dvorkin, Kfir Blum, and Marc Kamionkowski, “Constraining Dark Matter-Baryon Scattering with Linear Cosmology,” Phys. Rev. D 89, 023519 (2014), arXiv:1311.2937 [astro-ph.CO] .
  * Xu  _et al._ (2018) Weishuang Linda Xu, Cora Dvorkin, and Andrew Chael, “Probing sub-GeV Dark Matter-Baryon Scattering with Cosmological Observables,” Phys. Rev. D 97, 103530 (2018), arXiv:1802.06788 [astro-ph.CO] .
  * Gluscevic and Boddy (2018) Vera Gluscevic and Kimberly K. Boddy, “Constraints on Scattering of keV–TeV Dark Matter with Protons in the Early Universe,” Phys. Rev. Lett. 121, 081301 (2018), arXiv:1712.07133 [astro-ph.CO] .
  * Fialkov  _et al._ (2014a) Anastasia Fialkov, Rennan Barkana, Arazi Pinhas, and Eli Visbal, “Complete history of the observable 21-cm signal from the first stars during the pre-reionization era,” Mon. Not. Roy. Astron. Soc. 437, 36 (2014a), arXiv:1306.2354 [astro-ph.CO] .
  * Chen and Miralda-Escude (2004) Xue-Lei Chen and Jordi Miralda-Escude, “The spin - kinetic temperature coupling and the heating rate due to Lyman - alpha scattering before reionization: Predictions for 21cm emission and absorption,” Astrophys. J. 602, 1–11 (2004), arXiv:astro-ph/0303395 .
  * Chuzhoy and Shapiro (2006) Leonid Chuzhoy and Paul R. Shapiro, “Uv pumping of hyperfine transitions in the light elements, with application to 21-cm hydrogen and 92-cm deuterium lines from the early universe,” Astrophys. J. 651, 1–7 (2006), arXiv:astro-ph/0512206 .
  * Furlanetto and Pritchard (2006) Steven R. Furlanetto and Jonathan R. Pritchard, “The scattering of Lyman-series photons in the intergalactic medium,” Mon. Not. Roy. Astron. Soc. 372, 1093–1103 (2006), arXiv:astro-ph/0605680 [astro-ph] .
  * Venumadhav  _et al._ (2018) Tejaswi Venumadhav, Liang Dai, Alexander Kaurov, and Matias Zaldarriaga, “Heating of the intergalactic medium by the cosmic microwave background during cosmic dawn,” Phys. Rev. D 98, 103513 (2018), arXiv:1804.02406 [astro-ph.CO] .
  * Ali-Haïmoud  _et al._ (2014) Yacine Ali-Haïmoud, P. Daniel Meerburg, and Sihan Yuan, “New light on 21 cm intensity fluctuations from the dark ages,” Phys. Rev. D89, 083506 (2014), arXiv:1312.4948 [astro-ph.CO] .
  * Mertens  _et al._ (2021) F. G. Mertens, B. Semelin, and L. V. E. Koopmans, “Exploring the Cosmic Dawn with NenuFAR,” in  _Semaine de l’astrophysique française 2021_ (2021) arXiv:2109.10055 [astro-ph.CO] .
  * Koopmans  _et al._ (2015) L. Koopmans, J. Pritchard, G. Mellema, J. Aguirre, K. Ahn, R. Barkana, I. van Bemmel, G. Bernardi, A. Bonaldi, F. Briggs, A. G. de Bruyn, T. C. Chang, E. Chapman, X. Chen, B. Ciardi, P. Dayal, A. Ferrara, A. Fialkov, F. Fiore, K. Ichiki, I. T. Illiev, S. Inoue, V. Jelic, M. Jones, J. Lazio, U. Maio, S. Majumdar, K. J. Mack, A. Mesinger, M. F. Morales, A. Parsons, U. L. Pen, M. Santos, R. Schneider, B. Semelin, R. S. de Souza, R. Subrahmanyan, T. Takeuchi, H. Vedantham, J. Wagg, R. Webster, S. Wyithe, K. K. Datta, and C. Trott, “The Cosmic Dawn and Epoch of Reionisation with SKA,” in  _Advancing Astrophysics with the Square Kilometre Array (AASKA14)_ (2015) p. 1, arXiv:1505.07568 [astro-ph.CO] .
  * Perez and Granger (2007) Fernando Perez and Brian E. Granger, “IPython: A System for Interactive Scientific Computing,” Computing in Science and Engineering 9, 21–29 (2007).
  * Kluyver  _et al._ (2016) Thomas Kluyver _et al._ , “Jupyter notebooks - a publishing format for reproducible computational workflows,” in  _ELPUB_ (2016).
  * Hunter (2007) J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing In Science & Engineering 9, 90–95 (2007).
  * Harris  _et al._ (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant, “Array programming with NumPy,” Nature 585, 357–362 (2020).
  * Talman (1978) James D Talman, “Numerical fourier and bessel transforms in logarithmic variables,” Journal of Computational Physics 29, 35–48 (1978).
  * Hamilton (2000) A. J. S. Hamilton, “Uncorrelated modes of the nonlinear power spectrum,” Mon. Not. Roy. Astron. Soc. 312, 257–284 (2000), arXiv:astro-ph/9905191 .
  * Werthmüller (2020) Dieter Werthmüller, “prisae/pyfftlog: First packaged release,”  (2020).
  * Waskom  _et al._ (2017) Michael Waskom _et al._ , “mwaskom/seaborn: v0.8.1 (september 2017),”  (2017).
  * McKinney (2010) Wes McKinney, “Data structures for statistical computing in python,” in  _Proceedings of the 9th Python in Science Conference_ , edited by Stéfan van der Walt and Jarrod Millman (2010) pp. 51 – 56.
  * Virtanen  _et al._ (2020) Pauli Virtanen _et al._ , “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods (2020), https://doi.org/10.1038/s41592-019-0686-2.
  * da Costa-Luis (2019) Casper O da Costa-Luis, “tqdm: A fast, extensible progress meter for python and cli,” JOSS 4, 1277 (2019).
  * Field (1959) George B. Field, “The Time Relaxation of a Resonance-Line Profile.” Astrophys. J. 129, 551 (1959).
  * Pritchard and Loeb (2012) Jonathan R. Pritchard and Abraham Loeb, “21-cm cosmology,” Rept. Prog. Phys. 75, 086901 (2012), arXiv:1109.6012 [astro-ph.CO] .
  * Zygelman (2005) B. Zygelman, “Hyperfine Level-changing Collisions of Hydrogen Atoms and Tomography of the Dark Age Universe,” Astrophys. J. 622, 1356–1362 (2005).
  * Furlanetto and Furlanetto (2007a) Steven Furlanetto and Michael Furlanetto, “Spin Exchange Rates in Electron-Hydrogen Collisions,” Mon. Not. Roy. Astron. Soc. 374, 547–555 (2007a), arXiv:astro-ph/0608067 .
  * Furlanetto and Furlanetto (2007b) Steven Furlanetto and Michael Furlanetto, “Spin Exchange Rates in Proton-Hydrogen Collisions,” Mon. Not. Roy. Astron. Soc. 379, 130–134 (2007b), arXiv:astro-ph/0702487 .
  * Fialkov  _et al._ (2014b) Anastasia Fialkov, Rennan Barkana, and Eli Visbal, “The Observable Signature of Late Heating of the Universe during Cosmic Reionization,” Nature 506, 197 (2014b), arXiv:1402.0940 [astro-ph.CO] .
  * Fialkov  _et al._ (2013) Anastasia Fialkov, Rennan Barkana, Eli Visbal, Dmitriy Tseliakhovich, and Christopher M. Hirata, “The 21-cm signature of the first stars during the Lyman-Werner feedback era,” Mon. Not. Roy. Astron. Soc. 432, 2909 (2013), arXiv:1212.0513 [astro-ph.CO] .
  * Barkana and Loeb (2004) Rennan Barkana and Abraham Loeb, “Unusually large fluctuations in the statistics of galaxy formation at high redshift,” Astrophys. J. 609, 474–481 (2004), arXiv:astro-ph/0310338 .
  * Press and Schechter (1974) William H. Press and Paul Schechter, “Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation,” Astrophys. J. 187, 425–438 (1974).
  * Bond  _et al._ (1991) J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, “Excursion Set Mass Functions for Hierarchical Gaussian Fluctuations,” Astrophys. J. 379, 440 (1991).
  * Sheth and Tormen (1999) Ravi K. Sheth and Giuseppe Tormen, “Large scale bias and the peak background split,” Mon. Not. Roy. Astron. Soc. 308, 119 (1999), arXiv:astro-ph/9901122 .
  * Cohen  _et al._ (2019) Aviad Cohen, Anastasia Fialkov, Rennan Barkana, and Raul Monsalve, “Emulating the Global 21-Cm Signal from Cosmic Dawn and Reionization,”  (2019), 10.1093/mnras/staa1530, arXiv:1910.06274 [astro-ph.CO] .
  * Cohen  _et al._ (2018) Aviad Cohen, Anastasia Fialkov, and Rennan Barkana, “Charting the Parameter Space of the 21-cm Power Spectrum,” Mon. Not. Roy. Astron. Soc. 478, 2193–2217 (2018), arXiv:1709.02122 [astro-ph.CO] .

Observable 21-cm Fluctuations from Millicharged Dark Matter

Supplemental Material
Rennan Barkana, Anastasia Fialkov, Hongwan Liu and Nadav Outmezguine

In this supplemental material, we first review 21-cm cosmology for high-energy
physicists that are new to this area. We next provide all of the details
necessary for our power spectrum calculation, first writing down a general
expression for the two-point correlation function, before examining the small-
separation and large-separation limits of the function. We also provide some
useful expressions for the bulk relative velocity correlation functions.
Finally, we end with an elaboration of the simulations that we ran as part of
our work in the main body.

##  Appendix A 21-cm Cosmology for High-Energy Physicists

In this section, we briefly review the physics of the 21-cm emission in a
language familiar to high-energy physicists, and detail the full calculation
for obtaining the brightness temperature $T_{21}$, paying special attention to
deriving expressions that are valid even when the spin temperature is
comparable to the hyperfine splitting, a scenario that is possible due to the
significant cooling of baryons in the two-fluid dark sector model. For an
extensive review, we refer the readers to Ref. Barkana (2016).

Consider a line-of-sight through a large region of space, parametrized by
comoving coordinates $\bm{x}$ and redshift $z$. Radiation from the cosmic
microwave background (CMB) with temperature $T_{\gamma}(z)$ is incident on
neutral hydrogen or HI gas at $\bm{x}$ with number density
$n_{\text{HI}}(\bm{x},z)$. Since the 21-cm line is very narrow (it has a decay
width of $A_{10}=$2.85\text{\times}{10}^{-15}\text{\,}{\mathrm{s}}^{-1}$$,
compared to the energy level separation corresponding to a frequency of
$\nu_{10}=$1.42\text{\,}\mathrm{GHz}$$), we can safely treat all interactions
with neutral hydrogen as occurring only when the energy of the photon is
exactly given by the hyperfine splitting $\omega_{10}\equiv 2\pi\nu_{10}$ at
each point in $z$. Photons approaching the point with frequency just above
$\nu_{10}$ enter with the CMB blackbody intensity
$I_{\nu,\text{BB}}(\nu=\nu_{10},T_{\gamma}(z))$, with intensity being defined
with respect to frequency $\nu$.111This assumes that the only 21-cm photons
passing through this point originates only from the CMB, and that the CMB is a
perfect blackbody with no other sources of distortion. At this point, 21-cm
photons can be absorbed or emitted by the gas, leading to a change in its
intensity $\Delta I_{\nu_{10}}(z)$. This change in the intensity with respect
to the blackbody is the observable in 21-cm cosmology.

There are three processes at redshift $z$ that contribute to $\Delta
I_{\nu_{10}}(z)$: spontaneous emission, absorption and stimulated emission. In
the study of radiative transfer, it is conventional to treat absorption and
stimulated emission together, while spontaneous emission is regarded as a
source-term at each point. Under the narrow-width approximation, we can write
the absorption cross section as $\sigma_{0\to 1}\equiv
S\delta(\omega-\omega_{10})$, where $S$ is a constant, so that the number of
photons absorbed per volume per time is simply
$n_{0}(dn_{\gamma}/d\omega)|_{\omega_{10}}S$, where
$(dn_{\gamma}/d\omega|_{\omega_{10}})$ is the number density of photons per
unit energy at frequency $\omega_{10}$, where $n_{0}$ is the number density of
neutral hydrogen atoms in the ground state of the hyperfine splitting.
Similarly, stimulated emission can be encapsulated in an effective cross
section $\sigma_{1\to 0}=R\delta(\omega-\omega_{10})$, and the number of
photons emitted by stimulated emission per volume per time written as
$n_{1}(dn_{\gamma}/d\omega)|_{\omega_{10}}R$ for another constant $R$, where
$n_{1}$ is the number density of neutral hydrogen atoms in the excited state
of the hyperfine splitting. Then the usual detailed balance argument used in
deriving the Einstein coefficients tells us that in equilibrium at any
temperature $T$,

$\displaystyle
n_{1}^{\text{eq}}A_{10}+n_{1}^{\text{eq}}R\left.\frac{dn_{\gamma}^{\text{eq}}}{d\omega}\right|_{\omega_{10}}=n_{0}^{\text{eq}}S\left.\frac{dn_{\gamma}^{\text{eq}}}{d\omega}\right|_{\omega_{10}}\,,$

where “eq” denotes equilibrium quantities. The various number densities are
given by $n_{1}^{\text{eq}}/n_{0}^{\text{eq}}=3e^{-\omega_{10}/T}$, and

$\displaystyle\left.\frac{dn_{\gamma}^{\text{eq}}}{d\omega}\right|_{\omega_{10}}=\frac{\omega_{10}^{2}}{\pi^{2}}\frac{1}{e^{\omega_{10}/T}-1}$

Using the fact that detailed balance applies at any temperature, we find the
following relationship:

$\displaystyle R=\frac{S}{3}=\frac{\pi^{2}}{\omega_{10}^{2}}A_{10}\,.$ (7)

We now define the optical depth of a photon passing through the point at
redshift $z$ to be

$\displaystyle\tau(z)=\int
dt\,[n_{0}(z)S-n_{1}(z)R]\delta(\omega-\omega_{10})=\int
dt\,n_{0}(z)\left[1-e^{-\omega_{10}/T_{\text{S}}(z)}\right]S\,\delta(\omega-\omega_{10})\,,$
(8)

where we have used the definition of the spin temperature,

$\displaystyle n_{1}/n_{0}\equiv 3\exp(-\omega_{10}/T_{\text{S}}).$ (9)

Notice that the optical depth is defined using the net rate of absorption and
stimulated emission. As the photon travels through the point at redshift $z$,
it comes into resonance with $\dot{\omega}=H(z)\omega$. We can therefore
rewrite the integral in terms of $d\omega=H(z)\omega\,dt$, and perform the
integral to obtain

$\displaystyle\tau(z)=\frac{3\pi^{2}A_{10}n_{\text{HI}}(z)}{H(z)\omega_{10}^{3}}\frac{1-e^{-\omega_{10}/T_{\text{S}}(z)}}{1+3e^{-\omega_{10}/T_{\text{S}}(z)}}\,,$

with $n_{\text{HI}}=n_{0}+n_{1}$. In this expression, we have avoided the
common approximation $\omega_{10}\ll T_{\text{S}}$, since baryonic cooling in
our model can lead to very small values of $T_{\text{b}}$ and hence
$T_{\text{S}}$. In the limit $\omega_{10}\ll T_{\text{S}}$, owing to the
highly populated excited triplet state, the optical depth is small, and is
numerically given by Barkana (2016)

$\displaystyle\tau(z)\simeq 9.85\times
10^{-3}\left(\frac{T_{\gamma}(z)}{T_{\text{S}}(z)}\right)\left(\frac{\Omega_{\text{b}}h}{0.0327}\right)\left(\frac{\Omega_{\text{m}}}{0.307}\right)^{-1/2}\left(\frac{1+z}{10}\right)^{1/2}\,.$
(10)

After passing through the gas, the 21-cm intensity changes because of a
combination of absorption and stimulated emission—encapsulated by
$\tau(z)$—and spontaneous emission. We can write the change in intensity as

$\displaystyle\Delta
I_{\nu_{10}}(z)=-I_{\nu,\text{BB}}(\nu_{10},T_{\gamma}(z))\left[1-e^{-\tau(z)}\right]+S(z)\,,$
(11)

where $I_{\nu,\text{BB}}(\nu_{10},z)$ is the incoming intensity of 21-cm CMB
blackbody radiation, while $S(z)$ is the contribution from spontaneous
emission, which only depends on the properties of neutral hydrogen atoms. We
know that in thermal equilibrium, i.e. if the incoming intensity were a
perfect blackbody with temperature $T_{\text{S}}$, then we must have $\Delta
I_{\nu_{10}}(z)=0$. From this, we can conclude that

$\displaystyle
S(z)=I_{\nu,\text{BB}}(\nu_{10},T_{\text{S}}(z))\left[1-e^{-\tau(z)}\right]\,.$

In radio astronomy, the intensity $I_{\nu}$ at a particular frequency $\nu$ is
often expressed as a brightness temperature $\Theta$ instead, with
$\Theta(\omega)\equiv 2\pi^{2}I_{\nu}/\omega^{2}$. For a blackbody, the
relation between $\Theta$ and the thermodynamic temperature $T$ is

$\displaystyle\Theta(\omega)=T\times\frac{\xi}{e^{\xi}-1}\,,$ (12)

where $\xi\equiv\omega/T$, with $\Theta=T$ in the limit $\xi\to 0$. The
expression in Eq. (11) can therefore be written as

$\displaystyle\Delta\Theta(z)=\left[\frac{\xi(z)}{e^{\xi(z)}-1}T_{\text{S}}(z)-T_{\gamma}(z)\right]\left[1-e^{-\tau(z)}\right]\,,$

where $\xi(z)\equiv\omega_{10}/T_{\text{S}}$, and we have taken
$\omega_{10}\ll T_{\gamma}(z)$. Finally, the observed 21-cm brightness
temperature is precisely this absorption or emission relative to the
background $T_{\gamma}(z)$, redshifted to the present day, i.e.

$\displaystyle
T_{21}(z)=\frac{1}{1+z}\left[\frac{\xi(z)}{e^{\xi(z)}-1}T_{\text{S}}(z)-T_{\gamma}(z)\right]\left[1-e^{-\tau(z)}\right]\,.$
(13)

Once again, we have not taken the usual approximation $\omega_{10}\ll
T_{\text{S}}(z)$ or equivalently $\xi(z)\ll 1$, since this assumption can be
violated with baryonic cooling. Adopting this limit allows one to drop the
$\xi(z)/(e^{\xi(z)}-1)$ term, recovering the more usual expression Barkana
(2016).

Other than $\Lambda$CDM parameters, the only remaining unknown parameter that
determines $T_{21}(z)$ is the spin temperature. $T_{\text{S}}$ is determined
by a competition between 1) scattering of HI atoms with the CMB, which causes
$T_{\text{S}}\to T_{\gamma}$; 2) collisions between HI atoms, which causes
$T_{\text{S}}\to T_{\text{b}}$, and 3) Ly$\alpha$ scattering, which couples
$T_{\text{S}}$ to the color temperature of the Ly$\alpha$ photons
$T_{\text{S}}\to T_{\text{C}}$ through the Wouthuysen-Field (WF) effect
Wouthuysen (1952); Field (1958, 1959). The spin temperature can be expressed
as a weighted mean Field (1958)

$\displaystyle
T_{\text{S}}^{-1}=\frac{T_{\gamma}^{-1}+x_{c}T_{\text{b}}^{-1}+x_{\alpha}T_{\text{C}}^{-1}}{1+x_{c}+x_{\alpha}}\,,$
(14)

where $x_{c}$ and $x_{\alpha}$ represent coupling coefficients through
collisions and Ly$\alpha$ scattering respectively. The collisional coupling
coefficient can be written as Pritchard and Loeb (2012)

$\displaystyle
x_{c}=\frac{\omega_{10}}{A_{10}T_{\gamma}}\left[\kappa_{1-0}^{\text{HH}}(T_{\text{b}})n_{\text{HI}}+\kappa_{1-0}^{e\text{H}}(T_{\text{b}})n_{e}+\kappa_{1-0}^{p\text{H}}(T_{\text{b}})n_{p}\right]\,,$
(15)

where $n_{p}$ and $n_{e}$ are the number densities of free protons and
electrons; the rates $\kappa$ are calculated and tabulated in Refs. Zygelman
(2005); Furlanetto and Furlanetto (2007a, b). The Ly$\alpha$ coefficient is
Madau  _et al._ (1997); Barkana (2016)

$\displaystyle
x_{\alpha}=\frac{16\pi^{2}\alpha\omega_{10}}{27A_{10}m_{e}T_{\gamma}}J_{\alpha}\,,$
(16)

where $\alpha\approx 1/137$ is the fine-structure constant, and $J_{\alpha}$
is defined as

$\displaystyle
J_{\alpha}=\int\frac{d\Omega}{4\pi}\frac{I_{\nu_{\alpha}}}{\omega_{\alpha}}\,,$

where $I_{\nu_{\alpha}}$ is the photon intensity at the Ly$\alpha$ frequency,
and $\omega_{\alpha}\approx$10.2\text{\,}\mathrm{e\kern-0.50003ptV}$$ is the
Ly$\alpha$ transition energy. Including the effect of atomic recoil during
Ly$\alpha$ scattering as well as the possibility of multiple scatterings
allows us to write Eq. (14) as

$\displaystyle
T_{\text{S}}^{-1}=\frac{T_{\gamma}^{-1}+x_{\text{tot,eff}}T_{\text{b}}^{-1}}{1+x_{\text{tot,eff}}}\,,$
(17)

where $x_{\text{tot,eff}}=x_{c}+x_{\alpha,\text{eff}}$, and

$\displaystyle x_{\alpha,\text{eff}}\equiv
x_{\alpha}\left(1+\frac{T_{\text{se}}}{T_{\text{b}}}\right)^{-1}\exp\left[-2.06\left(\frac{\Omega_{\text{b}}h}{0.0327}\right)^{1/3}\left(\frac{\Omega_{\text{m}}}{0.307}\right)^{-1/6}\left(\frac{1+z}{10}\right)^{1/2}\left(\frac{T_{\text{b}}}{T_{\text{se}}}\right)^{-2/3}\right]\,.$
(18)

In this expression, $T_{\text{se}}\equiv
m_{\text{H}}(\omega_{10}/\omega_{\alpha})^{2}\approx$0.402\text{\,}\mathrm{K}$$,
with $m_{\text{H}}$ being the mass of the hydrogen atom Barkana (2016). For
the results presented in the Letter, we adopted 140 phenomenologically viable
models of $x_{\alpha}(z)$ from Ref. Reis  _et al._ (2021), chosen for
particular large values of $x_{\alpha}(z)$ so that $T_{\text{S}}$ is coupled
strongly to the $T_{\text{b}}$, leading to optimistically but realistically
large values of $T_{21}$ that can be potentially probed in near-future 21-cm
experiments.

##  Appendix B Details of the Power Spectrum Calculation

As we discussed in the main Letter, the statistics of the baryon-CDM bulk
relative velocity $v_{\text{bC}}$ fully specifies the statistics of $T_{21}$
in the two-fluid dark sector model. $v_{\text{bC}}(\bm{x})$ in this model is a
Gaussian random field, as it is in $\Lambda$CDM cosmology; it is fully
specified by its two-point function $P_{v}(k)$, defined as
$\langle\tilde{v}_{\text{bC}}(\bm{k})\tilde{v}_{\text{bC}}(\bm{k}^{\prime})\rangle=(2\pi)^{3}\delta^{(3)}(\bm{k}+\bm{k}^{\prime})P_{v}(k)$,
where
$\hat{k}\tilde{v}_{\text{bC}}(\bm{k})\equiv\tilde{\bm{v}}_{\text{bC}}(\bm{k})$,
the Fourier transform of $\bm{v}_{\text{bC}}(\bm{x})$. In particular, the one-
point probability density function (PDF) $f(\bm{v}_{\text{bC}})$ is

$\displaystyle f(\bm{v}_{\text{bC}})=\left(\frac{3}{2\pi
v_{\text{rms}}^{2}}\right)^{3/2}\exp\left(-\frac{3v_{\text{bC}}^{2}}{2v_{\text{rms}}^{2}}\right)\,,$
(19)

where $v_{\text{rms}}\equiv\langle
v_{\text{bC}}^{2}(\bm{x})\rangle^{1/2}\approx$29\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$
at $z=1010$. With this expression, we can easily obtain the mean value of
$T_{21}$ at redshift $z$, $\langle T_{21}\rangle(z)$, given any particular
astrophysics model or new physics parameters, by integrating over
$\bm{v}_{\text{bC}}$, i.e.

$\displaystyle\langle T_{21}\rangle(z)=\int
d^{3}\bm{v}_{\text{bC}}f(\bm{v}_{\text{bC}})T_{21}(v_{\text{bC}},z)\,.$ (20)

Similarly, the two-point correlation function (2PCF) is given by

$\displaystyle\xi_{T_{21}}(\bm{x})=\int d^{3}\bm{v}_{\text{bC},a}\int
d^{3}\bm{v}_{\text{bC},b}\mathcal{P}(\bm{v}_{\text{bC},a},\bm{v}_{\text{bC},b};\bm{x})\delta_{T_{21}}(v_{\text{bC},a})\delta_{T_{21}}(v_{\text{bC},b})\,,$
(21)

where $\delta_{T_{21}}\equiv(T_{21}-\langle T_{21}\rangle)/\langle
T_{21}\rangle$ is the spatial fluctuation of $T_{21}$, and $\mathcal{P}$ is
the two-point PDF, i.e. the joint probability density function of 3D
velocities at two different points, separated by the vector $\bm{x}$, given by

$\displaystyle\mathcal{P}(\bm{v}_{\text{bC},a},\bm{v}_{\text{bC},b};\bm{x})$
$\displaystyle=\frac{1}{(2\pi)^{3}\sqrt{|\mathsf{C}(\bm{x})|}}\exp\left[-\frac{1}{2}\vec{U}^{\mathsf{T}}\mathsf{C}^{-1}(\bm{x})\vec{U}\right]\,,$
(22)

where $\vec{U}\equiv(\bm{v}_{\text{bC},a},\bm{v}_{\text{bC},b})$ is a 6D
vector, and $\mathsf{C}$ is the covariance matrix for this multivariate
Gaussian. To differentiate between different types of quantities, we use a
boldface letters to represent 3D vectors, an arrow to indicate 6D vectors,
sans-serif font for $6\times 6$ matrices, and underscores for $3\times 3$
matrices. Note that ultimately, after performing the velocity integrals,
$\xi_{T_{21}}$ is only a function of $x\equiv|\bm{x}|$ in a homogeneous and
isotropic Universe.

The covariance matrix of the 6D Gaussian is given by Dalal  _et al._ (2010);
Ali-Haïmoud  _et al._ (2014)

$\frac{\mathsf{C}(\bm{x})}{\sigma_{v}^{2}}=\mathbb{1}_{6\times
6}+\left(\begin{array}[]{cc}0&\underline{c}(\bm{x})\\\
\underline{c}(\bm{x})&0\end{array}\right).$ (23)

Here, $\sigma_{v}^{2}\equiv v_{\rm rms}^{2}/3=\int d\log k\,\Delta^{2}_{v}(k)$
is the 1D velocity dispersion, with $\Delta^{2}_{v}(k)\equiv
k^{3}P_{v}(k)/(2\pi^{2})$. The elements of the $3\times 3$ matrix
$\underline{c}$ are given by

$\displaystyle\underline{c}(\bm{x})_{ij}\equiv\frac{\langle
v_{\text{bC}}^{i}v_{\text{bC}}^{j}\rangle(\bm{x})}{\sigma_{v}^{2}}=c_{\|}(x)\hat{\bm{x}}^{i}\hat{\bm{x}}^{j}+c_{\perp}(x)\left(\delta^{ij}-\hat{\bm{x}}^{i}\hat{\bm{x}}^{j}\right)\,,$
$\displaystyle
c_{\|}(x)=\frac{1}{v_{\mathrm{rms}}^{2}}\int\frac{dk}{k}\Delta_{v}^{2}(k)\left[j_{0}(kx)-2j_{2}(kx)\right]\,,$
(24) $\displaystyle
c_{\perp}(x)=\frac{1}{v_{\mathrm{rms}}^{2}}\int\frac{dk}{k}\Delta_{v}^{2}(k)\left[j_{0}(kx)+j_{2}(kx)\right]\,,$

where $j_{\ell}$ are the spherical Bessel functions of order $\ell$, $i$ and
$j$ denote spatial components, $\delta^{ij}$ is the Kronecker delta symbol,
and $\hat{\bm{x}}$ is a unit vector in the direction of $\bm{x}$.
$c_{\parallel}$ and $c_{\perp}$ give the correlation of the velocity component
parallel and perpendicular to the separation vector $\bm{x}$ respectively. For
later uses, in a coordinate system where $\bm{x}=x\hat{z}$, the 3D matrix
takes the form

$\underline{c}(x\hat{z})=\left(\begin{array}[]{ccc}c_{\|}&0&0\\\ 0&c_{\|}&0\\\
0&0&c_{\perp}\end{array}\right).$ (25)

Throughout this appendix, we find it convenient to use
$\bm{u}\sigma_{v}=\bm{v}_{\text{bC}}$ as our integration variable. Using the
Fourier transform of a 6D Gaussian, the 2PCF of any function $f(v)$ can be
expressed as

$\xi_{f}(x)=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\int\frac{d^{3}\bm{\omega}_{a}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\omega}_{b}}{(2\pi)^{3}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{a}}e^{i\bm{\omega}_{b}\cdot\bm{u}_{b}}\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\tilde{\cal
P}(\bm{\omega}_{a},\bm{\omega}_{b};\bm{x}),$ (26)

where $\delta_{f}\equiv f/\left\langle f\right\rangle-1$, ($\bm{\omega}_{a}$,
$\bm{\omega}_{b}$) are the Fourier transforms of ($\bm{u}_{a}$,
$\bm{\omega}_{b}$), and

$\tilde{\cal
P}(\bm{\omega}_{a},\bm{\omega}_{b};\bm{x})=\exp\left[-\frac{1}{2}\left(\omega_{a}^{2}+\omega_{b}^{2}\right)-\bm{\omega}_{a}^{T}\underline{c}(\bm{x})\bm{\omega}_{b}\right].$
(27)

Since $T_{21}$ is an isotropic function of velocity, integrating over the
angles of both $d^{3}u$ can be performed easily

$d^{3}\bm{u}_{a}\frac{d^{3}\bm{\omega}_{a}}{(2\pi)^{3}}d^{3}\bm{u}_{b}\frac{d^{3}\bm{\omega}_{b}}{(2\pi)^{3}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{a}}e^{i\bm{\omega}_{b}\cdot\bm{u}_{b}}=u_{a}^{2}du_{a}\frac{d^{3}\bm{\omega}_{a}}{2\pi^{2}}u_{b}^{2}du_{b}\frac{d^{3}\bm{\omega}_{b}}{2\pi^{2}}\frac{\sin(u_{a}\omega_{a})}{u_{a}\omega_{a}}\frac{\sin(u_{b}\omega_{b})}{u_{b}\omega_{b}}\,.$
(28)

We can simplify the non-diagonal part of the Gaussian in Eq. (27) by noting
that

$\bm{\omega}_{a}^{T}\underline{c}(\bm{x})\bm{\omega}_{b}=\omega_{a}\left|\underline{c}(\bm{x})\bm{\omega}_{b}\right|\cos\theta_{a}=\omega_{a}\omega_{b}\cos\theta_{a}\cdot
R(x,\cos\theta_{b})\,,$ (29)

where $\theta_{a}$ is the angle between $\bm{\omega}_{a}$ and
$\underline{c}\bm{\omega}_{b}$, and $\theta_{b}$ is the angle between
$\bm{\omega}_{b}$ and $\bm{x}$. Where through Eq. (25) we have identified

$\displaystyle
R(x,\cos\theta_{b})\equiv|\underline{c}(\bm{x})\bm{\omega}_{b}|=\sqrt{c_{\perp}^{2}(x)+\cos^{2}\theta_{b}\left[c_{\parallel}^{2}(x)-c_{\perp}^{2}(x)\right]}\,.$
(30)

This result can be directly computed using the expressions in Eq. (B).

At this point, we can integrate over the two remaining azimuthal angles and
$\theta_{a}$; after some algebra, we arrive at

$\xi_{f}(x)=\frac{2}{\pi}\int du_{1}\,u_{1}^{2}\int
du_{2}\,u_{2}^{2}\cdot\delta_{f}(u_{1}\sigma_{v})\delta_{f}(u_{2}\sigma_{v})W(u_{1},u_{2},x),$
(31)

with

$W(u_{1},u_{2},x)=\frac{1}{u_{1}u_{2}}\int_{-1}^{1}\frac{d\cos\theta_{b}}{2}\frac{\exp\left[-\frac{1}{2}\frac{u_{1}^{2}+u_{2}^{2}}{1-R^{2}}\right]\sinh\left[\frac{Ru_{1}u_{2}}{1-R^{2}}\right]}{R\sqrt{1-R^{2}}}\,.$
(32)

Eq. (31) can now be integrated numerically for a given
$\delta_{T_{21}}(v_{\text{bC}})$ to obtain the 2PCF, $\xi_{T_{21}}(x)$. Note
that the numerical integration that needs to be performed is a 3D integral,
which is numerically much simpler to evaluate than the equivalent 4D integral
expressions found in Refs. Dalal  _et al._ (2010); Ali-Haïmoud  _et al._
(2014).

Implementing Eq. (31) numerically poses some mild numerical challenges. For
separations $x$ much shorter than the sound horizon, which sets the typical
scale of spatial fluctuations, $c_{\perp},c_{\parallel}\to 1$ and $R\to 1$,
and the 3D integral becomes extremely peaked and difficult to evaluate
numerically with a regular mesh. On the other hand, for distances much larger
than the sound horizon, the correlation between the two points becomes weak; a
series expansion in the limit of large $x$ produces a simple expression for
$\xi_{T_{21}}$ that is highly accurate, allowing us to avoid performing the
relatively expensive numerical integration. We will now discuss the small $x$
and large $x$ limits in turn.

###  II.I The small separation limit

At distances much smaller than the sound horizon, we expect
$c_{\parallel},c_{\perp}\to 1$. To understand what happens in this limit, we
first write Eq. (26) as

$\displaystyle\xi_{f}(x)=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\int\frac{d^{3}\bm{\omega}_{a}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\omega}_{b}}{(2\pi)^{3}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{a}}e^{i\bm{\omega}_{b}\cdot\bm{u}_{b}}\exp\left[-\frac{1}{2}(\bm{\omega}_{a}+\bm{\omega}_{b})^{2}+\bm{\omega}_{a}^{\mathsf{T}}\underline{d}(x)\bm{\omega}_{b}\right]\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\,,$
(33)

with $\underline{d}\equiv 1-\underline{c}$, which is a small parameter in this
limit. We now rewrite
$\exp(\bm{\omega}_{a}^{\intercal}\underline{d}(x)\bm{\omega}_{b})$ as
derivatives acting on
$\exp(i\bm{\omega}_{a}\cdot\bm{u}_{a})\exp(i\bm{\omega}_{b}\cdot\bm{u}_{b})$,
to obtain

$\displaystyle\xi_{f}(x)$ $\displaystyle=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\exp\left[-\partial_{a}^{i}\partial_{b}^{j}\underline{d}_{ij}\right]\int\frac{d^{3}\bm{\omega}_{a}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\omega}_{b}}{(2\pi)^{3}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{a}}e^{i\bm{\omega}_{b}\cdot\bm{u}_{b}}\exp\left[-\frac{1}{2}(\bm{\omega}_{a}+\bm{\omega_{b}})^{2}\right]\,,$
(34)

where we introduce the notation $\partial_{a}^{i}$ is the partial derivative
with respect to the component $u_{a}^{i}$ (and likewise for
$\partial_{b}^{j}$), and we adopt Einstein summation convention from here on.
The last two integrals are now inverse Fourier transforms that we can
evaluate, to obtain

$\displaystyle\xi_{f}(x)$ $\displaystyle=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\exp\left[-\partial_{a}^{i}\partial_{b}^{j}\underline{d}_{ij}\right]\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\delta_{D}^{(3)}(\bm{u}_{a}-\bm{u}_{b})$
$\displaystyle=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{D}^{(3)}(\bm{u}_{a}-\bm{u}_{b})\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\exp\left[-\partial_{a}^{i}\partial_{b}^{j}\underline{d}_{ij}\right]\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\,,$
(35)

where in the last line we have performed an integration by parts to move the
derivatives over.

We are now ready to expand in terms of the small parameter
$\underline{d}_{ij}$. Expanding the exponential to second order, we have

$\displaystyle\exp\left[-\partial_{a}^{i}\partial_{b}^{j}\underline{d}_{ij}\right]\simeq
1-\partial_{a}^{i}\partial_{b}^{j}\underline{d}_{ij}+\frac{1}{2}\partial_{a}^{i}\partial_{b}^{j}\partial_{a}^{k}\partial_{b}^{l}\underline{d}_{ij}\underline{d}_{kl}\,.$
(36)

At leading order, we have

$\displaystyle\xi_{f}^{(0)}(x\to 0)=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{D}^{(3)}(\bm{u}_{a}-\bm{u}_{b})\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})=\langle\delta_{f}^{2}(v)\rangle\,,$
(37)

since the PDF of 3D bulk relative velocities is
$g(v)d^{3}\bm{v}=(2\pi\sigma_{v}^{2})^{-3/2}\exp[-v^{2}/(2\sigma_{v}^{2})]d^{3}\bm{v}$,
with $v=\sigma_{v}u$.

At the next order, evaluating the derivatives, we find

$\displaystyle\xi_{f}^{(1)}(x\to 0)$ $\displaystyle=-d_{ij}\int
d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{D}^{(3)}(\bm{u}_{a}-\bm{u}_{b})\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\frac{u_{a}^{i}u_{b}^{j}}{u_{a}u_{b}}\sigma_{v}^{2}\delta_{f}^{\prime}(u_{a}\sigma_{v})\delta_{f}^{\prime}(u_{b}\sigma_{v})$
$\displaystyle=-d_{ij}\int
d^{3}\bm{u}_{a}\,\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\frac{u_{a}^{i}u_{a}^{j}}{u_{a}^{2}}\sigma_{v}^{2}\delta_{f}^{\prime
2}(u_{a}\sigma_{v})\,.$ (38)

At this point, we make use of the identity

$\displaystyle\int d^{3}\bm{u}\,u^{i}u^{j}f(u)=\frac{\delta^{ij}}{3}\int
d^{3}\bm{u}\,u^{2}f(u)\,,$ (39)

which can be deduced from the fact that the symmetry of the integral makes it
proportional to the rank-2 isotropic tensor, i.e. the Kronecker delta
function. Putting everything together, we finally obtain

$\displaystyle\xi_{f}^{(1)}(x\to
0)=-\frac{1}{3}\underline{d}\indices{{}^{i}_{i}}\sigma_{v}^{2}\langle\delta_{f}^{\prime
2}\rangle=-\frac{1}{3}{\rm
Tr}\left[\underline{d}\right]\sigma_{v}^{2}\langle\delta_{f}^{\prime
2}\rangle\,$

where $\rm Tr$ is the trace of a matrix.

Finally, at the second order, we first expand the derivatives carefully:

$\displaystyle\partial_{a}^{i}\partial_{a}^{k}=\frac{\delta^{ik}}{u_{a}}\partial_{a}+\frac{u_{a}^{i}u_{a}^{k}}{u_{a}^{2}}\left(\partial_{a}^{2}-\frac{1}{u_{a}}\partial_{a}\right)\,,$
(40)

which can be applied to the second order expression to give

$\displaystyle\xi_{f}^{(2)}(x\to
0)=\frac{1}{2}\underline{d}_{ij}\underline{d}_{kl}\int
d^{3}\bm{u}_{b}\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\left[\left(\frac{\delta^{ik}}{u_{b}}-\frac{u_{b}^{i}u_{b}^{k}}{u_{b}^{3}}\right)\sigma_{v}\delta_{f}^{\prime}(u_{a}\sigma_{v})+\frac{u_{b}^{i}u_{b}^{k}}{u_{b}^{2}}\sigma_{v}^{2}\delta_{f}^{\prime\prime}(u_{a}\sigma_{v})\right]\times(i\to
j,k\to l)\,.$ (41)

To simplify the expression, we apply Eq. (39), as well as the analogous result
at rank-4:

$\displaystyle\int
d^{3}\bm{u}u^{i}u^{j}u^{k}u^{l}f(u)=\frac{1}{15}(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\int
d^{3}\bm{u}\,u^{4}f(u)\equiv\frac{1}{15}S_{ijkl}\int
d^{3}\bm{u}\,u^{4}f(u)\,,$ (42)

to obtain

$\displaystyle\xi_{f}^{(2)}(x\to 0)$
$\displaystyle=\frac{1}{2}\underline{d}_{ij}\underline{d}_{kl}\int
d^{3}\bm{u}_{b}\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\left(\left[\frac{\delta^{ik}\delta^{jl}}{3u_{b}^{2}}+\frac{S_{ijkl}}{15u_{b}^{2}}\right]\sigma_{v}^{2}\delta_{f}^{\prime
2}+2\left[\frac{\delta^{ik}\delta^{jl}}{3u_{b}}-\frac{S_{ijkl}}{15u_{b}}\right]\sigma_{v}^{3}\delta_{f}^{\prime}\delta_{f}^{\prime\prime}+\frac{S_{ijkl}}{15}\sigma_{v}^{4}\delta_{f}^{\prime\prime
2}\right)$ $\displaystyle=\frac{1}{2}\underline{d}_{ij}\underline{d}_{kl}\int
d^{3}\bm{u}_{b}\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}\left(\left[\frac{1}{3}\delta^{ik}\delta^{jl}+\frac{2S_{ijkl}}{15u_{b}^{2}}-\frac{S_{ijkl}}{15}\right]\sigma_{v}^{2}\delta_{f}^{\prime
2}+\frac{S_{ijkl}}{15}\sigma_{v}^{4}\delta_{f}^{\prime\prime 2}\right)\,,$
(43)

where in the last line we have performed an integration by parts. Contracting
the indices gives finally

$\displaystyle\xi_{f}^{(2)}(x\to
0)=\left[3\underline{d}^{ik}\underline{d}_{ik}-\underline{d}\indices{{}^{i}_{i}}\underline{d}\indices{{}^{k}_{k}}\right]\frac{\sigma_{v}^{2}\langle\delta_{f}^{\prime
2}\rangle}{30}+\frac{1}{30}\left(2\underline{d}^{ik}\underline{d}_{ik}+\underline{d}\indices{{}^{i}_{i}}\underline{d}\indices{{}^{k}_{k}}\right)\left(2\sigma_{v}^{4}\left\langle\frac{\delta_{f}^{\prime
2}}{v_{\text{bC}}^{2}}\right\rangle+\sigma_{v}^{4}\langle\delta_{f}^{\prime\prime
2}\rangle\right)\,.$ (44)

The combined expression can then be written as

$\langle f\rangle^{2}\xi_{f}(x\to 0)\simeq\langle f^{2}\rangle-\langle
f\rangle^{2}\\\ -\frac{1}{9}\left\\{{\rm Tr}[\underline{d}]-\frac{1}{2}{\rm
Tr}\left[\underline{d}^{2}\right]\right\\}\langle(v_{\text{rms}}f^{\prime})^{2}\rangle+\frac{1}{270}\left\\{{\rm
Tr}\left[\underline{d}\right]^{2}+2{\rm
Tr}\left[\underline{d}^{2}\right]\right\\}\left<\left(2\frac{v_{\text{rms}}^{2}}{v_{\text{bC}}^{2}}-3\right)(v_{\text{rms}}f^{\prime})^{2}+(v_{\text{rms}}^{2}f^{\prime\prime})^{2}\right>\,.$
(47)

###  II.II Large separation limit

At large separations, both $c_{\parallel}$ and $c_{\perp}$ approach $0$. Once
again, we can start from Eq. (26) and write

$\displaystyle\xi_{f}(x)=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\exp\left[\partial_{a}^{i}\partial_{b}^{j}\underline{c}_{ij}\right]\int\frac{d^{3}\bm{\omega}_{a}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\omega}_{b}}{(2\pi)^{3}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{a}}e^{i\bm{\omega}_{a}\cdot\bm{u}_{b}}\exp\left[-\frac{1}{2}(\omega_{a}^{2}+\omega_{b}^{2})\right]\,.$
(48)

We can immediately take the Fourier transform of the Gaussian to find

$\displaystyle\xi_{f}(x)=\int d^{3}\bm{u}_{a}\int
d^{3}\bm{u}_{b}\,\delta_{f}(u_{a}\sigma_{v})\delta_{f}(u_{b}\sigma_{v})\exp\left[\partial_{a}^{i}\partial_{b}^{j}\underline{c}_{ij}\right]\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\frac{e^{-u_{b}^{2}/2}}{(2\pi)^{3/2}}.$
(49)

At this point, we can expand the exponential in the limit of small
$\underline{c}$, noting that any term with odd derivatives in $u_{a}$ is zero,
since the resulting expression is odd under $\bm{u}_{a}\to-\bm{u}_{a}$.
Because $\langle\delta_{f}\rangle$ is zero, the leading order result goes as
two powers of $\underline{c}$:

$\displaystyle\xi_{f}^{(1)}(x\to\infty)=\frac{1}{2}\underline{c}_{ij}\underline{c}_{kl}\left(\int
d^{3}\bm{u}_{a}\,\delta_{f}(u_{a}\sigma_{v})\partial_{a}^{i}\partial_{a}^{k}\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\right)\times(a\to
b,i\to j,k\to l).$ (50)

As in the small separation limit, we can exploit the symmetry of the integral
to rewrite $\partial_{a}^{i}\partial_{a}^{k}\to\partial_{a}^{2}\delta^{ik}/3$,
giving

$\displaystyle\xi_{f}^{(1)}(x\to\infty)$
$\displaystyle=\frac{1}{18}\underline{c}^{ik}\underline{c}_{ik}\left(\int
d^{3}\bm{u}_{a}\,\delta_{f}(u_{a}\sigma_{v})\partial_{a}^{2}\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\right)^{2}\,,$
(51)

which we can evaluate to obtain

$\displaystyle\xi_{f}^{(1)}(x\to\infty)=\frac{1}{2}{\rm
Tr}\left[\underline{c}^{2}\right]\frac{\langle
v^{2}\delta_{f}\rangle^{2}}{v_{\text{rms}}^{4}}\,.$ (52)

At the next order, we have

$\displaystyle\xi_{f}^{(2)}(x\to\infty)=\frac{1}{24}\underline{c}_{ij}\underline{c}_{kl}\underline{c}_{mn}\underline{c}_{pq}\left(\int
d^{3}\bm{u}_{a}\,\delta_{f}(u_{a}\sigma_{v})\partial_{a}^{i}\partial_{a}^{k}\partial_{a}^{m}\partial_{a}^{p}\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\right)\times(a\to
b,i\to j,k\to l,m\to n,p\to q)\,.$ (53)

Once again, we can replace
$\partial_{a}^{i}\partial_{a}^{k}\partial_{a}^{m}\partial_{a}^{p}\to
S^{ikmp}\bm{\nabla}_{a}^{2}\bm{\nabla}_{a}^{2}/15$, where
$\bm{\nabla}_{a}^{2}\equiv\partial_{a}^{i}\partial_{a,i}$, leading to

$\displaystyle\xi_{f}^{(2)}(x\to\infty)=\frac{S^{ikmp}S^{jlnq}}{24(15)^{2}}\underline{c}_{ij}\underline{c}_{kl}\underline{c}_{mn}\underline{c}_{pq}\left(\int
d^{3}\bm{u}_{a}\,\delta_{f}(u_{a}\sigma_{v})\bm{\nabla}_{a}^{2}\bm{\nabla}_{a}^{2}\frac{e^{-u_{a}^{2}/2}}{(2\pi)^{3/2}}\right)^{2}\,.$
(54)

Contracting the tensor indices gives

$\displaystyle
S^{ikmp}S^{jlnq}\underline{c}_{ij}\underline{c}_{kl}\underline{c}_{mn}\underline{c}_{pq}$
$\displaystyle=S^{ikmp}\left(\underline{c}^{2}_{ik}\underline{c}^{2}_{mp}+\underline{c}^{2}_{im}\underline{c}^{2}_{kp}+\underline{c}^{2}_{ip}\underline{c}^{2}_{km}\right)$
$\displaystyle=3\left[(\underline{c}^{ik}\underline{c}_{ik})^{2}+2(\underline{c}^{2})^{ik}(\underline{c}^{2})_{ik}\right]=3\left\\{{\rm
Tr}\left[\underline{c}^{2}\right]^{2}+2{\rm
Tr}\left[\underline{c}^{4}\right]\right\\}\,.$ (55)

Evaluating the integral and simplifying leads us to the final expression,

$\displaystyle\xi_{f}^{(2)}(x\to\infty)=\frac{1}{8}\left\\{{\rm
Tr}\left[\underline{c}^{2}\right]^{2}+2{\rm
Tr}\left[\underline{c}^{4}\right]\right\\}\left(\frac{3}{5}\frac{\langle
v^{4}\delta_{f}\rangle}{v_{\text{rms}}^{4}}-2\frac{\langle
v^{2}\delta_{f}\rangle}{v_{\text{rms}}^{2}}\right)^{2}\,.$ (56)

The final combined result is

$\displaystyle\xi_{f}(x\to\infty)\simeq\frac{1}{2}{\rm
Tr}\left[\underline{c}^{2}\right]\frac{\langle
v^{2}\delta_{f}\rangle^{2}}{v_{\text{rms}}^{4}}+\frac{1}{8}\left\\{{\rm
Tr}\left[\underline{c}^{2}\right]^{2}+2{\rm
Tr}\left[\underline{c}^{4}\right]\right\\}\left(\frac{3}{5}\frac{\langle
v^{4}\delta_{f}\rangle}{v_{\text{rms}}^{4}}-2\frac{\langle
v^{2}\delta_{f}\rangle}{v_{\text{rms}}^{2}}\right)^{2}\,,$ (57)

which to leading order is in agreement with Ref. Ali-Haïmoud  _et al._ (2014).

###  II.III Velocity correlation functions

Finally, in this section, we compute the correlation functions for powers of
$v_{\text{bC}}$, which we denote as $v$ for simplicity in this section. We
also write $\xi_{f}$ in terms of various bias factors multiplying velocity
correlation functions, which is a common approximation scheme used in the
literature.

The correlation function for $v^{2}$ and $v^{4}$ are both simple Gaussian
integrals, resulting with

$\xi_{v^{2}}(x)=\frac{2}{9}{\rm
Tr}\left[\underline{c}^{2}\right]\,,\qquad\xi_{v^{4}}(x)=4\xi_{v^{2}}+\frac{8}{225}\left\\{{\rm
Tr}\left[\underline{c}^{2}\right]^{2}+2{\rm
Tr}\left[\underline{c}^{4}\right]\right\\}.$ (58)

Alternatively, these expressions could have been obtained by noting that the
large separation expansion is exact to first order for $\xi_{v^{2}}$ and to
second order for $\xi_{v^{4}}$. At large separation, we see that
$\xi_{v^{4}}\simeq 4\xi_{v^{2}}$ to leading order, in agreement with Eq. (57).
It is common to write the correlation function $\xi_{f}$ as a bias parameter
multiplied by $\xi_{v^{2}}$ and $\xi_{v^{4}}$; comparing our expressions here
and Eq. (57), we find in the large separation limit

$\xi_{f}(x\to\infty)\simeq\frac{9}{4}\left(\frac{\left\langle
v^{2}\delta_{f}\right\rangle}{\left\langle
v^{2}\right\rangle}\right)^{2}\xi_{v^{2}}(x)-\frac{225}{64}\left(2\frac{\left\langle
v^{2}\delta_{f}\right\rangle}{\left\langle
v^{2}\right\rangle}-\frac{\left\langle
v^{4}\delta_{f}\right\rangle}{\left\langle
v^{4}\right\rangle}\right)^{2}\left[4\xi_{v^{2}}(x)-\xi_{v^{4}}(x)\right].$
(59)

In the small separation limit, we write $\underline{c}=1-\underline{d}$ as
before, and noting that
$\underline{c}^{2}=1-2\underline{d}+\underline{d}^{2}$, we obtain

$\xi_{v^{2}}(x)=\frac{2}{3}-\frac{4}{9}{\rm
Tr}\left[\underline{d}\right]+\frac{2}{9}{\rm
Tr}\left[\underline{d}^{2}\right],$ (60)

which is an exact expression, while keeping terms up to order
$\underline{d}^{2}$, we find

$\xi_{v^{4}}(x\to 0)\simeq\frac{16}{5}-\frac{112}{45}\left({\rm
Tr}\left[\underline{d}\right]-\frac{1}{2}{\rm
Tr}\left[\underline{d}^{2}\right]\right)+\frac{32}{225}\left({\rm
Tr}\left[\underline{d}\right]^{2}+2{\rm
Tr}\left[\underline{d}^{2}\right]\right)\,.$ (61)

The above two expressions are in agreement with Eq. (47), which we can now
rewrite as

$\displaystyle\bar{\xi}_{f}(x\to 0)$
$\displaystyle\simeq\left(\frac{A}{4}-\frac{7B}{48}\right)\bar{\xi}_{v^{2}}(x)+\frac{5B}{192}\bar{\xi}_{v^{4}}(x)\,,$
$\displaystyle\bar{\xi}(x)=\xi(x)-\xi(0)\,,\qquad A$
$\displaystyle=\left\langle\left[v_{\rm
rms}f^{\prime}(v)\right]^{2}\right\rangle\,,\qquad
B=\left\langle\left(2\frac{v_{\rm rms}^{2}}{v^{2}}-3\right)\left[v_{\rm
rms}f^{\prime}(v)\right]^{2}+\left[v_{\rm
rms}^{2}f^{\prime\prime}(v)\right]^{2}\right\rangle\,.$ (62)

##  Appendix C Details of the Simulations

We rely on a large-scale, semi-numerical 21-cm code based on Refs. Visbal  _et
al._ (2012); Fialkov  _et al._ (2014a); Reis  _et al._ (2021); Fialkov  _et
al._ (2014b). Driven by the specifications of radio telescopes such as the
Square Kilometer Array, the simulation models large cosmic volumes ($384^{3}$
comoving $\text{\,}{\mathrm{Mpc}}^{3}$) with a resolution of 3 comoving
$\text{\,}\mathrm{Mpc}$. The initial conditions for density fields, bulk
relative velocities between dark matter and baryons Tseliakhovich and Hirata
(2010); Fialkov  _et al._ (2013); Visbal  _et al._ (2012), and IGM temperature
are generated at $z=60$. The halo abundance is calculated within each
resolution element using the approach of Ref. Barkana and Loeb (2004), which
is based on Refs. Press and Schechter (1974); Bond  _et al._ (1991); Sheth and
Tormen (1999). The resulting number of halos is biased by the local values of
the large-scale density and velocity fields. Subsequently, star formation is
derived assuming that every halo with a mass higher than the star-formation
threshold will form stars at a given star-formation efficiency Cohen  _et al._
(2019), which is a function of halo mass as well as the local value of the
Lyman-Werner (LW) radiative background, which suppresses star formation via
the molecular-cooling channel Fialkov  _et al._ (2013).

Radiation produced by stars and stellar remnants is propagated taking into
account redshifting and absorption in the IGM. We follow the evolution and
spatial fluctuations in several key radiative backgrounds: X-rays heat up and
mildly ionize the IGM, Ly$\alpha$ photons are responsible for the WF coupling
and contribute to IGM heating Reis  _et al._ (2021), LW photons affect the
efficiency of star formation, and ionizing photons drive the process of
reionization. We then compute the 21-cm signal of neutral hydrogen affected by
all the sources of light within the light-cone. The simulation produces three
dimensional cubes of the fluctuating 21-cm signal at a selection of redshifts,
which can be used to calculate $\langle T_{21}\rangle$ and $P_{T_{21}}$.

Since the Universe at the time of primordial star formation is practically
observationally unconstrained, we perform multiple simulations, varying
several free astrophysical parameters within their allowed ranges Cohen  _et
al._ (2018); Reis  _et al._ (2021). The relevant astrophysical parameters
include the minimum circular velocity of star forming halos $V_{c}$, the star
formation efficiency $f_{*}$, the X-ray spectral energy distribution (SED),
X-ray heating efficiency $f_{X}$, the ionizing efficiency of sources, and the
mean free path of ionizing photons. In this Letter, we aim to highlight the
discovery potential of the two-fluid dark sector model in the 21-cm power
spectrum; we therefore take an ensemble of simulations of 140 astrophysical
models from this realistic parameter range, chosen to minimize astrophysical
heating and maximize Ly$\alpha$ coupling, both of which result in a stronger
absorption for $T_{21}$. Specifically, we choose models where
$V_{c}>$16.5\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$$, $f_{X}<0.001$
and $f_{\text{radio}}<1$ Reis  _et al._ (2020). Although it is included in the
model, the process of reionization has a subdominant impact on the high
redshift signals discussed in this paper.