#  CMB and 21cm bounds on early structure formation boosted by primordial
black hole entropy fluctuations

Hiroyuki Tashiro1 and Kenji Kadota2
1 Department of physics and astrophysics, Nagoya University, Nagoya 464-8602,
Japan
2 Center for Theoretical Physics of the Universe, Institute for Basic Science
(IBS), Daejeon, 34051, Korea

###### Abstract

The dark matter (DM) can consist of the primordial black holes (PBHs) in
addition to the conventional weakly interacting massive particles (WIMPs). The
Poisson fluctuations of the PBH number density produce the isocurvature
perturbations which can dominate the matter power spectrum at small scales and
enhance the early structure formation. We study how the WIMP annihilation from
those early formed structures can affect the CMB (in particular the E-mode
polarization anisotropies and $y$-type spectral distortions) and global 21cm
signals. Our studies would be of particular interest for the light (sub-GeV)
WIMP scenarios which have been less explored compared with the mixed DM
scenarios consisting of PBHs and heavy ($\gtrsim 1$ GeV) WIMPs. For instance,
for the self-annihilating DM mass $m_{\chi}=1$ MeV and the thermally averaged
annihilation cross section $\langle\sigma v\rangle\sim 10^{-30}\rm cm^{3}/s$,
the latest Planck CMB data requires the PBH fraction with respect to the whole
DM to be at most ${\cal O}(10^{-3})$ for the sub-solar mass PBHs and an even
tighter bound (by a factor $\sim 5$) can be obtained from the global 21-cm
measurements.

##  I introduction

While the properties of the dark matter (DM) still remain unknown, there has
been a revived interest in the primordial black hole (PBH) DM in view of the
advancement of gravitational wave experiments Abbott et al. (2016). While the
allowed parameter range for the PBH to be the dominant DM component has been
narrowed down, the PBH being a partial component of DM still remains an
intriguing possibility Carr and Kuhnel (2020); Green and Kavanagh (2020); Carr
et al. (2020a).

We study the effects of the DM annihilation on the cosmic microwave background
(CMB) and 21-cm signals in the mixed DM scenarios consisting of the self-
annihilating DM and PBH. We focus on the boosted DM annihilation due to the
enhanced early structure formation in the presence of the Poisson noise
sourced by the PBHs Afshordi et al. (2003).

The randomly distributed PBHs can add the Poisson noise to the matter power
spectrum which can dominate the conventional adiabatic perturbations at small
scales. The possibility for the PBHs to be (partial) DM and the consequent
early structure formation due to such PBH sourced isocurvature perturbations
have been investigated to seek the potential signals of the PBHs, and the
astrophysical probes sensitive to those enhanced small scale structures such
as the gravitational lensing and Ly-$\alpha$ observations have been explored
Afshordi et al. (2003); Ali-Haïmoud (2018); Inman and Ali-Haïmoud (2019);
Desjacques and Riotto (2018); Oguri and Takahashi (2020); Kadota and Silk
(2021); Gong and Kitajima (2017); Mena et al. (2019).

When there exists self-annihilating DM (a typical example is WIMP), in
addition to PBH as a partial DM, the thermal evolution of baryons in the dark
ages is also affected due to the enhanced DM annihilation from those abundant
early formed halos. The additional energy injection from DM annihilation may
produce observable CMB distortions Zeldovich and Sunyaev (1969); Sunyaev and
Zeldovich (1970); Burigana et al. (1991); Hu and Silk (1993). The early
reionization due to the enhancement of DM annihilation affects the Thomson
optical depth and modifies the CMB angular power spectrum Aghanim et al.
(2020). The measurement of global 21-cm line signal is also expected to be a
useful tool to probe the thermal history of the Universe Madau et al. (1997);
Shaver et al. (1999); Pritchard and Loeb (2010). The recent result reported by
the EDGES Bowman et al. (2018) motivates to study the constraint on extra heat
sources Clark et al. (2018); Hektor et al. (2018); Mitridate and Podo (2018);
Cheng et al. (2018); D’Amico et al. (2018); Safarzadeh et al. (2018); Minoda
et al. (2019) and the nature of DM Tashiro et al. (2014); Barkana (2018). We
in this paper quantify the effects of the boosted DM annihilation on the CMB
anisotropy (in particular the E-mode polarization power spectrum), CMB
spectral distortion and global 21cm signals.

The early structure formation can result in abundant small halos some of which
can potentially survive the tidal disruptions due to the tightly bound
structures. We calculate the effects of DM annihilation from the substructures
due to survived small halos as well as those from the halos without
considering the substructures. For instance, for $m_{\chi}=1$MeV and
$\langle\sigma v\rangle\sim 10^{-30}~{}\rm cm^{2}$, the Planck CMB data can
give the bound on the PBH fraction to the total DM $f_{\rm PBH}\lesssim
10^{-2.7}$ for the sub-solar mass PBH assuming the survival of the halos if
produced at $z>250$ (note the CMB bound in the absence of PBHs is
$\langle\sigma v\rangle/m_{\chi}\lesssim 3\times 10^{-27}\rm cm^{3}/s/GeV$ for
the DM with s-wave annihilation Slatyer (2016)). We for simplicity assume the
monochromatic mass function for the PBHs throughout this paper. We also find
that the effects of the survived minihalos are negligible for $f_{\rm
PBH}\lesssim 10^{-2.6}$, and the aforementioned bound $f_{\rm PBH}\lesssim
10^{-2.7}$ hence can be considered as the robust bound arising only from
isolated halos without including the effects of substructures. We also show
that, using the redshift dependence of the EDGES data (without using their
larger-than-expected absorption signal amplitude information), the 21cm global
signal can lead to the tighter bounds on $f_{\rm PBH}$ by a factor 5 than that
from the CMB.

We mention that, independently from such ’Poisson effects’ (PBHs’ collective
effects on the large scale structure), one can also study the ’seed effects’
(the DM accretion into an individual PBH) which can be complimentary to the
studies in this paper. The mixed DM scenarios of PBHs and self-annihilating
WIMPs indeed have been actively investigated for the scenarios where the PBH
can be a ’seed’ for the ultracompact minihalo (UCMH) with a steep DM density
profile Gondolo and Silk (1999); Lacki and Beacom (2010); Boucenna et al.
(2018); Adamek et al. (2019); Eroshenko (2016); Carr et al. (2020b); Cai et
al. (2020); Delos et al. (2018); Kohri et al. (2014); Bertone et al. (2019);
Carr and Silk (2018); Carr and Rees (1984); Hertzberg et al. (2020);
Kashlinsky (2016); Tashiro and Kadota (2021). No observation of possible
enhanced DM annihilation signals from such a DM density spike around a PBH can
lead to the severe bounds on the allowed abundance of PBHs, and the
incompatibility of co-existence of the PBH and WIMP has been pointed out. Such
’seed effects’ have been mainly discussed for the heavy ($>1$ GeV) DM and the
conventional cross sections of order $\langle\sigma v\rangle\sim{\cal
O}(10^{-26})~{}\rm cm^{3}/s$. For instance, the tight bounds $f_{\rm
PBH}\lesssim{\cal O}(10^{-9})$ can be obtained for the typical WIMP parameters
(e.g. $m_{\chi}\sim 100~{}{\rm GeV},\langle\sigma v\rangle=3\times
10^{-26}~{}{\rm cm^{3}/s}$) due to no detection of such enhanced DM
annihilation signals in the gamma ray or CMB data Adamek et al. (2019);
Tashiro and Kadota (2021). Those tight bounds on $f_{\rm PBH}$ are usually
derived based on the characteristic steep DM density profile $\rho\propto
r^{-\gamma}$ (e.g. $\gamma\sim 9/4$), and those bounds cannot be
straightforwardly applied when the DM kinetic energy cannot be ignored
compared with its potential energy in estimating the halo profile around a PBH
Gondolo and Silk (1999); Lacki and Beacom (2010); Boucenna et al. (2018);
Adamek et al. (2019); Eroshenko (2016); Carr et al. (2020b); Cai et al.
(2020); Delos et al. (2018); Kohri et al. (2014); Bertone et al. (2019). Such
a UCMH profile in the presence of the PBH for the non-trivial initial
conditions including those for the light (e.g. $\lesssim 1$ MeV) DM still has
not been fully explored by the numerical simulations yet without definite
answers for the bounds on the PBH parameters. Fortunately, our Poisson bounds
are independent from those bounds involving the non-trivial DM accretion to
the PBHs and hence can offer complementary bounds to those seed effects. Our
studies would be of particular interest for the light (sub-GeV) DM for which
the UCMH profile is less steep and the seed effect bounds are less severe. The
bounds from the DM accretion onto the PBH are also relaxed for the lighter PBH
mass as well as for the lighter DM mass, and we focus our discussions on the
light (sub-GeV) DM mass and the light PBH mass (sub-solar mass) even though
our study can be straightforwardly extended to the larger DM and PBH masses
for which the PBH seed effects however would give more stringent bounds on
$f_{\rm PBH}$ than the Poisson effects.

Our paper is organized as follows. Sec. II reviews the matter power spectrum
in the presence of the PBH isocurvature perturbations and specifies the PBH
parameter range to be explored in our study. Sec. III outlines the formalism
to calculate the DM annihilation boost in the presence of the early formed
structures. Sec IV discusses how the evolutions of ionization fraction and
baryon temperature are affected due to the DM annihilation in the presence of
the PBHs. Such changes in the early history of the Universe are implemented
for the calculations of the bounds on $f_{\rm PBH}$ from the CMB (fluctuation
anisotropy and spectral distortion) and 21cm signals in Sec V.

##  II The Poisson fluctuations due to PBHs

The PBHs are assumed to be randomly distributed111This is a reasonable
assumption because the typical separation between the PBHs which can be formed
before the matter-radiation equality epoch would be larger than the horizon
scale at the formation of PBHs., and the PBHs can contribute to the matter
power spectrum as the Poisson noise

$\displaystyle P_{\rm PBH}=\frac{1}{n_{\rm PBH}},~{}\quad n_{\rm
PBH}=\frac{\Omega_{\rm DM}\rho_{\rm cri}f_{\rm PBH}}{M_{\rm PBH}},$ (1)

where $f_{\rm PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM}$ represents the
fraction of PBH contribution to the total DM. We focus on the scales larger
than the mean separation of the PBHs, $k\lesssim k_{*}=n_{\rm PBH}^{1/3}$, so
that we treat the PBHs as the ideal pressureless fluid. Such PBH fluctuations
appear only in the PBH component and independent from the adiabatic
perturbations, and the PBH Poisson fluctuations can be treated as the
isocurvature contribution to the total power spectrum Afshordi et al. (2003);
Ali-Haïmoud (2018); Inman and Ali-Haïmoud (2019); Desjacques and Riotto (2018)

$\displaystyle P(k,z)=D^{2}(z)\left(T^{2}_{ad}(k)P_{\rm adi}(k)+T^{2}_{\rm
iso}P_{\rm iso}(k)\right),$ (2)

where $P_{\rm adi}$ is the conventional adiabatic power spectrum and $P_{\rm
iso}=f_{\rm PBH}^{2}P_{\rm PBH}$. $D(z)$ is the growth function normalized by
$D(0)=1$ and $T$ represents the transfer function. The isocurvature transfer
function $T_{\rm iso}$ reads Peacock (1999)

$\displaystyle T_{\rm iso}(k)$ $\displaystyle=$
$\displaystyle\frac{3}{2}(1+z_{\rm eq})\quad\mbox{ for }k_{\rm eq}<k<k_{*},$
(3) $\displaystyle T_{\rm iso}(k)$ $\displaystyle=$ $\displaystyle
0\quad\mbox{ otherwise.}$ (4)

The dimensionless power spectra $\Delta^{2}(k)=P(k)k^{3}/2\pi^{2}$ in the
presence of PBH are shown in Fig. 1. The dimensionless power spectrum
$\Delta^{2}(k)$ is smaller with a smaller $M_{\rm PBH}$ for a given $k$, but
the cutoff $k_{*}$ is bigger for a smaller $M_{\rm PBH}$ leading to the common
peak height with a different $M_{\rm PBH}$ for a given $f_{\rm PBH}$.

---
Figure 1:  The linear matter power spectrum in the mixed DM scenarios at
$z=0$. The PBH contributions are parameterized by its fraction and mass
$f_{\rm PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM},M_{\rm PBH}$. The scenario
with the PBH fraction smaller than $10^{-4}$ is not distinguishable from the
standard $\Lambda$CDM with no PBH in this figure.  Figure 2:  The mass
variance $\sigma(M,z=0)$ in the mixed DM scenarios. The PBH mass dependence
becomes small for a small halo mass scale $M\lesssim 10^{-6}M_{\odot}$ in the
parameter range of our interest $f_{\rm PBH}\lesssim 0.1$,
$10^{-6}[M_{\odot}/h]\lesssim M_{\rm PBH}\lesssim 1[M_{\odot}/h]$.

This independence of the peak amplitudes on $M_{\rm PBH}$ is also reflected in
the insensitivity of $\sigma$ to $M_{\rm PBH}$ for a small halo mass as
illustrated in Fig. 2. The variance of the mass fluctuations smoothed on the
scale of a sphere containing a mass $M$ is defined by

$\displaystyle\sigma^{2}(M,z)=\int d\ln
k\frac{k^{3}P(k,z)}{2\pi^{2}}\left|W(kR)\right|^{2},$ (5)

where $W(kR)=3\left[\sin(kR)-(kR)\cos(kR)\right]/(kR)^{3}$ is the Fourier
transform of the real-space spherical top-hat window function which contains
mass $M$.

We study the effects of early formed minihalos on the CMB and 21cm signals,
and the parameters, besides the conventional $\Lambda$CDM parameters, which
are relevant in our study are the minimal halo mass $M_{\rm min}$, PBH
fraction and mass ($f_{\rm PBH}$ and $M_{\rm PBH}$) and the annihilating DM
mass and the thermally averaged DM annihilation cross section ($m_{\chi}$ and
$\langle\sigma v\rangle$). The minimum halo mass can vary in a wide range (for
instance $10^{-4}M_{\odot}/h$ to $10^{-12}M_{\odot}/h$) which is heavily
dependent on the nature of DM kinetic decoupling affecting the DM free
streaming and acoustic damping. We for concreteness use the conventional value
of $M_{\rm min}=10^{-6}M_{\odot}/h$ in our calculations Bertschinger (2006);
Loeb and Zaldarriaga (2005); Gondolo et al. (2012); Profumo et al. (2006);
Gondolo and Kadota (2016); Green et al. (2004, 2005); Bringmann and Hofmann
(2007); Diemand et al. (2005). The mass scale of $10^{-6}M_{\odot}/h$
corresponds to the comoving scale of order $k\sim 10^{6}h/\rm Mpc$, and we
focus on the PBH parameters whose Poisson noise cutoff $k_{*}<10^{6}h/\rm Mpc$
so that the PBH fluctuation cutoff shows up above the DM minimum halos mass
scale. For instance, $k_{*}\gtrsim 4\times 10^{6}h/\rm Mpc$ for $M_{\rm
PBH}\lesssim 10^{-7}M_{\odot}/h$ with $f_{\rm PBH}=0.1$. We hence for
concreteness focus on $M_{\rm PBH}>10^{-6}M_{\odot}$ in our discussions unless
stated otherwise. We also note, for this parameter range of our interest,
$\sigma(M_{\rm min})$ is not dependent on $M_{\rm PBH}$ for a given $f_{PBH}$.
Due to this insensitivity of $\sigma(M_{\rm min})$ on $M_{\rm PBH}$, our
constraints on $f_{\rm PBH}$ turn out to be insensitive to $M_{\rm PBH}$ as
discussed in the following analysis.

##  III Annihilation boosting by the early structure formation due to PBHs

The annihilation of self-annihilating DM is proportional to the density
squared and the ratio of the fluctuation contribution to the smooth background
contribution is the so-called boost factor $B$

$\displaystyle\langle\rho_{\rm DM}^{2}(z)\rangle=\bar{\rho}_{\rm
DM}^{2}(z)\langle(1+\delta(z))^{2}\rangle=\bar{\rho}_{\rm
DM}^{2}(z)(1+\langle\delta^{2}(z)\rangle)\equiv\bar{\rho}_{\rm
DM}^{2}(z)(1+B(z))$ (6)

where $\bar{\rho}_{\rm DM}$ is the homogeneous DM background density and
$\delta=(\rho-\bar{\rho})/\bar{\rho}$ is the density contrast. We adopt the
halo model approach in calculating the boost factor to take account of the DM
structure formation. The energy injection rate per volume due to the DM
annihilation (we assume the self-annihilating Majorana DM $\chi$) into the
cosmic plasma reads Cooray and Sheth (2002); Poulin et al. (2015)

$\frac{d^{2}E}{dVdt}(z)=\left.\frac{d^{2}E}{dVdt}\right|_{\rm
sm}+\left.\frac{d^{2}E}{dVdt}\right|_{\rm str},$ (7)

where the subscript sm denotes the smooth background contributions and str
represents the DM halo contributions to account for the effects from the
structure formation.

The smooth background contribution is obtained from

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm sm}=f_{\rm ann}\frac{\langle\sigma
v\rangle}{m_{\chi}}\rho_{\rm cri}^{2}\Omega^{2}_{\chi}(1+z)^{6},$ (8)

where $f_{\rm ann}$ is the fraction of the energy absorbed into the plasma at
a redshift $z$. Although it depends on the annihilating DM mass and the
redshift, we assume for simplicity that $f_{\rm ann}$ is constant Galli et al.
(2013). $\Omega_{\chi}$ is the energy density parameter of the self-
annihilating DM $\chi$. Since we also consider PBHs as DM, $\Omega_{\rm
DM}=\Omega_{\rm PBH}+\Omega_{\chi}$.

The halo contribution is given as Giesen et al. (2012)

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm halo}=f_{\rm ann}\frac{\langle\sigma
v\rangle}{m_{\chi}}(1+z)^{3}\int^{\infty}_{M_{\rm
min}}dM~{}\frac{dn}{dM}(M,z)\left(\int^{r_{200}}_{0}dr~{}4\pi r^{2}\rho_{\rm
halo}^{2}(r)\right),$ (9)

where $M_{\rm min}$ is the minimum mass of DM halos, $dn/dM$ is the mass
function of DM halos and $r_{200}$ is the typical radius of a DM halo with
mass $M$ in which an averaged density equals to 200 times the background
density and $\rho_{\rm halo}$ is the radial density profile of self-
annihilating DM.

When we use the NFW profile for the halo density $\rho_{\rm halo}$ and take
into account that only a fraction $\Omega_{\chi}/\Omega_{\rm M}$ of the total
matter can annihilate, the integral of the self-annihilating DM density
profile can be written as

$\displaystyle\int^{r_{200}}_{0}dr4\pi r^{2}\rho_{\rm
halo}^{2}(r)=\frac{M\bar{\rho}(z_{\rm
f})}{3}\left(\frac{\Omega_{\chi}}{\Omega_{\rm M}}\right)^{2},$ (10)
$\displaystyle\bar{\rho}(z_{\rm f})=200\rho_{\rm cri}\Omega_{\rm M}(1+z_{\rm
f})^{3}F(c(z_{\rm f})),~{}\quad F(c(z_{\rm
f}))=\frac{c^{3}}{3}\frac{1-(1+c)^{-3}}{\left(\ln(1+c)-c(1+c)^{-1}\right)^{2}},$
(11)

where $z_{\rm f}$ is the formation redshift of DM halos. $c(z_{\rm f})$ is the
concentration parameter of the NFW profile for DM halos forming at $z_{\rm
f}$. We conservatively take $c(z_{\rm f})=2$ for all the early formed halos
because it is the minimal value of the concentration parameters found in the
simulations performed for a wide range of halo masses Wang et al. (2020);
Ishiyama et al. (2020); Sánchez-Conde and Prada (2014); Moline et al. (2017);
Zhao et al. (2009); Prada et al. (2012); Ackermann et al. (2015). A bigger
concentration parameter results in a bigger annihilation rate.

Assuming the Press-Schechter (PS) mass function can simplify the integral of
the mass function and the mass fraction in collapsed objects, $f_{\rm coll}$,
becomes, according to the PS formalism,

$f_{\rm coll}(z)=\int^{\infty}_{M_{\rm min}}dM~{}\frac{M}{\Omega_{\rm
M}\rho_{\rm cri}}\frac{dn}{dM}(M,z)={\rm erfc}(u_{\rm min}(z)),\quad u_{\rm
min}(z)=\frac{\delta_{c}}{\sqrt{2}\sigma(M_{\rm min},z)},$ (12)

where $\delta_{\rm c}\approx 1.686$ is the critical density contrast for
collapse. Therefore, we can rewrite the DM halo contribution as

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm halo}=\frac{\bar{\rho}(z_{\rm
f})}{3}\left(\frac{\Omega_{\chi}}{\Omega_{\rm M}}\right)^{2}f_{\rm
ann}\frac{\langle\sigma v\rangle}{m_{\chi}}(1+z)^{3}f_{\rm coll}(z).$ (13)

which represents the energy injection at the redshift $z$ from the halos
formed at $z_{\rm f}$.

In a simple treatment where the substructure contribution is ignored (such as
in the conventional PS formalism which does not take account of subhalo
abundance once they merge into the larger halos), $z_{\rm f}$ is set to
$z_{\rm f}=z$ in Eq. (13). The notable feature of our PBH scenarios arises due
to the enhancement of early structure formation caused by the PBH isocurvature
fluctuations. Those early formed halos are shown to be dense enough (the
density $\propto(1+z_{\rm f})^{3}$) to survive the potential tidal disruptions
for $z_{\rm f}\gtrsim 250$ Blinov et al. (2020); Kadota and Silk (2021)222The
tidal disruptions include the stellar encounters, encountering among small
halos and tidal stripping at the core of a host halo after the infall due to
the dynamical friction. We refer the readers to Refs. Zhao et al. (2007); van
den Bosch et al. (2018); Green and van den Bosch (2019); Kadota and Silk
(2021); Blinov et al. (2020); Berezinsky et al. (2003); Arvanitaki et al.
(2020); Berezinsky et al. (2006, 2003); Dai and Miralda-Escudé (2020); Xiao et
al. (2021) for the further discussions on the survival of the early formed
halos. One can analytically estimate the time-scale of tidal disruptions due
to stellar encounters (which are the most significant tidal disruption
processes), based on the impulse approximation, as Kadota and Silk (2021)
$\displaystyle t\sim
70Gyr\left(\frac{1+z_{f}}{100}\right)^{3/2}\left(\frac{10^{6}M_{\odot}/kpc^{3}}{\bar{n}_{*}M_{*}}\right)$
(14) where $M_{*}$ is the mass of a stellar object (such as a star) and
$\bar{{n}}_{*}$ is the corresponding mean number density of stars in a host
halo. We can hence infer that the minihalos in a galactic host halo could well
survive if they were formed at $z_{f}\geq 100$. Such early formed dense
structures resilient to the tidal disruptions make our scenarios strikingly
different from the $\Lambda$CDM for which the typical halo formation occurs at
$z\lesssim 20$.

We introduce a characteristic redshift $z_{*}$ (for concreteness, we use
$z_{*}=250$ in our calculations Blinov et al. (2020); Kadota and Silk (2021)).
We call the small halos produced at $z_{\rm f}>z_{*}$ ’minihalos’ and assume
those minihalos can survive until now while keeping the density
$200\bar{\rho}(z_{\rm f})$ at the forming epoch for simplicity. On the other
hand, we assume the halos forming at $z_{\rm f}<z_{*}$ cannot survive and lose
their identities in the hierarchical structure formation. In other words, for
the halos formed at $z_{\rm f}<z_{*}$, we apply the conventional PS formalism
(we hence do not account for the substructure contribution and ignore their
survivals at the later epoch as minihalos). The density of such halos at a
redshift $z$ is simply proportional to the background density
$200\bar{\rho}(z)$. The actual DM annihilation effects could be bigger than
our estimations for those halos produced at $z<z_{*}$ because we completely
ignore their survival. For the minihalos collapsed at $z>z_{*}$, we assume the
DM annihilation continues in those dense surviving minihalos even after their
collapse epochs. The effects of the survived minihalos compared with those
produced at $z<z_{*}$ become more prominent at a lower redshift because of a
bigger relative difference in $(1+z)^{3}$ factor in their densities. This can
be seen, for instance, in Fig. 3 where the baryon temperature in the presence
of the survived minihalos exceeds that without substructures at a low redshift
$z\lesssim 50$.

In order to take into account those early formed minihalos, we divide the
contribution from the DM structure formation in Eq. (7) into two parts,

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm
str}=\left.\frac{d^{2}E}{dVdt}\right|_{\rm
mh}+\left.\frac{d^{2}E}{dVdt}\right|_{\rm halo}.$ (15)

Here the term with the subscript mh represents the contribution from the
survived minihalos forming before $z_{*}$. Analogously to Eq. (9), the
minihalo contribution at a redshift $z$ can be given as

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm mh}=f_{\rm ann}\frac{\langle\sigma
v\rangle}{m_{\chi}}(1+z)^{3}\int_{z_{\rm min}}^{z_{\rm max}}dz_{\rm
f}\int^{\infty}_{M_{\rm min}}dM~{}\frac{d}{dz_{\rm f}}\left[\frac{dn_{\rm
mh}}{dM}(M,z_{\rm f})\right]\left(\int^{r_{200}}_{0}dr~{}4\pi r^{2}\rho_{\rm
halo}^{2}(r)\right),$ (16)

where ${dn_{\rm mh}}/{dM}(M,z_{\rm f})$ is the mass function of minihalos
forming at $z_{\rm f}$. In the above equation, $z_{\rm min}={\rm
max}[z_{*},z]$. The redshift $z_{\rm max}$ represents the maximum redshift for
the structure formation. Although the DM density fluctuation can grow
logarithmically in the radiation dominated epoch, it cannot collapse until the
matter dominated epoch. In our fiducial model, we set $z_{\rm max}=2000$ and
we conservatively assume that, even if its density contrast exceeds the
critical density contrast for collapse before $z_{\rm max}$, the overdensity
region collapses to a minihalo at $z_{\rm max}$.

In estimating Eq. (16), assuming that DM halos forming before $z_{*}=250$ can
survive as minihalos, we can consider the fraction of matter collapsed into
minihalos $f_{\rm coll,mh}$ as

$f_{\rm coll,mh}(z)=f_{\rm coll}({\rm max}[z_{*},z]),$ (17)

which can also be written as

$f_{\rm coll,mh}(z)=\int_{z}^{z_{\rm max}}dz_{\rm f}~{}\frac{df_{\rm
coll,mh}}{dz_{\rm f}}(z,z_{\rm f})=\frac{1}{\rho_{\rm cri}\Omega_{\rm
M}}\int_{z}dz_{\rm f}~{}\int_{M_{\rm min}}^{\infty}dMM\frac{d}{dz_{\rm
f}}\left[\frac{dn_{\rm mh}}{dM}(M,z_{\rm f})\right].$ (18)

Note that our formalism does not require the explicit expression for
differential fraction $df_{\rm coll,mh}/dz_{\rm f}$ or minihalo mass function
$dn_{\rm mh}/dM(M,z_{\rm f})$ because $f_{\rm coll,mh}$ can be obtained from
Eqs. (12) and (17) (see Ref. Delos et al. (2018) for more details on the
derivation of differential mass fraction $df/dz$). Accordingly, the minihalo
contribution can be rewritten with $f_{\rm coll,mh}$

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm
mh}=\left(\frac{\Omega_{\chi}}{\Omega_{\rm M}}\right)^{2}f_{\rm
ann}\frac{\langle\sigma v\rangle}{m_{\chi}}(1+z)^{3}\int dz_{\rm
f}\frac{\bar{\rho}(z_{\rm f})}{3}\frac{df_{\rm coll,mh}}{dz_{\rm f}}(z,z_{\rm
f}).$ (19)

In Eq. (15), ${d^{2}E}/{dVdt}|_{\rm halo}$ represents the contribution from DM
which is not included in a survived minihalo formed at $z>z_{*}$ but resides
in a halo formed at $z<z_{*}$. The collapse fraction of such DM is given by
$f_{\rm coll,halo}=f_{\rm coll}-f_{\rm coll,mh}$. Avoiding the double counting
of DM in minihalos in Eq. (13), ${d^{2}E}/{dVdt}|_{\rm halo}$ is given as

$\left.\frac{d^{2}E}{dVdt}\right|_{\rm halo}=\frac{\bar{\rho}(z_{\rm
f})}{3}\left(\frac{\Omega_{\chi}}{\Omega_{\rm M}}\right)^{2}f_{\rm
ann}\frac{\langle\sigma v\rangle}{m_{\chi}}(1+z)^{3}f_{\rm coll,halo}(z).$
(20)

The boost factor also consists of two parts,

$B(z)=B_{\rm mh}(z)+B_{\rm halo}(z).$ (21)

Each boost factor is given by the ratio of Eq. (19, 20) to the smooth
background contribution, Eq. (8). Therefore, we can calculate boost factors
from

$\displaystyle B_{\rm mh}(z)$ $\displaystyle=$
$\displaystyle\frac{200}{3}\int^{z_{\rm mh}}_{z}dz_{\rm f}\frac{(1+z_{\rm
f})^{3}}{(1+z)^{3}}F(c(z_{\rm f}))\frac{df_{\rm coll,mh}}{dz_{\rm f}}(z,z_{\rm
f}),$ (22) $\displaystyle B_{\rm halo}(z)$ $\displaystyle=$
$\displaystyle\frac{200}{3}F(c(z))f_{\rm coll,halo}(z).$ (23)

##  IV Evolution of the ionization fraction and the baryon temperature

The injected energy due to the WIMP annihilation can affect the evolutions of
the ionization fraction and temperature in the cosmic plasma. Including the
WIMP annihilation effects, these evolutions are calculated from

$\displaystyle(1+z)\frac{dx_{e}}{dz}$
$\displaystyle=\frac{1}{H(z)}\left[R_{s}(z)-I_{s}(z)-I_{x}(z)\right],$ (24)
$\displaystyle(1+z)\frac{dT_{k}}{dz}$
$\displaystyle=2T_{k}+\frac{8\sigma_{T}a_{R}T_{\gamma}^{4}}{3m_{e}cH(z)}\frac{x_{e}(T_{k}-T_{\gamma})}{(1+f_{\rm
He}+x_{e})}-\frac{2}{3k_{B}H(z)}\frac{K_{h}(z)}{(1+f_{\rm He}+x_{e})}.$ (25)

Eq. (24) gives the evolution of the ionization fraction, where $R_{s}$ and
$I_{s}$ are the standard primordial hydrogen recombination rate and ionization
rate respectively Peebles (1968); Zeldovich et al. (1968). Eq. (25) provides
the evolution of the baryon temperature. In these equations, the contributions
of the WIMP annihilation are represented in the terms $I_{x}$ and $K_{h}$,

$\displaystyle I_{x}=\frac{\chi_{x}}{{n_{\rm H}(z)}}\frac{d^{2}E}{dVdt}(z),$
(26) $\displaystyle K_{h}=\frac{\chi_{h}}{{n_{\rm
H}(z)}E_{i}}\frac{d^{2}E}{dVdt}(z).$ (27)

where ${d^{2}E(z)}/{dVdt}$ is the energy injection rate of the WIMP
annihilation given in Eq (7), $E_{i}$ is the ionization energy of hydrogen,
$\chi_{x}$ and $\chi_{h}$ are the fractions of energy used for the ionization
and heating of the cosmic plasma, respectively. In Eqs. (26) and (27), for
simplicity, we take the on-the-spot approximation in which the injected energy
is assumed to be absorbed into the cosmic plasma and, instantaneously,
consumed to ionize and heat the plasma Padmanabhan and Finkbeiner (2005). The
functions $\chi_{x}$ and $\chi_{h}$ mainly depend on the ionization fraction
of the Universe, $x_{e}$, and we use the fitting formulae given in Ref. Poulin
et al. (2017) which is based on the result in Ref. Galli et al. (2013).

The CMB bounds depend on the combination of parameters $f_{\rm
ann}\langle\sigma v\rangle/m_{\chi}$ rather than separately on each of the DM
mass and annihilation cross section, and the effects of DM annihilation are
conventionally parameterized by a quantity $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}$. The current CMB bound is of order $f_{\rm
ann}\langle\sigma v\rangle/m_{\chi}\lesssim 4\times 10^{-28}\rm cm^{3}/s/GeV$,
and the canonical thermal WIMP relic value $\langle\sigma v\rangle=3\times
10^{-26}\rm cm^{3}/s$ for the DM with s-wave annihilation, for instance, can
be already excluded for $m_{\chi}\lesssim 10$ GeV for a typical range of
$f_{\rm ann}\sim{\cal O}(0.1)$ Slatyer (2016). We for the illustration purpose
use $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times 10^{-28}\rm
cm^{3}/s/GeV$ as our fiducial value in the figures unless stated
otherwise.333Using for instance the value $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}=3\times 10^{-29}\rm cm^{3}/s/GeV$, instead of $3\times
10^{-28}\rm cm^{3}/s/GeV$, does not affect our following discussions, even
though the constraint on $f_{PBH}$ becomes less tight by a factor $\sim$ 5 as
mentioned in the next section. We solve Eqs. (24) and (25) by using the public
code HyRec Ali-Haïmoud and Hirata (2011) in which the primordial helium
contribution is also included. In Fig. 3, we show the effect of the WIMP
annihilation in the mixed DM scenario on the ionization fraction (left panel)
and the baryon temperature (right panel). In this figure using $f_{\rm
ann}\langle\sigma v\rangle/m_{\chi}=3\times 10^{-28}\rm cm^{3}/s/GeV$, the
blue line represents the evolution with $\log_{10}f_{\rm PBH}=-2.5$, and, for
comparison, we also plot the results for no PBH case (a conventional adiabatic
perturbation scenario without PBH isocurvature perturbations) and for no DM
annihilation case with the thin red line and black dashed line, respectively.
In order to illustrate the impact of minihalo contributions, we also plot the
evolution for the case with no survived minihalos (i.e., $f_{\rm coll,mh}=0$)
in the orange line in Fig. 3. The figure tells us that the survived minihalo
contribution is significant in the late universe. When $\log_{10}f_{\rm
PBH}=-2.5$, the mass dispersion is small $\sigma(M_{\rm
min},z_{*})=0.4~{}(<\delta_{c}\sim 1.686)$ at $z_{*}=250$ and the collapsed
matter fraction locked in the survived minihalos would be consequently small.
As a result, the survived minihalo contribution is not significant until a
later epoch $z\sim 50$. As the universe expands further, the core density in
survived minihalos relatively becomes significantly large, compared with the
background density or the DM halo density. Around $z\sim 50$, the boost factor
of survived minihalos cannot be negligible and finally dominates. The
ionization and heating from the annihilation in survived minihalos can
significantly affect the Thomson optical depth of the CMB and 21-cm signals in
the dark ages. We also show the evolution of the ionization fraction and
baryon temperature for different $f_{\rm PBH}$ values in Fig. 4. A large
$f_{\rm PBH}$ induces the early structure formation and the effect of the halo
contribution start to become important at higher redshifts. The survived
minihalo contribution becomes smaller when $f_{\rm PBH}$ becomes smaller, and
we found that, when $f_{\rm PBH}\lesssim 10^{-2.6}$, the contribution from the
survived minihalos, the first term in the right-hand side of Eq. (15), is
negligible in the thermal evolution of the Universe.

---
Figure 3: The ionization fraction (left panel) and the baryon temperature
(right panel) as functions of a redshift. We set $f_{\rm PBH}=10^{-2.5}$. The
blue line represents the total contribution from DM, while the orange line is
for the case without survived minihalos ($f_{\rm coll,mh}=0$). For reference,
we plot the evolution without PBHs and without annihilation.
---
Figure 4: The dependence on the PBH abundance $f_{\rm PBH}$ in the ionization
fraction (left panel) and the baryon temperature (right panel). From the top
to the bottom, the solid lines represent the cases with $f_{\rm
PBH}=10^{-2.4}$, $f_{\rm PBH}=10^{-2.5}$ and $f_{\rm PBH}=10^{-2.6}$.

##  V cosmological constraints on the PBH abundance

###  V.1 CMB anisotropy

As shown in the last section, the WIMP annihilation can be enhanced in the
mixed DM scenario and hence induce the early reionization. One of cosmological
probes on such early reionization is the CMB anisotropy measurement. In
particular, the CMB polarization is created due to the scattering of free
electrons in the propagation of CMB photons from the last scattering surface.
The free electron density indeed increases for a bigger $f_{\rm PBH}$ as shown
in Fig. 4.

Taking into account the enhancement of the ionization fraction, we evaluate
the CMB polarization anisotropy by the public Boltzmann code CLASS Blas et al.
(2011). Fig. 5 shows the angular power spectra of the CMB E-mode polarization
for different $f_{\rm PBH}$ values. We adopt the ”tanh”-shape reionization
history with $z_{\rm reio}=7.68$ Aghanim et al. (2018) to include the standard
reionization scenario driven by first stars and galaxies. The figure tells us
that the early reionization due to a large $f_{\rm PBH}$ increases the CMB
polarization signal represented by the change in the reionization bump. The
detailed CMB anisotropy measurement hence can provide a constraint on the PBH
abundance. In order to obtain the constraint from the current CMB measurement,
we perform the MCMC analysis using Monte Python Audren et al. (2013). In this
analysis we use the baseline likelihood (TTTEEE-lowl-lowE) from the Planck
2018 data release Aghanim et al. (2020) and adopt the standard six
cosmological parameters with the PBH abundance
$\\{\omega_{b},\omega_{d},100\theta_{s},\ln(10^{10}A_{s}),n_{s},z_{\rm
reio},f_{\rm PBH}\\}$.

Fig 6 shows 2-D contour plot for $z_{\rm reio}$ and $f_{\rm PBH}$. The
annihilation enhancement occurs when the fluctuations grow enough to produce
halos and induces the abrupt early reionization. PBH isocurvature fluctuations
lead to the early structure formation and significantly enhance such effects.
The MCMC analysis is sensitive to the peak height of the acoustic part in the
temperature anisotropy, which is proportional to $A_{s}e^{-2\tau}$ where
$\tau$ is the optical depth to the last scattering surface. This optical depth
$\tau$ is a sum of contributions from the epoch soon after the last scattering
surface ($200\lesssim z\lesssim 1100$), the early reionization epoch
($20\lesssim z\lesssim 200$) and the late reionization epoch ($z\lesssim 20$)
(this late reionization is parameterized by ”tanh”-shape reionization model
centered at $z=z_{\rm reio}$). The DM annihilation in halo structures largely
affect the ionization fraction in the early reionization part by the abrupt
increase of the reionization fraction even if $f_{\rm PBH}=0$ (because of the
conventional adiabatic perturbations). Not to affect the total $\tau$, a
smaller $z_{\rm reio}$ is preferable. We also note the best fit value of
$A_{s}$ does not change dramatically in the presence of PBH in our analysis
partly because the uncertainty in the polarization data is larger than that of
the temperature anisotropy. Hence the MCMC tends to vary $z_{\rm reio}$ rather
than $A_{s}$ to keep the amplitude of $A_{s}e^{-2\tau}$. In other words, to
maintain the peak height in the CMB temperature anisotropy, the smaller
contribution of the conventional tanh reionization model is preferable to
compensate the additional contribution from WIMP annihilation reionization.
This hence leads to the smaller $z_{\rm reio}$ compared with $z_{\rm
reio}=7.67\pm 0.73$ in the standard $\Lambda$CDM case Aghanim et al. (2018)
(without including DM annihilation).

Our MCMC analysis provides the constraint on the PBH abundance, $f_{\rm
PBH}<10^{-2.66}$ for $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times
10^{-28}\rm cm^{3}/s/GeV$, at the 95% confidence level. For the DM
annihilating into the SM particle pairs, a typical range of $f_{\rm ann}$ is
of order $f_{\rm ann}\approx 0.2-0.6$ Slatyer (2016) and we can obtain the
constraint

$\displaystyle f_{\rm PBH}<10^{-2.66}\quad{\rm for~{}}\langle\sigma
v\rangle/m_{\chi}\sim 10^{-27}\rm cm^{3}/s/GeV,$ (28)

at the 95% confidence level. As mentioned in the last section, the survived
minihalo contribution is negligible when $\log_{10}f_{\rm PBH}<-2.6$.
Therefore, this constraint given by Eq. (28) does not contain the model
uncertainty related to the abundance of survived minihalos. On the other hand,
the constraint strongly depends on the annihilation rate $f_{\rm
ann}\langle\sigma v\rangle/m_{\chi}$. When we take $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}=3\times 10^{-29}\rm cm^{3}/s/GeV$, the constraint is relaxed
to $f_{\rm PBH}<10^{-2.38}$ at the 95% confidence level.

Fig. 7 is the summary of our constraint on $f_{\rm PBH}$. The colored region
represents the excluded parameter region for each $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}$ by our analysis. As mentioned in the previous section, when
$f_{\rm PBH}\gtrsim 10^{-2.6}$, the contribution from the survived minihalos,
the first term in the right-hand side of Eq. (15), becomes significant. The
minihalo contribution becomes too big for $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}\gtrsim 10^{-28}\rm cm^{3}/s/GeV$ if $f_{\rm PBH}\gtrsim
10^{-2.6}$, and the excluded region for $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}\gtrsim 10^{-28}\rm cm^{3}/s/GeV$ represents the parameter
sets where the survived minihalo contribution cannot be ignored compared with
the smooth background contribution. For $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}\lesssim 10^{-28}\rm cm^{3}/s/GeV$, there is room for the
minihalo contribution to dominate the smooth background contribution, and the
slope of the constraint becomes smaller compared with that for a higher
$f_{\rm ann}\langle\sigma v\rangle/m_{\chi}$ for which there is little room
for the minihalo contribution.

Figure 5: Angular power spectrum of CMB E-mode polarization in the mixed DM
scenario of PBHs and self-annihilating DM. From the top to the bottom, the
solid lines are the power spectra for $f_{\rm PBH}=10^{-2.4}$ and $f_{\rm
PBH}=10^{-2.5}$. For reference, we show the angular power spectrum without the
DM annihilation in the dashed line. Figure 6: Two-dimensional contour plot for
$z_{\rm reio}$ and $f_{\rm PBH}$ from the MCMC analysis with Planck 2018 data.
Here we set $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times 10^{-28}\rm
cm^{3}/s/GeV$.  Figure 7: The constraints on $f_{\rm PBH}$ for different
values of $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}$. The colored region is
excluded by the CMB anisotropy measurement.

###  V.2 CMB $y$-type distortion

The CMB distortion can offer another powerful probe on the thermal history of
the Universe Tashiro (2014); Chluba and Jeong (2014); Chluba (2014). The
injected energy from the DM annihilation creates the deviation from the black-
body spectrum in the CMB energy spectrum. In particular, the enhancement of
the structure formation due to the PBHs can result in the enhanced $y$-type
distortion. The CMB $y$-type distortion can be estimated by Chluba (2016)

$y\approx\frac{1}{4}\int^{z_{\mu y}}_{z_{0}}\frac{dz}{(1+z)}\frac{\dot{Q}^{\rm
heat}}{H(z)\rho_{\gamma}(z)},$ (29)

where $\rho_{\gamma}(z)$ is the CMB energy density at $z$ and we set $z_{\mu
y}=5\times 10^{4}$. Before $z_{\mu y}$, $\mu$-type distortion is created
instead of $y$-type distortion. We take $z_{0}=200$ because the energy
transfer from baryons to CMB becomes inefficient after $z_{0}$ Chluba (2016).
Fig. 8 shows the created $y$-type distortion as a function of the PBH
fraction, $f_{\rm PBH}$. The blue line represents our fiducial case with
$f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times 10^{-28}\rm cm^{3}/s/GeV$.
As $f_{\rm PBH}$ increases, the larger PBH isocurvature fluctuations can lead
to the larger $y$-type distortion. However, when $f_{\rm PBH}>0.1$, $y$-type
distortion decreases with an increasing $f_{\rm PBH}$ because the amount of
annihilating DM is not sufficient to generate the large CMB distortion (we
remind readers that $\Omega_{\rm DM}=\Omega_{\rm PBH}+\Omega_{\chi}$).

One of the model parameters, which strongly affects $y$-type distortion, is
$z_{\rm max}$ which represents the maximum redshift for halo formations.
Although we conservatively set $z_{\rm max}=2000$ in the fiducial model, halos
can form when the density fluctuations grow sufficiently after the epoch of
the matter-radiation equality. In order to evaluate the impact of $z_{\rm
max}$ on $y$-type distortion, we calculate $y$-type distortion with $z_{\rm
max}=3000$ and plot the result in the thin orange line in Fig. 8. When $f_{\rm
PBH}$ is large, $y$-type distortion strongly depends on $z_{\rm max}$. A
larger $f_{\rm PBH}$ induces an earlier formation of minihalos whose DM
density $(\propto(1+z_{\rm f})^{3})$ can be big, and a larger $y$-type
distortion can be created from the DM annihilation with a higher $z_{\rm
max}$. However, when $\log_{10}f_{\rm PBH}<-1.5$, the halo formation at
$z>2000$ is not so efficient and the effect on $y$-type distortion is small
even if $z_{\rm max}=2000$ is changed to $z_{\rm max}=3000$ in our analysis.
The current limit of $y$-type distortion is given by COBE/FIRAS, $y<1.5\times
10^{-5}$ Fixsen et al. (1996), while the projected 1-$\sigma$ detection
sensitivity for the next-generation PIXIE-like experiment Kogut et al. (2011)
is $\sim 3.4\times 10^{-9}$ Chluba et al. (2019). As discussed for the CMB
anisotropy, $f_{\rm PBH}$ is constrained as $f_{\rm PBH}<10^{-2.66}$ from
Planck 2018 data. The PBH abundance which can create the PIXIE detection level
of $y$-type distortion is already ruled out by Planck observation. This is
partly because the spectral distortion is affected by the DM annihilation
energy going into the abundant CMB photons while the CMB anisotropy bounds
mainly come from the ionization fraction which concerns the DM annihilation
energy going into the gas. It is hence reasonable to expect the latter can
give the tighter bounds on $f_{\rm PBH}$ considering the much bigger number
density of photons than baryons, in addition to the current poor precision on
the CMB spectral distortion measurements compared with the CMB anisotropy
data.

Figure 8: CMB $y$-type distortion as a function of $f_{\rm PBH}$. Our fiducial
case is represented in a thick blue line with $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}=3\times 10^{-28}~{}\rm cm^{3}/s/GeV$ and $z_{\rm max}=2000$.
To represent the dependence on $z_{\rm max}$, we plot the result with $z_{\rm
max}=3000$ in a thin orange line.

###  V.3 Global 21-cm signal

Redshifted 21-cm signals can probe the IGM evolution of the thermal and
ionization history. The multifrequency radio observation can measure the 21-cm
signals from different redshifts corresponding to different observation
frequencies. The measured quantity in 21-cm observations is the so-called
differential brightness temperature, which is the deviation of the redshifted
21-cm brightness temperature from the CMB temperature. Similarly to CMB, both
global (sky-averaged) signal and spatial fluctuations of differential
brightness temperature can provide valuable information about our Universe
(for a comprehensive review, see Ref. Furlanetto et al. (2006)). Here we focus
on the global signal.

The global differential brightness temperature from a redshift $z$ is provided
by Madau et al. (1997); Pritchard and Loeb (2010)

$\delta
T_{b}(z)=\frac{3}{32\pi}\frac{hc^{3}A_{10}}{k_{B}\nu_{0}^{2}}\frac{x_{\rm
HI}n_{H}}{(1+z)^{2}(dv_{||}/dr_{||})}\left(1-\frac{T_{\gamma}}{T_{S}}\right),$
(30)

where $dv_{||}/dr_{||}$ is the gradient of the proper velocity along the line
of sight, which is the Hubble expansion contribution in the global signal.
$T_{S}$ is the spin temperature of the neutral hydrogen hyperfine structure.
The spin temperature is determined by the balance in the processes of the
hyperfine excitation and de-excitation,

$T_{S}=\frac{T_{\gamma}+y_{\rm kin}T_{k}}{1+y_{\rm kin}},$ (31)

where $y_{\rm kin}$ represents the efficiency ratio between the absorption of
CMB photons and the thermal collisions in the hyperfine transition. We adopt
the approximated analytical form of $y_{\rm kin}$ in Ref. (Kuhlen et al.,
2006). In Eq. (30), we do not include the contribution from Ly-$\alpha$
coupling, because it is not important until the beginning of the first star
formation. In the left panel of Fig. 9, we plot the evolution of spin and
baryon temperatures with $f_{\rm PBH}=10^{-3}$ in solid lines. For comparison,
we also plot the evolutions in the standard cosmology case (”no annihilation”
case). Here we do not consider other heating sources including stars and
galaxies. For $f_{\rm PBH}=10^{-3}$, the annihilation from DM halos becomes
effective below $z\sim 100$ when baryons are heated up and the thermal
evolution of baryons deviate from the adiabatic evolution,
$T_{k}\propto(1+z)^{2}$, in ”no annihilation” case. Due to heated baryons, the
spin temperature also becomes large with $f_{\rm PBH}=10^{-3}$, compared with
”no annihilation” case. The spin temperature lies between the CMB and baryon
temperatures. Therefore, after the baryon temperature becomes larger than the
CMB temperature, the spin temperature also becomes higher than the CMB
temperature for $f_{\rm PBH}=10^{-3}$. On the other hand, the spin temperature
never exceeds the CMB temperature in ”no annihilation” case.

The right panel of Fig. 9 shows the global differential brightness temperature
as a function of a redshift. The sign of the differential brightness
temperature depends on the spin temperature. When the spin temperature is
smaller than the CMB temperature, the differential brightness temperature is
negative and observed as absorption signals on the CMB frequency spectrum. As
$f_{\rm PBH}$ decreases, the effects of the DM annihilation becomes small and
the redshift, $z_{\rm tr}$, at which the global differential brightness
temperature shifts from the absorption to the emission signal becomes smaller.
In the case of ”no annihilation”, the differential brightness temperature
cannot turn to the positive side because of the lack of heating sources and
the spin temperature stays lower than the CMB temperature. Therefore,
identifying $z_{\rm tr}$ observationally can provide the constraint on the
heating source including the DM annihilation.

---
Figure 9: Left panel: The spin and baryon temperatures as functions of
redshifts. The solid blue and orange lines represent the spin and baryon
temperatures for $f_{\rm PBH}=10^{-3.0}$. For comparison, we plot their
evolutions in the no annihilation DM case in the dashed lines. Right panel:
The evolution of the global differential brightness temperature with different
$f_{\rm PBH}$. From top to bottom, the solid lines are for $f_{\rm
PBH}=10^{-3.0}$, $f_{\rm PBH}=10^{-3.2}$ and $f_{\rm PBH}=10^{-3.4}$. For
reference, the dashed line shows the evolution in the case of no annihilation
DM.

In Fig. 10, we show the relation between $z_{\rm tr}$ and $f_{\rm PBH}$. A
small $f_{\rm PBH}$ provides a low $z_{\rm tr}$. Recently the EDGES reported
the detection of the global 21-cm absorption in the redshift range between
$z\sim 21$ and $z\sim 15$ Bowman et al. (2018). This result suggests that the
baryon temperature is lower than the CMB temperature until $z\sim 15$ and can
limit the heating source causing the early reionization. As shown in Fig. 10,
$f_{\rm PBH}>10^{-3.4}$ leads to the emission signal above $z\sim 15$. The
EDGES absorption signal at $z\sim 15$ hence requires $f_{\rm PBH}<10^{-3.4}$
with $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times 10^{-28}{\rm
cm^{3}/s/GeV}$. Therefore, for $f_{\rm ann}\sim{\cal O}(0.1)$ (a typical range
for the SM particle annihilation final states Slatyer (2016)), our constraint
is

$f_{\rm PBH}<10^{-3.4}\quad{\rm with~{}}\langle\sigma v\rangle/m_{\chi}\sim
10^{-27}{\rm cm^{3}/s/GeV}.$ (32)

As mentioned before, when $f_{\rm PBH}<10^{-2.6}$, the contribution from
survived minihalos is negligible. Therefore, this constraint given by Eq. (32)
does not include the substructure contributions and hence does not suffer from
the model parameter uncertainties related to $z_{\rm max}$ and $z_{\rm f}$. On
the other hand, $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}$ strongly affects
the constraint. If $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times
10^{-29}{\rm cm^{3}/s/GeV}$, the EDGES provides the constraint $f_{\rm
PBH}\lesssim 10^{-2.5}$. We also note that the baryon temperature larger than
the CMB temperature can turn the absorption into the emission to result in the
tight bound on $f_{\rm PBH}$, but it is not high enough to sufficiently ionize
the IGM. This is a reason why the bounds on $f_{\rm PBH}$ from the global 21cm
signals can become tighter than those from the CMB.

Figure 10: The dependence of the transition redshift, $z_{\rm tr}$, from the
absorption to emission on the PBH fraction $f_{\rm PBH}$.

##  VI conclusion

In this paper, we studied the bounds on the PBH abundance in the mixed DM
scenarios consisting of the self-annihilating DM and PBHs. The Poisson
distribution of PBHs can lead to the isocurvature perturbations and
consequently to the early structure formation. The enhancement of the DM
annihilation occurs in those early formed dense halos and modifies the
ionization and temperature evolutions of baryons. Such modifications affect
early Universe observables such as the CMB and 21cm signals. In order to
obtain the constraint on the PBH abundance from CMB observations, we performed
the MCMC analysis with the latest Planck data. The obtained constraint is
$f_{\rm PBH}<10^{-2.66}$ for $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}=3\times 10^{-28}\rm cm^{3}/s/GeV$. The constraint depends on
the annihilation rate and, for instance, a weaker bound $f_{\rm
PBH}<10^{-2.38}$ arises for a smaller $f_{\rm ann}\langle\sigma
v\rangle/m_{\chi}=3\times 10^{-29}\rm cm^{3}/s/GeV$. The energy injection from
DM annihilation can also create the CMB spectral distortion, and we found that
the PBH abundance which can generate the CMB distortion observable by the
PIXIE-like future observation has been already ruled out by the Planck data.
The redshifted 21-cm observation is also a promising probe on the heating
sources such as the DM annihilation in the dark ages. The observational
determination of the redshift at which the global 21-cm signal shifts from the
absorption to the emission can provide the constraint on the efficiency of the
heating source in the dark ages because, at this transition redshift, the
heated baryon temperature becomes larger than the CMB temperature. We
discussed how this transition redshift is affected in the presence of the PBH
isocurvature perturbations. Recent EDGES observation reported that they
detected the absorption signals of global redshifted 21-cm lines between
$15\lesssim z\lesssim 22$. While the justification of their large absorption
signal amplitude is under an active debate, we solely focused on the redshift
dependence of the signal and obtained the PBH abundance constraint $f_{\rm
PBH}<10^{-3.4}$ for $f_{\rm ann}\langle\sigma v\rangle/m_{\chi}=3\times
10^{-28}{\rm cm^{3}/s/GeV}$. We also argued that our bounds on $f_{\rm PBH}$
are insensitive to $M_{\rm PBH}$ for the parameter range of interest in our
analysis ($M_{\rm PBH}\gtrsim 10^{-6}M_{\odot}$ when the minimum dark matter
halo mass $M_{\rm min}=10^{-6}M_{\odot}$).

Our Poisson effect constraints would be of particular interest for the light
(sub-GeV) WIMP which has been less explored compared with the heavier ($>1$
GeV) WIMP in the presence of PBHs. For instance, the CMB can give the bounds
$f_{\rm PBH}\lesssim{\cal O}(10^{-3})$ for $m_{\chi}=1$ MeV and $f_{\rm
ann}\langle\sigma v\rangle=3\times 10^{-31}\rm cm^{3}/s$. Our bounds are
independent from and complementary to the seed effects which consider the DM
accretion into individual PBHs. The DM accretion into PBHs can form the steep
profile $\rho(r)\propto r^{-9/4}$ when the DM kinetic energy is negligible
compared with their potential energy, and the consequent enhanced DM
annihilation can lead to the tight bounds on PBHs using the observation data
such as the Fermi gamma ray and Planck CMB data Boucenna et al. (2018); Adamek
et al. (2019); Eroshenko (2016); Carr et al. (2020b); Cai et al. (2020); Delos
et al. (2018); Kohri et al. (2014); Bertone et al. (2019); Carr and Silk
(2018); Carr and Rees (1984); Hertzberg et al. (2020); Kashlinsky (2016);
Tashiro and Kadota (2021). Even though the detailed numerical analysis has not
been performed yet for such seed effects when the dark matter is light such
that the kinetic energy cannot be ignored in estimating the dark matter
profile around a PBH, the analytical estimation indicates the bounds on
$f_{\rm PBH}$ would be significantly weakened for small $m_{\chi}$ and $M_{\rm
PBH}$ due to a less steep profile around a PBH Eroshenko (2020); Boucenna et
al. (2018); Carr et al. (2020b). Ref. Carr et al. (2020b), for instance,
analytically estimated the upper bound on $f_{\rm PBH}$ assuming the Maxwell-
Boltzmann distribution for the DM velocity and, for example, $f_{\rm
PBH}\lesssim 0.1$ for $m_{\chi}\sim 1$ MeV and $M_{\rm PBH}\sim 0.01M_{\odot}$
(the upper bound of $f_{\rm PBH}$ scales as $\propto m_{\chi}^{-3.7}M_{\rm
PBH}^{-1.5}$ assuming the DM kinetic decoupling dependence on the DM mass as
$T_{\rm KD}\propto m_{\chi}^{5/4}$ typical for a bino-like WIMP Bringmann and
Hofmann (2007)). The bounds on $f_{\rm PBH}$ when the dark matter kinetic
energy is not negligible would require a more detailed numerical analysis, and
such numerical studies which also should take account of a DM model dependence
such as the nature of DM kinetic decoupling Loeb and Zaldarriaga (2005);
Bertschinger (2006); Gondolo et al. (2012); Profumo et al. (2006); Gondolo and
Kadota (2016); Green et al. (2004, 2005); Bringmann and Hofmann (2007) would
further quantitatively clarify the complementary between the Poisson and seed
effects of PBHs in the presence of the WIMPs. The bounds on the PBH parameters
from DM annihilation when a halo includes multiple PBHs (which is likely when
$f_{\rm PBH}>10^{-4}$ Inman and Ali-Haïmoud (2019)) and when a PBH interacts
with another (e.g. forming the binary) would be also worth pursuing.

This work was supported by the Institute for Basic Science (IBS-R018-D1) and
Grants-in-Aid for Scientific Research from JSPS (21K03533). KK thanks the
Kobayashi-Maskawa Institute at Nagoya University for the hospitality.

## References

  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), eprint 1602.03837.
  * Carr and Kuhnel (2020) B. Carr and F. Kuhnel (2020), eprint 2006.02838.
  * Green and Kavanagh (2020) A. M. Green and B. J. Kavanagh (2020), eprint 2007.10722.
  * Carr et al. (2020a) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama (2020a), eprint 2002.12778.
  * Afshordi et al. (2003) N. Afshordi, P. McDonald, and D. Spergel, Astrophys. J. Lett. 594, L71 (2003), eprint astro-ph/0302035.
  * Ali-Haïmoud (2018) Y. Ali-Haïmoud, Phys. Rev. Lett. 121, 081304 (2018), eprint 1805.05912.
  * Inman and Ali-Haïmoud (2019) D. Inman and Y. Ali-Haïmoud, Phys. Rev. D 100, 083528 (2019), eprint 1907.08129.
  * Desjacques and Riotto (2018) V. Desjacques and A. Riotto, Phys. Rev. D 98, 123533 (2018), eprint 1806.10414.
  * Oguri and Takahashi (2020) M. Oguri and R. Takahashi (2020), eprint 2007.01936.
  * Kadota and Silk (2021) K. Kadota and J. Silk, Phys. Rev. D 103, 043530 (2021), eprint 2012.03698.
  * Gong and Kitajima (2017) J.-O. Gong and N. Kitajima, JCAP 08, 017 (2017), eprint 1704.04132.
  * Mena et al. (2019) O. Mena, S. Palomares-Ruiz, P. Villanueva-Domingo, and S. J. Witte, Phys. Rev. D 100, 043540 (2019), eprint 1906.07735.
  * Zeldovich and Sunyaev (1969) Y. B. Zeldovich and R. A. Sunyaev, Astrophysics & Space Science 4, 301 (1969).
  * Sunyaev and Zeldovich (1970) R. A. Sunyaev and Y. B. Zeldovich, Astrophysics & Space Science 7, 20 (1970).
  * Burigana et al. (1991) C. Burigana, L. Danese, and G. de Zotti, Astronomy & Astrophysics 246, 49 (1991).
  * Hu and Silk (1993) W. Hu and J. Silk, Phys. Rev. Lett. 70, 2661 (1993).
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A5 (2020), eprint 1907.12875.
  * Madau et al. (1997) P. Madau, A. Meiksin, and M. J. Rees, Astrophys. J. 475, 429 (1997), eprint astro-ph/9608010.
  * Shaver et al. (1999) P. A. Shaver, R. A. Windhorst, P. Madau, and A. G. de Bruyn, Astronomy & Astrophysics 345, 380 (1999), eprint astro-ph/9901320.
  * Pritchard and Loeb (2010) J. R. Pritchard and A. Loeb, Phys. Rev. D 82, 023006 (2010), eprint 1005.4057.
  * Bowman et al. (2018) J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen, and N. Mahesh, Nature (London) 555, 67 (2018), eprint 1810.05912.
  * Clark et al. (2018) S. J. Clark, B. Dutta, Y. Gao, Y.-Z. Ma, and L. E. Strigari, Phys. Rev. D 98, 043006 (2018), eprint 1803.09390.
  * Hektor et al. (2018) A. Hektor, G. Hütsi, L. Marzola, M. Raidal, V. Vaskonen, and H. Veermäe, Phys. Rev. D 98, 023503 (2018), eprint 1803.09697.
  * Mitridate and Podo (2018) A. Mitridate and A. Podo, JCAP 5, 069 (2018), eprint 1803.11169.
  * Cheng et al. (2018) H.-C. Cheng, L. Li, and R. Zheng, Journal of High Energy Physics 9, 98 (2018), eprint 1805.12139.
  * D’Amico et al. (2018) G. D’Amico, P. Panci, and A. Strumia, Phys. Rev. Lett. 121, 011103 (2018), eprint 1803.03629.
  * Safarzadeh et al. (2018) M. Safarzadeh, E. Scannapieco, and A. Babul, ApJL 859, L18 (2018), eprint 1803.08039.
  * Minoda et al. (2019) T. Minoda, H. Tashiro, and T. Takahashi, Mon. Not. Roy. Astron. Soc. 488, 2001 (2019), eprint 1812.00730.
  * Tashiro et al. (2014) H. Tashiro, K. Kadota, and J. Silk, Phys. Rev. D 90, 083522 (2014), eprint 1408.2571.
  * Barkana (2018) R. Barkana, Nature 555, 71 (2018), eprint 1803.06698.
  * Slatyer (2016) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016), eprint 1506.03811.
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999), eprint astro-ph/9906391.
  * Lacki and Beacom (2010) B. C. Lacki and J. F. Beacom, Astrophys. J. Lett. 720, L67 (2010), eprint 1003.3466.
  * Boucenna et al. (2018) S. M. Boucenna, F. Kuhnel, T. Ohlsson, and L. Visinelli, JCAP 07, 003 (2018), eprint 1712.06383.
  * Adamek et al. (2019) J. Adamek, C. T. Byrnes, M. Gosenca, and S. Hotchkiss, Phys. Rev. D 100, 023506 (2019), eprint 1901.08528.
  * Eroshenko (2016) Y. Eroshenko, Astron. Lett. 42, 347 (2016), eprint 1607.00612.
  * Carr et al. (2020b) B. Carr, F. Kuhnel, and L. Visinelli (2020b), eprint 2011.01930.
  * Cai et al. (2020) R.-G. Cai, X.-Y. Yang, and Y.-F. Zhou (2020), eprint 2007.11804.
  * Delos et al. (2018) M. S. Delos, A. L. Erickcek, A. P. Bailey, and M. A. Alvarez, Phys. Rev. D 98, 063527 (2018), eprint 1806.07389.
  * Kohri et al. (2014) K. Kohri, T. Nakama, and T. Suyama, Phys. Rev. D 90, 083514 (2014), eprint 1405.5999.
  * Bertone et al. (2019) G. Bertone, A. M. Coogan, D. Gaggero, B. J. Kavanagh, and C. Weniger, Phys. Rev. D 100, 123013 (2019), eprint 1905.01238.
  * Carr and Silk (2018) B. Carr and J. Silk, Mon. Not. Roy. Astron. Soc. 478, 3756 (2018), eprint 1801.00672.
  * Carr and Rees (1984) B. J. Carr and M. J. Rees, Monthly Notices of the RAS 206, 801 (1984).
  * Hertzberg et al. (2020) M. P. Hertzberg, S. Nurmi, E. D. Schiappacasse, and T. T. Yanagida (2020), eprint 2011.05922.
  * Kashlinsky (2016) A. Kashlinsky, Astrophys. J. Lett. 823, L25 (2016), eprint 1605.04023.
  * Tashiro and Kadota (2021) H. Tashiro and K. Kadota (2021), eprint 2104.09738.
  * Peacock (1999) J. A. Peacock, _Cosmological Physics_ (1999).
  * Bertschinger (2006) E. Bertschinger, Phys. Rev. D 74, 063509 (2006), eprint astro-ph/0607319.
  * Loeb and Zaldarriaga (2005) A. Loeb and M. Zaldarriaga, Phys. Rev. D 71, 103520 (2005), eprint astro-ph/0504112.
  * Gondolo et al. (2012) P. Gondolo, J. Hisano, and K. Kadota, Phys. Rev. D 86, 083523 (2012), eprint 1205.1914.
  * Profumo et al. (2006) S. Profumo, K. Sigurdson, and M. Kamionkowski, Phys. Rev. Lett. 97, 031301 (2006), eprint astro-ph/0603373.
  * Gondolo and Kadota (2016) P. Gondolo and K. Kadota, JCAP 06, 012 (2016), eprint 1603.05783.
  * Green et al. (2004) A. M. Green, S. Hofmann, and D. J. Schwarz, Mon. Not. Roy. Astron. Soc. 353, L23 (2004), eprint astro-ph/0309621.
  * Green et al. (2005) A. M. Green, S. Hofmann, and D. J. Schwarz, JCAP 08, 003 (2005), eprint astro-ph/0503387.
  * Bringmann and Hofmann (2007) T. Bringmann and S. Hofmann, JCAP 04, 016 (2007), [Erratum: JCAP 03, E02 (2016)], eprint hep-ph/0612238.
  * Diemand et al. (2005) J. Diemand, B. Moore, and J. Stadel, Nature 433, 389 (2005), eprint astro-ph/0501589.
  * Cooray and Sheth (2002) A. Cooray and R. K. Sheth, Phys. Rept. 372, 1 (2002), eprint astro-ph/0206508.
  * Poulin et al. (2015) V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP 12, 041 (2015), eprint 1508.01370.
  * Galli et al. (2013) S. Galli, T. R. Slatyer, M. Valdes, and F. Iocco, Phys. Rev. D 88, 063502 (2013), eprint 1306.0563.
  * Giesen et al. (2012) G. Giesen, J. Lesgourgues, B. Audren, and Y. Ali-Haimoud, JCAP 12, 008 (2012), eprint 1209.0247.
  * Wang et al. (2020) J. Wang, S. Bose, C. S. Frenk, L. Gao, A. Jenkins, V. Springel, and S. D. White, Nature 585, 39 (2020), eprint 1911.09720.
  * Ishiyama et al. (2020) T. Ishiyama et al. (2020), eprint 2007.14720.
  * Sánchez-Conde and Prada (2014) M. A. Sánchez-Conde and F. Prada, Mon. Not. Roy. Astron. Soc. 442, 2271 (2014), eprint 1312.1729.
  * Moline et al. (2017) A. Moline, M. A. Sanchez-Conde, S. Palomares-Ruiz, and F. Prada, Mon. Not. Roy. Astron. Soc. 466, 4974 (2017), eprint 1603.04057.
  * Zhao et al. (2009) D. Zhao, Y. Jing, H. Mo, and G. Boerner, Astrophys. J. 707, 354 (2009), eprint 0811.0828.
  * Prada et al. (2012) F. Prada, A. A. Klypin, A. J. Cuesta, J. E. Betancort-Rijo, and J. Primack, Mon. Not. Roy. Astron. Soc. 423, 3018 (2012), eprint 1104.5130.
  * Ackermann et al. (2015) M. Ackermann et al. (Fermi-LAT), JCAP 09, 008 (2015), eprint 1501.05464.
  * Blinov et al. (2020) N. Blinov, M. J. Dolan, and P. Draper, Phys. Rev. D 101, 035002 (2020), eprint 1911.07853.
  * Zhao et al. (2007) H. Zhao, J. E. Taylor, J. Silk, and D. Hooper, Astrophys. J. 654, 697 (2007), eprint astro-ph/0508215.
  * van den Bosch et al. (2018) F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, Mon. Not. Roy. Astron. Soc. 474, 3043 (2018), eprint 1711.05276.
  * Green and van den Bosch (2019) S. B. Green and F. C. van den Bosch, Mon. Not. Roy. Astron. Soc. 490, 2091 (2019), eprint 1908.08537.
  * Berezinsky et al. (2003) V. Berezinsky, V. Dokuchaev, and Y. Eroshenko, Phys. Rev. D 68, 103003 (2003), eprint astro-ph/0301551.
  * Arvanitaki et al. (2020) A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, Phys. Rev. D 101, 083014 (2020), eprint 1909.11665.
  * Berezinsky et al. (2006) V. Berezinsky, V. Dokuchaev, and Y. Eroshenko, Phys. Rev. D 73, 063504 (2006), eprint astro-ph/0511494.
  * Dai and Miralda-Escudé (2020) L. Dai and J. Miralda-Escudé, Astron. J. 159, 49 (2020), eprint 1908.01773.
  * Xiao et al. (2021) H. Xiao, I. Williams, and M. McQuinn (2021), eprint 2101.04177.
  * Peebles (1968) P. J. E. Peebles, Astrophys. J. 153, 1 (1968).
  * Zeldovich et al. (1968) Y. B. Zeldovich, V. G. Kurt, and R. A. Syunyaev, Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki 55, 278 (1968).
  * Padmanabhan and Finkbeiner (2005) N. Padmanabhan and D. P. Finkbeiner, Phys. Rev. D 72, 023508 (2005), eprint astro-ph/0503486.
  * Poulin et al. (2017) V. Poulin, J. Lesgourgues, and P. D. Serpico, JCAP 2017, 043 (2017), eprint 1610.10051.
  * Ali-Haïmoud and Hirata (2011) Y. Ali-Haïmoud and C. M. Hirata, Phys. Rev. D 83, 043513 (2011), eprint 1011.3758.
  * Blas et al. (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 2011, 034 (2011), eprint 1104.2933.
  * Aghanim et al. (2018) N. Aghanim et al. (Planck) (2018), eprint 1807.06209.
  * Audren et al. (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), eprint 1210.7183.
  * Tashiro (2014) H. Tashiro, Progress of Theoretical and Experimental Physics 2014, 06B107 (2014).
  * Chluba and Jeong (2014) J. Chluba and D. Jeong, Monthly Notices of the RAS 438, 2065 (2014), eprint 1306.5751.
  * Chluba (2014) J. Chluba, arXiv e-prints arXiv:1405.6938 (2014), eprint 1405.6938.
  * Chluba (2016) J. Chluba, Monthly Notices of the RAS 460, 227 (2016), eprint 1603.02496.
  * Fixsen et al. (1996) D. Fixsen, E. Cheng, J. Gales, J. C. Mather, R. Shafer, and E. Wright, Astrophys. J. 473, 576 (1996), eprint astro-ph/9605054.
  * Kogut et al. (2011) A. Kogut, D. J. Fixsen, D. T. Chuss, J. Dotson, E. Dwek, M. Halpern, G. F. Hinshaw, S. M. Meyer, S. H. Moseley, M. D. Seiffert, et al., JCAP 2011, 025 (2011), eprint 1105.2044.
  * Chluba et al. (2019) J. Chluba, M. H. Abitbol, N. Aghanim, Y. Ali-Haimoud, M. Alvarez, K. Basu, B. Bolliet, C. Burigana, P. de Bernardis, J. Delabrouille, et al., arXiv e-prints arXiv:1909.01593 (2019), eprint 1909.01593.
  * Furlanetto et al. (2006) S. R. Furlanetto, S. P. Oh, and F. H. Briggs, Phys. Rept. 433, 181 (2006), eprint astro-ph/0608032.
  * Kuhlen et al. (2006) M. Kuhlen, P. Madau, and R. Montgomery, ApJL 637, L1 (2006), eprint astro-ph/0510814.
  * Eroshenko (2020) Y. Eroshenko, Int. J. Mod. Phys. A 35, 2040046 (2020), eprint 1910.01564.