# Revisiting CMB constraints on dark matter annihilation

Masahiro Kawasaki  ICRR, University of Tokyo, Kashiwa, 277-8582, Japan  Kavli
IPMU (WPI), UTIAS, University of Tokyo, Kashiwa, 277-8583, Japan  Hiromasa
Nakatsuka  ICRR, University of Tokyo, Kashiwa, 277-8582, Japan  Kazunori
Nakayama  Department of Physics, Faculty of Science, The University of Tokyo,
Bunkyo-ku, Tokyo 113-0033, Japan  Toyokazu Sekiguchi  Theory Center, IPNS,
KEK, Tsukuba, Ibaraki 305-0801, Japan

###### Abstract

The precision measurements of the cosmic microwave background power spectrum
put a strong constraint on the dark matter annihilation cross section since
the electromagnetic energy injection by the dark matter annihilation affects
the ionization history of the universe. In this paper, we update our previous
simulation code for calculating the ionization history with the effect of dark
matter annihilation by including Helium interactions and improving the
precision of calculations. We give an updated constraint on the annihilation
cross section and mass of dark matter using the modified RECFAST code with the
Planck 2018 datasets.

###### Contents

  1. I Introduction
  2. II Radiative processes
    1. II.1 Electron Interactions
    2. II.2 Photon Interactions
    3. II.3 Treatment of positron
  3. III Energy injection as shower process
  4. IV Effect of energy injection on CMB
  5. V Constraints on the Dark matter annihilation
  6. VI Conclusion
  7. A Numerical scheme of shower code

##  I Introduction

The existence of dark matter (DM), contributing to $24$% of the present energy
density of the universe, is one of the most important mysteries in cosmology
and particle physics. DM is likely composed of particles beyond the standard
model of particle physics, and many candidates have been proposed so far.
Among them, WIMPs (weakly interacting particles) have been most intensively
studied and searched for directly and indirectly.

The so-called indirect searches look for the DM annihilation signals. Several
observations can put constraints on the DM annihilation cross section
$\braket{\sigma v}$: CMB (Cosmic microwave background) Padmanabhan and
Finkbeiner (2005); Mapelli  _et al._ (2006); Zhang  _et al._ (2006); Kanzaki
and Kawasaki (2008); Kanzaki  _et al._ (2010); Slatyer  _et al._ (2009);
Slatyer (2013); Galli  _et al._ (2013); Slatyer (2016a, b); Liu  _et al._
(2016); Acharya and Khatri (2019a); Cang  _et al._ (2020), $\gamma$-rays
Ackermann  _et al._ (2015); Ahnen  _et al._ (2016), cosmic rays Aguilar  _et
al._ (2013); Adriani  _et al._ (2009); Aharonian  _et al._ (2009); Albert  _et
al._ (2020); Gelmini (2015) and light elements synthesized in BBN (Big Bang
Nucleosynthesis) Reno and Seckel (1988); Jedamzik (2004); Hisano  _et al._
(2009, 2011); Kawasaki  _et al._ (2015); Depta  _et al._ (2019). For the DM
mass about $10$–$1000$ GeV, the $\gamma$-ray observations of dwarf galaxies
may give the most stringent constraint Acciari  _et al._ (2010); Abramowski
_et al._ (2014); Ackermann  _et al._ (2015); Ahnen  _et al._ (2016). However,
the constraint may contain sizable astronomical uncertainties because it is
sensitive to the DM profiles in dwarf galaxies which are difficult to
determine precisely. On the other hand, CMB can give a robust constraint,
which does not suffer from such astronomical uncertainties, since the dominant
effect of DM annihilation on the CMB comes from the epoch around the
recombination when the DM spatial distribution is almost homogeneous. In the
early universe, the DM annihilation into the standard model particles injects
energy into the background plasma and it changes the ionization history of the
universe around/after the recombination epoch ($z<3000$ in terms of the
redshift) Adams  _et al._ (1998); Chen and Kamionkowski (2004); Poulin  _et
al._ (2017). This then affects the CMB angular power spectrum. The Planck
observation reaches the $\mathcal{O}(1)\%$ precision in determining the
various cosmological parameters, which strongly constrains the energy
injection from the DM annihilation in the early universe.

Let us focus on the CMB constraint on DM annihilation. In order to obtain a
reliable constraint, we have to estimate the ionization fractions of hydrogen
and helium taking all relevant radiative processes into account. In the
literature, the efficiency factor $f_{\text{eff}}(z)$ is often used for
parameterizing the energy injection of DM annihilation. $f_{\text{eff}}(z)$ is
the ratio of absorbed energy to the injected one at a given redshift Aghanim
_et al._ (2020a) and it depends on the interaction between the background
plasma and the high energy particles injected by the DM annihilation. The
high-energy particles lose their energies by causing electromagnetic shower,
in which injected particles produce a lot of scattered particles with energy
large enough to ionize hydrogen and helium. The electromagnetic interactions
that take place in the showers are usually faster than the cosmic expansion.
However, the cosmic plasma is transparent against photons with high energy
($10^{8}$ eV -$10^{10}$ eV at $z\sim 10^{3}$) and they lose their energy
slowly in a longer time-scale than the Hubble time Chen and Kamionkowski
(2004); Slatyer  _et al._ (2009), for which $f_{\text{eff}}(z)$ cannot
describe the energy injection correctly. Thus we need to solve the shower
process over the whole energy and time ranges, by piling up every interaction
from low to high energy and from low to high redshift Kanzaki and Kawasaki
(2008); Kanzaki  _et al._ (2010); Slatyer  _et al._ (2009). The high-energy
particles finally lose their energy by heating the plasma, ionizing hydrogen
and helium, and the cosmic expansion. The modified ionization history changes
the CMB power spectrum. Using the current CMB data, we can constrain the DM
parameter space $(m_{\text{DM}},\braket{\sigma v})$ where $m_{\text{DM}}$ is
the DM mass. We use the Markov Chain Monte Carlo methods (MCMC) to estimate
the allowed parameter space.

Particle production and scattering in electromagnetic showers have been
investigated by several works. Among them, Refs. Galli  _et al._ (2013);
Slatyer (2016a, b) separated the shower process into the low and high energy
parts, where the Monte Carlo method is used for cooling of electrons through
atomic processes with energy smaller than $3\,$keV and the recursive method is
used with simplified interactions with energy larger than $3\,$keV. In the
high energy regime, for example, the absorbed energy by photoionization is
neglected and the collisional ionization on helium is treated as that of
hydrogen. They calculated the shower process with very fine time steps
$\Delta(\log z)=10^{-3}$ and made the table of the ratio of the injected
energy to that used for various processes (e.g. ionization) with 40 log-spaced
bins for energy between 1 keV and 10 TeV and 63 log-spaced bins for redshift
from $1+z=10$ to 3000. On the other hand, in our previous works Kanzaki and
Kawasaki (2008); Kanzaki  _et al._ (2010); Kawasaki  _et al._ (2016), we
efficiently calculated the shower process over the whole energy range by
recursive method with the “discrete” and “continuous” interactions, where the
scattered particle loses larger (smaller) energy than the size of the energy
bin $\Delta E$ in discrete (continuous) interaction. We can consistently
calculate shower process in both high and low energy regions in this method
without unnecessary assumptions. Some other works Acharya and Khatri (2019b,
a) also adopted this method and calculated the constraints on the decaying
dark matter.

In this paper, we calculate the electromagnetic energy injection due to DM
annihilation and the subsequent radiative processes based on the previously
developed method Kanzaki and Kawasaki (2008); Kanzaki  _et al._ (2010), and
obtain constraints on the DM annihilation cross section from the CMB power
spectrum. We update the simulation code Kanzaki and Kawasaki (2008); Kanzaki
_et al._ (2010); Kawasaki  _et al._ (2016) by (1) including the helium
interactions, (2) updating the Planck data sets Aghanim  _et al._ (2020b), (3)
fixing the incorrect simplification on some interactions and (4) increasing
the time steps from 50 bins to 200 bins for $z\in[0,10^{4}-1]$, which are
summarized in section II. We confirm that our simulation code reproduces the
results from the Monte Carlo computation of the electron injection into the
interstellar gas Shull and van Steenberg (1985). The injected energy increases
the ionization rate of the hydrogen and helium, and we modified the RECFAST
code Seager  _et al._ (2000, 1999); Wong  _et al._ (2008) to take account of
the effect of extra ionization. Then we can calculate the CMB power spectrum
by the CAMB code Lewis  _et al._ (2000) with the modified RECFAST code.
Finally, we search for the allowed parameter space of the annihilation cross
section and mass of DM using the CosmoMC code Lewis and Bridle (2002) with the
Planck 2018 datasets Aghanim  _et al._ (2020b).

We introduce various radiative processes related to our calculation in section
II. In section III, we show a formula which describes how the high-energy
particles lose their energy in the early universe. We then numerically
calculate the fractions of the injected energy used for heating, ionization,
and excitation. In section IV, we discuss the modification of the RECFAST to
include the effect of extra energy injection. The results of the parameter
search, using the CosmoMC code, are described in section V.

We use the following notation in this paper. $n_{\text{HI}}$ and
$n_{\text{HII}}$ describe the number densities of the neutral and ionized
hydrogen. The total number of the hydrogen atoms is $n_{\text{H}}\equiv
n_{\text{HI}}+n_{\text{HII}}$. For helium, we define $n_{\text{HeI}}$,
$n_{\text{HeII}}$, $n_{\text{HeIII}}$ as number densities of neutral,
partially ionized and fully ionized helium and $n_{\text{He}}\equiv
n_{\text{HeI}}+n_{\text{HeII}}+n_{\text{HeIII}}$. Since we focus on the
thermal history near and after the last scattering epoch, we approximate
$n_{\text{HeIII}}=0$. We denote $n_{e}$ as the free electron number density.
The ionization fractions for hydrogen and helium are given by
$x_{\text{H}}\equiv n_{\text{HII}}/n_{\text{H}}$ and $x_{\text{He}}\equiv
n_{\text{HeII}}/n_{\text{He}}$, respectively. The fraction of free electrons
is given by $x_{e}\equiv
n_{e}/n_{\text{H}}=x_{\text{H}}+f_{\text{He}}x_{\text{He}}$ with the number
fraction of helium atom to hydrogen atom
$f_{\text{He}}\equiv{n_{\text{He}}}/{n_{\text{H}}}=Y_{\rm He}/(4(1-Y_{\rm
He}))\simeq 0.081$ and mass fraction $Y_{\rm He}\simeq 0.245$.

##  II Radiative processes

In this section, we summarize newly added and updated processes in our
calculation of electromagnetic shower. The annihilation of DM injects high-
energy particles, which immediately decay into photons, electrons, positrons,
and neutrinos. The injected photon, electron and positron lose their energy
through various radiative processes. For energetic electrons, we consider the
excitation and ionization of hydrogen and helium atoms, Coulomb scattering off
the background electrons, the momentum transfer with neutral hydrogen,and
inverse Compton scattering off the background photons. For energetic photons,
we consider photoionization of the hydrogen and helium atoms, pair creation,
Compton scattering off the background electrons, photon-photon scattering, and
double-photon pair creation off the background photons. For energetic
positrons, we consider the same scattering processes as those of electrons and
the annihilation with the background electrons. The details of the above
processes are given in Kanzaki and Kawasaki (2008), and here we describe our
improvements.

In the present work, we update our previous code by newly including the effect
of helium. We add the following new interactions : the excitation and
collisional ionization of helium by high energy electrons, the photo-
ionization of helium by high energy photons, and the pair creation of
electrons and positrons induced by the photon-helium collisions.

We also update some naive assumptions in our previous code. First, we used the
common cross section for the pair creation in nuclei in the previous code.
However, in the thermal plasma, the pair creation is induced by the free
electrons, neutral and ionized hydrogen, and neutral and ionized of helium,
and cross sections depend on those target particles. We separately treat these
pair creation processes. Second, the cross section of the double-photon pair
creation is sensitive to the energy of the background photons. In the previous
code, however, it was assumed that the background photons have a monochromatic
energy distribution. We include the energy distribution of the CMB photons to
improve the calculation. Third, the time resolution of the previous
calculation was 50 log-spaced bins over the $z\in[0,10^{4}-1]$, which was not
sufficient. In the revised calculation we take 200 bins for
$z\in[0,10^{4}-1]$, which are enough to achieve the $\mathcal{O}(1)$%
precision of results.

We include the above modification in our calculation. The explicit formulas of
calculation are shown in the following, and we also comment on the treatment
of positron, which is not explicitly explained in our previous work.

###  II.1 Electron Interactions

We include electron interactions with helium following Refs. Slatyer  _et al._
(2009); Acharya and Khatri (2019a).

For the excitation of the helium atom, we only include the excitation to 2p
states. For low energy electrons, we use the table of the cross section in
Ref. Stone  _et al._ (2002). For high energy electrons, we use the following
formula Stone  _et al._ (2002):

$\displaystyle\sigma(T)_{\text{He~{}exc}}=\frac{4\pi
a_{0}^{2}\text{Ry}}{T+B+E}\left[a\ln(T/\text{Ry})+b+c\,\text{Ry}/T\right]\frac{f_{\mathrm{acc}}}{f_{\mathrm{sc}}},$
(1)

with the Bohr radius $a_{0}$, the Rydberg constant Ry, the kinetic energy of
incident electron $T$, the excitation energy $E=21.218~{}\mathrm{eV}$, and the
binding energy $B=24.5874~{}\mathrm{eV}$. The dimensionless constants are
given by $a=0.165601,~{}b=-0.076942,~{}c=0.033306$ and
$(f_{\mathrm{acc}}/f_{\mathrm{sc}})=0.2762/0.2583$. The final state electron
has the energy $E_{f}=T-E$, and the differential cross section is given by
$\frac{\text{d}\sigma}{\text{d}E_{f}}(T,E_{f})=\sigma(T)_{\text{He~{}exc}}~{}\delta(E_{f}-(T-E))$.

For the collisional ionization of the helium atom, we use the table of cross
section in the CCC database 111http://atom.curtin.edu.au/CCC-WWW/ for a low
energy collision. For a high energy collision, we use the following formula:
Arnaud and Rothenflug (1985); Slatyer  _et al._ (2009)

$\displaystyle\sigma(E_{i})=10^{-14}\text{cm}^{2}\frac{1}{u(I/\mathrm{eV})^{2}}\left[A(1-u^{-1})+B(1-u^{-1})^{2}+C\ln
u+Du^{-1}\ln u\right],$ (2)

with the ionization threshold $I$ and the normalized energy of incident
electron $u=E_{i}/I$. For the neutral helium, the constants are $A=17.8$,
$B=-11.0$, $C=7.0$, $D=-23.2$, $I=24.6~{}\mathrm{eV}$. The differential cross
section is approximately given by

$\displaystyle\frac{\text{d}\sigma(E,\epsilon)}{\text{d}\epsilon}=\frac{A_{\sigma}(E)}{1+(\epsilon/\bar{\epsilon})^{2}}\quad\text{for}\quad
0<\epsilon<\frac{1}{2}(E-I),$ (3)

where $E$ is the energy of a incident electron, $\epsilon$ is that of a
ejected electron, $\bar{\epsilon}=15.8~{}\mathrm{eV}$ and $A_{\sigma}(E)$ is a
constant independent of $\epsilon$ Opal  _et al._ (1971); Shull (1979);
Furlanetto and Johnson Stoever (2010).

###  II.2 Photon Interactions

Incident photons interact with hydrogen and helium atoms through the photo-
ionization and the pair-creation processes. We obtain the helium photo-
ionization cross section by interpolating the values in the table of Reilman
and Manson (1979).

We update the cross section of pair creation based on Ref. Slatyer  _et al._
(2009). We separately calculate the cross sections of pair creation with
singly ionized hydrogen (HII) and helium (HeII) Motz  _et al._ (1969), neutral
hydrogen (HI) and helium (HeI) Zdziarski and Svensson (1989) and free
electrons Joseph and Rohrlich (1958); they are given by

$\displaystyle\sigma_{\rm HII}$ $\displaystyle=\sigma_{\rm HeII}=\alpha
r_{0}^{2}\left(\frac{28}{9}\ln\left(\frac{2E_{i}}{m_{e}}\right)-\frac{218}{27}\right),$
(4) $\displaystyle\sigma_{\rm HI}$ $\displaystyle=5.4\alpha
r_{0}^{2}\ln\left(\frac{513E_{i}}{E_{i}+825m_{e}}\right),$ (5)
$\displaystyle\sigma_{\rm HeI}$ $\displaystyle=8.76\alpha
r_{0}^{2}\ln\left(\frac{513E_{i}}{E_{i}+825m_{e}}\right),$ (6)
$\displaystyle\sigma_{{\rm free}\ e}$ $\displaystyle=\alpha
r_{0}^{2}\left(\frac{28}{9}\ln\left(\frac{2E_{i}}{m_{e}}\right)-\frac{100}{9}\right),$
(7)

where $E_{i}$ is the initial energy of a photon, $\alpha$ is the fine
structure constant, and $r_{0}=\alpha^{2}a_{0}$ is the classical electron
radius. The differential cross section and the spectrum of produced particles
is estimated by the Bethe-Heitler formula Motz  _et al._ (1969). (See Kanzaki
and Kawasaki (2008) for the used formula in our code. )

The incident photon may cause the double photon pair creation by interacting
with the CMB photons once its energy overcomes the threshold $E_{i}\geq
m_{e}^{2}/E_{\rm cmb}$ where $E_{\rm cmb}$ is the CMB photon energy. We
average the cross section over the energy distribution of the CMB photon to
precisely estimate the effect of the double photon pair creation, following
Ref. Slatyer  _et al._ (2009).

###  II.3 Treatment of positron

For the shower process, positrons are produced through the pair-creation and
double-photon-pair-creation processes. We assume that positrons lose their
energy in the same way as electrons since both behave similarly for
interactions at high energy Slatyer  _et al._ (2009). The produced positrons
quickly lose their kinetic energy within one Hubble time since the background
plasma is not transparent for positrons as well as electrons. Then, they
annihilate the background electrons and produce photons with $E_{\gamma}=511$
keV. In our code, the kinetic energy loss processes of the positrons are
treated in the same way as the electrons and the positron rest mass energy is
converted to two photons with $E_{\gamma}=511$ keV.

Both contributions are included in our shower code.

##  III Energy injection as shower process

We briefly summarize the calculation of electromagnetic showers induced by DM
annihilation based on Kanzaki and Kawasaki (2008); Kanzaki  _et al._ (2010);
Kawasaki  _et al._ (2016). We treat the energy injection as a perturbation on
the standard thermal history of the $\Lambda$CDM model. We focus on the
redshift up to $z=10^{4}$, since the most important effects on the CMB come
from the energy injection near and after the last scattering epoch. In our
calculation, we discretize the redshift and energy of particles and use
logarithmic bins for them.

Let us describe energy injection and the induced electromagnetic shower in the
thermal plasma. Since the electromagnetic processes mainly produce high-energy
photons, electrons, and positrons, we assume that the initial particle “$a$”
is an electron or photon ($a=e$ or $\gamma$). We consider an injected particle
“$a$” with energy $E_{i}$ at redshift $z=z_{i}$, which collides with another
particle “$b$”, produces a cascade shower, and finally its energy is used for
heating the background plasma, ionization, and excitation of the background
hydrogen and helium. Here we denote those processes by $\alpha=$ heat,
ion${}_{\text{H}}$, ion${}_{\text{He}}$, exc, respectively. We refer
$\text{d}\epsilon^{a}_{\alpha}(E_{i},z_{i},z_{f})$ to the partial energy used
for the process $\alpha$ during $[z_{f},z_{f}-\text{d}z_{f}]$. Note that the
redshift $z_{f}$ can be different from $z_{i}$ since the early universe can be
transparent for the high-energy photons. We want to evaluate the following
quantity:

$\displaystyle
Q^{a}_{\alpha}(E_{i},z_{i},z_{f})=\frac{\text{d}\epsilon^{a}_{\alpha}(E_{i},z_{i},z_{f})}{\text{d}z_{f}}.$
(8)

When a particle is injected at slightly earlier time, $z_{i}\to z_{i}+\Delta
z_{i}$, $Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$ is changed by scattering process
during $\Delta z_{i}$. Suppose that the injected particle “$a$” scatters off a
particle “$b$” with the cross section $\sigma(E_{i})$, where the particle
“$b$” is a photon, electron, hydrogen, or helium in the background plasma. The
probability that the particle “$a$” loses its energy to $E_{f}$ for $\Delta
z_{i}$ is given by

$\displaystyle\Delta p^{a}(E_{i},E_{f})$
$\displaystyle=\frac{\text{d}\nu^{a}(E_{i},E_{f})}{\text{d}E_{f}}\frac{\text{d}t}{\text{d}z}\Delta
z_{i},$ (9)

where $\text{d}\nu^{a}/\text{d}E_{f}$ is the differential scattering frequency
of the injected particle,

$\displaystyle\frac{\text{d}\nu^{a}(E_{i},E_{f})}{\text{d}E_{f}}=n_{t}v_{p}\frac{\text{d}\sigma^{a}(E_{i},E_{f})}{\text{d}E_{f}}.$
(10)

Here $n_{t}$ is the number of target particles,
$\text{d}\sigma^{a}/\text{d}E_{f}$ is the differential cross section and
$v_{p}$ is the relative velocity. The differential cross section is evaluated
for the various radiative processes summarized in Section II. We denote
$E^{\text{abs}}_{\alpha}(E_{i},E_{f})$ as the partial energy consumed in the
process $\alpha$, e.g. the binding energy of a hydrogen atom for the
ionization process. When the scattering process occurs at $z_{i}$,
$E^{\text{abs}}_{\alpha}(E_{i},E_{f})$ contributes to
$Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$ as $\Delta
p^{a}(E_{i},E_{f})E^{\text{abs}}_{\alpha}(E_{i},E_{f})\delta(z_{i}-z_{f})$. In
the shower process, the multiple scatterings take place, and
$Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$ depends on that of scattered particle,
$Q^{a}_{\alpha}(E_{f},z_{i},z_{f})$. As a result, the $Q$ changes as

$\displaystyle\Delta Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$
$\displaystyle=\delta(z_{i}-z_{f})\int\text{d}E_{f}~{}\Delta
p^{a}(E_{i},E_{f})E^{\text{abs}}_{\alpha}$
$\displaystyle+\int\text{d}E_{f}~{}\Delta
p^{a}(E_{i},E_{f})\left[Q^{a}_{\alpha}(E_{f},z_{i},z_{f})+Q^{b}_{\alpha}(E_{b},z_{i},z_{f})\right],$
(11)

where $E_{b}(E_{i},E_{f})(\equiv E_{i}-E_{f}-E_{\alpha}^{abs})$ is the final
state energy of the particle $b$.

We rewrite $\Delta Q$ by the difference of $Q$ as

$\displaystyle\Delta Q^{a}_{\alpha}(E_{i},z_{i},z_{f})=$ $\displaystyle
Q^{a}_{\alpha}\left(E_{i}+\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}\Delta
z_{i},z_{i}+\Delta z_{i},z_{f}\right)-Q^{a}_{\alpha}(E_{i},z_{i},z_{f}),$ (12)

where $\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}$ is the energy shift
by the cosmic expansion and is given by
$\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}=E_{i}/(z_{i}+1)$ for a
relativistic particle and
$\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}=2E_{i}/(z_{i}+1)$ for a
non-relativistic particle. The injected energy at $z_{i}+\Delta z_{i}$ should
be shifted to evaluate the difference due to radiative processes correctly.
Using Eqs.(11) and (12), the time evolution of $Q$ is described by

$\displaystyle\frac{\partial Q^{a}_{\alpha}(E_{i},z_{i},z_{f})}{\partial
z_{i}}$
$\displaystyle+\left[\frac{\text{d}E_{i}}{\text{d}z}\right]_{H}\frac{\partial
Q^{a}_{\alpha}(E_{i},z_{i},z_{f})}{\partial E_{i}}$
$\displaystyle=\delta(z_{i}-z_{f})\int\text{d}E_{f}~{}\frac{\text{d}\nu^{a}(E_{i},E_{f})}{\text{d}E_{f}}\frac{\text{d}t(z_{i})}{\text{d}z}E^{\text{abs}}_{\alpha}$
$\displaystyle+\int\text{d}E_{f}~{}\frac{\text{d}\nu^{a}(E_{i},E_{f})}{\text{d}E_{f}}\frac{\text{d}t(z_{i})}{\text{d}z}\left[Q^{a}_{\alpha}(E_{f},z_{i},z_{f})+Q^{b}_{\alpha}(E_{b},z_{i},z_{f})\right].$
(13)

In the numerical calculation, we discretize the energy range
$E_{i}\in[1,10^{13}]\,\mathrm{eV}$ into $5200$ log-spaced bins and the
redshift range $z_{i},z_{f}\in[0,10^{4}-1]$ into $200$ log-spaced bins,
respectively. We solve Eq.(36) based on Eq.(13) (see Appendix A for details)
by using the radiative processes summarized in Section II and evaluate the
$Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$ for wide ranges of injected energy,
injected time, and consumed time.

We check our calculation for low injected energy ($3$ keV, $5$ keV) and high
injected energy ($6$ TeV). The results are shown in Figs. 1-3, where we plot
the ratios of the initial energy to the energy used for heating and ionization
during $z_{f}\in[9,z_{i}]$ as a function of the injected time $z_{i}$ for a
given injected energy $E_{i}$, which is given by

$\displaystyle\frac{\epsilon^{a}_{\alpha}(E_{i},z_{i})}{E_{i}}\equiv\int_{9}^{z_{i}}\text{d}z_{f}\frac{Q^{a}_{\alpha}(E_{i},z_{i},z_{f})}{E_{i}}.$
(14)

Except for a high energy photon, the
$\epsilon^{e}_{\alpha}(E_{i},z_{i})/E_{i}$ is mainly sourced by the
$Q^{a}_{\alpha}(E_{i},z_{i},z_{f})$ at $z_{f}\sim z_{i}$. Following Slatyer
(2016b), we choose the lower bound of the integration at $z_{f}=9$ to avoid
the uncertainty of the cosmic reionization, which is neglected in our code.

In Fig. 1, we show $\epsilon^{e}_{\alpha}(E_{i},z_{i})/E_{i}$ for the injected
electron with $E_{i}\simeq 3\,\mathrm{keV}$ (solid lines). The red, blue, and
green lines represent the energy ratios used for heating, ionization of
hydrogen and helium, respectively. We compare our results with the absorption
rate calculated by Shull and Steenberg Shull and van Steenberg (1985) to check
whether our calculation properly includes the effect of helium. Shull and
Steenberg performed the Monte Carlo computation of the electron injection into
the interstellar gas, with two assumptions on the interstellar gas, $n_{\rm
He}/n_{\text{H}}=0.1$ and $x\equiv x_{\text{H}}=x_{\rm He}$. The derived
fitting formula for the energy ratios to the injected electron energy
$E_{i}\gg 100\,\mathrm{eV}$ as a function of ionization fraction $x$ are given
by

$\displaystyle\frac{\epsilon^{e}_{\rm heat(SS)}}{E_{i}}$
$\displaystyle=C\left[1-(1-x^{a})^{b}\right]\quad,\quad$
$\displaystyle(C,a,b)=(0.9971,0.2663,1.3163),$ (15)
$\displaystyle\frac{\epsilon^{e}_{\rm ionH(SS)}}{E_{i}}$ $\displaystyle=C_{\rm
H}(1-x^{a_{\rm H}})^{b_{\rm H}}\quad,\quad$ $\displaystyle(C_{\rm H},a_{\rm
H},b_{\rm H})=(0.3908,0.4092,1.7592),$ (16)
$\displaystyle\frac{\epsilon^{e}_{\rm ionHe(SS)}}{E_{i}}$
$\displaystyle=C_{\rm He}(1-x^{a_{\rm He}})^{b_{\rm He}}\quad,\quad$
$\displaystyle(C_{\rm He},a_{\rm He},b_{\rm He})=(0.0554,0.4614,1.6660).$ (17)

We plot Eqs. (15)-(17) by dotted lines in Fig. 1 taking $x=x_{\text{H}}$.
Here, we fix $n_{\rm He}/n_{\text{H}}=0.1$ in our code to compare the effect
of helium in the electromagnetic shower. It is seen that our calculation well
reproduces the ionization ratio of He in Ref. Shull and van Steenberg (1985).
Thus, we confirm that our code properly includes the helium interaction around
$E_{i}\leq 3$ keV. We also comment on the slight difference of the heating and
ionization ratios, where our results predict the larger ratios than those of
Shull and Steenberg, which might be due to the different treatment of
interactions. A similar difference of the ionization ratio is reported in Fig.
1 of Slatyer (2016b). For $z\gtrsim 10^{3}$, our energy ratio used for the
helium ionization is much larger than that in Ref. Shull and van Steenberg
(1985), because the ionization fraction of helium is much smaller than that of
hydrogen, i.e, $x_{\rm He}\ll x_{\rm H}$.

Figure 1:  The energy ratios used for heating and ionization [Eq. (14)] in
case of the injected electron $E_{i}\simeq 3\,{\rm keV}$ and injected time
$z_{i}$ integrated over $z_{f}\in[9,z_{i}]$ with $f_{\mathrm{He}}=0.1$. The
solid and dotted lines represent results of our simulation code and the Monte
Carlo computation of the electron injection into the interstellar gas [Eqs.
(15)-(17)] by Shull and Steenberg Shull and van Steenberg (1985),
respectively. The red, blue, and green lines represent the energy ratio used
for heating, ionization of hydrogen, and ionization of helium, respectively.

Figure 2:  The energy ratio integrated over $z_{f}\in[9,z_{i}]$ for the
injected photon with $E_{i}=5\,\mathrm{keV}$ and $6\,$TeV [Eq. (14)]. The red,
blue, and green lines represent the energy ratios used for heating, the total
ionization of hydrogen and helium, and ionization of only helium,
respectively. The dotted lines are taken from Fig. 3 in Slatyer (2016b), where
we rename the sum of “heating” and “sub-10.2 eV continuum photons” in Slatyer
(2016b) to “heating” (red dotted line) since our code does not distinguish
them.

Figure 3:  Same as Fig. 2 except for the injected electron with its mass
energy.

Next, we compare our results with that of the work by Slatyer Slatyer (2016b),
where the numerical table of results over
$E_{i}\in[5\,\mathrm{keV},6\,\mathrm{TeV}]$ are provided. The results at
$E_{i}=6\,\mathrm{TeV}$ include all information of interactions since the
shower process at high energy includes the effects of the low energy
processes. Note that the interaction with highest energy scale is the double
photon pair creation, which occurs at $E_{i}>m_{e}^{2}/E_{\mathrm{cmb}}\sim 1$
TeV at $z=10^{3}$. We also compare our results and the Slatyer’s data at the
lowest energy scale $5$ keV to check the consistency of results on a low
energy scale. The energy absorption ratio for photon injection is shown in
Fig. 2. The dotted lines represent the results of the work by Slatyer Slatyer
(2016b). The blue lines represent the total ionization of hydrogen and helium.
We consider the sum of the “heating” and “sub-10.2 eV continuum photons” in
Ref. Slatyer (2016b) as “heating” in this paper since our code does not
distinguish them. Compared to the photon injection, the definition of energy
injection by an electron is different between our result and Slatyer’s work.
Ref. Slatyer (2016b) considers the $e^{+}e^{-}$ pair injection and the total
injected energy includes the rest mass energy of the electron and positron, so
the energy absorption ratio is written by using our results as

$\displaystyle\frac{E_{i}}{E_{i}+m_{e}}\frac{\epsilon^{e}(E_{i})}{E_{i}}+\frac{m_{e}}{E_{i}+m_{e}}\frac{\epsilon^{\gamma}(m_{e})}{m_{e}},$
(18)

where $E_{i}$ represents the kinetic energy of the electron and we assume that
the rest mass energy of the positron is converted to two $511$ keV photons. We
plot the energy absorption ratio of electron injection in Fig. 3.

One can see that our energy ratio used for the heating and total ionization
(red and blue lines) are similar to those in Ref. Slatyer (2016b). On the
other hand, our ratio used for the ionization of helium (green lines) is
larger because of different treatment of the ionization process. We separately
treat ionization of hydrogen and helium, while Ref. Slatyer (2016b) uses the
approximation that collisional ionization of hydrogen and helium are not
distinguished for high injected energy $E_{i}>3$ keV. The helium atoms are
mostly ionized by the energetic photons with energy
$E\sim(100-1000)\,\mathrm{eV}$. Such photons are produced through inverse
Compton scatterings by electrons with $E\gtrsim 10^{7}\,\mathrm{eV}$ at
$z=10^{3}$ and lose their energy mainly by photo-ionization of helium.
Nevertheless, our calculation for the total energy used for ionization of
hydrogen and helium atoms (blue line) is roughly the same as that in Ref.
Slatyer (2016b), and hence their effects on the CMB also do not change much.

We also mention the difference between our new calculation and the previous
one Kanzaki  _et al._ (2010). While the previous work assumes that all the
baryon number is carried by the hydrogen, we correctly distribute the baryon
number to hydrogen and helium. Thus, the total number of bounded electrons is
almost the same between the new and previous calculations when the helium and
hydrogen are mostly neutral, which leads to a similar energy ratio for the
total ionization. The differences between our previous and new code are also
induced by other modifications summarized in Sec. II and the updated
background cosmological parameters.

##  IV Effect of energy injection on CMB

In the previous section, we calculate the energy fraction
$Q^{a}_{\alpha}(E_{i},z,z^{\prime})$ for $a=e$ and $\gamma$. In general, DM
may annihilate into various particles, e.g. $\mu^{+}\mu^{-}$ and $W^{+}W^{-}$,
which finally decay into stable particles. When DM has an annihilation channel
$F$, we calculate the net energy injection by integrating the number of
electrons, positrons, and photons of the final decay products as

$\displaystyle Q^{(F)}_{\alpha}(E,z,z^{\prime})$
$\displaystyle=\int\text{d}E^{\prime}\bigg{[}\left(\frac{\text{d}N_{F}^{(e^{+})}(E,E^{\prime})}{\text{d}E^{\prime}}+\frac{\text{d}N_{F}^{(e^{-})}(E,E^{\prime})}{\text{d}E^{\prime}}\right)Q^{(e)}_{\alpha}(E^{\prime},z,z^{\prime})$
(19)
$\displaystyle+\frac{\text{d}N_{F}^{(\gamma)}(E,E^{\prime})}{\text{d}E^{\prime}}Q^{(\gamma)}_{\alpha}(E^{\prime},z,z^{\prime})\bigg{]},$
(20)

where $dN_{f}^{a}(E,E^{\prime})/dE^{\prime}$ is the number of
$a(=e^{\pm},\gamma)$ produced through the decay channel $F$. In this paper, we
consider $F=e^{+}e^{-},~{}\gamma\gamma,~{}\mu^{+}\mu^{-}$ and $W^{+}W^{-}$ as
a typical decay channel. (Notice that
$dN_{f}^{\gamma}(E,E^{\prime})/dE^{\prime}=\delta(E-E^{\prime})$ for
$F=\gamma\gamma$ and
$dN_{f}^{e^{+}}(E,E^{\prime})/dE^{\prime}=dN_{f}^{e^{-}}(E,E^{\prime})/dE^{\prime}=\delta(E-E^{\prime})/2$
for $F=e^{+}e^{-}$.) The decay fractions into $\mu^{+}\mu^{-}$ and
$W^{+}W^{-}$ are calculated by the PYTHIA package Sjostrand  _et al._ (2006).

We calculate the thermal history adding the following three terms in the
RECFAST code:

$\displaystyle-\left[\frac{\text{d}x_{\text{H}}}{\text{d}z}\right]_{\text{DM}}$
$\displaystyle=\sum_{F}\int_{z}\frac{\text{d}z^{\prime}}{H(z^{\prime})(1+z^{\prime})}\frac{n_{\text{DM}}^{2}(z^{\prime})\braket{\sigma
v}_{F}}{n_{\text{H}}(z^{\prime})}\frac{Q_{\text{ion}_{\text{H}}}^{(F)}(m_{\text{DM}},z^{\prime},z)}{B_{\text{H}}},$
(21)
$\displaystyle-\left[\frac{\text{d}x_{\text{He}}}{\text{d}z}\right]_{\text{DM}}$
$\displaystyle=\sum_{F}\int_{z}\frac{\text{d}z^{\prime}}{H(z^{\prime})(1+z^{\prime})}\frac{n^{2}_{DM}(z^{\prime})\braket{\sigma
v}_{F}}{n_{\text{He}}(z^{\prime})}\frac{Q_{\text{ion}_{\text{He}}}^{(F)}(m_{\text{DM}},z^{\prime},z)}{B_{\text{He}}},$
(22) $\displaystyle-\left[\frac{\text{d}T_{M}}{\text{d}z}\right]_{\text{DM}}$
$\displaystyle=\sum_{F}\int_{z}\frac{\text{d}z^{\prime}}{H(z^{\prime})(1+z^{\prime})}\frac{2}{3}\frac{n^{2}_{DM}(z^{\prime})\braket{\sigma
v}_{F}}{(1+f_{\text{He}}+x_{e})n_{H}(z^{\prime})}Q_{\text{heat}}^{(F)}(m_{\text{DM}},z^{\prime},z),$
(23)

where $n_{\text{DM}}$ is the DM number density and the binding energies of HI
and HeII are $B_{\text{H}}=13.6057$ eV and $B_{\text{He}}=24.5874$ eV. For the
fixed DM energy density $\rho_{\text{DM}}$, the above modifications explicitly
depends on $n_{\text{DM}}^{2}\braket{\sigma
v}m_{\text{DM}}=\rho_{\text{DM}}^{2}\braket{\sigma v}/m_{\text{DM}}$. Thus, we
expect the constraints on $\braket{\sigma v}/m_{\text{DM}}$ by the CMB, which
slightly depends on $m_{\text{DM}}$ through
$Q^{(F)}_{\alpha}(m_{\text{DM}},z^{\prime},z)/m_{\text{DM}}$.

The modified thermal history changes the CMB power spectrum mainly through the
enhanced Thomson scattering by the free electron. The Thomson scattering
suppresses the CMB temperature mode by increasing the optical depth, and at
the same time, the CMB polarization mode is produced by the Thomson
scattering. The latter effect mainly contributes to the constraint on the
annihilation cross section. This modification depends on the total number of
ionization but not on the species of the ionized atom, e.g. hydrogen or helium
Galli  _et al._ (2013); Aghanim  _et al._ (2020c).

In this paper, we calculate the CMB power spectrum using the CAMB Lewis  _et
al._ (2000) code with the modified RECFAST code Seager  _et al._ (2000, 1999);
Wong  _et al._ (2008). The energy injection changes the power spectrum of the
CMB. As an example, in Fig. 4 we calculate the difference of TT- and EE- mode
power spectrum defined by

$\displaystyle\Delta_{C_{l}}=\frac{C^{\rm X}_{l}-C^{\rm X,(0)}_{l}}{C^{\rm
X,(0)}_{l}},$ (24)

where $X=TT$ or $EE$, and $C^{\rm X,(0)}_{l}$ is the power spectrum of CMB
without energy injection by DM and $C^{\rm X}_{l}$ is the spectrum with energy
injection by $e^{+}e^{-}$ channel for $m_{\rm DM}=100$ GeV and $\braket{\sigma
v}/m_{\rm DM}=1.0\times 10^{-26}{\rm cm}^{3}{\rm sec}^{-1}{\rm GeV}^{-1}$. In
Fig. 4, the blue, green, red and black lines represent the modified power
spectra including the effects of ionization of hydrogen [Eq. (21)], ionization
of helium [Eq. (22)], heating [Eq. (23)], and all three effects, respectively.
It is seen that the dominant effect on the CMB power spectrum comes from the
ionization of hydrogen, while the effect of helium is about 20% of hydrogen
and heating hardly affects the power spectrum.

Figure 4:  The difference of TT- and EE-mode power spectrum [Eq. (24)] by the
energy injection of DM annihilation into $e^{+}e^{-}$ channel for $m_{\rm
DM}=100$ GeV and $\braket{\sigma v}/m_{\rm DM}=1.0\times 10^{-26}{\rm
cm}^{3}{\rm sec}^{-1}{\rm GeV}^{-1}$. The blue, green, red and black lines
represent the modified power spectrum including the effects of ionization of
hydrogen in Eq. (21), ionization of helium in Eq. (22), heating in Eq. (23),
and all three contributions in Eqs. (21)-(23), respectively.

##  V Constraints on the Dark matter annihilation

In this section, we derive constraints on the DM annihilation cross section
from the CMB observational data. We run the CosmoMC code Lewis and Bridle
(2002) with the modified CAMB code as described in the previous sections and
search for the DM parameter space $(m_{\text{DM}},\braket{\sigma
v}/m_{\text{DM}})$. The range of DM mass is
$m_{\text{DM}}\in[1,10^{4}]\,\mathrm{GeV}$ for $2\gamma,~{}e^{+}e^{-}$, and
$\mu^{+}\mu^{-}$ channel, and $m_{\text{DM}}\in[80,10^{4}]\,\mathrm{GeV}$ for
$W^{+}W^{-}$ channel. The range of the cross section divided by mass is given
by $\braket{\sigma v}/m_{\text{DM}}\in[0,3\times 10^{-27}]\text{
cm}^{3}\sec^{-1}\mathrm{GeV}^{-1}$ for $e^{+}e^{-}$ and $2\gamma$ channels,
and $\braket{\sigma v}/m_{\text{DM}}\in[0,6\times 10^{-27}]\text{
cm}^{3}\sec^{-1}\mathrm{GeV}^{-1}$ for $\mu^{+}\mu^{-}$ and $W^{+}W^{-}$
channels. We use the top-hat prior for $\log(m_{\text{DM}})$ and
$\braket{\sigma v}/m_{\text{DM}}$. As for the cosmological model, we use the
base-$\Lambda$CDM model in the Planck collaboration Aghanim  _et al._ (2020c).
Except for the DM mass and annihilation cross section, we vary the six
cosmological parameters
$\\{\Omega_{b}h^{2},~{}\Omega_{c}h^{2},~{}100\theta_{MC},~{}\tau,~{}n_{s},~{}{\rm{ln}}(10^{10}A_{s})\\}$
and neutrinos with normal mass hierarchy and $\sum m_{\nu}=0.06$ eV.

We use two combinations of data sets to constrain $\braket{\sigma
v}/m_{\text{DM}}$. The first one is the CMB anisotropy from the Planck 2018
data set (“Planck only”) including low-l ($l\leq 29$) TT- and EE-mode (“lowl”
and “simall_EE”) , and high-l ($l\geq 30$) TT-, TE- and EE-mode
(“plik_rd12_HM_v22_TTTEEE”) Aghanim  _et al._ (2020b). The other one,
“Planck+ext”, consists of the Planck CMB data and other cosmological data on
galaxy correlation functions including the baryon acoustic oscillation (BAO)
York  _et al._ (2000); Anderson  _et al._ (2014); Ross  _et al._ (2015);
Beutler  _et al._ (2011), Dark Energy Survey (DES) Abbott  _et al._ (2018),
and the CMB lensing power spectrum. The latter choice is expected to put
stronger constraints than the former, since it can solve the parameter
degeneracy.

We show the constraints on the annihilation cross section in the
$(\braket{\sigma v}/m_{\rm DM},m_{\rm DM})$ plane in Fig. 5 and in the
$(\braket{\sigma v},m_{\rm DM})$ plane in Fig. 6, respectively. The contours
plotted by the GetDist code Lewis (2019) describe the $95\%$ exclusion region.
Comparing “Planck only” and “Planck+ext” in Fig. 5, one can see that the extra
data of BAO, DES and CMB lensing power spectrum tighten the constraints about
$O(10\%)$. The relative improvements are slightly dependent on the channel. It
is not easy to figure out the physical reasons of this dependence because the
role of the extra data (BAO, DES and CMB lensing) to constrain the DM model
parameters is indirect and complicated. Note also that our numerical precision
may be close to the limitation to resolve this slight channel dependence. The
constraints for $W^{+}W^{-}$ and $\mu^{+}\mu^{-}$ modes are about twice weaker
than those for $e^{+}e^{-}$ and $2\gamma$ modes. These features are consistent
with our previous paper Kawasaki  _et al._ (2016). Since we change the prior
of the Monte Carlo calculation of CosmoMC from the flat prior on
$\braket{\sigma v}$ to the flat prior on $\braket{\sigma v}/m_{\rm DM}$, the
upper bound on $\braket{\sigma v}$ is changed about $\mathcal{O}(10\%)$.

Figure 5:  CMB constraints on the annihilation cross section for four decay
channels of dark matter. We plot the 95% exclusion contours of $\braket{\sigma
v}/m_{\rm DM}$ (vertical axes) for mass of dark matter (horizontal axes). We
calculate the posterior distribution using two sets of data: “Planck only”
uses the Planck CMB data for TT-, TE- and EE-modes Aghanim  _et al._ (2020b)
and “Planck+ext” uses BAO York  _et al._ (2000); Anderson  _et al._ (2014);
Ross  _et al._ (2015); Beutler  _et al._ (2011), DES Abbott  _et al._ (2018),
and the CMB lensing power spectrum data in addition to the Planck CMB data.
Figure 6:  CMB constraints on the annihilation cross section for four DM
annihilation channels. We plot the 95% exclusion contours of $\braket{\sigma
v}$, which represent the same contours in Fig. 5.

Constraints on the DM annihilation cross section has also been studied in many
previous works Padmanabhan and Finkbeiner (2005); Kanzaki and Kawasaki (2008);
Kanzaki  _et al._ (2010); Slatyer  _et al._ (2009); Slatyer (2013); Galli  _et
al._ (2013); Slatyer (2016a, b); Liu  _et al._ (2016); Acharya and Khatri
(2019a); Cang  _et al._ (2020). Among them, Refs. Slatyer (2016a, b) include
the effects of electromagnetic shower by injected particles based on Refs.
Slatyer  _et al._ (2009); Galli  _et al._ (2013), where electromagnetic shower
is calculated differently above and below the threshold energy $3$ keV. We
confirm that our calculation reproduces similar results on the electromagnetic
shower to their calculation in Sec. III. Our “Planck only” result for the
$e^{+}e^{-}$ channel in Eq. (6) is consistent with Fig. 4 of Ref. Slatyer
(2016a) within an accuracy of 10%.

##  VI Conclusion

DM candidates like WIMPs are once in thermal equilibrium in the early universe
and decouple from the thermal bath when their annihilation rate becomes slower
than the cosmic expansion. The abundance of such thermal relics is determined
by the annihilation cross section. Although the DM annihilation rate is
smaller than the expansion rate, the annihilation of DM particles continues to
take place after decoupling and may cause significant effects on cosmology and
astrophysics. We revisited the CMB constraints on the DM annihilation cross
section.

In this paper, we calculated the electromagnetic energy injection due to DM
annihilation and the subsequent radiative processes based on the previously
developed method Kanzaki and Kawasaki (2008); Kanzaki  _et al._ (2010), and
obtain constraints on the annihilation cross section from the CMB power
spectrum. We have improved our previous calculation Kawasaki  _et al._ (2016)
on the following points: we have updated the observational data set including
the Planck 2018 result Aghanim  _et al._ (2020b), carefully included the
helium interactions, fixed the incorrect simplification on some interactions,
and improve the precision of each calculation. Although it was pointed out in
Galli  _et al._ (2013); Aghanim  _et al._ (2020c) that the helium interaction
in the shower process does not much change the constraints on DM annihilation
cross section, the introduction of helium interaction improves the reliability
of our code, and the consistent treatment of the shower process will become
important for the other application of the energy injection into the thermal
plasma.

Our calculation of the electromagnetic shower process is comparable to that in
Ref. Slatyer (2016a) and our method is simpler. While in Ref. Slatyer (2016a)
the shower process is separately calculated by different codes for the
low/high energy region, our code can consistently calculate the shower process
in both high and low energy regions without unnecessary assumptions. We use
the recursive method by separating the “discrete” and “continuous”
interactions, where the scattered particle loses larger (smaller) energy than
the size of energy bin $\Delta E$ in discrete (continuous) interactions. Our
results are fairly consistent with the results of Ref. Slatyer (2016a) for the
$e^{+}e^{-}$ channel. We also confirm that our calculation reproduces the
Monte Carlo calculation on the electromagnetic shower of particle scattering
given in Ref. Shull and van Steenberg (1985). Our results are summarized in
Fig. 6, where the 95% exclusion limits are shown for four decay channels:
$2\gamma$, $e^{+}e^{-}$, $\mu^{+}\mu^{-}$ and $W^{+}W^{-}$.

##  Appendix A Numerical scheme of shower code

In this Appendix, we describe our scheme of the numerical calculation. We
discretize the energy and redshift in log-scale. First, we define the
normalized energy ratio $\chi^{a}_{\alpha}$ as

$\chi^{a}_{\alpha}(E_{i},z_{m},z_{n})\equiv\Delta
z_{n}~{}Q^{a}_{\alpha}(E_{i},z_{m},z_{n})/E_{i},$ (25)

with a size of bin $\Delta z_{n}$. Here $z_{m}$ and $z_{n}$ represent the
redshifts at injection and absorption, respectively. Using $\chi_{\alpha}^{a}$
we rewrite Eq.(12) as

$\displaystyle\Delta Q^{a}_{\alpha}(E_{i},z_{m},z_{n})$
$\displaystyle=\frac{E_{i}}{\Delta
z_{n}}\left[\chi^{a}_{\alpha}(E_{i},z_{m},z_{n})-\frac{E_{i}-\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}\Delta
z_{m}}{E_{i}}\chi^{a}_{\alpha}\left(E_{i}-\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}\Delta
z_{m},z_{m}-\Delta z_{m},z_{n}\right)\right],$ (26)

where we slightly change the time variable from $z_{i}$ to $z_{m}-\Delta
z_{m}$. In the same way, we rewrite Eq.(11) as

$\displaystyle\Delta Q^{a}_{\alpha}(E_{i},z_{m},z_{n})$
$\displaystyle=\frac{E_{i}}{\Delta z_{n}}\bigg{[}\Delta
z_{n}\delta(z_{n}-z_{m})\sum_{j}\Delta E_{j}~{}\Delta
p^{a}(E_{i},E_{j})\frac{E^{\text{abs}}_{\alpha}}{E_{i}}$
$\displaystyle+\sum_{j}\Delta E_{j}~{}\Delta
p^{a}(E_{i},E_{j})\left(\frac{E_{j}}{E_{i}}\chi^{a}_{\alpha}(E_{j},z_{m},z_{n})+\frac{E_{b}}{E_{i}}\chi^{b}_{\alpha}(E_{b},z_{m},z_{n})\right)\bigg{]}$
(27) $\displaystyle=\frac{E_{i}}{\Delta
z_{n}}\bigg{[}\delta_{m,n}\sum_{j}\Delta
P^{a}(E_{i},E_{j},z_{m})\frac{E^{\text{abs}}_{\alpha}}{E_{i}}$
$\displaystyle+\sum_{j}\Delta
P^{a}(E_{i},E_{j},z_{m})\left(\frac{E_{j}}{E_{i}}\chi^{a}_{\alpha}(E_{j},z_{m},z_{n})+\frac{E_{b}}{E_{i}}\chi^{b}_{\alpha}(E_{b},z_{m},z_{n})\right)\bigg{]},$
(28)

where we use $\Delta z_{n}\delta(z_{n}-z_{m})\to\delta_{m,n}$ and define the
differential scattering probability as

$\displaystyle\Delta P^{a}(E_{i},E_{j},z_{m})$
$\displaystyle\equiv\Delta\nu^{a}(E_{i},E_{j})\frac{\text{d}t(z_{i})}{\text{d}z}\Delta
z_{n},$ (29)

with the differential collisional frequency:

$\displaystyle\Delta\nu^{a}(E_{i},E_{j})$ $\displaystyle\equiv\Delta
E_{j}~{}\frac{\text{d}\nu^{a}(E_{i},E_{j})}{\text{d}E_{j}}=n_{t}v_{p}\frac{\text{d}\sigma^{a}(E_{i},E_{j})}{\text{d}E_{j}}\Delta
E_{j}.$ (30)

We should be careful when the difference of energy before and after scattering
is small. At first, when some fraction of the scattered particles belong to
the same energy bin as the initial particles do, the effective scattering
probability should be smaller than that calculated by Eq. (30). We call such
scattering “continuous loss”, and its collisional frequency is given by the
ratio of losing energy to $\Delta E_{i}$ as

$\displaystyle\Delta\nu^{a}(E_{i},E_{i-1})=$ $\displaystyle\frac{1}{\Delta
E_{i}}\left[\frac{\text{d}E}{\text{d}t}\right]_{c},$ (31)

where $\left[\tfrac{\text{d}E}{\text{d}t}\right]_{c}$ is the energy loss by
the continuous loss process. Second, when the plasma is optically thick for an
injected particle, the sum of the scattering probability during $\Delta z_{m}$
can exceed $1$. In this case, the scattering probability at one channel is
given by the ratio of its differential scattering frequency to the total as

$\displaystyle\Delta P^{a}(E_{i},E_{j},z_{m})$
$\displaystyle=\frac{\Delta\nu^{a}(E_{i},E_{j})}{\sum_{k}\Delta\nu^{a}(E_{i},E_{k})+\nu_{H}^{a}(E_{i})}.$
(32)

Here $\nu_{H}^{a}(E_{i})$ is the effective frequency of the Hubble friction
defined by 222 Note that Eq.(32) coincides with Eq.(29) in the continuous
limit of time step since we can derive  $\displaystyle\nu_{H}(E_{i})$
$\displaystyle=\frac{H(z)E_{i}}{\Delta
z_{m}\frac{\text{d}E_{i}}{\text{d}z}}=\frac{H(z)E_{i}}{\Delta
t\frac{\text{d}E_{i}}{\text{d}t}}=\frac{H(z)E_{i}}{\Delta tHE_{i}}=(\Delta
t)^{-1}.$ (33)

$\displaystyle\nu_{H}(E_{i})\equiv\frac{H(z)E_{i}}{\Delta_{z}E_{i}},$ (34)

where $\Delta_{z}E_{i}=E_{i}\,\Delta z_{m}/(1+z_{m})$. We also normalize the
ratio of redshifted particles in the second term of Eq.(26) using the
following factor:

$\displaystyle\Delta P^{a}_{H}(E_{i},E_{j},z_{m})$
$\displaystyle=\frac{\nu_{H}^{a}(E_{i})}{\sum_{k}\Delta\nu^{a}(E_{i},E_{k})+\nu_{H}^{a}(E_{i})}.$
(35)

Combining these formulas, we can perform the numerical calculation of the
shower process based on the following formula:

$\displaystyle\chi^{a}_{\alpha}(E_{i},z_{m},z_{n})$
$\displaystyle=\delta_{mn}\sum_{j}\Delta
P^{a}(E_{i},E_{j},z_{m})\frac{E^{\text{abs}}_{\alpha}}{E_{i}}$
$\displaystyle\quad+\sum_{j<i}\Delta
P^{a}(E_{i},E_{j},z_{m})\left[\frac{E_{j}}{E_{i}}\chi^{a}_{\alpha}(E_{j},z_{m},z_{n})+\frac{E_{b}}{E_{i}}\chi^{b}_{\alpha}(E_{b},z_{m},z_{n})\right]$
$\displaystyle\quad+\delta_{m\neq
n}P_{H}^{a}(E_{i},z_{m})\chi^{a}_{\alpha}(E_{i}-\Delta_{z}E_{i},z_{m-1},z_{n})\frac{E_{i}-\Delta_{z}E_{i}}{E_{i}},$
(36)

where $\delta_{m\neq n}$ is 0 for $m=n$ and 1 for $m\neq n$.

To check the energy conservation, we also calculate the energy used for the
cosmic expansion (we denote this as process “$z$”, i.e. $\alpha=z$). Then we
can write $\chi^{a}_{z}$ as

$\displaystyle\chi_{z}^{a}(E_{i},z_{m},z_{n})$
$\displaystyle=\delta_{mn}P_{H}^{a}(E_{i},z_{m})\frac{\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}\Delta
z_{m}}{E_{i}}$ $\displaystyle\quad+\sum_{j<i}\Delta
P^{a}(E_{i},E_{j},z_{m})\left[\frac{E_{j}}{E_{i}}\chi^{a}_{z}(E_{j},z_{m},z_{n})+\frac{E_{b}}{E_{i}}\chi^{b}_{z}(E_{b},z_{m},z_{n})\right]$
$\displaystyle\quad+\delta_{m\neq
n}P_{H}^{a}(E_{i},z_{m})\chi^{a}_{z}(E_{i}-\Delta_{z}E_{i},z_{m-1},z_{n})\frac{E_{i}-\Delta_{z}E_{i}}{E_{i}}.$
(37)

We also calculate the energy fraction $\chi^{a}_{r}$ which is not used for any
process until $z_{m}$ and $\chi^{a}_{r}$ is written as

$\displaystyle\chi^{a}_{r}(E_{i},z_{m},z_{n})$
$\displaystyle=\delta_{mn}P_{H}^{a}(E_{i},z_{m})\frac{E_{i}-\left[\tfrac{\text{d}E_{i}}{\text{d}z}\right]_{H}\Delta
z_{m}}{E_{i}}$ $\displaystyle\quad+\sum_{j<i}\Delta
P^{a}(E_{i},E_{j},z_{m})\left[\frac{E_{j}}{E_{i}}\chi^{a}_{r}(E_{j},z_{m},z_{n})+\frac{E_{b}}{E_{i}}\chi^{b}_{r}(E_{b},z_{m},z_{n})\right]$
$\displaystyle\quad+\delta_{m\neq
n}P_{H}^{a}(E_{i},z_{m})\chi^{a}_{r}(E_{i}-\Delta_{z}E_{i},z_{m-1},z_{n})\frac{E_{i}-\Delta_{z}E_{i}}{E_{i}}.$
(38)

We then calculate $\chi^{a}_{\alpha}(E_{i},z_{m},z_{n})$ ($\alpha=$ heat,
ion${}_{\text{H}}$, ion${}_{\text{He}}$, exc, z, r) from low to high energy
and from low to high redshift. We have confirmed the energy conservation,
$\sum_{\alpha}\chi_{\alpha}^{a}(E_{i},z_{m},z_{m})=1$ for $z_{m}=z_{n}$ and
$\sum_{\alpha}\chi_{\alpha}^{a}(E_{i},z_{m},z_{n})=\chi_{r}^{a}(E_{i},z_{m},z_{n+1})$
for $z_{m}>z_{n}$.

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Nos. 17H01131 (M.K.), 17K05434
(M.K.), 20H05851(M.K.), 21K03567(M.K.), JP19J21974 (H.N.), 18K03609 (K.N.),
17H06359 (K.N.), Advanced Leading Graduate Course for Photon Science (H.N.),
and World Premier International Research Center Initiative (WPI Initiative),
MEXT, Japan (M.K.).

## References

  * Padmanabhan and Finkbeiner (2005) N. Padmanabhan and D. P. Finkbeiner, Phys. Rev. D 72, 023508 (2005), arXiv:astro-ph/0503486 .
  * Mapelli  _et al._ (2006) M. Mapelli, A. Ferrara, and E. Pierpaoli, Mon. Not. Roy. Astron. Soc. 369, 1719 (2006), arXiv:astro-ph/0603237 .
  * Zhang  _et al._ (2006) L. Zhang, X.-L. Chen, Y.-A. Lei, and Z.-G. Si, Phys. Rev. D 74, 103519 (2006), arXiv:astro-ph/0603425 .
  * Kanzaki and Kawasaki (2008) T. Kanzaki and M. Kawasaki, Phys. Rev. D78, 103004 (2008), arXiv:0805.3969 [astro-ph] .
  * Kanzaki  _et al._ (2010) T. Kanzaki, M. Kawasaki, and K. Nakayama, Prog. Theor. Phys. 123, 853 (2010), arXiv:0907.3985 [astro-ph.CO] .
  * Slatyer  _et al._ (2009) T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, Phys. Rev. D80, 043526 (2009), arXiv:0906.1197 [astro-ph.CO] .
  * Slatyer (2013) T. R. Slatyer, Phys. Rev. D 87, 123513 (2013), arXiv:1211.0283 [astro-ph.CO] .
  * Galli  _et al._ (2013) S. Galli, T. R. Slatyer, M. Valdes, and F. Iocco, Phys. Rev. D 88, 063502 (2013), arXiv:1306.0563 [astro-ph.CO] .
  * Slatyer (2016a) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016a), arXiv:1506.03811 [hep-ph] .
  * Slatyer (2016b) T. R. Slatyer, Phys. Rev. D 93, 023521 (2016b), arXiv:1506.03812 [astro-ph.CO] .
  * Liu  _et al._ (2016) H. Liu, T. R. Slatyer, and J. Zavala, Phys. Rev. D 94, 063507 (2016), arXiv:1604.02457 [astro-ph.CO] .
  * Acharya and Khatri (2019a) S. K. Acharya and R. Khatri, JCAP 1912, 046 (2019a), arXiv:1910.06272 [astro-ph.CO] .
  * Cang  _et al._ (2020) J. Cang, Y. Gao, and Y.-Z. Ma, Phys. Rev. D 102, 103005 (2020), arXiv:2002.03380 [astro-ph.CO] .
  * Ackermann  _et al._ (2015) M. Ackermann _et al._ (Fermi-LAT), Phys. Rev. Lett. 115, 231301 (2015), arXiv:1503.02641 [astro-ph.HE] .
  * Ahnen  _et al._ (2016) M. L. Ahnen _et al._ (MAGIC, Fermi-LAT), JCAP 02, 039 (2016), arXiv:1601.06590 [astro-ph.HE] .
  * Aguilar  _et al._ (2013) M. Aguilar _et al._ (AMS), Phys. Rev. Lett. 110, 141102 (2013).
  * Adriani  _et al._ (2009) O. Adriani _et al._ (PAMELA), Nature 458, 607 (2009), arXiv:0810.4995 [astro-ph] .
  * Aharonian  _et al._ (2009) F. Aharonian _et al._ (H.E.S.S.), Astron. Astrophys. 508, 561 (2009), arXiv:0905.0105 [astro-ph.HE] .
  * Albert  _et al._ (2020) A. Albert _et al._ (ANTARES), Phys. Lett. B 805, 135439 (2020), arXiv:1912.05296 [astro-ph.HE] .
  * Gelmini (2015) G. B. Gelmini, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders_ (2015) arXiv:1502.01320 [hep-ph] .
  * Reno and Seckel (1988) M. H. Reno and D. Seckel, Phys. Rev. D 37, 3441 (1988).
  * Jedamzik (2004) K. Jedamzik, Phys. Rev. D 70, 083510 (2004), arXiv:astro-ph/0405583 .
  * Hisano  _et al._ (2009) J. Hisano, M. Kawasaki, K. Kohri, T. Moroi, and K. Nakayama, Phys. Rev. D 79, 083522 (2009), arXiv:0901.3582 [hep-ph] .
  * Hisano  _et al._ (2011) J. Hisano, M. Kawasaki, K. Kohri, T. Moroi, K. Nakayama, and T. Sekiguchi, Phys. Rev. D 83, 123511 (2011), arXiv:1102.4658 [hep-ph] .
  * Kawasaki  _et al._ (2015) M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, Phys. Lett. B 751, 246 (2015), arXiv:1509.03665 [hep-ph] .
  * Depta  _et al._ (2019) P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, JCAP 04, 029 (2019), arXiv:1901.06944 [hep-ph] .
  * Acciari  _et al._ (2010) V. A. Acciari _et al._ (VERITAS), Astrophys. J. 720, 1174 (2010), arXiv:1006.5955 [astro-ph.CO] .
  * Abramowski  _et al._ (2014) A. Abramowski _et al._ (H.E.S.S.), Phys. Rev. D 90, 112012 (2014), arXiv:1410.2589 [astro-ph.HE] .
  * Adams  _et al._ (1998) J. A. Adams, S. Sarkar, and D. W. Sciama, Mon. Not. Roy. Astron. Soc. 301, 210 (1998), arXiv:astro-ph/9805108 .
  * Chen and Kamionkowski (2004) X.-L. Chen and M. Kamionkowski, Phys. Rev. D 70, 043502 (2004), arXiv:astro-ph/0310473 .
  * Poulin  _et al._ (2017) V. Poulin, J. Lesgourgues, and P. D. Serpico, JCAP 03, 043 (2017), arXiv:1610.10051 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 [astro-ph.CO] .
  * Kawasaki  _et al._ (2016) M. Kawasaki, K. Nakayama, and T. Sekiguchi, Phys. Lett. B756, 212 (2016), arXiv:1512.08015 [astro-ph.CO] .
  * Acharya and Khatri (2019b) S. K. Acharya and R. Khatri, Phys. Rev. D 99, 043520 (2019b), arXiv:1808.02897 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 [astro-ph.CO] .
  * Shull and van Steenberg (1985) J. M. Shull and M. E. van Steenberg, The Astrophysical Journal 298, 268 (1985).
  * Seager  _et al._ (2000) S. Seager, D. D. Sasselov, and D. Scott, Astrophys. J. Suppl. 128, 407 (2000), arXiv:astro-ph/9912182 .
  * Seager  _et al._ (1999) S. Seager, D. D. Sasselov, and D. Scott, Astrophys. J. Lett. 523, L1 (1999), arXiv:astro-ph/9909275 .
  * Wong  _et al._ (2008) W. Y. Wong, A. Moss, and D. Scott, Mon. Not. Roy. Astron. Soc. 386, 1023 (2008), arXiv:0711.1357 [astro-ph] .
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000), arXiv:astro-ph/9911177 .
  * Lewis and Bridle (2002) A. Lewis and S. Bridle, Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 .
  * Stone  _et al._ (2002) P. M. Stone, Y.-K. Kim, and J. Desclaux, Journal of research of the National Institute of Standards and Technology 107, 327 (2002).
  * Arnaud and Rothenflug (1985) M. Arnaud and R. Rothenflug, Astronomy and Astrophysics Supplement Series 60, 425 (1985).
  * Opal  _et al._ (1971) C. Opal, W. Peterson, and E. Beaty, The Journal of Chemical Physics 55, 4100 (1971).
  * Shull (1979) J. M. Shull, Astrophys. J. 234, 761 (1979).
  * Furlanetto and Johnson Stoever (2010) S. R. Furlanetto and S. Johnson Stoever, Monthly Notices of the Royal Astronomical Society (2010), 10.1111/j.1365-2966.2010.16401.x.
  * Reilman and Manson (1979) R. F. Reilman and S. T. Manson, ApJS 40, 815 (1979).
  * Motz  _et al._ (1969) J. Motz, H. Olsen, and H. Koch, Rev. Mod. Phys. 41, 581 (1969).
  * Zdziarski and Svensson (1989) A. A. Zdziarski and R. Svensson, The Astrophysical Journal 344, 551 (1989).
  * Joseph and Rohrlich (1958) J. Joseph and F. Rohrlich, Rev. Mod. Phys. 30, 354 (1958).
  * Sjostrand  _et al._ (2006) T. Sjostrand, S. Mrenna, and P. Z. Skands, JHEP 05, 026 (2006), arXiv:hep-ph/0603175 .
  * Aghanim  _et al._ (2020c) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020c), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * York  _et al._ (2000) D. G. York _et al._ (SDSS), Astron. J. 120, 1579 (2000), arXiv:astro-ph/0006396 .
  * Anderson  _et al._ (2014) L. Anderson _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 441, 24 (2014), arXiv:1312.4877 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Monthly Notices of the Royal Astronomical Society 416, 3017 (2011).
  * Abbott  _et al._ (2018) T. M. C. Abbott _et al._ (DES), Phys. Rev. D 98, 043526 (2018), arXiv:1708.01530 [astro-ph.CO] .
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .