DarkHistory: A code package for calculating modified cosmic ionization and thermal histories with dark matter and other exotic energy injections

TL;DR

DarkHistory is a Python package that models how dark matter and exotic energy injections influence the early universe's ionization and temperature history, enabling detailed studies of their cosmological impacts.

Contribution

It introduces a self-consistent, flexible code for calculating the effects of dark matter annihilation and decay on cosmic ionization and thermal evolution, including backreaction effects.

Findings

Demonstrates the code's ability to handle arbitrary dark matter parameters.

Shows how to incorporate astrophysical and exotic energy sources.

Provides examples of impact on 21cm signals and reionization.

Abstract

We present a new public Python package, DarkHistory, for computing the effects of dark matter annihilation and decay on the temperature and ionization history of the early universe. DarkHistory simultaneously solves for the evolution of the free electron fraction and gas temperature, and for the cooling of annihilation/decay products and the secondary particles produced in the process. Consequently, we can self-consistently include the effects of both astrophysical and exotic sources of heating and ionization, and automatically take into account backreaction, where modifications to the ionization/temperature history in turn modify the energy-loss processes for injected particles. We present a number of worked examples, demonstrating how to use the code in a range of different configurations, in particular for arbitrary dark matter masses and annihilation/decay final states. Possible applications of DarkHistory include mapping out the effects of dark matter annihilation/decay on the global 21cm signal and the epoch of reionization, as well as the effects of exotic energy injections other than dark matter annihilation/decay. The code is available at https://github.com/hongwanliu/DarkHistory with documentation at https://darkhistory.readthedocs.io . Data files required to run the code can be downloaded at https://doi.org/10.7910/DVN/DUOUWA .

Tables

Table 1. Table 1: List of q 𝑞 q -coefficients for use in Eqs. ( 41 )–( 43 ). The variable q 𝑞 q is defined in Eq. ( 44 ).

Method	$q_{\text{H}}^{\gamma }$	$q_{\text{H}}^e$	$q_{\text{He}}^{\gamma }$	$q_{\text{He,a}}^e$	$q_{\text{He,b}}^e$
1	1	1	0	0	0
2	1	$q$	0	$1-q$	$1-q$
3	$q$	$q$	$1-q$	$1-q$	0

Table 2. Table 2: Series expressions for the relevant indefinite integrals of the form shown in Eq. ( 78 ). Here, y 𝑦 y is the integration variable, and x 𝑥 x denotes the integration limit of interest. These expressions are used for x < 2 𝑥 2 x<2 .

$f\left(y\right)$
$y^n,n\ge 1$	$\sum _{k=0}^{\infty }\frac{B_k{x}^{k+n}}{k!\left(k+n\right)}$
$y\log y$	$x\log x-x+\sum _{k=1}^{\infty }\frac{B_k{x}^{k+1}}{k!\left(k+1\right)}\left[\log x-\frac{1}{k+1}\right]$
$y\log \left(y+a\right),a>-x$	$\sum _{k=0}^{\infty }\frac{B_k{x}^{k+1}}{k!\left(k+1\right)}\left[\log \left(x+a\right)+\frac{R(k,-x/a)}{k+1}-\frac{1}{k+1}\right]$
$1$	$\log \left(1-e^{-x}\right)$
$\log y$	$\frac{1}{2}{\log }^{2}x+\sum _{k=1}^{\infty }\frac{B_k{x}^k}{k!k}\left[\log x-\frac{1}{k}\right]$
$\log \left(y+a\right),a>0$	$\log x\log a-{\text{Li}}_2\left(-x/a\right)+\sum _{k=1}^{\infty }\frac{B_k{x}^k}{k!k}\left[\log \left(x+a\right)-\frac{x}{a\left(k+1\right)}R(k,-x/a)\right]$
	$\log \left(-x/a\right)\log \left(x+a\right)-{\text{Li}}_2\left(1+x/a\right)+\sum _{k=1}^{\infty }\frac{B_k{x}^k}{k!k}\left[\log \left(x+a\right)-\frac{x}{a\left(k+1\right)}R(k,-x/a)\right]$
$\frac{1}{y+a},a>-x$	$\frac{1}{a}\log \left(\frac{x}{x+a}\right)+\sum _{k=1}^{\infty }\frac{B_k{x}^k}{k!k}\left[\frac{1}{a}-\frac{kx}{\left(n+1\right)a^2}R(k,-x/a)\right]$
$y^{-n},n\ge 1$	$\sum _{k=0}^{n-1}\frac{B_k}{k!}\frac{x^{k-n}}{k-n}+\frac{B_n}{n!}\log x+\sum _{k=1}^{\infty }\frac{B_{k+n}}{\left(k+n\right)!}\frac{x^k}{k}$

Table 3. Table 3: Series expressions for the relevant indefinite integrals of the form shown in Eq. ( 78 ). Here, y 𝑦 y is the integration variable, and x 𝑥 x denotes the integration limit of interest. These expressions are used for x ≥ 2 𝑥 2 x\geq 2 .

$f\left(y\right)$	$P_f,x\ge 2$
${\left(y+a\right)}^n,\forall n\in \mathbb{Z},a>-x$	$\sum _{k=1}^{\infty }\frac{e^{ka}\Gamma (n+1,k\left(x+a\right))}{k^{n+1}}=\sum _{k=1}^{\infty }\frac{e^{ka}{E}_{-n}\left(k\left(x+a\right)\right)}{{\left(x+a\right)}^{-k-1}}$
$y\log \left(y+a\right),a>-x$	$\sum _{k=1}^{\infty }\frac{e^{ka}}{k^2}\left[\left(1+kx\right)e^{-k\left(x+a\right)}\log \left(x+a\right)+\left(1+kx\right)E_1\left(k\left(x+a\right)\right)+E_2\left(k\left(x+a\right)\right)\right]$
$\log \left(y+a\right),a>-x$	$\sum _{k=1}^{\infty }\frac{e^{ka}}{k}\left[e^{-k\left(x+a\right)}\log \left(x+a\right)+E_1\left(k\left(x+a\right)\right)\right]$

Table 4. Table 4: List of coefficients C n , j subscript 𝐶 𝑛 𝑗 C_{n,j} for use in Eq. ( 93 ).

$C_{n,j}$	$x_1$	$x_1^2$	$x_1^3$	$x_1^5$	$x_1^7$	$x_1^9$	$x_1^{11}$
${\beta }^0$	32/3	-16/3	8/9	-2/135	1/2835	-1/113400	1/4490640
${\beta }^2$	-64/9	0	32/27	-4/45	8/1701	-1/4860	1/124740
${\beta }^4$	-256/225	0	32/27	-296/1125	1208/42525	-64/30375	389/3118500
${\beta }^6$	-832/1575	0	32/27	-1828/3375	31352/297675	-10669/850500	10267/9355500

Table 5. Table 5: List of coefficients D n , j subscript 𝐷 𝑛 𝑗 D_{n,j} for use in Eq. ( 101 ).

$D_{n,j}$	$A$	$A^2$	$A^3$	$A^5$	$A^7$	$A^9$	$A^{11}$
${\beta }^2$	-8	16/3	-10/9	7/270	-1/1260	11/453600	-13/17962560
${\beta }^4$	-164/5	256/15	-134/45	161/2700	-19/9450	359/4536000	-289/89812800
${\beta }^6$	-40676/315	208/3	-1312/105	4651/18900	-416/59535	989/4536000	-173/22453200

Table 6. Table 6: A list of the important definitions used in DarkHistory . In this table, all spectra are discretized spectra as described in Sec.  III.2 . Spectra without overlines are normalized so that their entries contain number (per baryon) of particles produced in a redshift step. A primed energy denotes the energy of an injected particle, and by energy we mean kinetic energy. In this table, when we refer to electrons we will always mean electrons plus positrons.

Category	Symbol	Definition
General	$y$	log-redshift, $y\equiv \log \left(1+z\right)$.
	$\Delta y$, $\Delta t$	log-redshift step size and associated time step size.
	$𝐱$	Ionization levels: $𝐱\equiv (x_{\text{HII}},x_{\text{HeII}},x_{\text{HeIII}})\equiv (n_{\text{HII}}/n_{\text{H}},n_{\text{HeII}}/n_{\text{H}},n_{\text{HeIII}}/n_{\text{H}})$ i.e. the fractional abundance of ionized hydrogen atoms, singly-ionized helium atoms and doubly-ionized helium atoms with respect to the number of hydrogen atoms (both neutral and ionized).
	${\zeta }_i$	$\text{arctanh}\left[\left(2/{\chi }_i\right)\left(n_i/n_{\text{H}}-{\chi }_i/2\right)\right]$ where $i\in \{\text{HII},\text{HeII},\text{HeIII}\}$, convenient reparametrization of $𝐱$ introduced for numerical purposes.
	$T_m$	Temperature of the IGM.
	$T_m^{\left(0\right)}$, $x_{\text{HII}}^{\left(0\right)}\left(z\right)$	Baseline temperature and ionization histories, obtained from Eq. (1).
	$m_{\chi }$, $\tau$, $⟨\sigma v⟩$	Dark matter mass, lifetime, and velocity-averaged annihilation cross section.
	${\left(\frac{dE}{dVdt}\right)}_{\text{inj}}$	Energy injection rate per volume for exotic forms of energy injection, given for dark matter annihilation/decay in Eq. (3).
Spectra	$G\left(z\right)$	Conversion factor between the rate of injected events per volume to the number of injected events per baryon within a log-redshift step, as defined in Eq. (10).
	${\overline{𝐍}}_{\text{inj}}^{\alpha }\left[E_{\alpha ,i}^{\prime }\right]$	Spectrum containing number of particles of type $\alpha \in \{\gamma ,e\}$ injected into energy bin $E_{\alpha ,i}^{\prime }$ per annihilation event.
	${𝐍}_{\text{inj}}^{\alpha }[E_{\alpha ,i}^{\prime },y^{\prime }]$	Spectrum containing the number of particles per baryon in a log-redshift step of type $\alpha$ injected into energy bin $E_{\alpha ,i}^{\prime }$ at log-redshift $y^{\prime }$, as defined in Eq. (9).
	${\overline{𝐍}}_{\text{pos}}^{\gamma }\left[E_{\gamma ,i}^{\prime }\right]$	Spectrum of photons produced from a single positronium annihilation event.
	${𝐍}_{\text{new}}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$	Sum of the spectra of primary injected photons, and secondary photons produced by the cooling of electrons, as defined in Eq. (23).
	${𝐍}_{\text{prop}}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$	Spectrum of propagating photons with energies greater than $13.6$ eV that do not photoionize or get otherwise deposited into low-energy photons.
	${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$	${𝐍}_{\text{prop}}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]+{𝐍}_{\text{new}}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$, as defined in Eq. (28).
	${𝐍}_{\text{low}}^{\alpha }[E_{\alpha ,i},y]$	Low-energy photons ($\alpha =\gamma$) or electrons ($\alpha =e$) at log-redshift $y$.
Photon Cooling	${\overline{𝖯}}^{\gamma }[E_{\gamma ,i}^{\prime },E_{\gamma ,j},y^{\prime },\Delta y,𝐱]$	Transfer function for propagating photons, which multiplies ${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$ and produces ${𝐍}_{\text{prop}}^{\gamma }[E_{\gamma ,j},y^{\prime }-\Delta y]$, as defined in Eq. (27).
	${\overline{𝖣}}^e[E_{\gamma ,i}^{\prime },E_{e,j},y^{\prime },\Delta y,𝐱]$	Low-energy electron deposition transfer function, which multiplies ${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$ and produces ${𝐍}_{\text{low}}^e[E_{e,j},y^{\prime }-\Delta y]$, as defined in Eq. (29).
	${\overline{𝖣}}^{\gamma }[E_{\gamma ,i}^{\prime },E_{\gamma ,j},y^{\prime },\Delta y,𝐱]$	Low-energy photon deposition transfer function, which multiplies ${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$ and produces ${𝐍}_{\text{low}}^{\gamma }[E_{\gamma ,j},y^{\prime }-\Delta y]$, as defined in Eq. (31).
	${\overline{𝖣}}_c^{\text{high}}[E_{\gamma ,i}^{\prime },y^{\prime },\Delta y,𝐱]$	High-energy deposition transfer matrix, which multiplies ${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$ and returns the total energy that greater than $3\mathrm{keV}$ electrons produce during the cooling process deposit into channel $c\in \{$‘ion’, ‘exc’, ‘heat’$\}$, as defined in Eq. (32) in the next log-redshift step at $y^{\prime }-\Delta y$.
	${\left({\overline{𝖯}}_{1/2}^{\gamma }\right)}^n$	Coarsened propagating photon transfer function with a coarsening factor of $n$, as defined in Eq. (33), which multiplies ${𝐍}^{\gamma }[E_{\gamma ,i}^{\prime },y^{\prime }]$ and produces ${𝐍}_{\text{prop}}^{\gamma }[E_{\gamma ,j},y^{\prime }-n\Delta y]$.
Electron Cooling	$\overline{𝖭}[E_{e,i}^{\prime },E_{e,j}]$	Spectrum of secondary electrons produced due to the cooling of a single injected electron with initial energy $E_{e,i}^{\prime }$.
	${\overline{𝐑}}_c\left[E_{e,i}^{\prime }\right]$	High-energy deposition vector containing the total energy deposited into channel $c\in${‘ion’, ‘exc’, ‘heat’} by a single injected electron with kinetic energy $E_{e,i}^{\prime }$, as defined in Eq. (13).
	${\overline{𝐑}}_{\text{CMB}}\left[E_{e,i}^{\prime }\right]$	Total initial energy of CMB photons that are upscattered via ICS due to the cooling of a single electron of energy $E_{e,i}^{\prime }$.
	${\overline{𝖳}}_{\text{ICS,0}}[E_{e,i}^{\prime },E_{\gamma ,j}]$	Spectrum of photons produced with energy $E_{\gamma ,j}$ due to the cooling of a single electron of energy $E_{e,i}^{\prime }$, as defined in Eq. (14).
	${\overline{𝖳}}_{\text{ICS}}[E_{e,i}^{\prime },E_{\gamma ,j}]$	The same as ${\overline{T}}_{\text{ICS,0}}$, but with the pre-scattering spectrum of upscattered CMB photons subtracted out, as defined in Eq. (16).
	${\overline{𝖳}}_e[E_{e,i}^{\prime },E_{e,j}]$	Low-energy electron spectrum produced due to the cooling of a single electron of energy $E_{e,i}^{\prime }$, as defined in Eq. (17).
Low-energy deposition	$f_c(z,𝐱)$	Ratio of deposited to injected energy, as a function of redshift $z$ and the ionization level $𝐱$, into channels $c\in \{$‘H ion’, ‘He ion’, ‘exc’, ‘heat’, ‘cont$\}$, as defined in Eq. (4).
	${\left(\frac{d{E}^{\alpha }}{dVdt}\right)}_c$	Energy deposited per volume and time by low-energy photons ($\alpha =\gamma$) or electrons ($\alpha =e$) into channel $c$.
	$E_c^{\text{high}}\left[y\right]$	Total amount of high-energy deposition into channels $c\in${‘ion’, ‘exc’, ‘heat’} at log-redshift $y$.

Equations

$${\dot{T}}_m^{(0)}$$

$${\dot{x}}_{\text{HII}}^{(0)}$$

$${\Gamma }_C=\frac{x_e}{1+{\mathcal{F}}_{\text{He}}+x_e}\frac{8{\sigma }_T{a}_r{T}_{\text{CMB}}^4}{3m_e},$$

$${\left(\frac{dE}{dVdt}\right)}^{\text{inj}}=\{\begin{array}{ll}{\rho }_{\chi ,0}^2(1+z{)}^6⟨\sigma v⟩/m_{\chi },& \text{annihilation},\\ {\rho }_{\chi ,0}(1+z{)}^3/\tau ,& \text{decay},\end{array}$$

$${\left(\frac{dE}{dVdt}\right)}_c^{\text{dep}}=f_c(z,\mathbf{x}){\left(\frac{dE}{dVdt}\right)}^{\text{inj}},$$

$${\dot{T}}_m^{\text{inj}}$$

$${\dot{x}}_{\text{HII}}^{\text{inj}}$$

$${\dot{T}}_m$$

$${\dot{x}}_{\text{HII}}$$

$$\frac{dN}{dE}(E_i,A_j,B_k,\cdots )\approx \mathsf{S}[E_i,A_j,B_k,\cdots ].$$

$$\mathsf{N}[E_i,A_j,B_k,\cdots ]\equiv \mathsf{S}[E_i,A_j,B_k,\cdots ]\times E_i\Delta \log E_i,$$

$${\mathbf{N}}_{\text{inj}}^{\alpha }[E_i^{\prime },z]={\overline{\mathbf{N}}}_{\text{inj}}^{\alpha }[E_i^{\prime }]{\left(\frac{dN}{dVdt}\right)}^{\text{inj}}G(z),$$

$$G(z)\equiv \frac{\Delta \log (1+z)}{n_B(z)H(z)},$$

$$y\equiv \log (1+z),$$

$$R_c(E^{\prime })=\int dE{R}_c(E)\frac{dN}{dE}+P_c(E^{\prime }).$$

$${\overline{\mathbf{R}}}_c[E_i^{\prime }]=\sum _j\overline{\mathsf{N}}[E_i^{\prime },E_j]{\overline{\mathbf{R}}}_c[E_j]+{\overline{\mathbf{P}}}_c[E_i^{\prime }],$$

$${\overline{\mathsf{T}}}_{\text{ICS},0}[E_{e,i}^{\prime },E_{\gamma ,j}]$$

$$+{\overline{\mathsf{N}}}_{\text{ICS}}[E_{e,i}^{\prime },E_{\gamma ,j}],$$

$${\overline{\mathbf{R}}}_{\text{CMB}}[E_i^{\prime }]=\sum _j\overline{\mathsf{N}}[E_i^{\prime },E_j]{\overline{\mathbf{R}}}_{\text{CMB}}[E_j]+{\overline{\mathbf{P}}}_{\text{CMB}}[E_i^{\prime }],$$

$${\overline{\mathsf{T}}}_{\text{ICS}}[E_{e,i}^{\prime },E_{\gamma ,j}]$$

$$-{\overline{\mathbf{R}}}_{\text{CMB}}[E_{e,i}^{\prime }]{\overline{\mathbf{N}}}_{\text{CMB}}[E_{\gamma ,j}],$$

$${\overline{\mathsf{T}}}_e[E_{e,i}^{\prime },E_{e,j}]$$

$$+{\overline{\mathsf{N}}}_{\text{low}}[E_{e,i}^{\prime },E_{e,j}],$$

$$\frac{\tilde{E}[E_i^{\prime }]}{E_i^{\prime }}{\overline{\mathbf{R}}}_c[E_i^{\prime }]=\sum _j\tilde{\mathsf{N}}[E_i^{\prime },E_j]{\overline{\mathbf{R}}}_c[E_j]+{\overline{\mathbf{P}}}_c[E_i^{\prime }],$$

$$\tilde{\mathsf{N}}[E_i^{\prime },E_j]$$

$$\tilde{E}[E_i^{\prime }]$$

$$+\sum _j{\overline{\mathsf{T}}}_{\text{ICS}}[E_i^{\prime },E_{\gamma ,j}]E_{\gamma ,j}.$$

$${\mathbf{N}}_{\text{pos}}^{\gamma }[E_i]=\frac{1}{2}{\overline{\mathbf{N}}}_{\text{pos}}^{\gamma }[E_i]\sum _j{\mathbf{N}}_{\text{inj}}^e[E_j^{\prime }],$$

$$\frac{d{N}_{\text{new}}^{\gamma }}{d{E}_j^{\prime }dy}(E_j^{\prime })\times E_j^{\prime }\log \Delta E_j^{\prime }\times \Delta y\approx {\mathbf{N}}_{\text{new}}^{\gamma }[E_j^{\prime }],$$

$${\mathbf{N}}_{\text{new}}^{\gamma }[E_j^{\prime }]$$

$$+\sum _i{\overline{\mathsf{T}}}_{\text{ICS}}[E_{e,i}^{\prime },E_j^{\prime }]{\mathbf{N}}_{\text{inj}}^e[E_{e,i}^{\prime }].$$

$${\frac{d{N}_{\text{prop}}^{\gamma }}{dE}|}_y$$

$$+\int d{E}^{\prime }{\int }_y^{y^{\prime }}d\eta {\overline{P}}^{\gamma }(E^{\prime },E,\eta ,y){\frac{d{N}_{\text{new}}^{\gamma }}{d{E}^{\prime }d\eta }|}_{\eta }.$$

$${\overline{\mathsf{P}}}^{\gamma }[E_i^{\prime },E_j,y^{\prime },\Delta y]E_i^{\prime }\Delta \log E_i^{\prime }$$

Full text

††preprint: MIT-CTP/5113

I Introduction

Dark matter annihilation or decay and other exotic sources of energy injection can significantly alter the ionization and temperature histories of the universe. In this paper we describe a new public code package, DarkHistory, that allows fast and accurate computation of these possible effects of exotic energy injection on astrophysical and cosmological observables.

In particular, we will focus on interactions that allow dark matter (DM) to decay or annihilate into electromagnetically interacting Standard Model particles. This case has been studied extensively in the literature: stringent constraints on the dark matter annihilation cross section and decay lifetime have been derived from the way these Standard Model products would distort the anisotropies of the cosmic microwave background (CMB) Slatyer et al. (2009); Slatyer (2013); Slatyer and Wu (2017); Kanzaki et al. (2010), or increase the temperature of the Inter-Galactic Medium (IGM), consequently affecting 21-cm and Lyman-$\alpha$ line emission Liu et al. (2016); Lopez-Honorez et al. (2016); Liu and Slatyer (2018); Diamanti et al. (2014).

DarkHistory facilitates the calculation of these observables and the resulting constraints. In particular, DarkHistory makes the temperature constraint calculations significantly more streamlined, self-consistent, and accurate. It has a modular structure, allowing users to easily adjust individual inputs to the calculation – e.g. by changing the reionization model, or the spectrum of particles produced by dark matter annihilation/decay. Compared to past codes developed for such analyses Stöcker et al. (2018), DarkHistory has a number of important new features:

• •

the first fully self-consistent treatment of exotic energy injection. Exotic energy injections can modify the evolution of the IGM temperature $T_{\text{IGM}}$ and free electron fraction $x_{e}$, and previously this modification has been treated perturbatively, assuming the backreaction effect due to these modifications on the cooling of injected particles is negligible. This assumption can break down toward the end of the cosmic dark ages for models that are not yet excluded Liu et al. (2016). DarkHistory solves simultaneously for the temperature and ionization evolution and the cooling of the injected particles, avoiding this assumption;

• •

a self-contained treatment of astrophysical sources of heating and reionization, allowing the study of the interplay between exotic and conventional sources of energy injection;

• •

a large speed-up factor for computation of the full cooling cascade for high-energy injected particles (compared to the code employed in e.g. Liu et al. (2016)), via pre-computation of the relevant transfer functions as a function of particle energy, redshift and ionization level;

• •

support for treating helium ionization and recombination, including the effects of exotic energy injections; and

• •

a new and more correct treatment of inverse Compton scattering (ICS) for mildly relativistic and non-relativistic electrons; previous work in the literature has relied on approximate rates which are not always accurate.

Due to these improvements, DarkHistory allows for rapid scans over many different prescriptions for reionization, either in the form of photoheating and photoionization rates, or a hard-coded background evolution for $x_{e}$. The epoch of reionization is currently rather poorly constrained, making it important to understand the observational signatures of different scenarios, and the degree to which exotic energy injections might be separable from uncertainties in the reionization model. Previous attempts to model the effects of DM annihilation and decay into the reionization epoch have typically either assumed a fixed ionization history Stöcker et al. (2018) – requiring a slow re-computation of the cooling cascade if that history is changed Liu et al. (2016) – or made an approximation for the effect of a modified ionization fraction on the cooling of high-energy particles Lopez-Honorez et al. (2013); Diamanti et al. (2014); Poulin et al. (2015, 2017); Lopez-Honorez et al. (2016).

Despite our emphasis on dark matter annihilation and decay, DarkHistory can be used to explore the effect of other forms of exotic particle injection. Other such possible sources include Hawking radiation from black holes Poulin et al. (2017); Clark et al. (2018), radiation from accretion onto black holes Hektor et al. (2018a), and processes from new physics such as de-excitation of dark matter or decay of meta-stable species Hektor et al. (2018b).

In Section II we review the physics of the ionization and temperature evolution, in the context of the three-level-atom (TLA) approximation, including the possibility of exotic energy injections. In Section III we discuss the overall structure of DarkHistory, which self-consistently combines the TLA evolution of the ionization and gas temperature with the cooling of particles injected by exotic processes. This section also describes the implementation of various physical processes in the code, in particular the treatment of cooling and production of secondaries by electrons and photons. In Section IV we relate these processes to the various modules of DarkHistory, before providing a number of worked examples in Section V. We present our conclusions and discuss some future directions in Section VI. We discuss our improved treatment of ICS in detail in Appendix A, provide the photon spectra from positronium annihilation in Appendix B, discuss a series of cross checks in Appendix C, and provide a table of definitions used throughout this paper in Appendix D.

II Ionization and Thermal Histories

DarkHistory computes the ionization and temperature evolution of the universe in the presence of an exotic source of energy injection, such as dark matter annihilation or decay, using a modified version of the three-level atom (TLA) model for both hydrogen and helium, based on RECFAST Seager et al. (2000, 1999). The reader may refer to Ref. Ali-Haimoud and Hirata (2011) for a detailed derivation of the unmodified TLA equations with hydrogen only, and Refs. Seager et al. (2000, 1999); Wong et al. (2008) for the treatment of helium recombination in RECFAST. In this section, we will neglect the evolution of helium for simplicity, leaving a detailed discussion of our treatment of helium to Sec. III.7.1.

In the absence of any source of energy injection, the TLA model, first derived in Peebles (1968); Zeldovich et al. (1969), provides a pair of coupled differential equations for the matter temperature in the IGM and the hydrogen ionization fraction:

[TABLE]

where $H$ is the Hubble parameter, $n_{\text{H}}$ is the total number density of hydrogen (both neutral and ionized), $x_{\text{HII}}\equiv n_{\text{HII}}/n_{\text{H}}$ where $n_{\text{HII}}$ is the number density of free protons, $x_{e}\equiv n_{e}/n_{\text{H}}$ is the free electron fraction with $n_{e}$ being the free electron density, and $E_{21}=$10.2\text{,}\mathrm{e\kern-0.50003ptV}$$ is the Lyman-$\alpha$ transition energy. $T_{m}$ and $T_{\text{CMB}}$ are the temperatures of the IGM and the CMB respectively.111We follow the standard astrophysical convention in which H and H+ are denoted HI and HII, while He, He+ and He2+ are denoted HeI, HeII and HeIII respectively. $\alpha_{\text{H}}$ and $\beta_{\text{H}}$ are case-B recombination and photoionization coefficients for hydrogen respectively,222The value of $\beta_{\text{H}}$ used in DarkHistory includes the constant and gaussian fudge factors used by version 1.5.2 of RECFAST. and $\mathcal{C}$ is the Peebles-C factor that represents the probability of a hydrogen atom in the $n=2$ state decaying to the ground state before photoionization can occur Peebles (1968); Ali-Haimoud and Hirata (2011). The photoionization coefficient is evaluated at the radiation temperature, $T_{\text{CMB}}$, in agreement with Ref. Chluba et al. (2015). $\Gamma_{C}$ is the Compton scattering rate, given by

[TABLE]

where $\sigma_{T}$ is the Thomson cross section, $a_{r}$ is the radiation constant, $m_{e}$ is the electron mass, and $\mathcal{F}_{\text{He}}\equiv n_{\text{He}}/n_{\text{H}}$ is the relative abundance of helium nuclei by number. In the absence of helium, note that $x_{e}=x_{\text{HII}}$. The solutions to Eq. (1) — i.e. without any sources of energy injection — define what we will call the baseline temperature and ionization histories, $T_{m}^{(0)}(z)$ and $x_{\text{HII}}^{(0)}(z)$.

Exotic sources may inject additional energy into the universe, altering the thermal and ionization evolution shown in Eq. (1). For example, the rate of energy injection from DM annihilating with some velocity averaged cross section $\langle\sigma v\rangle$, or decaying with some lifetime $\tau$ much longer than the age of the universe, is given by

[TABLE]

where $\rho_{\chi,0}$ is the mass density of DM today, and $m_{\chi}$ is the DM mass. This injected energy, however, does not in general manifest itself instantaneously as ionization, excitation, or heating of the gas. Instead, the primary particles injected into the universe may cool over timescales significantly larger than the Hubble time, producing secondary photons that may redshift significantly before depositing their energy into the gas.

Although the primary particles injected into the universe may be any type of Standard Model particle, we will only need to consider the cooling of photons and electron/positron pairs Slatyer et al. (2009). This simplification occurs because either the primaries are stable particles like photons, electrons and positrons, neutrinos, protons and anti-protons, and heavier nuclei, or are unstable particles that resolve into these particles on time scales much shorter than the cosmological time scales under consideration. For typical sources of energy injection we can neglect heavier nuclei because they are produced in negligible amounts, and neutrinos because they are very ineffective at depositing their energy. Protons and antiprotons generally form a subdominant component of stable electromagnetic particles across all possible Standard Model primaries Cirelli et al. (2011), and deposit energy less effectively than electrons, positrons, and photons (although their effects are not completely negligible Weniger et al. (2013)). We therefore only decompose the injection of any primary into an effective injection of photons, electrons, and positrons, in accordance with Ref. Slatyer et al. (2009) and subsequent works. Adding the contribution from protons and antiprotons may strengthen these constraints by a small amount.

A significant amount of work has been done on computing the cooling of high energy photons, electrons, and positrons Slatyer et al. (2009); Furlanetto and Stoever (2010); Valdes et al. (2010); Slatyer (2013); Evoli et al. (2012); Galli et al. (2013); Evoli et al. (2014); Slatyer (2016); Kanzaki and Kawasaki (2008); Kawasaki et al. (2016). Once the cooling of injected primary particles is determined, the energy deposited into channel $c$ (hydrogen ionization, excitation, or heating) can be simply parametrized as

[TABLE]

with all of the complicated physics condensed into a single numerical factor that is dependent on the redshift and the ionization fractions of all of the relevant species in the gas, which we denote $\mathbf{x}\equiv(x_{\text{HII}},x_{\text{HeII}},x_{\text{HeIII}})$. When helium is neglected, the ionization dependence of these $f_{c}$ functions simplifies to a dependence on $x_{\text{HII}}=x_{e}$. These $f_{c}$ functions also depend on the energies and species of the injected particles, but for simplicity of notation we will not write these arguments explicitly.

The effect of energy injection on the thermal and ionization history can now be captured by additional source terms,

[TABLE]

where $\mathcal{R}=$13.6\text{,}\mathrm{e\kern-0.50003ptV}$$ is the ionization potential of hydrogen.

Prior to this work, $f_{c}(z,\mathbf{x})$ has largely been computed assuming the standard ionization history computed by recombination codes $\mathbf{x}_{\text{std}}(z)$, essentially making $z$ the only independent variable of $f_{c}$ as a function. These calculations are therefore applicable only so long as any perturbations to the assumed ionization history (e.g. by additional sources of energy injection) are sufficiently small. This is generally a good approximation near recombination: at these redshifts, the ionization history is well-constrained by CMB power spectrum measurements, and therefore large perturbations to $x_{e}$ are highly disfavored. For $z\lesssim 100$, however, ionization levels that exceed the standard value of $x_{e}\sim 2\times 10^{-4}$ by several orders of magnitude are experimentally allowed Liu et al. (2016). Moreover, star formation during the process of reionization rapidly ionizes and heats the universe at $z\lesssim 20$, causing the ionization and thermal history to diverge from the baseline histories.

The primary effect of an increase in ionization levels is to decrease the number of neutral hydrogen and helium atoms available to ionize, decreasing the fraction of injected power that goes into ionization of these species; on the other hand, increasing $x_{e}$ increases the number of charged particles available for low-energy electrons to scatter off and heat the IGM, increasing the fraction of power going into heating. Since energy injection processes generally increase $x_{e}$ with time, the power into heating increases at an accelerated rate at late times, making a proper calculation of $f_{c}(z,\mathbf{x})$ crucial for an accurate computation of the temperature history.

Computing the full $\mathbf{x}$-dependence of $f_{c}(z,\mathbf{x})$ also allows us to perform, for the first time, a consistent calculation of the temperature and ionization histories with both exotic energy injection processes and reionization. At the onset of reionization, stars begin to form, and the ionizing radiation emitted by these objects injects a large amount of energy into the IGM. There remains a large degree of uncertainty regarding how reionization proceeds, but given some model for the photoionization and photoheating rates, and including other important energy transfer processes such as collisional ionization and excitation, additional terms $\dot{T}_{m}^{\text{re}}$ and $\dot{x}_{\text{HII}}^{\text{re}}$ (as well as the corresponding terms for helium) can be included in Eq. (1) to model reionization. These terms are discussed in much greater detail in Sec. III.7.

To summarize, DarkHistory computes the ionization and thermal history in the presence of exotic sources of energy injection, with the evolution equations in the absence of helium given by

[TABLE]

In the rest of the paper, we will describe how we calculate the inputs required to integrate these equations, i.e. $f_{c}(z,\mathbf{x})$, $\dot{T}_{m}^{\text{re}}$, $\dot{x}_{\text{HII}}^{\text{re}}$ and the modifications necessary to include helium.

III Code Structure and Content

In this section we discuss the structure and physics content of the DarkHistory package.

III.1 Overview

Fig. 1 shows a flowchart depicting the overall structure of DarkHistory. The overall goal of the code is to take in some injected spectrum of photons and electron/positron pairs at a given redshift, and partition the energy into several categories as they lose their energy over a small redshift step:

High-energy deposition. This is the total amount of energy deposited into ionization, excitation and heating by any high-energy (above $3\text{\,}\mathrm{ke\kern-0.50003ptV}$) electron generated during any of the cooling processes; 2. 2.

Low-energy electrons. These are electrons that have kinetic energy below $3\text{\,}\mathrm{ke\kern-0.50003ptV}$ where atomic cooling processes typically dominate over ICS after recombination. These electrons are separated out at each step in order to treat their energy deposition (which occurs in a timescale much shorter than the time step) more carefully; 3. 3.

Low-energy photons. These are photons with energies below $3\text{\,}\mathrm{ke\kern-0.50003ptV}$ that either photoionize within the redshift step, or lie below $13.6\text{\,}\mathrm{e\kern-0.50003ptV}$. Such photons either lose all their energy within the redshift step, or cool only through redshifting, and thus can be treated in a simplified manner; and 4. 4.

Propagating photons. These are photons that are present at the end of the redshift step and are not included in the low-energy photons category.

Throughout the paper, we use the word “electrons” to refer to both electrons and positrons. Although the interactions of electrons and positrons with the gas differ, the ICS cross-sections are identical, and ICS dominates the energy losses down to energy scales where the positron is nonrelativistic Ferrara et al. (2012). For nonrelativistic positrons, their mass energy is converted into photons through annihilation with electrons. Since the positron mass is much larger than the kinetic energy in this regime, neglecting differences in kinetic energy loss between electrons and positrons is unlikely to be important. In a future version of DarkHistory we plan to include a more sophisticated treatment of low energy electrons and positrons.

The outputs in the first three categories are used to compute the evolution of the ionization and temperature history at this redshift step, before the code moves on to the next step and performs the same calculation again. A brief description of a step in this loop is as follows:

Input. Before the code begins, the user specifies a DM energy injection model or some other redshift-dependent energy injection rate, as well as the photon and $e^{+}e^{-}$ spectra produced per energy injection event. By default, DarkHistory starts from an initial redshift of $1+z=3000$, ensuring that the spectra of particles present at and after recombination (at $z\sim 1000$) are accurate. Details are provided in Sec. III.3. Inputs to the code are provided to the function evolve() found in the module main; some tools for obtaining spectra from an arbitrary injection of Standard Model particles can be found in the pppc module; 2. 2.

Injected electron cooling. Injected electrons (and positrons) cool through a combination of atomic processes and ICS. Transfer functions that map these injected electrons to high-energy deposition, secondary photons from ICS and positron annihilation, and low-energy electrons are computed and applied to the injected electrons. A discussion of these calculations can be found in Sec. III.4 and in the electrons module of the code.

The sum of the secondary photons produced by electron cooling, photons injected on this timestep, and propagating photons from the previous timestep are used as input to the photon cooling transfer functions, which we describe next; 3. 3.

Photon propagation and energy deposition. At this stage, we have a spectrum of photons that can undergo a range of cooling processes to lose their energy over this redshift step. The effect of these cooling processes on the photon spectrum can be reduced to three transfer functions that we will describe in detail in Sec. III.5. These transfer functions have been pre-computed separately and can be downloaded at https://doi.org/10.7910/DVN/DUOUWA, together with all the other data required to run the code. These transfer functions determine how photons in this redshift step turn into propagating photons that continue on to the next redshift step, and low-energy photons and low-energy electrons that undergo further processing. All of these computations occur in the main module; 4. 4.

Calculating $f_{c}(z,\mathbf{x})$. The low-energy photons and low-energy electrons from this redshift step deposit their energy into ionization, heating and excitation of atoms, and the value of $f_{c}(z,\mathbf{x})$ at this step is computed by comparing the energy deposited in each channel to the energy injection rate for this timestep. Details of this computation are given in Sec. III.6, and can be found in the low_energy module; 5. 5.

TLA integration and reionization. With $f_{c}(z,\mathbf{x})$ at this step, we can now integrate the TLA across this redshift step. We can also include a reionization model, or track helium ionization, both of which add more terms to the TLA, as detailed in Sec. III.7. We now know the $\mathbf{x}$ and $T_{m}$ that are reached at the end of this step. These calculations are done in the history module; and 6. 6.

Next step. The $\mathbf{x}$ and $T_{m}$ values computed above are passed to the next redshift step, so that all transfer functions at the next step can be computed at the appropriate ionization level. The propagating photons found above are also passed to the next step, and the loop repeats.

Because $f_{c}(z,\mathbf{x})$ is computed by integrating the TLA at each step, and all transfer functions are evaluated at the value of $\mathbf{x}$ in the step, the backreaction of increased ionization levels is now fully accounted for.

In the next several sections, we will describe both the physics and numerical methods that go into the loop.

III.2 Discretization

Before describing in detail each part of DarkHistory, we will first describe how discretization occurs in our code, and the notation we will use throughout this paper. Typically, we will deal with some smooth spectrum of particles $dN/dE(E,A,B,\cdots)$, which is a function of the energy abscissa $E$, and several other variables that we denote here as $A,B,\cdots$. Smooth functions that are derivatives will always use ‘$d$’ to denote differentiation, and parentheses to denote functional dependence. We shall always discretize such spectra as

[TABLE]

The discretized spectrum $\mathsf{S}$ is a matrix of dimension equal to the number of variables it depends on, where $i,j,k,...$ index discrete values of these variables. Throughout this paper, we will denote vectors (quantities which depend on a single variable) by a bold typeface and matrices (quantities that depend on multiple variables) by a sans-serif typeface. Discrete steps or changes are denoted by ‘$\Delta$’, and discrete functional dependencies are written in square brackets.

$\mathsf{S}$ times the bin width should always be regarded as a matrix of number of particles inside some bin, all with energy given by $E_{i}$. This matrix is mathematically defined as

[TABLE]

where $\Delta\log E_{i}$ is the log-energy bin width. We will always take $E_{i}\Delta\log E_{i}$ to be the bin width by convention. In DarkHistory, spectra are binned into energy values that are evenly log-spaced. $E_{i}$ should be regarded as the bin center, with the bin boundaries occurring at the geometric mean of adjacent energy values, and the boundaries of the first and last bin are taken to be symmetric (in log-space) about the bin centers.

III.3 Input

To initialize the loop described above, the user must specify the discretized photon and electron spectra produced per injection event, which we denote $\overline{\mathbf{N}}^{\gamma}_{\text{inj}}[E_{j}^{\prime}]$ and $\overline{\mathbf{N}}^{e}_{\text{inj}}[E_{j}^{\prime}]$. Bars denote spectra or transfer functions that have been normalized by some process or quantity, while spectra without any markings denote a number of particles per baryon from here on, unless otherwise specified.

Given the redshift-dependent rate of injection events per volume $\left(dN/dV\,dt\right)^{\text{inj}}$ we can determine the spectrum of particles $\mathbf{N}_{\text{inj}}^{\alpha}$ injected within a log-redshift step of width $\Delta\log(1+z)$ per baryon by

[TABLE]

where $\alpha$ take on values $\gamma$ or $e$, and

[TABLE]

where $n_{B}$ is the number density of baryons. $G(z)$ converts between the rate of injection events per volume and the number of injection events per baryon in the log-redshift step.

In the following sections, we will be mostly concerned with log-redshift steps, and so it is convenient to define

[TABLE]

and likewise $\Delta y\equiv\Delta\log(1+z)$.

III.4 Injected Electron Cooling

After specifying the injected spectra, the next step of the code is to resolve the injected electron/positron pairs, $\mathbf{N}^{e}_{\text{inj}}$. High-energy electrons and positrons cool through atomic processes (collisional ionization, collisional excitation and Coulomb heating), as well as ICS off CMB photons. After losing their kinetic energy to these processes, positrons ultimately annihilate with free electrons in the IGM, producing high-energy photons. All of these processes occur within a timescale much shorter than the timesteps considered in DarkHistory. Because of this, the code converts all input high-energy electrons into energy deposited into ionization, excitation, heating, scattered photons from ICS, and low-energy electrons (below $3\text{\,}\mathrm{ke\kern-0.50003ptV}$), which we treat separately. The photons produced from ICS are added to those that are injected promptly from the DM energy injection process, as well as propagating photons from the previous step.

We will first briefly discuss our calculation of the scattered photon and electron spectra from ICS, and then move on to describe the numerical method used to compute electron cooling.

III.4.1 Inverse Compton Scattering

ICS off CMB photons is an important energy loss mechanism for electrons/positrons over a large range of energies and redshifts. The efficiency of ICS as a cooling mechanism relative to atomic cooling processes has been the subject of some confusion in the literature, with some earlier studies Valdes et al. (2010); Evoli et al. (2012) underestimating the cooling rate of the electrons. ICS becomes more important relative to atomic processes as the electron energy increases, but a correct treatment shows that even nonrelativistic electrons can have ICS as the main cooling mechanism in the early universe; at $z\sim 600$, for example, it is the primary energy loss mechanism for electrons with kinetic energy $\gtrsim$10\text{,}\mathrm{ke\kern-0.50003ptV}$$ Galli et al. (2013); Slatyer (2016). Existing work on electron cooling has focused on the highly nonrelativistic regime (electron kinetic energy below $3\text{\,}\mathrm{ke\kern-0.50003ptV}$) Furlanetto et al. (2006), where ICS is unimportant compared to atomic cooling processes, or on the relativistic regime Valdes et al. (2010); Evoli et al. (2012); Hansen and Haiman (2004).

Earlier work by one of the authors Slatyer (2013, 2016) already incorporates ICS cooling for electrons across both the Thomson and the relativistic regimes. DarkHistory improves the accuracy of the calculation in the Thomson regime by using the full expression for the spectrum of scattered photons, with no further approximation. As a result, the code is able to accurately calculate the scattered photon spectrum and the energy loss spectrum of electrons. This means that we fully cover all relevant regimes for ICS for electrons of arbitrary energy scattering off the CMB at all redshifts $z\sim 10^{9}$ and below.333Above this redshift, photons have energies comparable to the electron mass $m_{e}$, and Klein-Nishina scattering can occur between photons and non-relativistic electrons, which falls outside of the two regimes considered here. These calculations are fast and numerically stable even for nonrelativistic electrons, where conventional numerical integration can be unreliable due to the presence of catastrophic cancellations between large terms.

We leave a full discussion of how DarkHistory treats ICS to Appendix A. In summary, the code is able to compute the scattered photon and electron spectra that are produced per unit time due to ICS off the CMB across all relevant kinematic regimes. These spectra are then taken as inputs for the numerical computation of how an electron cools taking into account all processes, which is described below.

III.4.2 Numerical Method

Consider an injected electron (or positron) with kinetic energy $E^{\prime}$ (all quantities associated with injected particles throughout this paper will be denoted with ′). Let $R_{c}(E^{\prime})$ be the energy eventually deposited into some channel $c$ by this electron, once it has lost all of its initial energy. Within a short time interval $\Delta t$ (taken to be $1\text{\,}\mathrm{s}$ in our calculation), the electron undergoes all possible cooling processes with some probability, producing the (averaged) secondary electron spectrum $dN/dE$. Within this same interval $\Delta t$, some portion of the energy $P_{c}(E^{\prime})$ is also deposited promptly into the channel under consideration. The secondary electron spectrum then deposits its energy according to $R_{c}$ for energies lower than $E^{\prime}$. We can thus write the following recursive equation:

[TABLE]

Note that $R_{c}(E^{\prime})$ does not include deposition to the channel $c$ via secondary photons from ICS or positron annihilation; because the cooling times of secondary photons can be much longer than a timestep, they must be treated separately. $R_{c}(E^{\prime})$ as defined here is the “high-energy deposition” from electrons within the timestep, as described in Section III.1. The relevant channels are $c=\{\text{`ion'},\text{`exc'},\text{`heat'}\}$ for deposition into collisional ionization, collisional excitation and heating respectively. The ‘ion’ and ‘exc’ channels include ionization and excitation off all species.

As long as the time step $\Delta t$ is much shorter than the characteristic interaction timescale of all of the interactions, $dN/dE$ is simply the sum of all of the scattered electron spectra due to each process within $\Delta t$, normalized to a single injected electron. A detailed accounting of the relevant cross sections and secondary spectra is provided in Ref. Slatyer et al. (2009), and these results can be used to calculate $dN/dE$ and $P_{c}$. We will denote the discretized version of the normalized scattered electron spectra by $\overline{\mathsf{N}}$, since it is normalized to one electron.

Numerically, we would like to compute $\overline{\mathbf{R}}_{c}$, a vector containing the energy deposited into channel $c$, with each entry corresponding to a single electron with initial kinetic energy $E^{\prime}$. The overline notation serves as a reminder that the quantity is normalized to one injected electron. The discretized version of Eq. (12) reads

[TABLE]

where $\overline{\mathbf{P}}_{c}$ is the vector of the prompt energy deposition in channel $c$ per electron. This is a linear system of equations, and we can solve for each $\overline{\mathbf{R}}_{c}$ given $\overline{\mathsf{N}}$ and $\overline{\mathbf{P}}_{c}$.

A similar procedure also works for finding the ICS photon spectrum after an electron completely cools. Let the discretized spectrum be $\overline{\mathsf{T}}_{\text{ICS},0}[E_{e,i}^{\prime},E_{\gamma,j}]$, where $E^{\prime}_{e}$ is the initial electron kinetic energy, and $E_{\gamma}$ is the photon energy. Then the ICS photon spectrum produced after complete cooling of a single electron satisfies

[TABLE]

with $\overline{\mathsf{N}}_{\text{ICS}}$ being the discretized version of the scattered photon spectrum defined in Eq. (58) within $\Delta t$, and indices $e$ and $\gamma$ have been inserted to clarify the difference between electron and photon energies. This spectrum consists of CMB photons that are upscattered by the injected electron; in order to be able to track energy conservation, we also need to keep track of the initial energy of the upscattered photons. We therefore also need to solve

[TABLE]

where $\overline{\mathbf{P}}_{\text{CMB}}$ is the total initial energy of photons upscattered in $\Delta t$.444We do not have to track the photon spectrum, since the initial CMB photon energy is only significant for nonrelativistic injected electrons, which are always in the Thomson regime and hence scatter in a frequency-independent manner. For relativistic electrons, the initial CMB photon energy is neglected, as the photon is overwhelmingly upscattered to a much higher final energy. At this point, we now define $\overline{T}_{\text{ICS}}$ to be the ICS photon spectrum with the upscattered CMB spectrum subtracted out, so that $\overline{T}_{\text{ICS}}$ now represents a distortion to the CMB spectrum:

[TABLE]

where $\overline{\mathbf{N}}_{\text{CMB}}$ is the CMB spectrum normalized to unit total energy. The total of energy of $\overline{T}_{\text{ICS}}$ for each $E_{e,i}^{\prime}$ therefore gives the energy lost by the incoming electron through ICS.

Finally, the low-energy electron spectrum produced is similarly given by

[TABLE]

where $\overline{\mathsf{N}}_{\text{high}}$ ($\overline{\mathsf{N}}_{\text{low}}$) is $\overline{\mathsf{N}}$ with only high-energy (low-energy) $E_{e,k}$ included.

In DarkHistory, we choose a square matrix $\overline{\mathsf{N}}$ with the same abscissa for both injected and scattered electron energies. As a result, $\overline{\mathsf{N}}$ has diagonal values that are very close to 1, since most particles do not scatter within $\Delta t$. Because of this, we find that it is numerically more stable to solve the equivalent equation

[TABLE]

where

[TABLE]

The variables $\widetilde{\mathsf{N}}$ and $\widetilde{E}$ are simply the number of electrons and total energy excluding electrons that remained in the same energy bin after $\Delta t$. Eqs. (16) and (17) can be similarly transformed in the same way as Eq. (18) and solved. Since $\widetilde{N}$ is a triangular matrix, the SciPy function solve_triangular() is used for maximum speed.555The upscattering of electrons during ICS is negligible: see Appendix A for more details.

Having calculated $\overline{\mathbf{R}}_{c}$, $\overline{\mathsf{T}}_{\text{ICS}}$ and $\overline{\mathsf{T}}_{e}$, all normalized to a single electron, the final result when an arbitrary electron spectrum $\mathbf{N}^{e}_{\text{inj}}[E^{\prime}_{e,i}]$ completely cools is simply given by contracting these quantities with $\mathbf{N}^{e}_{\text{inj}}$. Note that all of these quantities are also dependent on redshift: we have simply suppressed this dependence for notational simplicity in this section.

Finally, after positrons have lost all of their kinetic energy, they are assumed to form positronium and annihilate promptly, producing a gamma ray spectrum that also gets added to the propagating photon spectrum. The positronium spectrum is given simply by

[TABLE]

where $\overline{\mathbf{N}}^{\gamma}_{\text{pos}}$ is the positronium annihilation spectrum normalized to a single positron, shown in Appendix B. The factor of $1/2$ accounts for the fact that $\mathbf{N}^{e}_{\text{inj}}$ contains both electrons and positrons in equal number.

Since all calculated quantities depend on $z$ and $\mathbf{x}$, all quantities discussed in this section have to be computed at each redshift step. This allows us to properly capture the effect of changing ionization levels on the energy deposition process.

III.5 Photon Propagation and Energy Deposition

After resolving the injected electrons and obtaining the photons produced from their cooling, the spectrum of photons that have been newly injected per baryon per log-redshift can be discretized as

[TABLE]

where $\mathbf{N}^{\gamma}_{\text{new}}$ is the sum of photons injected directly by the injection event, and photons produced by the cooling of injected electrons, i.e.

[TABLE]

These photons can cool through a number of processes, including redshifting, pair production, Compton scattering and photoionization. Within a particular log-redshift step, low-energy photons and low-energy electrons are produced, and some high-energy deposition from high-energy electrons produced by $\mathbf{N}^{\gamma}_{\text{new}}$ occur. On the other hand, some part of the photon spectrum lies above $13.6\text{\,}\mathrm{e\kern-0.50003ptV}$ and does not photoionize within the log-redshift step; instead, these photons propagate forward to the next step.

The resulting deposition into low-energy photons and electrons was used to compute $f_{c}$ in Ref.  Slatyer (2016), assuming the fixed baseline ionization history. In order to capture the dependence on ionization history, however, we need to be able to calculate the propagation and deposition processes at any ionization level, redshift and injected particle energy.

One of the main ideas of DarkHistory is to capture the photon cooling processes as precomputed transfer functions with injection energy, redshift and ionization levels as the dependent variables. These transfer functions then act on some incoming spectrum and produce a spectrum of propagating particles, a spectrum of deposited particles or some amount of deposited energy within a log-redshift step. These transfer functions can be evaluated at various points in injection energy, redshift, and ionization levels, and interpolated at other points. With a given injection model, we can then string together these transfer functions to work out the propagation of photons and the deposition of energy, over an extended redshift range, given any exotic source of energy injection.

III.5.1 Propagating Photons

Consider a spectrum of photons per baryon denoted $dN^{\gamma}/dE^{\prime}$ that is present in the universe at some log-redshift $y$. As these photons propagate, various cooling processes result in these photons being scattered into energies below $13.6\text{\,}\mathrm{e\kern-0.50003ptV}$, or they may photoionize on an atom in the gas. Those particles that do not undergo either process within a redshift step are called “propagating photons”, and continue to propagate into the next redshift step.

We define the transfer function for propagating photons $\overline{P}^{\gamma}(E^{\prime},E,y^{\prime},y)$ through the following relation:

[TABLE]

$\overline{P}^{\gamma}$ takes a spectrum of photons that are present at $y^{\prime}$ and propagates them forward to a spectrum of propagating photons at $y$. $\overline{P}^{\gamma}(E^{\prime},E,y^{\prime},y)$ is exactly the number of propagating photons per unit energy that results from a single photon injected at log-redshift $y^{\prime}$ with energy $E^{\prime}$ cooling until log-redshift $y$. The $\overline{P}^{\gamma}$ functions are calculated separately using the code described in Ref. Slatyer et al. (2009); Slatyer (2016).

We distinguish between two different sources of photons between two redshifts $y^{\prime}$ and $y$ (with $y^{\prime}>y$): propagating photons at $y^{\prime}$, $dN_{\text{prop}}^{\gamma}/dE^{\prime}$, and the newly injected photons between the redshifts $y^{\prime}$ and $y$, defined in discretized form in Eq. (23). With these sources, we can write the spectrum of propagating photons at $y$ as

[TABLE]

We discretize this expression by defining the following discrete quantities according to the conventions set down in Eqs. (7) and (8):

[TABLE]

where we have chosen some fixed value of $\Delta y$, so that the final redshift is $y=y^{\prime}-\Delta y$. In DarkHistory, the default value is $\Delta y=10^{-3}$, although this can be adjusted by the process of coarsening, described in Sec. III.5.3. Dropping the dependence on $\Delta y$ for simplicity, the discretized version of Eq. (25) reads

[TABLE]

where we have defined

[TABLE]

III.5.2 Energy Deposition

Aside from $\overline{\mathsf{P}}^{\gamma}$, we also have three deposition transfer functions describing the energy losses of $\mathbf{N}^{\gamma}$ into high-energy deposition, low-energy electrons and low-energy photons, as defined in Sec. III.1. These transfer functions are defined by their action on the discretized photon spectrum, $\mathbf{N}^{\gamma}$, and are discretized in a similar manner.

The low-energy electron deposition transfer matrix, $\overline{\mathsf{D}}^{e}$, yields the low-energy electrons produced via cooling of $\mathbf{N}^{\gamma}$. Adding the low-energy electrons produced directly from the injected electrons $\mathbf{N}^{e}_{\text{inj}}$, we obtain the full low-energy electron spectrum $\mathbf{N}^{e}_{\text{low}}[E_{j},y]$ at a particular redshift step:

[TABLE]

where

[TABLE]

while the deposition transfer matrix $\overline{\mathsf{D}}^{\gamma}$ yields the low-energy photons,

[TABLE]

$\mathbf{N}^{\gamma}_{\text{low}}$ is computed as a distortion to the CMB spectrum, with $\overline{\mathsf{D}}^{\gamma}$ computed with the initial spectrum of upscattered CMB photons subtracted, in the same way as $\overline{T}_{\text{ICS}}$, as shown in Eq. (16).

As the propagating photons cool over a single log-redshift step, they generate high-energy electrons along the way. These are handled in a similar manner to injected high-energy electrons as described in Sec. III.4, but instead of performing the calculation at each step, we simply provide transfer functions $\overline{\mathbf{D}}_{\text{c}}^{\text{high}}$ that act on propagating photons and return the high-energy deposition into the channels $c=${‘ion’, ‘exc’, ‘heat’}.666For legacy reasons, DarkHistory actually computes the transfer function that returns the high-energy deposition per second; this is just a difference in convention. We can then combine this with the result from electron cooling to obtain the high-energy deposition per baryon within a log-redshift step into each channel $c$:

[TABLE]

To summarize, we have defined the following transfer functions: $\overline{\mathsf{P}}^{\gamma}$ for propagating photons, and $\overline{\mathsf{D}}^{\gamma}$, $\overline{\mathsf{D}}^{e}$ and $\overline{\mathbf{D}}^{\text{high}}_{c}$ for deposition into low-energy photons, low-energy electrons and high-energy deposition channels respectively. These transfer functions act on the spectrum of photons $\mathbf{N}^{\gamma}$ (from both the injection source and the cooling of injected electrons). Together with the transfer functions for the cooling of injected electrons, we have all the information needed to propagate injected particles and compute their energy deposition as a function of redshift.

III.5.3 Coarsening

The propagating photons transfer function $\overline{\mathsf{P}}^{\gamma}$ can always be evaluated with the same input and output energy abscissa, so that the 2D transfer matrix at each $y$ is square. If the transfer function $\overline{\mathsf{P}}^{\gamma}$ does not vary significantly over redshift, then in the interest of computational speed, we can make the following approximation of Eq. (27) for propagation transfer matrices:

[TABLE]

where repeated indices are summed. $i$ and $j$ index input and output energies, and $\overline{\mathsf{P}}^{\gamma}_{1/2}$ is $\overline{\mathsf{P}}^{\gamma}$ evaluated at log-redshift $y-n\Delta y/2$ to minimize interpolation error. When making this approximation, we also have to ensure that we redefine

[TABLE]

for both channels $\alpha=e$ and $\gamma$, so that we (approximately) include all of the particles injected between $y$ and $y-n\Delta y$.

Likewise, if both the deposition and propagation matrices do not vary significantly over redshift, we can approximate Eq. (31) as

[TABLE]

with repeated indices once again being summed over. $\overline{\mathsf{D}}^{\gamma}_{1/2}$ is defined in the same manner as $\overline{\mathsf{P}}^{\gamma}_{1/2}$. This equation essentially applies the deposition transfer matrix at $y-n\Delta y/2$ to all $n$ steps of the propagation of the spectrum $\mathbf{N}^{\gamma}$ from $y$ to $y-\Delta y$, which itself is approximated by $\overline{\mathsf{P}}^{\gamma}_{1/2}$. In our code, we call these approximations “coarsening”, and the number $n$ in both Eqs. (33) and (35) the “coarsening factor”.

III.5.4 Different Redshift Regimes

In DarkHistory we separate our transfer matrices into three redshift regimes: redshifts encompassing reionization ($z<50$), redshifts encompassing the times between recombination and reionization ($50\leq z\leq 1600$), and redshifts well before recombination ($z>1600$). During the redshifts encompassing reionization, we allow our transfer functions to be functions of $x_{\text{HII}}$ and $x_{\text{HeII}}$, enabling the use of reionization models that evolve hydrogen and helium ionization levels separately. We only consider singly-ionized helium in the current version of DarkHistory since we expect $x_{\text{HeIII}}$ not to play an important role until $z\sim 6$. We compute the transfer functions on a grid of $z^{k}$, $x^{m}_{\text{HII}}$, and $x^{n}_{\text{HeII}}$, and linearly interpolate over the grid of pre-computed transfer functions.

Between recombination and reionization, the helium ionization level lies at or below the hydrogen ionization level, since helium has a larger ionization potential at $24.6\text{\,}\mathrm{e\kern-0.50003ptV}$. After recombination, current experimental constraints typically forbid a large ionization fraction, i.e. we expect $x_{\text{HII}}\lesssim 0.1$ Liu et al. (2016). As such, setting $x_{\text{HeII}}=0$ is a good approximation for the photon propagation and deposition functions: since $\mathcal{F}_{\text{He}}\sim 8\%$, neglecting helium ionization only results in $\lesssim 8\%$ error to $x_{e}$, and $\lesssim 10\%$ error in the density of neutral helium. We therefore follow the same procedure as before, except we now calculate and interpolate the transfer functions over a grid of $z^{k}$ and $x_{\text{HII}}^{m}$ values while holding the helium ionization level fixed to zero.

Finally, well before recombination, we expect the universe to be close to 100% ionized and tightly coupled thermally to the CMB. Any extra source of exotic energy injection that is consistent with current experimental constraints will likely have a negligible effect on the ionization and thermal histories. We thus calculate and interpolate our transfer functions over a grid of $z^{k}$ values while holding the hydrogen and helium ionization levels to the baseline values provided by RECFAST Wong et al. (2008).

The actual grid values $z^{k}$, $x^{m}_{\text{HII}}$, and $x^{n}_{\text{HeII}}$ in each of these regimes can be found in the code, and have been chosen so that interpolation errors are at the sub-10% level when $f_{c}(z)$ is calculated using the same method detailed in Ref. Slatyer (2016). Our results for $f_{c}(z)$ without taking into account backreaction, including some improvements over Ref. Slatyer (2016), can be found in Appendix C.

III.6 Calculating $f_{c}(z)$

The low-energy photons $\mathbf{N}^{\gamma}_{\text{low}}[E_{i},z]$ and low-energy electrons $\mathbf{N}^{e}_{\text{low}}[E_{i},z]$, defined in Sec III.1, transfer their energy into ionization and excitation of atoms, heating of the IGM, and free-streaming photons to be added to the CMB continuum. In DarkHistory we keep track of how much energy low energy photons and electrons deposit into channels c $\in$ {‘H${}_{\text{ion}}$’, ‘He${}_{\text{ion}}$’, ‘exc’, ‘heat’, ‘cont’}, which represent hydrogen ionization, helium ionization, hydrogen excitation, heating of the IGM, and sub-$10.2\text{\,}\mathrm{e\kern-0.50003ptV}$ continuum photons respectively. The energy deposition fractions $f_{c}(z)$ are then found by normalizing the total energy deposited into channel c within a redshift step by the total energy injected within that step according to Eq (4). We closely follow the method for computing $f_{c}(z)$ described in Ref. Slatyer (2016).

Before calculating $f_{c}(z)$ for each channel, it is instructive to see how to calculate the total amount of energy deposited per unit time and volume, $\left(dE/dV\,dt\right)^{\text{dep}}$. The low-energy photon and electron spectra $\mathbf{N}^{\gamma}_{\text{low}}[E_{i}]$ and $\mathbf{N}^{e}_{\text{low}}[E_{i}]$ as defined above contain a number of particles per baryon deposited within each log-redshift bin (the $z$-dependence has been suppressed since all calculations in this section occur at the same redshift step). We can convert between these and spectra containing the number of particles produced per unit volume and unit time using the conversion factor $G(z)$ introduced in Eq. (10). For example, to obtain the total amount of energy deposited at a given redshift per unit time and volume, one simply sums over low-energy particle type and applies the conversion factor,

[TABLE]

To calculate the total amount of energy deposited we must also add the amount deposited by high energy electrons and photons, which we computed in Eq. (32):

[TABLE]

Then the total deposited energy summed over all channels is given by

[TABLE]

With this example in mind, we are now ready to understand how to split the energy deposition into the different channels.

III.6.1 Photons

We first compute $f_{c}(z)$ for low-energy photons, starting with energy deposition into continuum photons. These are photons with energy below $3\mathcal{R}/4=$10.2\text{,}\mathrm{e\kern-0.50003ptV}$$ that are unable to effectively transfer their energy to free electrons or atoms, so they just free stream. The energy of these photons constitutes deposition into the continuum channel, i.e.

[TABLE]

To calculate the total amount of energy deposited into hydrogen excitation, we make the approximation that all photons with energies between $3\mathcal{R}/4=$10.2\text{,}\mathrm{e\kern-0.50003ptV} and $\mathcal{R}=$13.6\text{\,}\mathrm{e\kern-0.50003ptV} deposit their energy instantaneously into hydrogen Lyman-$\alpha$ excitation, following Slatyer (2016):

[TABLE]

A more complete treatment of excitation would involve keeping track of sub-$13.6\text{\,}\mathrm{e\kern-0.50003ptV}$ energy photons as they redshift into the Lyman-$\alpha$ transition region at $10.2\text{\,}\mathrm{e\kern-0.50003ptV}$, and should also include two-photon excitation into the $2s$ state.777Two-photon $1s\to 2s$ transitions are in fact as important as Lyman-$\alpha$ transitions near recombination in determining the ionization history, due to the fact that the Lyman-$\alpha$ line is optically thick at this time. Finally, helium excitation has been neglected, since the de-excitation of helium atoms, which occurs quickly, produces photons that can eventually photoionize hydrogen. We therefore expect almost no net deposition of energy into helium excitation. Energy injection through helium excitation would mainly affect the process of helium recombination, when the probability of ionization after excitation to a higher state is significant due to the photon bath. However, we do not track this small effect, since the change to $x_{e}$ would be very small. We leave a more careful treatment of excitation that can correctly take into account all of these effects to future work.

We now move on to ionization. All photons above $\mathcal{R}=$13.6\text{,}\mathrm{e\kern-0.50003ptV} that are included in $N^{\gamma}_{\text{low}}$ have photoionized one of the atomic species (HI, HeI and HeII). However, after photoionizing a helium atom, the resulting ion may quickly recombine with an ambient free electron, producing an $\mathcal{R}_{\text{He}}=$24.6\text{\,}\mathrm{e\kern-0.50003ptV} or $4\mathcal{R}=$54.4\text{,}\mathrm{e\kern-0.50003ptV}$$ photon, which may then go on to photoionize hydrogen instead.888The photoionization rate on neutral hydrogen is much faster than the Hubble rate for $x_{\text{HII}}\lesssim 0.9999$ for $z>3$.

We can handle low-energy photons with energy $E_{\gamma}$ that photoionize neutral helium in one of the following three ways:

if helium is completely ignored, the photon is assumed to photoionize hydrogen, producing a low-energy electron with energy $E_{\gamma}-\mathcal{R}$ from photoionization and depositing $\mathcal{R}$ into hydrogen ionization. This is the approach used in previous calculations of $f_{c}(z)$ Slatyer (2016), but leaves us unable to self-consistently track $x_{\text{HeII}}$ if desired; 2. 2.

the photon produces a low-energy electron with energy $E_{\gamma}-\mathcal{R}_{\text{He}}$ from photoionization, depositing $\mathcal{R}$ into hydrogen ionization from the recombination photon (with energy $\mathcal{R}_{\text{He}}$) and producing an electron with energy $\mathcal{R}_{\text{He}}-\mathcal{R}$, which ultimately deposits energy into hydrogen excitation, heating and sub-$10.2\text{\,}\mathrm{e\kern-0.50003ptV}$ photons. This approach was previously discussed in Ref. Galli et al. (2013), and found to result in very little difference when compared to method (1); or 3. 3.

the photon produces a low-energy electron with energy $E_{\gamma}-\mathcal{R}_{\text{He}}$ from photoionization and deposits $\mathcal{R}_{\text{He}}$ into helium ionization.

The most accurate accounting of helium ionization lies somewhere between methods (2) and (3); however, either method will likely lead to very similar results in terms of $x_{e}$ and $T_{m}$, since the bulk of the energy is deposited by the electron from the initial photoionization for photon energies $E_{\gamma}\gg\mathcal{R}_{\text{He}}$, and the remaining energy always leads to one ionization event overall. DarkHistory offers the choice of these three options for implementing helium ionization.

We have checked that all three methods lead to similar ionization and temperature histories for DM models over a large range of masses decaying to both $e^{+}e^{-}$ and $\gamma\gamma$; these checks are shown in Appendix C. We recommend simply using method (1) with helium turned off if the user is interested in ionization and temperature histories well before reionization, and using both method (2) and (3) with helium turned on to bracket the uncertainties associated with energy deposition on helium if the user is interested in the epoch of reionization.

To summarize, the amount of deposited energy into hydrogen per unit time and volume is given by

[TABLE]

and into helium ionization by:

[TABLE]

producing a low-energy electron spectrum after photoionization of

[TABLE]

where $\delta[E_{i}-\mathcal{R}_{\text{He}}+\mathcal{R}]$ is one when the bin boundaries span the energy $\mathcal{R}_{\text{He}}-\mathcal{R}$ and is zero otherwise, and

[TABLE]

with the $\sigma$’s denoting the photoionization cross section of the appropriate species. $\mathbf{N}_{\text{ion}}^{e}$ is added to the low-energy electron spectrum, $\mathbf{N}^{e}_{\text{low}}$, which is then treated in the next section. The values of the $q$-coefficients depend on the method, and are shown in Table 1.

III.6.2 Electrons

To compute how low-energy electrons deposit their energy into the different channels, we use the results obtained by the MEDEA code Evoli et al. (2012); Valdes et al. (2010), following a similar treatment to Ref. Galli et al. (2013). Although DarkHistory also includes a calculation of electron energy deposition, which we discussed in Sec. III.4, the MEDEA results are more accurate in the sub-$3\text{\,}\mathrm{ke\kern-0.50003ptV}$ electron energy range, including a more detailed accounting of all possible atomic processes (such as $2s\to 1s$ deexcitations) and with more up-to-date cross sections. However, at mildly nonrelativistic to mildly relativistic regimes, our calculation of ICS is more accurate, as argued in Sec. III.4.1. Furthermore, the MEDEA results assume that hydrogen and helium are at similar ionization levels, which is not always a good assumption. In future versions of DarkHistory, an improved treatment of electrons may be a useful addition to the code.

The MEDEA code uses a Monte Carlo method to track high-energy electrons as they are injected into the IGM, and determines the fraction of the initial electron energy deposited into ionization, Lyman-$\alpha$ excitation, heating of the gas and sub-$10.2\text{\,}\mathrm{e\kern-0.50003ptV}$ photons. We use a table of these energy deposition fractions $p_{c}(E_{i},x_{e,j})$ Galli et al. (2013), where $c\in$ {‘H${}_{\text{ion}}$’, ‘He${}_{\text{ion}}$’, ‘exc’, ‘heat’, ‘cont’} as before, $x_{e,j}$ ranges between 0 and 1, and $E_{i}$ ranges between $14\text{\,}\mathrm{e\kern-0.50003ptV}$ and $3\text{\,}\mathrm{ke\kern-0.50003ptV}$, and perform an interpolation over these values. The energy deposition from electrons is then simply given by

[TABLE]

keeping in mind that $\mathbf{N}^{e}_{\text{ion}}$ has already been added to $\mathbf{N}^{e}_{\text{low}}$. Between energies of $10.2\text{\,}\mathrm{e\kern-0.50003ptV}$ and $13.6\text{\,}\mathrm{e}\mathrm{V}$, where collisional excitations of hydrogen are possible but not ionization, we use the result at $14\text{\,}\mathrm{e}\mathrm{V}$, but setting the component into hydrogen ionization to zero and normalizing to unit probability. Below $10.2\text{\,}\mathrm{e}\mathrm{V}$, electrons can only deposit energy through Coulomb heating.

III.6.3 High-Energy Deposition

Finally, the high-energy deposition component of the total energy deposited is given by:

[TABLE]

where $c\in$ { ‘ion’, ‘exc’, ‘heat’ }. Here, we add the high-energy excitation and ionization component to Lyman-$\alpha$ excitation and hydrogen ionization for simplicity, even though the high-energy deposition is computed for all atomic species. A more accurate computation of this together with a more consistent treatment of helium ionization will be a potential improvement in a future version of DarkHistory.

With the rate of energy deposition through both low-energy photons and low-energy electrons computed, the total energy deposition rate is then straightforwardly given by

[TABLE]

where $\alpha\in\{\gamma,e,\text{high}\}$.

III.7 TLA Integration and Reionization

DarkHistory offers several options for which set of assumptions should be used when integrating the ionization and thermal histories. In the simplest case, the user may integrate Eq. (6) at each redshift step based on the $f_{c}(z,\mathbf{x})$ calculated above, with the reionization terms switched off. As we have discussed, including this backreaction is already a significantly better treatment compared to calculations which assume a standard recombination history, i.e. using $f_{c}(z,\mathbf{x}_{\text{std}}(z))$ (although backreaction can also be switched off within DarkHistory).

The next significant improvement that is implemented within DarkHistory is the tracking of the neutral helium ionization fraction. Well before reionization, neglecting helium is a good approximation, since the number density of helium nuclei is only $\mathcal{F}_{\text{He}}\simeq 0.08$ of hydrogen, and we should expect only at most an 8% correction to $x_{e}$ if we include helium.

However, tracking helium allows us to accomplish a self-consistent modeling of exotic energy injection and the reionization of hydrogen and neutral helium. DarkHistory allows users to input a model of reionization, for the first time extending the validity of these energy injection calculations into a regime where hydrogen is fully ionized and helium is singly ionized.

III.7.1 Helium

The DarkHistory evolution equation governing helium without any energy injection is identical to the RECFAST model, and is given by Wong et al. (2008)

[TABLE]

The singlet and triplet ground states of helium must be treated separately, and terms relevant to the singlet or triplet state are represented with a superscript $s$ or $t$ respectively. Here, $\alpha_{\text{HeI}}$ and $\beta_{\text{HeI}}$ are the recombination and photoionization for HeI, $E_{21}^{\text{He}}$ represents the energy difference between the corresponding $n=1$ and $n=2$ states, and finally $\mathcal{C}_{\text{HeI}}$ is the analog to the Peebles-C coefficient found in Eq. (1), representing the probability of a helium atom in the $n=2$ state decaying to either the singlet or triplet ground state before photoionization can occur. The reader should refer to Refs. Wong et al. (2008); Kholupenko et al. (2008, 2007) for details on the numerical values of the coefficients, as well as how to compute $\mathcal{C}_{\text{HeI}}$.

We emphasize that although we have implemented all of the modifications to the standard TLA in Eq. (1), our code should not be used for high-precision cosmology, given that it has not been tested extensively, e.g. with different cosmological parameters from the central values used in DarkHistory. We find that our code agrees to within 3% of the RECFAST $x_{e}$ values for the cosmological parameters used here, which is sufficient for computing the effects of exotic energy injection at this stage.

In the presence of exotic sources of energy injection, low-energy photons and electrons can also change the helium ionization level. Once again, we express the energy injection source term as

[TABLE]

where $\mathcal{R}_{\text{He}}=$24.6\text{,}\mathrm{e\kern-0.50003ptV}$$ is the ionization potential of neutral helium. As we discussed in Sec. III.6.1, there are three different methods available to evaluate $f_{\text{He}_{\text{ion}}}$ which bracket the uncertainties involved in helium ionization.

To summarize, the user may opt to track the change in helium ionization levels. This means that in addition to Eq. (6), we also include

[TABLE]

where $\dot{x}_{\text{HeII}}^{\text{re}}$ is the contribution from processes that are active during reionization.

III.7.2 Reionization

The evolution equations shown in Eqs. (6) and (50) can be integrated with all reionization terms switched off if the user is primarily interested in temperatures or ionization levels well before reionization starts at $z\sim 20$. In this regime, turning off helium is also a reasonable approximation.

With reionization however, the helium ionization level should be solved as well for complete consistency. We solve the TLA differential equations shown in Eqs. (6) and (50) in two separate redshift regimes. Prior to some user-defined reionization redshift $1+z_{\text{re}}$ ($z_{\text{re}}\leq 50$), we set $\dot{T}_{m}^{\text{re}}$, $\dot{x}_{\text{HII}}^{\text{re}}$ and $\dot{x}_{\text{HeII}}^{\text{re}}$ to zero. Once reionization begins, we set $\dot{x}_{\text{HII}}^{(0)}$ and $\dot{x}_{\text{HeII}}^{(0)}$ to zero for $z<z_{\text{re}}$ instead, switching over to the specified reionization model with its own photoionization and recombination rates.999We do not set $\dot{T}_{m}^{(0)}=0$, since both adiabatic cooling and Compton scattering off the CMB remain active during reionization. We also begin tracking doubly-ionized helium $x_{\text{HeIII}}$, which is always assumed to be zero before reionization.

The $\dot{T}_{m}^{\text{re}}$, $\dot{x}_{\text{HII}}^{\text{re}}$ and $\dot{x}_{\text{HeII}}^{\text{re}}$ terms depend on the details of how reionization proceeds, which is still relatively uncertain. However, choosing a model for the formation of stars and active galactic nuclei (AGNs) and the associated photoionization and photoheating rates, these terms can be evaluated. DarkHistory by default includes the Puchwein+ model of Ref. Puchwein et al. (2018). We also demonstrate how to implement the older Madau and Haardt model Haardt and Madau (2012) in Example 8. Both models provide a photoionization rate $\Gamma_{\gamma X}^{\text{ion}}(z)$ and a photoheating rate $\mathcal{H}_{\gamma X}^{\text{ion}}(z)$ as a function of redshift and species $X$.

Along with these energy injection rates, we must also include other relevant processes that alter the ionization fraction of each species. Since these processes generally convert kinetic energy to atomic binding energy, cooling or heating of the gas due to these processes must also be included in $\dot{T}_{m}^{\text{re}}$. The processes we include are:

collisional ionization, occuring at a rate $\Gamma_{eX}^{\text{ion}}$ for each species $X$, and an associated cooling rate $-\mathcal{H}_{eX}^{\text{ion}}$; 2. 2.

case-A recombination, described by a rate coefficient $\alpha_{A,X}$ for each species $X$, and an associated cooling rate $-\mathcal{H}_{X}^{\text{rec}}$; 3. 3.

collisional excitation cooling, with a rate $-\mathcal{H}_{eX}^{\text{exc}}$; and 4. 4.

bremsstrahlung cooling, with a rate $-\mathcal{H}^{\text{br}}$.

The cooling rates here have been defined with a negative sign so that all quantities denoted by $\mathcal{H}$ contribute positively to any temperature change. Expressions for all of these rates can be found in Ref. Bolton and Haehnelt (2007). They are explicitly dependent on the ionization fraction of all three of the relevant species, namely $x_{\text{HI}}$, $x_{\text{HeI}}$ and $x_{\text{HeII}}$. The full expressions for the evolution of each of these fractions are as follows:

[TABLE]

with the temperature evolution given by

[TABLE]

Instead of specifying a full reionization model, the user may also choose the simpler alternative of fixing the value of $x_{\text{HII}}$ and $x_{\text{HeII}}$ as a function of redshift once reionization begins, and integrate the temperature evolution alone instead. We note that this approach is not self-consistent, since fixing the ionization levels forces us to neglect any additional contribution to ionization from exotic energy injection sources. However, if the contribution to ionization is known to be small, this can serve as a useful approximation.

III.7.3 Numerical Integration

To ensure that ionization fractions always remain appropriately bounded during integration, we introduce the variable

[TABLE]

where $\chi_{i}=1$ for HI and $\chi_{i}=\mathcal{F}_{\text{He}}$ for HeI and HeII. This transformed equation is then integrated using the standard odeint integrator provided by SciPy.

At early times, the equations we are integrating are very stiff, and solving them directly with numerical integration can often run into difficulties. We therefore assume that when $x_{\text{HII}}>0.99$ or $x_{\text{HeII}}>0.99\mathcal{F}_{\text{He}}$, either variable follows their Saha equilibrium values.

In Sec. V.4, we will show several thermal and ionization histories that showcase DarkHistory’s capabilities in tracking the helium ionization level, exotic energy injection and reionization all at the same time.

IV Modules

In this section we summarize the main modules in DarkHistory. We will pay particular attention to the modules shown in the flow chart in Fig. 1, and as far as possible provide links between the code and the text. Keep in mind that this is not a complete list and that it is subject to change in future versions of the code. There is more thorough documentation in DarkHistory itself that will be periodically updated at https://darkhistory.readthedocs.io, and will contain a more complete explanation of the code. In the interest of space, we only provide the full path of each module in the code when it is mentioned for the first time.

IV.1 Data

First, the user must download the data files found at https://doi.org/10.7910/DVN/DUOUWA. These files contain the photon propagation transfer function $\overline{\mathsf{P}}_{\gamma}$ and deposition transfer functions $\overline{\mathsf{D}}_{\gamma}$, $\overline{\mathsf{D}}_{e}$ and $\overline{\mathbf{D}}_{c}^{\text{high}}$, which have all been precomputed as discussed above. They also contain transfer functions for ICS calculations discussed in Appendix A, structure formation annihilation boost factors computed in Ref. Liu et al. (2016), the baseline thermal and ionization histories, data from pppc4dmid Cirelli et al. (2011) and $f_{c}(z)$ computed without backreaction for DM annihilation and decay, where photons and $e^{+}e^{-}$ are injected at a fixed set of energies.

IV.2 config

The config module contains the code required to access the downloaded data, and to store them in memory for use. Users should ensure that the variable data_path points to the directory containing the data files.

IV.3 main

The main module contains the function that implements the loop shown in Fig. 1, evolve(). The usage of this function will be discussed in great detail in Sec. V.

IV.4 darkhistory.physics

This module contains physical constants and useful functions found in cosmology, particle physics and atomic physics. We use units of $\text{\,}\mathrm{c}\mathrm{m}$ for length, $\text{\,}\mathrm{s}$ for time and $\text{\,}\mathrm{e\kern-0.50003ptV}$ for energy, mass and temperature. Some examples of functions that are included in this module include the Hubble parameter as a function of redshift, physics.hubble(), and the Peebles-C factor $\mathcal{C}$ found in Eq. (1), physics.peebles_C(). Constants provided in this module are taken from central values of the Planck 2018 TT,TE,EE+lowE results Aghanim et al. (2018) and the Particle Data Group review of particle physics Tanabashi et al. (2018).

IV.5 darkhistory.electrons

The electrons module contains all of the functions necessary to perform the electron cooling calculation. The positronium submodule contains functions that return the spectrum of photons obtained during positronium annihilation, which we denoted as $\overline{\mathbf{N}}_{\gamma}^{\text{pos}}$ in Eq. (21); Example 7 demonstrates how to use this module. The ics submodule contains all of the machinery necessary to compute the ICS scattered photon and electron spectra; for more details on how to use this submodule, refer to Example 4 in the code.

elec_cooling contains the code necessary to compute the transfer functions $\overline{\mathbf{R}}_{c}$, $\overline{\mathsf{T}}_{\text{ICS}}$ and $\overline{\mathsf{T}}_{e}$, as defined in Eqs. (13), (16) and (17) respectively; Example 6 shows how this module is used.

IV.6 darkhistory.history

This module contains our implementation of the TLA and reionization. The submodule tla corresponds to Sec. II where the function get_history implements the TLA, including all of the terms discussed in Eqs. (6) and Eqs. (50)–(III.7.2). The submodule reionization contains the Puchwein+ reionization model, and contains all of the coefficients found in Eqs. (III.7.2) and (III.7.2).

IV.7 darkhistory.low_energy

This module calculates $f_{c}(z)$. The lowE_photons and lowE_electrons submodules correspond to Sec. III.6.1 and Sec. III.6.2, respectively, implementing Eqs. (39)–(43) and Eq. (45) respectively. The lowE_deposition submodule then combines the energy deposited by photons, electrons (including photoionized electrons) and high-energy deposition to make $f_{c}(z,\mathbf{x})$.

IV.8 darkhistory.spec

This module contains functions for handling and generating spectra and transfer functions. All one dimensional spectra in the code can be handled using the class Spectrum, which stores not just the data of the spectrum, but also the abscissa, and other relevant information like redshift or the injection energy of the particle that produced the spectrum. This class includes many convenience functions, such as the ability to rebin the spectrum into a new binning while conserving total number and energy, or the ability to quickly obtain the total number of particles within some energy range. Example 1 in our code gives a quick introduction to this class.

The user may also want to store closely related spectra in one object. This may be desirable for spectra of the same particle type over different redshifts, or if they correspond to spectra from the same injected particle but at different injection energies. The class Spectra has been written to do exactly this. Example 2 provides a good overview of what this class can do.

IV.9 darkhistory.spec.pppc

Within the spec module, a dedicated submodule pppc has been written to calculate the electron and photon spectra from the injection of any arbitrary Standard Model particle, based on the pppc4dmid results. The function pppc.get_pppc_spec() is the main function to use for this end. See Example 4 for more information on how to use this function.

V Using the Code

We will now apply DarkHistory to perform a variety of calculations in order to highlight the key functionalities of the code. Each of the subsections corresponds to an example Jupyter notebook that has been provided as part of the code; the user should refer to these examples for a deeper look at the full capability of the code, as well as to the online documentation. In this paper, we will simply highlight capabilities and interesting physics results.

Within the code and in this section, the word “redshift” and variables that represent redshift (usually called rs in the code) refer to the quantity $1+z$, since this is the physically relevant quantity in many cosmological calculations.

V.1 A Simple Model: $\chi\chi\to b\bar{b}$

As a first example, we will demonstrate how to compute the ionization and thermal history of a simple annihilation model. Consider a $50\text{\,}\mathrm{Ge\kern-0.50003ptV}$ Majorana fermion DM particle that undergoes $s$-wave annihilation to a pair of $b\overline{b}$ quarks, with an annihilation cross section $\langle\sigma v\rangle=$2\text{\times}{10}^{-26}\text{,}{\mathrm{cm}}^{3}\text{,}{\mathrm{s}}^{-1}$$, close to the required thermal freezeout cross section for the correct relic abundance. Similar models have been considered as a possible dark matter explanation for the galactic center excess Calore et al. (2015) and the AMS-02 antiproton excess Cui et al. (2017); Cuoco et al. (2017). We perform the calculation in a relatively simplified setting: with no reionization, no backreaction included, but with a boost to the annihilation rate from structure formation. For more details, see Example 9 in the code.

The function that we use to compute histories is main.evolve(). There are many keyword parameters that can be used with this function, and the user should refer to the example notebooks and the online documentation for more information. To find the thermal history for this model, evolve() can be called in the following fashion:

⬇

import main

import darkhistory.physics as phys

bbbar_noBR = main.evolve(

    DM_process=’swave’, mDM=50e9,

    sigmav=2e-26,

    primary=’b’, start_rs=3000.,

    coarsen_factor=32, backreaction=False,

    struct_boost=phys.struct_boost_func()

)

The keyword parameters are as follows:

DM_process=’swave’ – specifies the DM process of interest. Currently, DarkHistory can handle $s$-wave annihilating and decaying DM models (DM_process=’decay’) with this keyword; 2. 2.

mDM=50e9 – specifies the DM mass, in $\text{\,}\mathrm{e\kern-0.50003ptV}$; 3. 3.

sigmav=2e-26 – specifies the velocity averaged annihilation cross section, in $\text{\,}{\mathrm{cm}}^{3}\text{\,}{\mathrm{s}}^{-1}$; 4. 4.

primary=’b’ – specifies the annihilation channel. The options include all of those offered by pppc4dmid Cirelli et al. (2011), and the spectra are extracted from the raw data provided by the cookbook. The $e^{+}e^{-}$ and photon spectra from the showering of a single $b\overline{b}$ pair are shown in Fig. 2. These are proportional to the injection spectra $\mathbf{N}_{\text{inj}}^{\alpha}$ defined in Sec. III.3, and can be generated using the function pppc.get_pppc_spec(); 5. 5.

start_rs=3000 – the redshift at which to start the evaluation. $1+z=3000$ is the highest redshift at which we have produced the photon cooling transfer functions, and represents the highest redshift that should be specified here. In this example, start_rs fixes the initial conditions of the TLA in Eq. (1) at the baseline ionization and temperature values at this redshift; 6. 6.

coarsen_factor=32 – the coarsening factor, defined in Sec. III.5.3. For a comparison between solutions with different coarsening factors, see Appendix C; 7. 7.

backreaction=False – this turns backreaction on and off; and 8. 8.

struct_boost=phys.struct_boost_func() – the structure formation prescription to use. Once dark matter halos start to collapse, the annihilation rate gets enhanced by the factor

[TABLE]

compared to the smooth annihilation rate shown in Eq. (3). Here, struct_boost is a function that takes redshift as the argument, and returns $1+\mathcal{B}(z)$. The user can make use of the structure formation boosts that are saved by default in DarkHistory in the physics module, which include the boost factors computed in Ref. Liu et al. (2016), and is used as the default boost factor by struct_boost_func().

By default, the solver integrates the equations down to $1+z=4$, and will not evolve the helium ionization levels. These choices can of course be changed with other keyword parameters. Note that the function is not limited to DM processes or pppc4dmid spectra; other keyword parameters allow the user to specify their own injection rates as a function of redshift (see the documentation for the keyword parameters rate_func_N and rate_func_eng), along with the spectra of photons and $e^{+}e^{-}$ injected (see the documentation for the keyword parameters in_spec_elec and in_spec_phot).

The output of evolve(), stored in bbbar_noBR, is a dictionary containing the redshift abscissa of the solutions, the ionization and temperature solutions, the propagating photon, low-energy photon and low-energy electron spectra, and the computed value of $f_{c}(z)$. To access the redshift, ionization and temperature, we can simply do:

⬇

Redshift abscissa.

rs_vec = bbbar_noBR[’rs’]

Matter temperature in eV.

Tm_vec = bbbar_noBR[’Tm’]

Ionization fraction.

Stored as 1+z by {xHII, xHeII, xHeIII}.

xHII_vec = bbbar_noBR[’x’][:,0]

In Fig. 3 we plot $T_{m}$ and $x_{\text{HII}}$ as a function of redshift for the $\chi\chi\to b\bar{b}$ model. For DM masses above $\gtrsim$10\text{,}\mathrm{Ge\kern-0.50003ptV}$$, values of $\langle\sigma v\rangle$ required for thermal freezeout are unconstrained by the CMB anisotropy power spectrum energy injection constraints: the ionization fraction, which changes by approximately 25% only at high redshifts, does not change enough to affect the power spectrum significantly. The sudden increase in ionization and temperature at $z\sim 30$ corresponds to an increase in the boost factor used (halos with an Einasto profile with halo substructure boost included Liu et al. (2016), found in physics.struct_boost_func()).

We also show in Fig. 3 for completeness the effect of turning on backreaction, i.e. including the effect of the increased ionization level on the evolution of the ionization and thermal histories. This is conveniently done by setting backreaction=True. In this particular example, the effect of backreaction is small, but we will show more scenarios where backreaction has large effect on $T_{m}$, and explain why this can be significant in the next example.

V.2 Backreaction

Let us explore the effects of backreaction a bit more using some of the code found in Example 10 of DarkHistory. As was described in Sec II, one of DarkHistory’s main improvements to ionization and temperature history calculations is its ability to include the effects of back-reaction. To see its importance, consider the example of $100\text{\,}\mathrm{Me\kern-0.50003ptV}$ dark matter decaying to a pair of $e^{+}e^{-}$, with a lifetime of $\tau=$3\text{\times}{10}^{25}\text{,}\mathrm{s}$$, a value that is close to the minimum lifetime allowed by constraints from the CMB power spectrum Slatyer and Wu (2017). The ionization and thermal histories can be evaluated in this way:

⬇

decay_BR = main.evolve(

DM_process=’decay’, mDM=1e8, lifetime=3e25,

primary=’elec_delta’, start_rs=3000.,

coarsen_factor=16, backreaction=True

)

The new keywords here are:

DM_process=’decay’ – specifies the DM process of interest to be decays; 2. 2.

lifetime=3e25 – specifies the decay lifetime in seconds; and 3. 3.

primary=’elec_delta’ – the primary channel options ’elec_delta’ and ’phot_delta’ can be used to inject an $e^{+}e^{-}$ and $\gamma\gamma$ pair respectively, with no electroweak corrections applied.

To do the calculation without backreaction, we can simply set backreaction=False. However, with primary=’elec_delta’ or ’phot_delta’, DarkHistory can instead rely on tabulated results of $f_{c}(z)$ for these two channels, using the same method based on results from Ref. Slatyer (2016), to calculate the ionization and thermal histories without evolving the input spectrum, leading to a significant speed-up. This can be done using the function tla.get_history():

⬇

import numpy as np

from darkhistory.tla import get_history

get_history takes a redshift vector:

rs_vec = np.flipud(np.arange(5, 3000, 0.1))

result = get_history(

 rs_vec, baseline_f=True, mDM=1e8,

 lifetime=3e25, DM_process=’decay’,

 inj_particle=’elec_delta’

)

with the following parameters:

rs_vec – the redshift vector, ordered from high to low, over which the temperature and ionization histories are to be evaluated; 2. 2.

baseline_f=True – this tells the code to use the baseline $f_{c}(z)$ computed by DarkHistory without backreaction. As we discussed in Sec. III.6, these $f_{c}(z)$ agree with those computed in Ref. Slatyer (2016) to within 10%, and 3. 3.

inj_particle=’elec_delta’ – used to specify one of two options ’elec_delta’ or ’phot_delta’.

The output result is an array of shape (len(rs_vec), 4), with the second dimension indexing $\{T_{m},x_{\text{HII}},x_{\text{HeII}},x_{\text{HeIII}}\}$. The temperature (in $\text{\,}\mathrm{e\kern-0.50003ptV}$) can be accessed through T_m = results[-1,0].

Although only the $f_{c}(z)$ values for the injection for an $e^{+}e^{-}$ and $\gamma\gamma$ pair have been saved for use with DarkHistory, the $f_{c}(z)$ for any arbitrary channel can be computed from a weighted average of the electron and photon results Slatyer (2016). We stress once again, however, that this can only be done assuming no backreaction.

The histories are shown in Fig 4, with and without backreaction turned on. First, even though the ionization level at $z\sim 10$ is three orders of magnitude larger than the baseline, such a scenario is actually still consistent with the CMB power spectrum constraints, owing to the fact that the ionization build-up occurs relatively late: the CMB constraints are sensitive to changes in $x_{e}$ near recombination, and become less sensitive at later times.

Comparing the temperature histories with and without backreaction, we see that the main effect of this increase in $x_{e}$ on the energy deposition processes is to increase energy deposition into heating. Ionization and excitation rates depend on the neutral fraction, which is still close to 100% even with energy deposition from DM. However, the energy rate into Coulomb heating is proportional to $x_{e}$, so taking into account the significantly elevated $x_{e}$ values leads to higher temperature levels. By about $z\sim 10$, $T_{m}$ with backreaction is larger than without backreaction by a factor of $\sim 4$, with the difference continuing to grow. Neglecting backreaction therefore leads to a severe underestimate of $T_{m}$, and including this effect consistently will certainly be important in understanding what measurements of $T_{m}$ at $z\simeq 20$ through the 21-cm signal or the Lyman-$\alpha$ power spectrum can tell us about exotic sources of energy injection.

We can perform the calculation over a range of DM masses by looping over values of mDM. For each value of $m_{\chi}$, we select the minimum lifetime $\tau$ which is consistent with the CMB power spectrum constraints, and compare the difference between the temperature history with backreaction ($T_{m,\text{BR}}$) and without ($T_{m,0}$) by computing the fractional change in temperature,

[TABLE]

In Fig 5 we plot this variable over a range of redshifts and dark matter masses for this particular channel ($\chi\to e^{+}e^{-}$), but also for decay and annihilation into $e^{+}e^{-}$ and $\gamma\gamma$, taking the maximum $\langle\sigma v\rangle$ again allowed by the CMB power spectrum constraints. At a redshift of $z\sim 17$ near the end of the cosmic dark ages, $\delta T_{m}/T_{m,0}\sim 100\%$ (i.e. $T_{m}$ with backreaction is a factor of 2 larger than without) or more can easily be obtained. Even larger deviations are possible at lower redshifts, depending on the channel under consideration.

V.3 21-cm Sensitivity

The global 21-cm signal is a measurement of the sky-averaged differential brightness temperature $T_{21}$ with respect to the background radiation. Measurements of this signal would open a window into the ionization and temperature histories of the universe at the cosmic dawn (see e.g. Ref. Pritchard and Loeb (2012) for a review of 21-cm cosmology). A first claim of such a measurement has already been made by the EDGES collaboration Bowman et al. (2018). The brightness temperature of the 21-cm hydrogen absorption line relative to the background radiation temperature is given by Pritchard and Loeb (2012):

[TABLE]

where $\Omega_{b}$ is the baryon energy density today as a fraction of the critical density, $h$ is the Hubble parameter today in $\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$, $T_{R}$ is the background radiation temperature (typically assumed to be the CMB temperature) and $T_{S}$ is the spin temperature of neutral hydrogen as a function of redshift, which determines the relative population of neutral hydrogen in the two hyperfine states. Due to the presence of an intense Lyman-$\alpha$ radiation field once stars begin to form, it is expected that $T_{S}\approx T_{m}$ at the cosmic dawn. This fact allows us to turn the 21-cm global signal into a limit on $T_{m}$ itself, assuming that $T_{R}=T_{\text{CMB}}$.

We will focus on $1+z\approx 18$, roughly the central value of the absorption trough measured by EDGES Bowman et al. (2018). At this redshift, almost all hydrogen is neutral, i.e. $x_{\text{HI}}\approx 1$, and we can invert Eq. (56) to find $T_{S}$ as a function of $T_{\text{21}}$. Since $T_{m}<T_{S}$, this yields the bound

[TABLE]

This temperature bound in turn puts a limit on the DM decay lifetime or cross-section because too much dark matter decay/annihilation would heat up $T_{m}$ past this point.

In contrast to the CMB power spectrum energy injection bounds, which is most sensitive to changes in $x_{e}$ around the time of recombination, the 21-cm global signal constraints are more sensitive to energy injection processes that are more active at late times, and are dependent primarily on $T_{m}$ instead. Since $T_{m}$ is significantly impacted by including the effects of backreaction, the calculation performed by DarkHistory becomes important for setting accurate constraints using the 21-cm global signal.

To illustrate this, we perform a simple sensitivity study by obtaining the constraints for a measured $T_{\text{21}}$ of $-50\text{\,}\mathrm{mK}$, and compare the constraints with and without backreaction taken into account. Although this value of $T_{21}$ is inconsistent with the EDGES experiment, it is impossible to interpret the EDGES result without proposing new physics that may be at play during the cosmic dark ages Liu and Slatyer (2018), which is a more complicated task and less relevant to helping users understand the code. The following analysis is worked out in more detail within the code in Example 11.

$T_{21}=$-50\text{,}\mathrm{mK} means that we require $T_{m}<$20.3\text{\,}\mathrm{K} according to Eq. (57). We once again scan over a grid of dark matter masses and lifetimes/cross-sections decaying/annihilating into $e^{+}e^{-}$ and $\gamma\gamma$, using get_history() for the case with no backreaction and evolve() for the case with backreaction, as explained in the previous section, to find where in parameter space dark matter energy injection leads to a violation of Eq. (57).

The resulting exclusion plots are shown in Fig. 6. We see that in each case the calculation with backreaction can be between 10%-50% stronger than without backreaction, which we would expect because backreaction leads to larger temperatures. We emphasize that this is the result for just one chosen value of $T_{21}$; for larger (less negative) $T_{21}$, we expect that the importance of backreaction will increase, since the energy injection is less constrained, allowing for larger values of $x_{e}$.

V.4 Helium, Dark Matter and Reionization

Finally, we will take a closer look at the different options one can use within the code to evaluate temperature and ionization histories. Throughout this section, we will demonstrate these different options mostly using get_history(), but similar options are also available in evolve(), which calls get_history() with all of the relevant options provided. We refer the reader to the online documentation and to Example 8 in the code for more details.

Without any exotic energy injection or any reionization, the function get_history() accepts a redshift vector, and simply returns the baseline ionization and temperature histories, obtained by solving Eq. (6):

⬇

import numpy as np

from darkhistory.tla import get_history

Redshift vector in decreasing order.

rs_vec = np.flipud(np.arange(1., 3000., 0.1))

soln_baseline = get_history(rs_vec)

Turning on helium evolution within get_history() is controlled by the flag helium_TLA, i.e.

⬇

soln_He = get_history(rs_vec,helium_TLA=True)

Fig. 7 shows the solution to Eq. (1) with just the “$(0)$” terms, i.e. without any energy injection or reionization, and compares that solution to one with Eq. (50) added as well. This is simply the standard ionization history with helium recombination ($z\sim 1800$) and hydrogen recombination ($z\sim 1100$), eventually leading to the residual ionization fraction at redshifts well below hydrogen recombination of about $x_{e}\sim 2\times 10^{-4}$. The inset of Fig. 7 shows that DarkHistory is able to correctly reproduce helium recombination; the entire ionization history agrees with RECFAST results at the central cosmological parameters used by DarkHistory to within $\sim 3\%$. We recommend that helium ionization levels are tracked when used in combination with reionization.

The next important option is whether to include the effects of reionization. This option is controlled by the flag reion_switch:

⬇

soln_default_reion = get_history(

rs_vec, helium_TLA=True, reion_switch=True

)

With no other options set, setting reion_switch to True causes DarkHistory to use the standard reionization model, which is based on the photoionization and photoheating rates provided in Puchwein et al. (2018). Fig. 8 shows the IGM temperature as well as the ionization levels of the different atomic species as a function of redshift. Both of these results agree well with the same result shown in Ref. Puchwein et al. (2018). Reionization of hydrogen and neutral helium is complete by about $z\sim 6$; soon after, HeII starts to become doubly ionized, leading to a decrease in $n_{\text{HeII}}$ and a corresponding increase in $n_{\text{HeIII}}$. Dips in $T_{m}$ correspond to a decrease in photoheating rates once a species becomes completely ionized and the production of high-energy electrons from photoionization off these species ceases.

Aside from the default reionization model, the user may also supply their own reionization models in two different ways: by either providing their own photoionization and photoheating rates on each atomic species (e.g. based on a model that is different from the default, e.g. Haardt and Madau (2012)), or by fixing the ionization history below a certain redshift, e.g. with a $\tanh$ model Lewis (2008); Adam et al. (2016). We leave a discussion of how to use these options to Example 8 in the code.

With the ability to include both helium and reionization, we can now add a new source of energy injection and compute the effects on ionization and temperature levels. We remind the reader that this means we are solving Eq. (6) together with Eq. (50). This is accomplished in the code with both reion_switch and helium_TLA set to True, and supplying the same keyword parameters used to inject energy from DM shown in Sec. V.2. We can add decaying DM with mass $100\text{\,}\mathrm{Me\kern-0.50003ptV}$ into an $e^{+}e^{-}$ pair with a lifetime of $3\text{\times}{10}^{25}\text{\,}\mathrm{s}$, like so (using evolve() in this example):

⬇

main.evolve(

DM_process=’decay’, mDM=1e8,

lifetime=3e25, primary=’elec_delta’,

start_rs=3000., coarsen_factor=1,

backreaction=True, helium_TLA=True,

reion_switch=True

)

By turning on and off the various flags backreaction, helium_TLA and reion_switch, we can produce histories including or excluding these various effects.

The results from different combinations of these switches are summarized in Fig. 9. The dashed lines shows the histories with DM decay only, and illustrates the significant difference that can arise after taking into account backreaction, which we have already seen in Fig. 4. Combining the DM energy injection with the reionization model gives the solid lines in Fig. 9. These curves should be compared to the default reionization model temperature and ionization histories, shown in the black, dashed lines. When computing the DM energy deposition without taking into account backreaction, we find that the amount of energy deposited into heating from DM is much smaller than heating from reionization processes once they begin in earnest, and so adding the DM decays on top of reionization produces only a small perturbation in the temperature history relative to $T_{m}$ for just the reionization model alone. In some cases, the addition of DM actually decreases $T_{m}$: this can happen due to reionization proceeding at a faster rate, leaving fewer atoms to photoionize and thus suppressing photoheating.

It is clear, however, that neglecting backreaction leads to a severe underestimation of the energy deposition into heating. Performing the full calculation with DM, reionization and backreaction correctly accounted for produces the line in red, which shows that the addition of DM significantly increases $T_{m}$ compared to both the reionization model and the case where DM energy deposition is added without backreaction. Reionization greatly enhances the energy deposition rate into heating of the IGM by increasing the number of free charged particles available for Coulomb heating, and properly accounting for backreaction using DarkHistory is critical to predicting the IGM temperature growth due to energy injection once reionization begins.

VI Conclusion

We have developed and made public a new code package for mapping out the effects of arbitrary exotic energy injections — including dark matter annihilation and decay to arbitrary Standard Model final states — on the temperature and ionization history of the early universe. DarkHistory is capable of self-consistently including the effects of conventional astrophysical sources of ionization and heating, and of including feedback effects that can significantly enhance the degree of heating. Additionally, the ICS module can be employed independently of the rest of the code, as an accurate and efficient numerical calculator of ICS across a very wide range of electron and photon energies. We have outlined here a number of worked examples, and provide more examples with the online code at https://github.com/hongwanliu/DarkHistory.

DarkHistory has a modular framework and can in the future be improved in several different directions, while keeping the same underlying structure. In this first version we have focused on the homogeneous signal, and neglected the possible effect of new radiation backgrounds and/or gas inhomogeneities on the cascade of secondaries produced by injected high-energy particles. Such effects may become important in the late cosmic dark ages and the epoch of reionization. The spectrum of low-energy photons produced by energy injection, and the resulting distortion to the spectrum of the CMB, is a possible observable in its own right; the current version of DarkHistory provides only a partial calculation of this spectral distortion, due to our approximate treatment of low-energy electrons, but we intend to improve this aspect in future work. The effects of other new physics on the temperature/ionization evolution – in particular, scattering between baryons and DM – can be incorporated within the same framework. We also intend to explore the possibility of interfacing DarkHistory with existing public codes for computing the recombination history, perturbations to the CMB, and 21cm signals.

The tools we have developed in this work can be used to understand the visible imprints of exotic energy injections that could appear in the CMB and the 21cm line of neutral hydrogen, and hence to place precise constraints on dark matter annihilation and decay. We hope they will help pave the way for a comprehensive description of the ways in which dark matter interactions, and other physics beyond the Standard Model, could reshape the early history of our cosmos.

VII Acknowledgement

The authors would like to thank Shi-Fan Stephen Chen for a fruitful early collaboration. We would also like to thank Yacine Ali-Haïmoud, Anastasia Fialkov, Katherine Mack, Vivian Poulin, Nicholas L. Rodd, Sarah Schön and Cannon Vogel for helpful discussions. We acknowledge the use of the following Python packages: NumPy Van Der Walt et al. (2011), SciPy jon (2001), Matplotlib Hunter (2007), as well as numpy_groupies for fast rebinning of spectra.

This work was supported by the Office of High Energy Physics of the U.S. Department of Energy under Grant No. DE-SC00012567 and DE-SC0013999. This work was also partly supported by the MIT Research Support Committee (NEC Corporation Fund for Research in Computers and Communications). GR is supported by a National Science Foundation Graduate Research Fellowship. TRS is partially supported by a John N. Bahcall Fellowship.

Appendix A Inverse Compton Scattering

In this appendix, we discuss in detail the methods used to compute the spectra of photons that are produced by the cooling of electrons through ICS. We restore $\hbar$, $c$ and $k_{B}$ in this appendix, since the exact numerical value of these spectra is important.

A.1 Scattered Spectra

We begin with some preliminaries that will be important in understanding our subsequent discussion of ICS. The goal is to determine the secondary photon spectrum produced on average by multiple scatterings of a single electron.

Consider an electron with energy $E_{e}$ and corresponding Lorentz factor $\gamma$ incident on some distribution of photons $n(\epsilon)$ with initial energy $\epsilon$ in the comoving frame. Since we are only interested in ICS off the CMB, we will only consider an isotropic photon bath in the co-moving frame, distributed as a blackbody with some temperature $T$. The electron has some probability per unit time of scattering the photons into some outgoing energy $\epsilon_{1}$, with some probability distribution $dN_{\gamma}/(d\epsilon\,d\epsilon_{1}\,dt)$, which we call the “differential scattered photon spectrum”. This quantity is proportional to the number density per unit energy of the photon bath $n(\epsilon)$, so that integrating over $\epsilon$ also integrates over the distribution of these photons. This can be interpreted as a normalized scattered photon spectrum for ICS by many electrons with the same energy. Integrating the differential scattered photon spectrum with respect to $\epsilon$ gives us the “scattered photon spectrum”,

[TABLE]

with $\epsilon_{\text{min}}$ and $\epsilon_{\text{max}}$ determined by the kinematics of ICS.

We further define the “scattered photon energy loss spectrum”,

[TABLE]

where $\Delta$ is the change in energy of a photon scattering by a single electron. This is simply the distribution of scattered photons as a function of the energy gained or lost by the photon during the scattering.

Now, consider some arbitrary injection spectrum of electrons $d\tilde{N}_{e}/dE_{1}$. The tilde serves to remind the reader that this is a distribution of electrons, and not a normalized quantity. From the definition of Eq. (59), we define the “scattered electron spectrum” as

[TABLE]

where $E_{1}$ is the energy of the scattered electron. However, this result allows some electrons to gain energy after scattering, significantly complicating our calculations. Intuitively, we expect electrons that upscatter from $E\to E_{1}$ to partially cancel with downscatters from $E_{1}-E$, justifying an approximate treatment where we simply cancel out photons that downscatter (and upscatters an electron) with photons that upscatter (and downscatters an electron). We leave a full justification of this to the end of this section, but for now, we will accordingly define the “scattered electron net energy loss spectrum”,

[TABLE]

with $\Delta\geq 0$ in the expression above. For relativistic electrons, the average energy lost due to an upscattering a photon is much larger than the average energy gained due to downscattering a photon, and it is therefore a good approximation to consider only scattering events where electrons lose their energy Blumenthal and Gould (1970). The upscattered photons also have outgoing energy $\epsilon_{1}\gg\epsilon$, and so a reasonable approximation to make in the relativistic limit is

[TABLE]

Now we can turn our attention to justifying the approximation laid out in Eq. (61). First, we split the exact integral in Eq. (60) into an integral from 0 to $E_{1}$, and from $E_{1}$ to $\infty$. The first integral can be rewritten as (dropping the $T$ dependence for clarity)

[TABLE]

where we have simply made the substitution $x=2E_{1}-E$. In this part of the integral, we are dealing with upscattered electrons and downscattered photons, and so we know that $dN_{\gamma}/(d\Delta\,dt)$ only has support when $E-E_{1}\sim T_{\text{CMB}}\ll E,E_{1}$, since ICS is only included for electrons with $E>$3\text{,}\mathrm{ke\kern-0.50003ptV}$$ Slatyer (2016). This implies that the integral only has support near $x=E_{1}$, and we can therefore make the following approximation:

[TABLE]

where in the last step we have trivially relabeled $x\to E$. We have therefore shown that

[TABLE]

and that is a good approximation due to the relatively low temperature of the CMB.

With these definitions in mind, we are now ready to understand how to compute these scattered spectra when the electron is in two limits. For $\gamma>20$, the spectra are computed in the relativistic limit, while below that, scattering with the CMB at all relevant redshifts lie well within the Thomson regime. Together, they cover all relevant kinematic regimes that we consider in our code.

A.1.1 Relativistic Electrons

The differential upscattered photon spectrum produced by ICS between an electron and the CMB blackbody spectrum in the relativistic regime ($\gamma\gg 1$) is given by Blumenthal and Gould (1970)

[TABLE]

where $r_{0}$ is the classical electron radius, $m_{e}$ is the electron mass, $\epsilon$ is the incident photon energy in the comoving frame, and $\epsilon_{1}$ is the scattered photon energy in the same frame, and we have defined

[TABLE]

We stress that Eq. (66) is strictly only correct when photons are upscattered by the incoming electron, which corresponds to the kinematic regime $\epsilon\leq\epsilon_{1}\leq 4\epsilon\gamma^{2}/(1+4\epsilon\gamma/m)$. In the opposite regime where $\epsilon/(4\gamma^{2})\leq\epsilon_{1}<\epsilon$ and photons get downscattered, we have Jones (1968)

[TABLE]

For ICS off CMB photons, the $n(\epsilon)$ is the number density of photons per unit energy; for a blackbody, this is

[TABLE]

where $T$ is the temperature of the CMB.

The complete upscattered photon spectrum for ICS off the CMB is therefore obtained by performing the integral in Eq. (66) over $\epsilon$, with the kinematic limits given by $1/4\gamma^{2}\leq q\leq 1$ Blumenthal and Gould (1970). Since the CMB photons at $z\lesssim 3000$ have energies less than $1\text{\,}\mathrm{e\kern-0.50003ptV}$, the amount of energy transferred by an electron is always completely dominated by Eq. (66). Furthermore, one can check that at $q=1/4\gamma^{2}$, $\epsilon\gg T$. We can therefore make the approximation that Eq. (66) gives the full ICS spectrum while neglecting Eq. (68), and take the integral limits to be $0\leq q\leq 1$ instead. This assumption is made in the ICS transfer functions provided as part of the downloaded data, but options are available in the ics module to turn these various assumptions off.

The quantity $\Gamma(\epsilon)$ separates the two kinematic regimes of Compton scattering: $\Gamma\gg 1$ for the Klein-Nishina regime, where Compton scattering in the electron rest frame is highly inelastic, and $\Gamma\ll 1$ for the Thomson regime, where it is almost elastic instead.101010Although the scattering process is almost elastic in the initial electron rest frame, it is certainly not elastic in the co-moving frame. In the co-moving frame, the electron loses a small fraction of its energy per collision, but each collision can upscatter a CMB photon by a significant factor. Eq. (66) applies to both regimes, with the only assumption being $\gamma\gg 1$.

To avoid computing the scattered photon spectrum repeatedly in the code, we use the following relation between spectra at different temperatures:

[TABLE]

for any real positive number $y$, even if $yE_{e}$ is unphysical.111111This trick can only be performed by integrating over $0\leq q\leq 1$, and is the key reason for making such an approximation. In DarkHistory, we evaluate the scattered photon spectrum at $1+z=400$, and use this relation to compute the subsequent spectra at lower redshifts by a straightforward interpolation.

A.1.2 Thomson Regime

In the Thomson regime, the rate at which photons are scattered is given by Blumenthal and Gould (1970)

[TABLE]

where $N_{\text{rad}}$ is the total number density of incident photons, with $\sigma_{T}=8\pi r_{0}^{2}/3$ being the Thomson cross section. Note that the scattering rate is independent on the incident photon energy. The energy loss rate of the electron is Blumenthal and Gould (1970)

[TABLE]

where $\beta$ is the velocity of the electron, with $U_{\text{rad}}$ being the total energy density of the incident photons.

While Eqs. (71) and (72) are well-known, the actual spectrum of scattered photons in the Thomson regime is much less so. The complete expression for the differential scattered photon spectrum with no further assumptions is, as far as the authors know, first given in Ref. Fargion et al. (1997), and we reproduce their final result here for completeness. For $(1-\beta)\epsilon_{1}/(1+\beta)<\epsilon<\epsilon_{1}$, we have

[TABLE]

and for $\epsilon_{1}<\epsilon<(1+\beta)\epsilon_{1}/(1-\beta)$,

[TABLE]

All other values of $\epsilon$ outside of the ranges specified are kinematically forbidden, and so to find the spectrum, we need to integrate over $\epsilon$ with $n(\epsilon)=n_{\text{BB}}(\epsilon)$ in the finite range specified above, i.e.

[TABLE]

The relationship between spectra at different temperatures is given by

[TABLE]

The scattered photon energy loss spectrum $dN_{\gamma}/(d\Delta\,dt)$ is similarly given by

[TABLE]

The relation shown in Eq. (76) between scattered photon spectra of different temperatures also holds for the energy loss spectrum, with $\epsilon_{1}\to\Delta$.

A.2 Numerical Methods

Computationally, to evaluate all of the scattered spectra, we need to perform numerical quadrature over a large range of electron and scattered photon energies; using a standard grid of $5000\times 5000$ energy values, the grid would take the standard numpy integrator over a day to populate. While a substantial speed-up may be obtained by using packages like Cython Behnel et al. (2011), numerical quadrature for ICS in the Thomson regime is also subject to significant numerical errors when the electron is nonrelativistic due to the existence of catastrophic cancellations. A semi-analytic approach provides both a faster method and a way to avoid such errors in a robust manner.

A.2.1 Thomson and Relativistic Regime: Large $\beta$

For $\beta\gtrsim 0.1$, we can obtain the scattered photon spectrum in Eq. (73) in the Thomson regime or Eq. (66) in the relativistic regime, as well as the scattered electron energy loss spectrum in the Thomson regime in Eq. (61), by direct integration.

The problem of integrating these expressions reduces to obtaining an expression for indefinite integrals over the Bose-Einstein distribution of the form

[TABLE]

Throughout this appendix, we ignore the constant of integration for such indefinite integrals, since we will ultimately be taking differences of such expressions to find definite integrals. For $f(y)\equiv y^{n}$ with integer $n\geq 0$, the indefinite integral is known explicitly:

[TABLE]

where $\text{Li}_{m}(z)$ is the polylogarithm function of order $m$ with argument $z$ (see Appendix A.4 for the definition). Note however that NumPy does not have a numerical function for the polylogarithm of order $m>2$, and so the semi-analytic method that we describe below is still necessary for $P_{y^{n}}$, $n\geq 2$ due to this limitation.

For other functions $f(y)$, closed-form solutions do not exist. However, an expression for the indefinite integral as an infinite series can be obtained Zdziarski and Pjanka (2013). Importantly, more than one series expression exists for all of the integrals $P_{f}(x)$ of interest in both the relativistic and nonrelativistic regimes, so that it is always possible to find a series expression that converges quickly for any integration limit. We tabulate the series expressions already found in Ref. Zdziarski and Pjanka (2013) for completeness, together with the many new series expressions derived in this paper required for the nonrelativistic limit in Appendix A.4.

A.2.2 Thomson Regime: Small $\beta$

In the Thomson regime for $\beta\lesssim 0.1$, catastrophic cancellations between terms in the integral make even the method described above insufficient. After integrating Eq. (75) over $\epsilon$ to get the scattered photon spectrum, for example, the final result must be $\mathcal{O}(\beta^{0})$, even though the prefactor in Eq. (73) is $\mathcal{O}(\beta^{-6})$. The integrals of all of the terms in the curly braces of Eq. (73) and their analog from Eq. (74) must therefore cancel among themselves to 1 part in $\beta^{-6}$; such a computation is impossible to perform for $\beta\lesssim 0.003$ due to floating point inaccuracy, even with double precision.

We avoid this problem by expanding the scattered photon spectrum in Eq. (75) and the mean electron energy loss spectrum in Eq. (61). Eq. (75) can be expanded straightforwardly in $\beta$, but Eq. (61) must be expanded in both $\beta$ and $\xi\equiv\Delta/T$, since catastrophic cancellations occur when either variable is small. In DarkHistory, we expand these expressions up to $\mathcal{O}(\beta^{6})$ and $\mathcal{O}(\xi^{6})$, but the precision of this calculation is systematically improvable by adding more terms to the code as desired. The exact expressions for the expansions, details of their derivations and several consistency checks for these expressions can be found in Appendix A.4.

A.3 Results

Figs. 10 and 11 show the scattered photon spectrum in the Thomson and relativistic regimes respectively as a function of electron energy, at a CMB temperature of $0.25\text{\,}\mathrm{e\kern-0.50003ptV}$, corresponding to a redshift of $z\approx 1065$ that is near recombination. By default, DarkHistory transitions between these two limits at $\gamma=20$. Fig. 12 shows the mean electron energy loss spectrum in the Thomson regime. Above $\gamma=20$, DarkHistory uses the approximation shown in Eq. (62). Finally, the computed secondary photon spectrum after completely cooling of all electrons and positrons through ICS is shown in Fig. 13.

All results shown here are computed using a $500\times 500$ grid of electron and photon energies/energy loss, and each can be completed under ten seconds on a typical personal computer.

A.4 Integrals and Series Expansions

We are now ready to detail the integrals and series expansions used in the numerical methods described above.

A.4.1 Bose Einstein Integrals

Each $P_{f}(x)$ that is of interest has a series that converges quickly for small values of $x$, and another that converges quickly for large values of $x$. DarkHistory by default chooses $x=2$ as the value to switch between the two expressions.

Suppose we approximate the indefinite integral in Eq. (78) $P_{f}(x)$ by the first $N$ terms of its series expression, which we denote $S_{N}(x)$. Let $S_{N}^{s}(x)$ and $S_{N}^{l}(x)$ be the series expressions we obtain for $x<2$ and $x\geq 2$ respectively. In all of the cases we are interested in, $S^{l}_{N\to\infty}(x\to\infty)=0$ (with the constant of integration taken to be zero) due to the exponential function in the denominator of the original integral, and so

[TABLE]

Then defining $\Delta S_{N}^{s,l}(a,b)=S_{N}^{s,l}(b)-S_{N}^{s,l}(a)$, the definite integral is evaluated as

[TABLE]

Terms are added sequentially until the next contribution to the full integral falls below a given relative tolerance; the default value for this tolerance used by DarkHistory is $10^{-10}$.

Before listing the series expressions, we must first introduce some notation that will be relevant. The numbers and analytic functions defined below are all well-known, but are often defined with different normalizations or given different names. We explicitly define all relevant functions used here for clarity.

$B_{n}$ are the Bernoulli numbers, defined through the following exponential generating function:

[TABLE]

with $B_{0,1,2,\cdots}=1,-1/2,1/6,\cdots$. Note that $B_{2j+1}=0$ for all integers $j>0$.

Next, we define the generalized exponential integrals

[TABLE]

and the closely related incomplete gamma function

[TABLE]

The polylogarithm of order $m$, denoted $\text{Li}_{m}(z)$, is defined as

[TABLE]

Finally, we define ${}_{2}F_{1}(a,b;c;z)$, the Gaussian hypergeometric function, as

[TABLE]

where $(x)_{n}\equiv x(x+1)\cdots(x+n-1)$ is the Pochhammer symbol. This function only appears in the form $R(n,x)\equiv\Re[_{2}F_{1}(1,n+1;n+2;x)]$, where $\Re$ denotes the real part; to avoid the slow evaluation of the hyp2f1 function in NumPy, we use instead the following relation:

[TABLE]

The list of all of the series expressions that we use, including those already derived in Zdziarski and Pjanka (2013), are shown in Tables 2 and 3 for $x<2$ and $x\geq 2$ respectively.

A.4.2 Nonrelativistic Thomson Limit: Small Parameter Expansion

The expression for the scattered photon spectrum in the Thomson limit, shown in Eq. (75), can be expanded in the small $\beta$ limit, to obtain

[TABLE]

where $x_{1}=\epsilon_{1}/T$, $A_{n}$ is a constant, and $P_{j,n}(x_{1})$ is some rational or polynomial function in $x_{1}$. These quantities are as follows. For $n=0$ i.e. $\mathcal{O}(\beta^{0})$,

[TABLE]

For $n=1$, i.e. $\mathcal{O}(\beta^{2})$,

[TABLE]

For $n=2$, i.e. $\mathcal{O}(\beta^{4})$,

[TABLE]

and finally for $n=3$, i.e. $\mathcal{O}(\beta^{6})$,

[TABLE]

Furthermore, when $x_{1}$ is small, it becomes numerically advantageous to expand Eq. (88) in $x_{1}$ as well, leaving a simple polynomial in $x_{1}$ and $\beta$, i.e.

[TABLE]

The values of $C_{n,j}$ are shown in Table 4.

Three checks can be performed to verify that this is indeed the correct expansion in $\beta$. First, taking $\beta\to 0$, the scattered photon spectrum simply becomes $dN_{\gamma}/(d\epsilon_{1}\,dt_{1})=n_{\text{BB}}(\epsilon_{1},T)\sigma_{T}c$, which is exactly the expected result for Thomson scattering in the rest frame of the electron: all photons simply scatter elastically at a rate governed by the Thomson scattering cross section, thus remaining in a blackbody distribution. Second, a more non-trivial check is to integrate Eq. (88) with respect to $\epsilon_{1}$, giving the total Thomson scattering rate given in Eq. (71). Since the scattering rate is independent of $\beta$, the $\mathcal{O}(\beta^{0})$ term in the series should integrate to exactly $\sigma_{T}cN_{\text{rad}}$ where $N_{\text{rad}}$ is the number density of the blackbody photons, while the other higher order terms should integrate to exactly zero. This is indeed the case for the series expansion shown here. Lastly, one can check that Eq. (88) agrees with the energy loss expression Eq. (72), by noting that

[TABLE]

where $u_{\text{BB}}(T)$ is the blackbody energy density with temperature $T$, i.e. the produced secondary photon spectrum must have the same energy as the upscattered CMB photons plus the energy lost from the scattering electron. This check has also been performed for the series expansions shown here.

For the scattered electron energy loss spectrum shown in Eq. (61), the small $\beta$ and $\xi$ expansion can be written as

[TABLE]

where $Q_{n,j}(x)$ is a polynomial, $A\equiv\Delta/(2\beta T)=\xi/(2\beta)$, and $R_{n}(A)$ is a sum of integrals of the form

[TABLE]

These integrals can be evaluated using the same methods detailed in Appendix A.4. The list of polynomials $Q_{n,j}$ and of $R_{n}(A)$ is given below. All expressions not listed should be taken to be zero. For $n=0$,

[TABLE]

For $n=1$,

[TABLE]

For $n=2$,

[TABLE]

And finally for $n=3$,

[TABLE]

These are all the terms necessary to work at order $\mathcal{O}(\beta^{6})$ and $\mathcal{O}(\xi^{6})$. As before, if $A$ becomes small, we should expand Eq. (95) as

[TABLE]

with the values of $D_{n,j}$ shown in Table 5.

These expressions are complicated, but can be checked in a similar fashion as the scattered photon spectrum by integrating over $\Delta\,d\Delta$ to obtain the mean energy loss rate of electrons scattering of a blackbody spectrum, given exactly in Eq. (72). Using the fact that

[TABLE]

one can verify that integrating the $\mathcal{O}(\beta^{6})$ expansion gives

[TABLE]

which is precisely the Taylor expansion of Eq. (72) in powers of $\beta$.

Appendix B Positronium Annihilation Spectra

The spin-triplet ${}^{3}S_{1}$ state of positronium annihilates to three photons, producing a photon spectrum per annihilation given by Ore and Powell (1949)

[TABLE]

where $x\equiv E_{\gamma}/m_{e}$. The kinematically allowed range is $0\leq x\leq 1$. Assuming that the formation of positronium by low energy positrons populates all of the degenerate ground states equally, the averaged photon spectrum per annihilation is

[TABLE]

Appendix C Cross Checks

C.1 Helium Deposition

In this section, we compare the various helium energy deposition methods discussed in Sec. III.5. We pick a dark matter candidate which decays to two photons with a lifetime of $3\text{\times}{10}^{24}\text{\,}\mathrm{s}$ as an example, but the results are similar across different dark matter masses and energy injection rates.

The user can switch between methods by passing the keyword parameter compute_fs_method to evolve() with the following strings for each method: (1) ’no_He’, (2) ’He_recomb’ and (3) ’He’, e.g.

⬇

helium_method_alt = main.evolve(

  DM_process=’decay’, mDM=1e8,

  lifetime=3e24, primary=’phot_delta’,

  start_rs=3000., backreaction=True,

  helium_TLA=True,

  compute_fs_method=’He_recomb’

)

Fig. 14 shows the helium ionization fraction $x_{\text{HeII}}$ as a function of redshift for each of the different methods. In method (1), $x_{\text{HeII}}$ is simply the baseline helium ionization fraction, which is almost entirely neutral once helium recombination is complete. No energy is assigned to helium iondization at all. Method (2) has no contribution to helium ionization from photons, since every ionized helium atom is assumed to recombine, producing a photon that photoionizes hydrogen instead (i.e. setting $q^{\gamma}_{\text{He}}=0$ in Eq. (42)). The helium ionization level therefore deviates from the baseline only from energy injection in the $\text{He}_{\text{ion}}$ channel from low-energy electrons. On the other hand, method (3) assumes that all helium atoms that get photoionized stay ionized, maximizing the amount of energy into $\text{He}_{\text{ion}}$ from photons (i.e. setting $q^{\gamma}_{\text{He}}=1-q$ in Eq. (42)). This explains the higher $x_{\text{HeII}}$ obtained.

Despite these differences in $x_{\text{HeII}}$, the evolution of $x_{e}$ remains almost identical, due to the fact that the total number of ionization events between both hydrogen and helium remains the same regardless of method used. This in turn ensures only a small difference in $T_{m}$ between the methods. The ionization and temperature histories for all three methods for the particular channel we have chosen are shown in Fig. 15. Users may bracket the uncertainty in the treatment of helium with methods (2) and (3).

C.2 Coarsening

In the absence of reionization, a coarsening factor of up to 32 has been found to yield a small relative difference of between 5–10% in the values of $f_{c}(z)$ across the full range of redshifts used in DarkHistory. With reionization, however, $T_{m}$ evolves more rapidly and attains larger values, and too much coarsening can lead to large absolute differences and somewhat larger relative differences in $T_{m}$, since we are averaging over the $T_{m}$ evolution over many redshift points. Fig. 16 shows the resultant temperature evolution as a function of redshift for the same model used in the previous section but with the default reionization model turned on, with coarsening factors of 1, 8 and 32. Once reionization starts, the difference in $T_{m}$ is $\sim 15\%$ for $n=32$ compared to the uncoarsened result at $z\sim 4$, corresponding to an absolute error of $\sim$5000\text{,}\mathrm{K}$$. Prior to reionization, the relative errors are slightly smaller at $\lesssim 10\%$.

We therefore recommend using a coarsening factor of up to 32 if no reionization models are used, depending on the level of precision desired, and to use coarsening with care once reionization is included. We also emphasize that when using coarsening, it is best to check for convergence by comparing the result with less coarsening.

C.3 $f_{c}(z)$ Contours

Fig. 17 show the computed $f_{c}(z)$ contours within DarkHistory for all channels of interest without any backreaction. The new $f_{c}(z)$ calculation by DarkHistory makes several small physics and numerical improvements over the previous calculation of these results Slatyer (2016); Liu et al. (2016), but still agree to within less than 10% when methodologies (cosmological parameters, methods of interpolation etc.) are standardized between the code used in Ref. Slatyer (2016) and DarkHistory. The new calculation also corrects a bug in earlier work in the treatment of prompt energy deposition from nonrelativistic and mildly relativistic injected electrons. This accounts for the bulk of the visible differences in Fig. 17 between the current contours and those of Refs. Slatyer (2016); Liu et al. (2016), which are most pronounced for DM annihilation/decay to electrons and low injected particle energies.

Appendix D Table of Definitions

Table LABEL:tab:variables shows a list of variables and their definitions for reference.