Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes
Abstract
We devise a new userfriendly tool interfaced with the Boltzmann code CLASS to deal with any kind of exotic electromagnetic energy injection in the universe and its impact on anisotropies of the Cosmic Microwave Background. It makes use of the results from standard electromagnetic cascade calculations develop in the context of WIMP annihilation, generalized to incorporate any injection history. We first validate it on a specific WIMP scenario, the Higgs Portal model, confirming that the standard effective onthespot treatment is accurate enough. We then analyze the more involved example of evaporating Primordial Black Holes (PBHs) with masses in the range g, for which the standard approximations break down. We derive robust CMB bounds on the relic density of evaporating PBHs, ruling out the possibility for PBHs with a monochromatic distribution of masses in the range g to represent all of the Dark Matter in our Universe. Remarkably, we confirm with an accurate study that the CMB bounds are several orders of magnitude stronger than those from the galactic gammaray background in the range g. A future CMB experiment like CORE+, or an experiment attempting at measuring the 21 cm signal from the Dark Ages could greatly improve the sensitivity to these models.
a]Patrick Stöcker,
a]Michael Krämer,
a]Julien Lesgourgues,
b]Vivian Poulin
\affiliation[a]Institute for Theoretical Particle Physics and Cosmology (TTK),
RWTH Aachen University, D52056 Aachen, Germany.
\affiliation[b]Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA.
\emailAddstoecker@physik.rwthaachen.de
\emailAddmkraemer@physik.rwthaachen.de
\emailAddvpoulin@jhu.edu
\emailAddlesgourg@physik.rwthaachen.de
TTKXX
1 Introduction
There are nowadays a wealth of observational evidence on a variety of scales in favor of the existence of 85% of matter in our universe in the form of a cold, noninteracting component called Dark Matter (DM) (1). Despite decades of experimental and theoretical efforts at unveiling this mystery, and many potential hints of its detection along the years, we are still lacking a clear nongravitational identification of the DM, which would potentially hold a lot of information about its nature. Nowadays, the only property commonly accepted is its relic density today (2). In the past, searches have focused on the WIMP paradigm, motivated by the theoretically appealing following observation: a single new particle whose mass and coupling to the standard model (SM) are close to those of the weak sector, will experience a standard chemical decoupling in the early universe (the freezeout process) naturally leading to such a relic density today. This has led to the development of three main channels of detection known as direct detection (where one tries to measure the recoil of a nucleon, eventually electron, being hit by a DM particle), indirect detection (where one tries to measure the antiparticles produced by residual annihilations of DM) and production at colliders (where one tries to measure the DM “missing energy” when reconstructing the invariant mass of particles produced during the collision). Remarkably, cosmological probes and in particular the Cosmic Microwave Background (CMB) can be used to perform both direct and indirect detection of DM. Indeed, nongravitational interactions between the DM and the SM may affect the growth of perturbations, leading to clear signatures on the CMB temperature and polarization power spectra (3); (4); (5); (6). On the otherhand, any energy injection from DM annihilation (or decay) affects the evolution of the free electron fraction and of the intergalactic medium temperature . For large enough energy injection rates, this mechanism could also have a strong impact on the CMB power spectra (7); (8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26).
It has been shown in the past that the phenomenology associated to exotic electromagnetic energy injection is extremely rich (see in particular Ref. (18); (27); (28) for a detailed analysis). Indeed, from the study of the CMB power spectra, it is possible to learn not only on the amount of electromagnetic (e.m.) energy injection but also on its time dependence, and thus to potentially distinguish between different scenarios leading to exotic e.m. energy injection.
In this paper, we develop and release a tool calculating the impact on the CMB of virtually any form of exotic electromagnetic energy injection. It basically “fills the gap” between existing public codes calculating, on the one hand, the development of the e.m. cascade in the plasma and its impact on the freeelectron fraction, and on the other hand, the CMB power spectra. It makes use of the results from Refs. (11); (29), developed in the context of WIMP annihilation, and it broadens their range of application to an arbitrary injection history, or in other words, to any particle injection spectrum and injection rate. With such a tool, it is possible to study for instance DM annihilation (including the impact of halo formation at low), DM decay (even when only a fraction of the total DM is unstable), or even Primordial Black Holes (PBHs) evaporation or accretion. We discuss in details two particular examples, which had only been studied using an effective parametrisation or some crude approximations: first, the scalar Higgs portal DM model, and second, lowmass evaporating PBHs with g. The first example allows us to validate our tool on a slightly nontrivial scenario and can be regarded as a “proofofprinciple”. It confirms previous results showing that in the context DM annihilation searches, one does not need to follow accurately the full redshift dependence of the energy deposition function (18). However, the second example is less trivial, given that in the case of an energy injection rate scaling like , any miscalculation of the energy deposited in the plasma leads potentially to large inaccuracies in CMB predictions, in particular when the peak of the energy injection is around recombination (28). Moreover, when PBHs evaporate, both the production rate and the spectrum of injected particles change with time, a complication which has always been neglected until now when calculating the CMB spectra.
The paper is structured as follows. In section 2, we first recall basic principles on e.m. energy injection and on its impact on the freeelectron fraction and IGM temperature. We then detail the content and the use of the new ExoCLASS branch of the Boltzmann code CLASS
2 The ExoCLASS branch of the Class code
ExoCLASS is a public branch of the CLASS code aimed at computing the CMB power spectra in the presence of an electromagnetic cascade (e.m. cascade) caused by any source of injection of electromagnetic energy. ExoCLASS implements several new features with respect to the main branch:

first, the user can specify several injection histories in the form of different expressions for the function
that will be defined in equation (6). In section 2.2 we will review different possibilities and their theoretical motivations. The main branch of CLASS only incorportates a simplified way to treat energy injection in the specific case of wave DM annihilation.

a new package written in python, called DarkAges, which role is to compute the efficiency function in different channels that will be defined later in section 2.2, based on the method developed in (11); (22); (29) and on the integral of equation (14). This module must be interfaced at least with the public tool kit of Ref. (29)
^{2} . It also needs some input for the spectrum of electrons, positrons, photons and other annihilation or decay products, that can be either computed from first principles (in this work we used the code PYTHIA^{3} v8.219 (31)), or taken as tabulated templates e.g. from the PPPC4DMID^{4} (32); (33). DarkAges can be used in two ways: as a standalone code, or as an internal CLASS module. We describe its structure and its basic usage in section 2.3.
A few days after the submission of this paper, the ExoCLASS branch will be available for cloning or download from the CLASS GitHub repository
2.1 Theoretical basics
The main impact of exotic energy injection on the CMB is related to a modification of the free electron fraction evolution. In particular, the modification of the Thompson optical depth and visibility function around and below will lead to very distinct features on the CMB power spectra, which just arise from modifications to the recombination equations.
We refer to electromagnetic energy injection as a generic word for the injection of photons and electromagnetically charged particles from the standard model (SM) following processes such as annihilations or decays. On the cosmological timescales of interest, all unstable particles will decay producing a primary spectrum of stable particles whose interaction with the cosmological plasma need to be accurately followed. Like in most of the literature, we restrict the analysis to the impact of injected positrons, electrons and photons. Indeed, neutrinos are basically invisible to the medium and simply carry away part of the energy. On the other hand, protons and antiprotons have been checked to loosen the bounds by about 10%, (36). Thus neglecting them only leads to slightly too conservative bounds, while permitting a significant reduction of the computing time.
Typically, the injected primary particles initiate an e.m. cascade by interacting with thermal photons, producing an increase in the number of nonthermal particles at the expense of a decrease in their average energy. When these extra particles cool down to energies of the order of a keV, they start interacting strongly with atoms of hydrogen (and subdominantly of helium (17)). To account for the ionization, excitation and heating of these atoms, we have to modify the equations governing the evolution of the fraction of free electrons, , taking into account both direct ionization and collisional excitation followed by photoionization by a CMB photon. At the same time, we must add to the equation for the evolution of the intergalactic medium (IGM) a term accounting for the associated heating, which has a feedback on the evolution of . Finally, at some point, the energy of the extra particles drops below the Lyman transition energy (10.2 eV). Then these particles are no longer able to interact with atoms and can be considered as “lost”. The threelevel atom approximation gives a good overall description of the processes at play, and can be fudged to achieve subpercent accuracy (37). In this approximation, the evolution equations of the free electron fraction and IGM temperature is governed by:
(1) 
where the and terms are the standard recombination and ionization rates given by
(2) 
The effective ionization rate can be decomposed as , where is the rate of direct ionization and that of excitation+ionization:
(3) 
while and are respectively the average ionization energy per baryon, and the Lyman energy. Finally, the rate at which the plasma is heated by DM decay or annihilation is defined as:
(4) 
We refer e.g. to the appendix of Ref. (28) for further definitions and more details on each of these coefficients. In CLASS, it is possible to use a fudged version of Recfast (38); (37) or HyRec code (34) to solve these recombination equations. The ExoCLASS branch proposes as a third possibility the use of CosmoRec (35).
All these processes have an impact on the CMB anisotropy angular power spectra through an enhancement of Thomson interactions between CMB photons and free electrons. In particular, the modified free electron fraction affects two very important thermodynamics quantities entering in the lineofsight solution of the CMB photon Boltzmann hierarchy (39), namely the Thompson optical depth and the visibility function ,
(5) 
In order to compute the three rates defined in equations (3, 4), we need to know the rate of energy density deposition in the plasma at redshift , , and how it is splitted between ionization, excitation of the Lyman transition, and heating of the IGM. It can be linked to the rate of energy injection per unit volume through three dimensionless functions :
(6) 
where the index runs over the three channels called (, , ) in equations (3, 4). The sum of the three components of should always be smaller than one, to account for the fraction of injected energy ending up in very low energy photons (10.2 eV) which are unable to interact.
In the standard formalism, the energy deposition functions can be obtained by convolving a given primary injection spectrum at redshift , , with a set of transfer functions defined as the fraction of the energy injected at that is deposited at in a given channel for a given particle . The main difficulty in solving the problem at hand is that, in order to correctly compute these transfer functions, it is necessary to follow the evolution of the daughter particle spectra over a very large range of energy and time scales. In particular, it has been shown in Ref. (11) that for WIMP annihilation, at redshift around and below recombination, the injected energy is usually not absorbed onthespot (at ), since at these times cooling processes are unefficient over a Hubble time scale (11). Technically, this means that the Boltzmann equations describing the evolution of the each particle spectrum cannot be solved in a stationary approximation, and thus require more involved resolution schemes such as Monte Carlo methods (e.g. (15); (16); (11)). Reference (28) extended this analysis to the case of decaying particles with arbitrary lifetimes. Through a detailed comparison of several energy deposition schemes, this work shows that an accurate determination of the energy deposition function is mandatory in order to get reliable predictions for the impact of energy injection on the CMB.
In this work, we use the most uptodate transfer functions computed in Ref. (23); (29). These functions were derived in the context of WIMP annihilation, but they actually have a wider range of application, as we shall discuss later. We refer the reader interested in understanding the resolution methods and the associated systematic uncertainties to (23); (29). We simply recall that one of the main assumptions in such calculations is that the altered free electron fraction will not significantly backreact onto the energy cascade evolution. One could improve over this approximation by computing some dependent transfer functions. We do not address such complications, given that no sign of exotic energy injection has been detected yet. Moreover, it has been shown in Ref. (26) that neglecting such a backreaction leads to conservative constraints.
2.2 Implementation of various injection histories in ExoCLASS
We have implemented in ExoCLASS most of the energy injection histories that have been studied so far:

DMannihilation, including the impact of halo formation at low: The impact of DM annihilation onto the CMB has been extensively studied in the literature, both in the smooth background (7); (8); (9); (10); (11); (14); (15); (16); (17); (18); (19); (20); (21); (22); (23) and within halos (8); (10); (12); (13); (14); (9); (21); (24); (25); (26). The rate of energy injection in the smooth background can be written as
(7) and involves the critical density today , the fractional DM density today , and the DM mass and annihilation crosssection . If DM is made of selfconjugated particles, such as Majorana fermions, one has , while if DM particles and antiparticles differ (as in the case of Dirac fermions) and are equally populated, . The main impact of structure formation is to enhance the average squared energy density with respect to the smooth background case, by an amount usually parametrized through a boost factor :
(8) Further details on this function can be found in the appendix C of Ref. (25). Currently, we have only implemented a simple parametrization based on the Halo model and on the PressSchechter formalism, but it would be trivial to generalise the expression of . The energy injection rate from annihilation in halos is thus:
(9) and the total energy injection rate is simply the sum of the two contributions. We describe how to make use of this injection history in ExoCLASS (and all other histories), the relevant key parameters and available options for the energy deposition in sec. B.
More details can also be found in the files explanatory.ini and README of the ExoCLASS branch.

DMdecay with arbitrary lifetimes from a fraction of the DM: The decay of exotic particles and its impact on the CMB has been scrutinized in the recent literature (40); (41); (42); (43); (22); (26); (44); (27); (28). The decay of a fraction of the DM with a constant rate and lifetime injects some energy density in the plasma at a rate
(10) Within CLASS, the user can pass the fraction of decaying DM and its lifetime in units of . More details on these parameters can be found in the file explanatory.ini. We follow the standard approach and define the deposited energy with respect to the injected energy of a longlived particle:
(11) This just means that we absorb the exponential factor within the definition of the function.

Accretion of matter onto heavy PBHs: Heavy () PBHs have received a lot of attention after the aLIGO discovery of binary black hole (BH) mergers of tens of solar masses (45); (46); (47); (48); (49). In particular, based on computation of the merging rate of PBHs today, it has been suggested that they could represent a large fraction (if not all) of the DM in our Universe (50); (51); (52); (53). Interestingly, such massive objects accrete matter, which heats up, gets eventually ionized and eventually emits highenergy radiation. This radiation can in turn leave distinct imprints on the CMB power spectra, eventually jeopardizing the success of CDM, and the CMB can thus put stringent bounds on the fraction of the DM made of PBHs. For black holes with relative density and mass , the total energy injection rate per unit volume is (54):
(13) where is the accretion luminosity, which is affected by large theoretical uncertainties. Usually, it is written as , where is the accretion rate and the radiative efficiency factor at which matter is converted to radiation. Moreover, it is necessary to specify an energy spectrum for the radiation, which is typically dominated by Bremsstrahlung emission up to hundreds of keV. All these quantities can be determined within a given accretion scenario, but the question of which scenario is most relevant is still unclear. We have implemented the spherical accretion scenario of Ref. (55), as well as its reevaluated version from Ref. (56), which corrected several mistakes. Recently, it has been noted that the hypothesis of spherical accretion might break down already around recombination. Accretion could thus happen in a disk, enhancing greatly the radiative efficiency and in turn the impact onto the CMB power spectra. We have also implemented a disk accretion scenario as advocated in Ref. (54), and we refer to this later article for details on these various scenarios.
2.3 The DarkAges module
The main role of the DarkAgespackage is to compute the energy deposition functions introduced above for any energy injection history, i.e. a spectra of injected particles at a given rate. These functions read in general
(14) 
where for an energy rate that scales like (e.g. for decays) and if it scales like (e.g. for wave annihilations). Note that we generalize the original formalism, by using the differential rate of the particle spectrum with time rather than a fixed reference spectrum. In general the time dependence can be factorized, but it is not always the case (e.g. for evaporating PBHs). The transfer functions have been tabulated for energy injection above 5 keV. However, for energy injection below that, we can safely extrapolate the transfer function down to 100 eV: It has been shown in Ref. (17) that the energy repartition functions are to an extremely good approximation independent of the initial particle energy in the range between 100 eV and a few keV. In fact, this behaviour is at the heart of the “low energy code” used by the authors of Ref. (29) to compute their transfer functions. Below 100 eV, the ionization efficiency starts to drop, and we conservatively cut the integral at the numerator at this energy.
The application of equation (14) to calculate from a given injected spectral rate and for a given scaling is the main purpose of the DarkAgesmodule. Apart from that, the module also includes its own interpolation class and routines to calculate the injected spectral rate for the scenarios given in section 2.2. For the sake of modularity and flexibility the module is intentionally written in Python, such that individual functions of the module can be easily used in a different context and to simplify the implementation of custom injection histories.
The DarkAgesmodule can be used, independently of ExoCLASS, as a standalone python routine in order to simply compute and print it in a file. Its usage is explained in details in App. A, and additional information can be found in the documentation of the module, which is located within the doc
folder. A basic example can also be found as a Jupyter notebook within the examples
folder.
Within the ExoCLASS branch, there are several ways to compute the impact of electromagnetic energy injection, typically leading to different level of accuracy and runtime. We detail the use of the ExoCLASS branch, the various new options available, and how to call the DarkAges module within the .inifile in App. B. Moreover, the ExoCLASS branch contains additional reionization parametrization, including an asymmetric reionisation as described in Ref. (57); (58) and a semianalytical resolution based on the model of star formation from Ref. (59) (introduced in Ref. (60) in the context of CMB studies). More details can be found in the file explanatory.ini within the CLASS folder of the ExoCLASS branch.
3 First application: the Higgs portal model
The scalar Higgs portal model (61); (62); (63) comprises the Standard Model and a real scalar field, , which is a singlet under all SM gauge groups and odd under a symmetry, . The scalar is thus stable and a viable dark matter candidate. The Lagrangian of the scalar Higgs portal model is given by
(15) 
After electroweak symmetry breaking, the mass and interaction terms of the scalar field become
(16) 
where , GeV. The physical mass of the scalar dark matter particle, , is given by . The quartic interaction affects the stability of the electroweak vacuum, but is not relevant for dark matter phenomenology. For our purposes, the model is therefore specified by only two free parameters, i.e. the dark matter mass, , and the coupling of dark matter to the SM Higgs, .
The scalar Higgs portal model is among the simplest UVcomplete dark matter models. While the model is arguably too simplistic in its minimal form, equation 16, a coupling between a gauge singlet dark sector and the Standard Model via a Higgs bilinear term is expected in a wide class of models, as is the only SM gauge singlet operator of mass dimension two.
The Higgs portal interaction provides an efficient dark matter annihilation mechanism. The Higgs portal model can thus describe the dark matter relic density, and it can be probed in indirect dark matter searches. Both the ray galactic center excess and a potential dark matter signal in cosmic rays antiprotons can be accommodated by the Higgs portal model (64). Furthermore, the Higgs portal interaction leads to invisible Higgs decays, , and a DMnucleon interaction through the exchange of a Higgs boson, which can be probed at the LHC and in direct detection experiments, respectively. For a recent discussion of the various constraints see e.g. (65); (66).
Within the Higgs portal model, dark matter annihilation into Standard Model final states proceeds through the exchange of a Higgs boson in the  and channel, and through the quartic interaction. Below the Higgs threshold, , only channel Higgs exchange is relevant, so that the relative weight of the different SM final states is determined by the SM Higgs branching ratios and independent of the portal coupling . For dark matter masses , final states contribute and the composition of the SM final state depends on the size of the . In figure 1 we display the relative size of the different SM annihilation channels as a function of the dark matter mass in the range 5 GeV 5 TeV. The Higgs portal coupling has been chosen such that for a given value of the correct relic density (2) is obtained. (We use micrOMEGAs (67) to calculate the relic density, see also (65).) To compute the annihilation spectrum, we combine the spectra of the various SM final states as provided in (68), weighted by their relative strength as displayed in figure 1. We also include the contributions of threebody final states from the annihilation into offshell gauge bosons, see (65).
The total energy deposition function (equation 14) of the scalar Higgs portal model is presented in figure 2. In the left panel we show a contour plot in the plane of the dark matter mass, , and the redshift, . The dependence of on the dark matter mass is small for the large redshifts, , which are probed by the CMB analysis. In figure 2, left panel, we also display the effective redshift, (dashed red line), defined by , where is the effective energy deposition factor as proposed in (23). The effective redshift is close to (solid white line) which has been used in most previous analyses to define an effective energy deposition. In the right panel of figure 2 we show the dependence of for various dark matter masses.
From our general analysis of the energy deposition function we can conclude that an effective energy deposition factor as proposed in Ref. (23) is an excellent approximation to derive CMB constraints on the Higgs portal model. We checked that this conclusion holds also for a future CMB experiments such as CORE+ (57). Moreover, the effective energy deposition is essentially identical for the different SM final states we consider, i.e. annihilation into quarks, gluons, gauge or Higgs bosons, and thus largely independent of the relative size of the various annihilation channels. Thus, the CMB constraints we derive for the Higgs portal model are applicable to a large class of models where direct annihilation into leptons or photons is suppressed.
To constrain the annihilation cross section using Planck data, we use MontePython (69) and perform a scan over the six CDM parameters
and two additional dark matter parameters
We take flat priors on these parameters, and restrict the dark matter mass to the range 5 GeV 5 TeV. Our data set consists of Planck high TT,TE,EE + simLOW (prior on ) and Planck CMB lensing (70). The results are shown in figure 3. Our analysis excludes large annihilation cross sections, as shown by the green shaded region. We also display the region of parameter space which corresponds to a overabundance of dark matter (red shaded area), which constrains the annihilation cross section from below.
4 Second application: evaporating Black Holes
4.1 Energy injection from Hawking radiation
We now consider a more involved case, still lacking an accurate treatment, namely the search for evaporating Primordial Black Holes (PBHs). PBHs are created by large density contrasts entering the horizon in the primordial Universe, when the fluid pressure is not able to counteract the gravitational force. PBHs can have a wide range of masses depending on the size of the horizon at collapse, and are often invoked as one of the most appealing alternatives to particle Dark Matter. Black Holes are expected to experience evaporation, i.e. loss of mass through the emission of the famous Hawking radiation, predicted by S. Hawking in 1975 for Schwarzschild black holes. The wellknown results by Hawking predict a Black Body radiation emitted at a temperature depending on the BH mass,
(17) 
The spectrum of particles with spin emitted by a black hole with temperature is given by
(18) 
(see (71)), where is the absorption crosssection. This function has in general a nontrivial form (see (72)), but it can be written in the optical (highenergy) limit as
(19) 
Given that this function quickly decrease in the opposite regime (73), we conservatively cut the spectrum at , just below the peak of the distributions for leptons and photons (71). We checked that slighlty varying this threshold doesn’t affect our conclusions. In the lowmass range, roughly below g, these BHs are expected to emit particles at energies and rates which can affect the CMB anisotropies, and can thus be looked for in CMB power spectra analysis. This has not been done before with good accuracy, mainly because both the spectrum and the injection rate evolve while the mass decreases. Hence, it is necessary to make use of numerical tools able to track the evolution of both the injected and deposited energy rates. This can be done fairly easily with the new ExoCLASS package. The PBH mass evolves according to (74)
(20) 
where the effective number of degrees of freedom is obtained by summing over all possible emitted particles,
(21) 
In this equation is the numbers of internal degrees of freedom of the particles , while the factor is an additional weighting factor taking into account the fact that particles of different spin and charge have slightly different blackbody spectra.
The last two factors in (21) are , the mass of a BH when its temperature is equal to the mass of the produced particles , and , which takes into account the shift between the peak of the blackbody distribution and the temperature.
In table 1 we give the value of these parameters for each of the particles that we consider.
Particles  

We show and the relative contribution of each particle species to the whole emission in figure (4).
We can see that within the mass range of interest, light quarks provide the leading contribution, even in the case where the black hole temperature is around the scale. This was overlooked in a previous analysis (28), leading to an overestimate of the bounds in the range to g.
One might wonder whether for temperature below the asymptotic freedom limit, BHs radiate fundamental particles or hadrons. We follow the argument in Refs. (72) and (74), and assume that a black hole only emit particles that are fundamental at the temperaturescale of the black hole. By this we mean that for temperatures below the QCDconfinement scale (roughly ), there are no quarks or gluons emitted, but rather pions that decay into photons and electrons. Above this threshold, we assume quarks and gluons to be emitted freely. In order to take this complication into account, we add to equation (21) an exponential suppression factor parametrized as
(22) 
where we choose and .
With this improved parametrization, primary electrons and photons are not strongly suppressed anymore, and the impact of BH evaporation on the CMB is stronger. The emission can mostly be described by electrons, photons, muons and pions (which subsequently decay into photons). Only for very low masses do quarks become the dominant contribution. The low mass quarks then hadronize and decay, eventually giving rise to photons and electrons. In this study for which particles are injected already at very low energies, it is not possible to make use of the PPPC4DMID spectra, valid only above 5 GeV. For the secondaries of the muons and pions, we thus make use of PYTHIA v8.219 (31), in which we turn off hadronization/radiative effects. We checked that neglecting these effects do not lead to any sizeable difference on the energy deposition functions even at high energy. This could be expected since these functions measure the calorimetric properties of the plasma. Hence, unlike for cosmic rays searches (33), we don’t need to compute the injection spectra in great details.
We are now ready to compute the relevant energy deposition functions per channel, which we plot in figure 7. As discussed in sec. 2, the injection rate of electromagnetic energy for low mass PBHs is given by
(23) 
Here is the fraction of the total emitted radiation energy that goes into electromagnetic particles. It is given by , where the sum giving runs over all particles excepted neutrinos. Equivalently, one could omit the factor in equation (23), but evaluate from equation (20) with neutrinos omitted in the sum of equation (21).
We plot on the left panel of figure 6 the residuals of the unlensed CMB TT and EE power spectra, taken between a universe with evaporating PBH and the vanilla CDM with . The density of PBH is set by the constraint at 95% C.L. obtained in this work and we show three different masses and . The effect of the PBH evaporation and energy injection is typical of any electromagnetic energy injection (see Ref. (28) for a review), although the exact details of the effect depend on the PBH mass. Basically, the increased freezeout fraction leads to additional Thomson scattering of photons off free electrons along the lineofsight, which manifests itself as a damping of anisotropies at high’s and an enhanced power in the low polarization spectrum. For the lowest masses which starts to evaporate earlier, the delayed recombination slightly shifts acoustic peaks and thus generates small wiggles (almost invisble for these tiny PBH fractions) at high multipoles in the residuals with respect to a standard CDM scenario.
On the right panel of the same figure, we plot the residuals of the CMB TT and EE power spectra now taken between a model where the impact of the energy deposition is computed with the effective “onthespot” approximation as advocated for decaying particles
4.2 Results
We perform our comparision with CMB data using ExoCLASS and MontePython (69). We run over the six CDM parameters plus (with a flat prior), for 10 different PBH masses distributed between g and g, as well as 4 masses between g and , where a greater accuracy is required. Our data set consists of Planck high TT,TE,EE + simLOW (prior on ) and Planck CMB lensing (70).
The results of the MCMC scan in the plane are shown in figure 8, together with constraints coming from the Extragalactic Gammaray Background (EGB) from Ref. (75). Remarkably, our constraints largely dominate in the range g to g. They are only slightly stronger than the approximate estimate from Ref. (28). On the other hand, they are a factor of a few weaker than the EGB ones in the range g to g. Note that, although a bit rough
We also perform the same analysis for a future CMB experiment whose specifications are identical to the recently proposed CORE+ experiment with a sky coverage (57). The results are shown as the blue band. One can see that, as found in similar analyses in the context of DM decay (28), up to a one order of magnitude improvement on the sensitivity to is expected. Interestingly, the CMB constraints would then be comparable to current EGB constraints even in the highest mass range.
5 Conclusions
In this paper, we have developed and released a tool allowing one to study the signatures and constraints associated to virtually any form of exotic electromagnetic energy injection. The key quantity to estimate in order to get a reliable prediction of the impact of the energy injection onto the CMB is the energy deposition function per channel. We made use of the results from Ref. (11); (29) obtained in the context of WIMP annihilation, and generalized them to study injection with arbitrary spectra of particles at a given rate. Together with this paper, we release the new ExoCLASS branch of the Boltzmann code CLASS. The new branch already incorporates DM annihilations (including the impact of halos (60)), decays from a fraction of DM (28), PBH evaporation and accretion (both spherical and disk (54)). The ExoCLASSbranch involves a new module called DarkAges, which is easy to modify in order to incorporate any alternative history that one would want to study. Hence, it makes it a very useful tool for the search of DM particle(s) beyond the WIMP paradigm.
We have applied our tools to two particular examples which had not been accurately analysed before: the scalar Higgs portal and low mass (g) evaporating PBHs. The first study allowed us to validate our tool against the standard “effective onthespot” approximation. We have derived the most uptodate cosmological constraints on the scalar mass and Higgs portal coupling to the SM: when requiring the relic density not to overshoot the current measurements, we can identify the surviving part of the parameter space which will partly be probed by the next generation of CMB experiments. Note that the limits derived from energy injection are rather generic and hold for a wide class of models where dark matter annihilates predominantly into gluons, quarks, gauge or Higgs bosons. In the second study concerning evaporating PBHs, we have shown that CMB constraints largely dominate in the range g to g, validating the results of Ref. (28). However, they are slightly weaker than EGB constraints above these masses. Still, this represents an independent cosmological constraints on the presence of evaporating PBHs, which rules out the possibility for them to comprise a large fraction of the DM for PBHs with monochromatic distribution of masses in the range g. A CMB nextgeneration experiment such as CORE+ on the other hand could provide constraints competitive with the EGB ones.
The ExoCLASS branch represents a step toward a more accurate characterisation of the impact of a broad range of DM models onto the reionization history and CMB power spectra. This is a key aspect in order to precisely measure properties of the particle or astrophysical DM candidate, if one day a signal was to be found. The next step in this characterization will be to improve the accuracy of the energy deposition tool and interface with the recombination history, today valid at the level (17). Also, it would be interesting and straightforward to use our tools to study the impact of energy injection on the 21cm line from the epoch of reionization and the Dark Ages, a channel known to be extremely sensitive to these scenarios (28). Moreover, it is important to get our tools ready given that numerous experiments such as PAPER 64
We thank warmly Pasquale D. Serpico for very useful discussions and comments. We thank Jan Heisig for providing us the standard model branching ratio of the scalar Higgs Portal model. This work has been partly done thanks to the facilities offered by the Université Savoie Mont Blanc MUST computing center. Parts of the simulations were performed with computing resources granted by RWTH Aachen University under project thes0264. PS is supported by the DFG Emmy Noether Grant No. KA 4662/11.
Appendix A Description and usage as a standalone package
Within the package folder, any call to the DarkAges script can be made by simply typing./bin/DarkAges
. To specify the spectra and the injection history, the script can be called with additional arguments. These are:

hist
This option is used to specify the injection history, as described in previous section. The valid options so far areannihilation
for wave annihilation,annihilation_halos
for wave annihilation with the impact of the halo boost ,decay
for a decaying dark matter component with constant lifetime,evaporating_PBH
for light primordial black holes, andaccreting_PBH
for heavy primordial black holes.
Depending on the history chosen, additional options need to be set.
For the case of decaying dark matter (
hist=decay
) the lifetime (in seconds) needs to be set withtdec
. 
For the case of halo boosted dark matter annihilation (
hist=annihilation_halos
) the parameters and for the parameterization of (equation (C.5) of (25)) need to be set withzh
andfh
) . 
For heavy primordial black holes (
hist=accreting_PBH
), the accretion recipe used to compute the energy injection needs to be specified withaccretion_recipe
. The options to choose arespherical_accretion
ordisk_accretion
.


m / mass
This is the mass of the dark matter candidate. For the scenarios included so far the units are for annihilating and decaying dark matter, for light (evaporating) primordial black holes, and for heavy primordial black holes.
Alternatively, it is also possible to specify the logarithm to base 10 of the mass with the optionlog10mass
. 
s / spectrum
In scenarios in which the injection spectra are not automatically set by the model,the injected spectrum of electrons, positrons and photons need to be given. As of now, the scenarios where this is needed are
annihilation
,annihilation_halos
, anddecay
. The code can work with multiple spectra which are given as a list, separated by blanks. There are three possibilities to pass the spectrum.
If the spectra of electrons/positrons and photons are stored in a file, the path to this file can be used as an input for the spectrum. In its default setting the code assumes the file to contain the following table,
… … … … … If more convenient, it is the possible to use rather than and rather than , by passing additional keyword arguments to the interpreter. For details we refer to the documentation and the examples provided in the directory of the DarkAges module. If the spectra in this table are given for multiple masses (masses need to be ordered in growing order), the spectra will be automatically interpolated for masses within the given range.

If the input is either
dirac_electron
ordirac_photon
, a dirac distribution at the energy given by the dark matter mass (for annihilation) or half the dark matter mass (for decay) will be used as the spectrum of and .^{18} .


b / branching
When more than one spectrum is used inspectrum
, the relative branching ratio between the channels need to be given separated by blanks. The list needs to have the same nuber of entries thanspectrum
and they need to add up exactly to one.
Per default, the spectra are sampled at the energies given in the tables of the transfer functions. The user has also the possibilty to give a custom energy range (equally spaced in logarithmic space). The options to be used there are:

log10Emin
The lower bound of the energy table (the energy is given in units of ). 
log10Emax
The upper bound of the custom energy table. 
nbins_table
Number of bins of the custom energy table.
In case the range asked by the user goes beyond the energy range sampled by the transfer functions, it is currently assumed that these are constant, equal to the very lowest and very highest point of the energy grid. As discussed in section 2, this is known to be a very good approximation for the lowest energies (up to eV, we conservatively set the transfer function to zero below). Inspection of the transfer functions at the highest energies show that it is also fairly good for these, although a dedicated study would be needed to safely go above the hundreds of TeV scale.
This concludes the use of the basic usage of the DarkAges module. However, when analyzing particular dark matter models, the preferable choice of inputs are often parameters of the model (e.g. masses of the dark matter particles and couplings) rather than the use of branching
and specfile
as inputs for every parameter point. For that reason we include the possibility to define custom physics models which translate the custom model parameters into given spectra and branching ratios. In that case, the user needs to provide the spectra and branching ratios on a grid of the model parameters which will later be interpolated.
In order to reduce the execution time of the code, this grid is read once in a preparation step at which the interpolation tables are initiated and saved. Later on, the execution step is performed and the spectra and branching ratios are sampled at the required point in parameter space.
This mode can be invoked via model=MODEL
. Here, MODEL
defines the name of the folder where the routines of the physics model are located (./models/MODEL/
). Inside this folder the code expects at least one Python script which defines the following routines

prepare()
This routine reads the branching ratios and spectra given on a grid of the model parameters and initiates the interpolations and saves them. This routine will only be run once, when DarkAges accesses the model for the first time. At a later point this step can only be accessed by forcing it with using the optionrebuildmodel
. 
run(*arguments, **DarkOptions)
This routine takes the interpolations of the spectra and branching ratios performed and saved by theprepare()
routine and samples them at the model parameters. Here,*arguments
is a list of the custom model parameters and**DarkOptions
are additional keyword arguements, such as precision parameters or handles for nondefault behaviour of the functions within DarkAges.
For more details we refer the reader to the documentation and to the “toy model” located at ./models/simple_mix
.
To pass additional options to the functions of the DarkAges module, most functions take additional keyword arguemts which are located within the globally accessible dictionary DarkOptions
. The option extraoptions=EXTRA.yaml
can be used to access and update this dictionary. Here, EXTRA.yaml
is a file in the YAMLDarkOptions
dictionary we refer to the documentation of DarkAges.
Appendix B Exotic energy injection within ExoCLASS and call to DarkAges
Within ExoCLASS, there are several ways to compute the impact of electromagnetic energy injection, typically leading to different level of accuracy/runtime. In order to call a specific energy injection history and treatment of the energy deposition, one needs to specify the following options in the standard input file.
Key parameters describing the injection history: The type of injection history (as described in section 2) will be set by key parameters of each history. Those are:

For annihilation, the DM mass and crosssection respectively specified with DM_mass (in GeV) and annihilation_cross_section (in cm/s). Alternatively, the user can choose to pass the usual annihilation parameter [m/s/kg], called annihilation within CLASS. Moreover, when calling the DarkAges module (its usage within ExoCLASS will be described below), one can specify the injection spectra and corresponding branching ratio thanks to the options injected_particle_spectra and injected_particle_branching_ratio respectively. The usage is the same as the spectrum and branching options described in sec. A when DarkAges is used as a standalone routine.

For annihilation in halos, one can additionnally pass the amplitude of the halo with the option annihilation_f_halo and the redshift at which halos start to form with annihilation_z_halo. Their exact definition can be found in Ref. (60), appendix C.

For decay, the DM mass, lifetime and fraction (defined with respect to the total DM). Those are specified respectively by DM_mass (in GeV), tau_dcdm (in s) and decay_fraction. Alike the annihilation case, one can additionnally pass the injected particle spectra and branching ratio when calling DarkAges.

For PBH evaporation, the PBH mass (in g) is specified via PBH_evaporating_mass and the PBH fraction of the total DM via PBH_fraction.

For PBH accretion, the PBH mass (in ) is specified via PBH_accreting_mass, the PBH fraction of the total DM via PBH_fraction. Additionnaly, on needs to specify the accretion recipe thanks to the parameter PBH_accretion_recipe, which can be one of spherical_accretion or disk_accretion.
Parameters setting the treatment of the energy deposition: In ExoCLASS, it is possible to pass several parameters adjusting the treatment of the energy deposition.

on the spot: First and foremost, just like in CLASS, one can decide to work in the “onthespot” approximation by setting the parameter on the spot = yes. As recalled in section 2, this approximation essentially consists in assuming that the deposition happens at the injection redshift. It is known to be a good approximation in the context of DM annihilation or longlived decay as long as the injected energy is weighted by a factor adequate for the model under study. There are several ways to specify this parameter. One can first set a constant parameter f_eff within the input file. However, given the degeneracy of this parameter with other parameters describing the energy injection, one can simply give the parameters relevant to the energy injection (annihilation_cross_section, decay_fraction, PBH_fraction) rescaled by the value of f_eff. In that case f_eff can be set to one, or ignored since it is the default value. This is for instance useful in the context of a MCMC runs with MontePython. Finally, one can follow the approach of Refs. (11) later called “3 keV prescription” in Ref. (29), which amounts in factorizing the energy deposition function into a redshiftdependent function (describing the highenergy cooling) and another freeelectronfraction dependent one (describing the repartition in each channel). To do so, one can specify the path to a file containing the redshiftdependent in the form
number of lines … … thanks to the command energy injection f_eff file=.... This approach has for instance been used in Refs. (76); (24); (60); (28).

energy_repartition_treatment
: When working in the effective on the spot picture, it is necessary to specify energy repartition functions often noted . In CLASS we have implemented an old prescription advocated by Shull, Van Steenberg, Chen and Kamionkowski (SSCK) (77); (40) in which the energy repartitions functions are and . In this prescription there is no “lost” photons and no differentiation between hydrogen and helium reionisation. There is also the possibility to use a more accurate description as computed by Galli, Slatyer, Valdes and Iocco (GSVI) (17). The user can also specified the path to a file containing a table of the formnumber of lines … … … … … … These prescriptions can be chosen by setting
energy_repartition_treatment =
SSCK, GSVI or from_file. The path to the file is specified within the input file via the parameterenergy repartition coefficient file=...
.When working beyond the onthespot approximation, it can be useful to bypass the energy repartition function, since those can be directly included in the energy repartition functions . This is done by setting
energy_repartition_treatment =
no_factorization
, in which case the ’s functions are effectively set to one within the code. 
energy_deposition_function:
This option can be used when working beyond the onthespot treatment, i.e. when one sets the option on the spot = no in the input file. The argument can be one of the following
from_file
: The user might find useful to give functions from an existing file. Its location needs to be given withenergy deposition function file=...
. The file needs to be in the formnumber of lines … … … … … … 
DarkAges
: Last but not least, it is possible to call the DarkAges package by settingenergy_deposition_function=DarkAges
in the input file. In that case, ExoCLASS automatically sets the mandatory parameterson the spot = no
andenergy_repartition_treatment = no_factorization
.
Calling the DarkAges module: Apart from its use as a standalone package to compute the functions, the DarkAgesmodule is designed to be used in conjuction with CLASS. It is called by setting energy_deposition
_function=DarkAges
in the input file.
In figure 9 we show the programflow of ExoCLASS: all relevant parameters, i.e. the cosmological parameters used by CLASS, as well as the input parameters relevant for the correct input of the DarkAgescomandline script, are passed as input parameters of CLASS in a standard .ini
file. The call to the module is done in the thermodynamics.c file. DarkAges then computes and return the functions to ExoCLASS which are used only within this module. If needed, the user can print the output of the computation in a file with the option print_energy_deposition_function = yes
.
The DarkAgesmodule can be used in several modes, which are specified via the option DarkAges_mode = MODE. Those were described in previous subsection and are called within ExoCLASS with the following keywords

built_in
: In this mode, ExoCLASS sets up an automatic call to the DarkAgesscript given the scenario of energy injection the user is interested in. From the key parameters of the injection history that are given to ExoCLASS, the code automatically detects what history the user is computing. 
user_command
: If the call to the DarkAgesscript cannot be performed with one of the builtin methods, this setting allows the user to have full control on the python commandcalled by ExoCLASS to get the table. The command is given via the optionext_fz_command=...
and there is the possibility to use up to five additional floating point parametersext_fz_par1
, …,ext_par_fz5
which are meant to be used as varying parameters when using ExoCLASS with MontePython (those can be for instance some relevant couplings). The full command which is executed by CLASS is then set up by the content ofext_fz_command
and the five additional valuesext_fz_par
are appended to it separated by blanks.
Thanks to clever preparation and interpolation routines, it is perfectly possible to perform intense MCMC runs using ExoCLASS and MontePython. ExoCLASS is especially fast in conjuction with the fuged version of Recfast implemented in CLASS, which has been shown to be accurate enough to describe the recombination history from detailed comparaison with CosmoRec and HyRec. This is illustrated on figure 10 in the context of PBH evaporation, where we plot the residuals of the CMB TT and EE power spectra computed using Recfast and HyRec for three different PBH masses and abundance set by the constraint at 95%C.L. obtained in this work. The difference is reassuringly well below cosmic variance.
Footnotes
 https://www.classcode.net
 http://nebel.rc.fas.harvard.edu/epsilon/
 http://home.thep.lu.se/~torbjorn/pythia8/pythia8219.tgz
 http://www.marcocirelli.net/PPPC4DMID.html
 http://cosmo.nyu.edu/yacine/hyrec/hyrec.html
 http://www.cita.utoronto.ca/~jchluba/Science_Jens/Recombination/CosmoRec.html
 https://github.com/lesgourg/class_public
 Those have a redshift dependence of the injection history close to that of evaporating PBH.
 we estimate g to be a more accurate cutoff, with a very fast loosening of the bound below g.
 http://eor.berkeley.edu
 http://21cma.bao.ac.cn
 http://www.mwatelescope.org
 http://www.lofar.org
 http://reionization.org
 http://www.skatelescope.org
 Within the CLASS folder the exact command will be ./DarkAgesModule/bin/DarkAges.
 The light quarks are treated as one particle
 Please note, when using a diraclike spectrum, it is not yet possible to combine them with spectra of the first two kinds.
 See http://www.yaml.org
References
 G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept., vol. 405, pp. 279–390, 2005.
 N. Aghanim et al., “Planck intermediate results. XLVI. Reduction of largescale systematic effects in HFI polarization maps and estimation of the reionization optical depth,” Astron. Astrophys., vol. 596, p. A107, 2016.
 C. Dvorkin, K. Blum, and M. Kamionkowski, “Constraining Dark MatterBaryon Scattering with Linear Cosmology,” Phys. Rev., vol. D89, no. 2, p. 023519, 2014.
 R. J. Wilkinson, J. Lesgourgues, and C. Boehm, “Using the CMB angular power spectrum to study Dark Matterphoton interactions,” JCAP, vol. 1404, p. 026, 2014.
 R. J. Wilkinson, C. Boehm, and J. Lesgourgues, “Constraining Dark MatterNeutrino Interactions using the CMB and LargeScale Structure,” JCAP, vol. 1405, p. 011, 2014.
 V. Gluscevic and K. K. Boddy, “Constraints on scattering of keV–TeV dark matter with protons in the early Universe,” 2017.
 N. Padmanabhan and D. P. Finkbeiner, “Detecting dark matter annihilation with CMB polarization: Signatures and experimental prospects,” Phys.Rev., vol. D72, p. 023508, 2005.
 A. V. Belikov and D. Hooper, “How Dark Matter Reionized The Universe,” Phys.Rev., vol. D80, p. 035007, 2009.
 M. Cirelli, F. Iocco, and P. Panci, “Constraints on Dark Matter annihilations from reionization and heating of the intergalactic gas,” JCAP, vol. 0910, p. 009, 2009.
 G. Huetsi, A. Hektor, and M. Raidal, “Constraints on leptonically annihilating Dark Matter from reionization and extragalactic gamma background,” Astron. Astrophys., vol. 505, pp. 999–1005, 2009.
 T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, “CMB Constraints on WIMP Annihilation: Energy Absorption During the Recombination Epoch,” Phys.Rev., vol. D80, p. 043526, 2009.
 A. Natarajan and D. J. Schwarz, “The effect of early dark matter halos on reionization,” Phys.Rev., vol. D78, p. 103524, 2008.
 A. Natarajan and D. J. Schwarz, “Dark matter annihilation and its effect on CMB and Hydrogen 21 cm observations,” Phys.Rev., vol. D80, p. 043529, 2009.
 A. Natarajan and D. J. Schwarz, “Distinguishing standard reionization from dark matter models,” Phys.Rev., vol. D81, p. 123510, 2010.
 M. Valdes, C. Evoli, and A. Ferrara, “Particle energy cascade in the Intergalactic Medium,” Mon. Not. Roy. Astron. Soc., vol. 404, pp. 1569–1582, 2010.
 C. Evoli, M. Valdes, A. Ferrara, and N. Yoshida, “Energy deposition by weakly interacting massive particles: a comprehensiv e study,” Mon. Not. Roy. Astron. Soc., vol. 422, pp. 420–433, 2012.
 S. Galli, T. R. Slatyer, M. Valdes, and F. Iocco, “Systematic Uncertainties In Constraining Dark Matter Annihilation From The Cosmic Microwave Background,” Phys.Rev., vol. D88, p. 063502, 2013.
 D. P. Finkbeiner, S. Galli, T. Lin, and T. R. Slatyer, “Searching for Dark Matter in the CMB: A Compact Parameterization of Energy Injection from New Physics,” Phys.Rev., vol. D85, p. 043522, 2012.
 G. Hutsi, J. Chluba, A. Hektor, and M. Raidal, “WMAP7 and future CMB constraints on annihilating dark matter: implications on GeVscale WIMPs,” Astron.Astrophys., vol. 535, p. A26, 2011.
 T. R. Slatyer, “Energy Injection And Absorption In The Cosmic Dark Ages,” Phys.Rev., vol. D87, no. 12, p. 123513, 2013.
 G. Giesen, J. Lesgourgues, B. Audren, and Y. AliHaïmoud, “CMB photons shedding light on dark matter,” JCAP, vol. 1212, p. 008, 2012.
 S. Galli, T. R. Slatyer, M. Valdes, and F. Iocco, “Systematic Uncertainties In Constraining Dark Matter Annihilation From The Cosmic Microwave Background,” Phys.Rev., vol. D88, p. 063502, 2013.
 T. R. Slatyer, “Indirect Dark Matter Signatures in the Cosmic Dark dAges I. Generalizing the Bound on swave Dark Matter Annihilation from Planck,” 2015.
 L. LopezHonorez, O. Mena, S. PalomaresRuiz, and A. C. Vincent, “Constraints on dark matter annihilation from CMB observationsbefore Planck,” JCAP, vol. 1307, p. 046, 2013.
 V. Poulin, P. D. Serpico, and J. Lesgourgues, “Dark Matter annihilations in halos and highredshift sources of reionization of the universe,” JCAP, vol. 1512, no. 12, p. 041, 2015.
 H. Liu, T. R. Slatyer, and J. Zavala, “The Darkest Hour Before Dawn: Contributions to Cosmic Reionization from Dark Matter Annihilation and Decay,” Phys. Rev., vol. D94, no. 6, p. 063507, 2016.
 T. R. Slatyer and C.L. Wu, “General Constraints on Dark Matter Decay from the Cosmic Microwave Background,” 2016.
 V. Poulin, J. Lesgourgues, and P. D. Serpico, “Cosmological constraints on exotic injection of electromagnetic energy,” JCAP, vol. 1703, no. 03, p. 043, 2017.
 T. R. Slatyer, “Indirect Dark Matter Signatures in the Cosmic Dark Ages II. Ionization, Heating and Photon Production from Arbitrary Energy Injections,” 2015.
 D. Blas, J. Lesgourgues, and T. Tram, “The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes,” JCAP, vol. 1107, p. 034, 2011.
 T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, “An Introduction to PYTHIA 8.2,” Comput. Phys. Commun., vol. 191, pp. 159–177, 2015.
 M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, “PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection,” JCAP, vol. 1103, p. 051, 2011. [Erratum: JCAP1210,E01(2012)].
 P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia, and A. Urbano, “Weak Corrections are Relevant for Dark Matter Indirect Detection,” JCAP, vol. 1103, p. 019, 2011.
 Y. AliHaïmoud and C. M. Hirata, “HyRec: A fast and highly accurate primordial hydrogen and helium recombination code,” Phys. Rev. D, vol. 83, p. 043513, Feb. 2011.
 J. Chluba and R. M. Thomas, “Towards a complete treatment of the cosmological recombination problem,” Mon. Not. Roy. Astron. Soc., vol. 412, p. 748, 2011.
 C. Weniger, P. D. Serpico, F. Iocco, and G. Bertone, “CMB bounds on dark matter annihilation: Nucleon energylosses after recombination,” Phys.Rev., vol. D87, no. 12, p. 123008, 2013.
 J. A. RubinoMartin, J. Chluba, W. A. Fendt, and B. D. Wandelt, “Estimating the impact of recombination uncertainties on the cosmological parameter constraints from cosmic microwave background experiments,” Mon. Not. Roy. Astron. Soc., vol. 403, p. 439, 2010.
 S. Seager, D. D. Sasselov, and D. Scott, “A new calculation of the recombination epoch,” Astrophys. J., vol. 523, pp. L1–L5, 1999.
 U. Seljak and M. Zaldarriaga, “A Line of sight integration approach to cosmic microwave background anisotropies,” Astrophys. J., vol. 469, pp. 437–444, 1996.
 X.L. Chen and M. Kamionkowski, “Particle decays during the cosmic dark ages,” Phys. Rev., vol. D70, p. 043502, 2004.
 S. Kasuya and M. Kawasaki, “Early reionization by decaying particles in light of three year WMAP data,” JCAP, vol. 0702, p. 010, 2007.
 L. Zhang, X. Chen, M. Kamionkowski, Z.g. Si, and Z. Zheng, “Constraints on radiative darkmatter decay from the cosmic microwave background,” Phys.Rev., vol. D76, p. 061301, 2007.
 S. Yeung, M. Chan, and M.C. Chu, “Cosmic Microwave Background constraints of decaying dark matter particle properties,” Astrophys.J., vol. 755, p. 108, 2012.
 I. M. Oldengott, D. Boriero, and D. J. Schwarz, “Reionization and dark matter decay,” JCAP, vol. 1608, no. 08, p. 054, 2016.
 B. P. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett., vol. 116, no. 6, p. 061102, 2016.
 B. P. Abbott et al., “GW151226: Observation of Gravitational Waves from a 22SolarMass Binary Black Hole Coalescence,” Phys. Rev. Lett., vol. 116, no. 24, p. 241103, 2016.
 B. P. Abbott et al., “GW170104: Observation of a 50SolarMass Binary Black Hole Coalescence at Redshift 0.2,” Phys. Rev. Lett., vol. 118, no. 22, p. 221101, 2017.
 B. P. Abbott et al., “GW170814: A ThreeDetector Observation of Gravitational Waves from a Binary Black Hole Coalescence,” Phys. Rev. Lett., vol. 119, no. 14, p. 141101, 2017.
 B. P. Abbott et al., “GW170608: Observation of a 19solarmass Binary Black Hole Coalescence,” 2017.
 S. Bird, I. Cholis, J. B. Muñoz, Y. AliHaïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, “Did LIGO detect dark matter?,” Phys. Rev. Lett., vol. 116, no. 20, p. 201301, 2016.
 M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, “Primordial Black Hole Scenario for the GravitationalWave Event GW150914,” Phys. Rev. Lett., vol. 117, no. 6, p. 061101, 2016.
 S. Clesse and J. GarcíaBellido, “The clustering of massive Primordial Black Holes as Dark Matter: measuring their mass distribution with Advanced LIGO,” Phys. Dark Univ., vol. 15, pp. 142–147, 2017.
 S. Clesse and J. GarcíaBellido, “Seven Hints for Primordial Black Hole Dark Matter,” 2017.
 V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, “CMB bounds on diskaccreting massive primordial black holes,” Phys. Rev., vol. D96, no. 8, p. 083524, 2017.
 M. Ricotti, J. P. Ostriker, and K. J. Mack, “Effect of Primordial Black Holes on the Cosmic Microwave Background and Cosmological Parameter Estimates,” Astrophys. J., vol. 680, p. 829, 2008.
 Y. AliHaïmoud and M. Kamionkowski, “Cosmic microwave background limits on accreting primordial black holes,” Phys. Rev., vol. D95, no. 4, p. 043534, 2017.
 E. Di Valentino et al., “Exploring Cosmic Origins with CORE: Cosmological Parameters,” 2016.
 R. Adam et al., “Planck intermediate results. XLVII. Planck constraints on reionization history,” Astron. Astrophys., vol. 596, p. A108, 2016.
 B. E. Robertson, R. S. Ellis, S. R. Furlanetto, and J. S. Dunlop, “Cosmic Reionization and Early Starforming Galaxies: a Joint Analysis of new Constraints From Planck and the Hubble Space Telescope,” Astrophys. J., vol. 802, no. 2, p. L19, 2015.
 V. Poulin, P. D. Serpico, and J. Lesgourgues, “Dark Matter annihilations in halos and highredshift sources of reionization of the universe,” JCAP, vol. 1512, no. 12, p. 041, 2015.
 V. Silveira and A. Zee, “SCALAR PHANTOMS,” Phys. Lett., vol. B161, p. 136, 1985.
 J. McDonald, “Gauge singlet scalars as cold dark matter,” Phys. Rev., vol. D50, pp. 3637–3649, 1994.
 C. P. Burgess, M. Pospelov, and T. ter Veldhuis, “The Minimal model of nonbaryonic dark matter: A Singlet scalar,” Nucl. Phys., vol. B619, pp. 709–728, 2001.
 A. Cuoco, J. Heisig, M. Korsmeier, and M. Krämer, “Probing dark matter annihilation in the Galaxy with antiprotons and gamma rays,” JCAP, vol. 1710, no. 10, p. 053, 2017.
 A. Cuoco, B. Eiteneuer, J. Heisig, and M. Krämer, “A global fit of the ray galactic center excess within the scalar singlet Higgs portal model,” JCAP, vol. 1606, no. 06, p. 050, 2016.
 P. Athron et al., “Status of the scalar singlet dark matter model,” Eur. Phys. J., vol. C77, no. 8, p. 568, 2017.
 G. Bélanger, F. Boudjema, A. Pukhov, and A. Semenov, “micrOMEGAs4.1: two dark matter candidates,” Comput. Phys. Commun., vol. 192, pp. 322–329, 2015.
 M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, “PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection,” JCAP, vol. 1103, p. 051, 2011. [Erratum: JCAP1210,E01(2012)].
 B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, “Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code,” JCAP, vol. 1302, p. 001, 2013.
 Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, and et al., “Planck 2015 results. XIII. Cosmological parameters,” Astron. Astrophys., vol. 594, p. A13, 2016.
 B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev., vol. D81, p. 104019, 2010.
 J. H. MacGibbon and B. R. Webber, “Quark and gluonjet emission from primordial black holes: The instantaneous spectra,” Phys. Rev. D, vol. 41, pp. 3052–3079, May 1990.
 D. N. Page, “Particle emission rates from a black hole: Massless particles from an uncharged, nonrotating hole,” Phys. Rev. D, vol. 13, pp. 198–206, Jan 1976.
 J. H. MacGibbon, “Quark and gluonjet emission from primordial black holes. II. The emission over the blackhole lifetime,” Phys. Rev. D, vol. 44, pp. 376–392, Jul 1991.
 B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev., vol. D81, p. 104019, 2010.
 R. Diamanti, L. LopezHonorez, O. Mena, S. PalomaresRuiz, and A. C. Vincent, “Constraining Dark Matter LateTime Energy Injection: Decays and PWave Annihilations,” JCAP, vol. 1402, p. 017, 2014.
 J. Shull and M. van Steenberg, “The ionization equilibrium of astrophysically abundant elements,” Astrophys.J.Suppl., vol. 48, pp. 95–107, 1982.