Cosmological Constraints on Unstable Particles: Numerical Bounds and Analytic Approximations

Many extensions of the Standard Model predict large numbers of additional unstable particles whose decays in the early universe are tightly constrained by observational data. For example, the decays of such particles can alter the ratios of light-element abundances, give rise to distortions in the cosmic microwave background, alter the ionization history of the universe, and contribute to the diffuse photon flux. Constraints on new physics from such considerations are typically derived for a single unstable particle species with a single well-defined mass and characteristic lifetime. In this paper, by contrast, we investigate the cosmological constraints on theories involving entire ensembles of decaying particles --- ensembles which span potentially broad ranges of masses and lifetimes. In addition to providing a detailed numerical analysis of these constraints, we also formulate a set of simple analytic approximations for these constraints which may be applied to generic ensembles of unstable particles which decay into electromagnetically-interacting final states. We then illustrate how these analytic approximations can be used to constrain a variety of toy scenarios for physics beyond the Standard Model. For ease of reference, we also compile our results in the form of a table which can be consulted independently of the rest of the paper. It is thus our hope that this work might serve as a useful reference for future model-builders concerned with cosmological constraints on decaying particles, regardless of the particular model under study.

Oct 2018 of observable consequences -especially if the final states into which these particles decay involve visible-sector particles. Indeed, electromagnetic or hadronic showers precipitated by unstable-particle decays within the recent cosmological past can alter the primordial abundances of light nuclei both during and after Big-Bang nucleosynthesis (BBN) [1][2][3][4], give rise to spectral distortions in the cosmic microwave background (CMB) [5,6], alter the ionization history of the universe [7][8][9][10], and give rise to characteristic features in the diffuse photon background. These considerations therefore place stringent constraints on models for new physics involving unstable particles.
Much previous work has focused on examining the cosmological consequences of a single particle species decaying in isolation, and the corresponding limits on the properties of such a particle species are now well established. Indeed, simple analytic approximations can be derived which accurately model the effects that the decays of such a particle can have on many of the relevant observables [2,5]. However, many theories for new physics involve not merely one or a few unstable particles, but rather a large -and potentially vast -number of such particles with a broad spectrum of masses, lifetimes, and cosmological abundances. For example, theories with additional spacetime dimensions give rise to infinite towers of Kaluza-Klein (KK) excitations for any field which propagates in the higher-dimensional bulk. Likewise, string theories generally predict large numbers of light moduli [11][12][13] or axion-like particles [14][15][16]. Collections of similar light fields also arise in supergravity [17]. There even exist approaches to the dark-matter problem such as the Dynamical Dark Matter (DDM) framework [18,19] which posit the existence of potentially vast ensembles of unstable dark-sector particles. It is therefore crucial to understand the cosmological consequences of entire ensembles of decaying particles in the early universe and, if possible, to formulate a corresponding set of analytic approximations which model the effects of these decays.
For a variety of reasons, assessing the effects of an entire ensemble of decaying particles with a broad range of masses, lifetimes, and cosmological abundances is not merely a matter of trivially generalizing the results obtained in the single-particle case. The decay of a given unstable particle amounts to an injection of additional electromagnetic radiation and/or other energetic particles into the evolution of the universe, and injection at different characteristic timescales during this evolution can have markedly different effects on the same observable. Moreover, since many of these observables evolve in time according to a complicated system of coupled equations, the effects of injection at any particular time t inj depend in a non-trivial way on the entire injection history prior to t inj through feedback effects.
In principle, the cosmological constraints on ensembles of decaying particles in the early universe can be evaluated numerically. Indeed, there are several publicly available codes [10,20,21] which can readily be modified in order to assess the effects of an arbitrary additional injection history on the relevant observables. Computational methods can certainly yield useful results in any particular individual case. However, another complementary approach which can provide additional physical insight into the underlying dynamics involves the formulation of approximate analytic expressions for the relevant observables -expressions analogous to those which already exist for a single particle species decaying in isolation. Our aim in this paper is to derive such a set of analytic expressions -expressions which are applicable to generic theories involving large numbers of unstable particles, but which nevertheless provide accurate approximations for the relevant cosmological observables. Thus, our results can serve as a useful reference for future model-builders concerned with cosmological constraints on decaying particles, regardless of the particular model under study.
In this paper, we shall focus primarily on the case in which electromagnetic injection dominates -i.e., the case in which the energy liberated by the decays of these particles is released primarily in the form of photons, electrons, and positrons rather than hadrons. We also emphasize that the approximations we shall derive in this paper are not ad-hoc in nature; in particular, they are not the results of empirical fits. Rather, as we shall see, they emerge organically from the underlying physics and thus carry direct information about the underlying processes involved.
This paper is organized as follows. In Sect. II, we begin by establishing the notation and conventions that we shall use throughout this paper. We also review the various scattering processes through which energetic photons injected by particle decay interact with other particles present in the radiation bath. In subsequent sections, we then turn our attention to entire ensembles of unstable particles, focusing on the electromagnetic injections arising from decays occurring after the BBN epoch. Each section is devoted to a different cosmological consideration arising from such injection, and in each case we ultimately obtain a simple, analytic approximation for the corresponding constraint. For example, in Sect. III we consider the constraints associated with the modification of the abundances of light nuclei after BBN, and in Sect. IV we consider limits on distortions of the CMB-photon spectrum. Likewise, in Sect. V we consider the constraints associated with the ionization history of universe and its impact on the CMB, and in Sect. VI we consider the constraints associated with additional contributions from unstable-particle decays to the diffuse photon background. Ultimately, the results from these sections furnish us with the tools needed to constrain decaying ensembles of various types. This is then illustrated in Sect. VII, where we consider how our results may be applied to two classes of ensembles whose constituents exhibit different representative mass spectra. Finally, in Sect. VIII, we conclude with a discussion of our main results and avenues for future work. For future reference, we also provide (in Table IV) a summary/compilation of our main results.

II. ELECTROMAGNETIC INJECTION: OVERVIEW AND CLASSIFICATION OF RELEVANT PROCESSES
Our aim in this paper is to assess the cosmological constraints on an ensemble consisting of a potentially large number N of unstable particle species χ i with masses m i and decay widths Γ i (or, equivalently, lifetimes τ i ≡ Γ −1 i ), where the index i = 0, 1, . . . , N − 1 labels these particle species in order of increasing mass m i . We shall characterize the cosmological abundance of each of the χ i in terms of a quantity Ω i which we call the "extrapolated abundance." This quantity represents the abundance that the species χ i would have had at present time, had it been absolutely stable. We shall assume that the total abundance Ω tot ≡ i Ω i of the ensemble is sufficiently small that the universe remains radiation-dominated until the time of matter-radiation equality t MRE ≈ 10 12 s. Moreover, we shall focus on the regime in which m i 1 GeV for all χ i and all of the ensemble constituents are non-relativistic by end of the BBN epoch. Within this regime, as we shall discuss in further detail below, the spectrum of energetic photons produced by electromagnetic injection takes a characteristic form which to a very good approximation depends only on the overall energy density injected [1]. By contrast, for much lighter decaying particles, the form of the resulting photon spectrum can differ from this characteristic form as a result of the immediate decay products lacking sufficient energy to induce the production of e + e − pairs by scattering off background photons [22,23]. The cosmological constraints on a single electromagnetically-decaying particle species with a mass below 1 GeV have recently been investigated in Ref. [24].
The considerations which place the most stringent constraints on the ensemble depend on the values of Ω i and τ i for the individual ensemble constituents. For ensembles of particles with lifetimes in the range 1 s τ i t now , where t now denotes the present age of the universe, the dominant constraints are those related to the abundances of light elements, to spectral distortions of the CMB, and to the ionization history of the universe. The effect of electromagnetic injection on the corresponding observables is sensitive to the overall energy density injected and to the timescales over which that energy is injected, but not to the details of the decay kinematics or the particular channels through which the χ i decay. Thus, in order to retain as much generality as possible in our analysis, we shall focus on ensembles for which all constituents with non-negligible Ω i have lifetimes within this range; moreover, we shall refrain from specifying any particular decay channel for the χ i when assessing the bounds on these ensembles due to these considerations. By contrast, the constraints on decaying ensembles that follow from limits on features in the diffuse photon background do depend on the particulars of the decay kinematics. Thus, when analyzing these constraints in Sect. VI, we not only present a general expression for the relevant observable -namely the contribution to the differential photon flux from the decaying ensemblebut also apply this result to a concrete example involving a particular decay topology.
Many of the constraints on electromagnetic injection are insensitive to the details of the decay kinematics because the injection of photons and other electromagnetically-interacting particles prior to CMB decoupling sets into motion a complicated chain of interactions which serve to redistribute the energies of these particles. In particular, the effects of electromagnetic injection on cosmological observables ultimately depend on the interplay between three broad classes of processes through which these photons interact with other particles in the background plasma. These are: • Class I: Cascade and cooling processes which rapidly redistribute the energy of the injected photons. Processes in this class include γ + γ BG → e + + e − and γ + γ BG → γ + γ, where γ BG denotes a background photon, as well as inverse-Compton scattering and e + e − pair production off nuclei. These processes occur on timescales far shorter than the timescales associated with other relevant processes, and thus may be considered to be effectively instantaneous. As we shall see, these processes serve to establish a non-thermal population of photons with a characteristic spectrum.
• Class II: Processes through which the non-thermal population of photons established by Class-I processes can have a direct effect on cosmological observables. These include the photoproduction and photodisintegration of light elements during or after BBN, as well as the photoionization of neutral hydrogen and helium after recombination.
• Class III: Processes which serve to bring the non-thermal population of photons established by Class-I processes into kinetic and/or thermal equilibrium with the radiation bath. Processes in this class include Compton scattering, bremsstrahlung, and e + e − pair production off nuclei.
We emphasize that these classes are not necessarily mutually exclusive, and that certain processes play different roles during different cosmological epochs. Any energy injected in the form of photons prior to last scattering is rapidly redistributed to lower energies due to the Class-I processes discussed above. The result is a non-thermal contribution to the photon spectrum at high energies with a normalization that depends on the total injected power and a generic shape which is essentially independent of the shape of the initial injection spectrum directly produced by χ i decays. This "reprocessed" photon spectrum serves as a source of for Class-II processes -processes which include, for example, reactions that alter the abundances of light nuclei and scattering processes which contribute to the ionization of neutral hydrogen and helium after recombination. Since all information about the detailed shape of the initial injection spectrum from decays is effectively washed out by Class-I processes in establishing this reprocessed photon spectrum, the results of our analysis are largely independent of the kinematics of χ i decay. This is ultimately why many of our results -including those pertaining to the alteration of light abundances after BBN, distortions in the CMB, and the ionization history of the universe -are likewise largely insensitive to the decay kinematics of the χ i .
The timescale over which injected photons can cause these alterations is controlled by the Class-III processes. Prior to CMB decoupling, these processes serve to "degrade" the reprocessed photon spectrum established by Class-I processes by bringing this non-thermal population of photons into kinetic or thermal equilibrium with the photons in the radiation bath. As this occurs, these Class-III interactions reduce the energies of the photons below the threshold for Class-II processes while also potentially altering the the shape of the CMB-photon spectrum. These Class-III processes eventually freeze out as well, after which point any photons injected by particle decays simply contribute to the diffuse extragalactic photon background.

III. IMPACT ON LIGHT-ELEMENT ABUNDANCES
We begin our analysis of the cosmological constraints on ensembles of unstable particles by considering the effect that the decays of these particles have on the abundances of light nuclei generated during BBN. We shall assume that these decays occur after BBN has concluded, i.e., after initial abundances for these nuclei are already established. We begin by reviewing the properties of the non-thermal photon spectrum which is established by the rapid reprocessing of injected photons from these decays by Class-I processes. We then review the corresponding constraints on a single unstable particle species [2]constraints derived from a numerical analysis of the coupled system of Boltzmann equations which govern the evolution of these abundances. We then set the stage for our eventual analysis by deriving a set of analytic approximations for the above constraints and demonstrating that the results obtained from these approximations are in excellent agreement with the results of a full numerical computation within our regime of interest. Finally, we apply our analytic approximations in order to constrain scenarios involving an entire ensemble of multiple decaying particles exhibiting a range of masses and lifetimes.

A. Reprocessed Injection Spectrum
As discussed in Sect. II, the initial spectrum of photons injected at time t inj is redistributed effectively instantaneously by Class-I processes. A detailed treatment of the Boltzmann equations governing these processes in a radiation-dominated epoch can be found, e.g., in Ref. [1]. For injection at times t inj 10 12 s, the resulting reprocessed photon spectrum turns out to take a characteristic form which we may parametrize as follows: The quantity dρ(t inj )/dt inj appearing in this expression, which specifies the overall normalization of the contribution to the reprocessed photon spectrum, represents the energy density injected by particle decays during the infinitesimal time interval from t inj to t inj + dt inj . The function K(E, t inj ), on the other hand, specifies the shape of the spectrum as a function of the photon energy E. This function is normalized such that It can be shown that for any process that injects energy primarily through electromagnetic (rather than hadronic) channels, the function K(E, t inj ) takes the universal form [1,25,26] K(E, t inj ) = 3) where K 0 is an overall normalization constant and where E C and E X are energy scales associated with specific Class-I processes whose interplay determines the shape of the reprocessed photon spectrum. The normalization convention in Eq. (3.2) implies that K 0 is given by Physically, the energy scales E C and E X appearing in Eq. (3.3) can be understood as follows. The scale E C represents the energy above which the photon spectrum is effectively extinguished by pair-production process γ + γ BG → e + +e − , in conjunction with interactions between the resulting electron and positron and other particles in the thermal bath. The energy scale E X represents the threshold above which γ + γ BG → γ + γ is the dominant process through which photons lose energy. By contrast, below this energy threshold, the dominant processes are Compton scattering and e + e − pair-production off nuclei. Note that while the normalization of the reprocessed photon spectrum is set by dρ(t inj )/dt inj , the shape of this spectrum is entirely controlled by the temperature T inj at injection. This temperature behaves like T inj ∝ t −1/2 inj in a radiation-dominated epoch. This implies that as t inj increases, the value of E C also increases. This reflects the fact that the thermal bath is colder at later injection times, and thus an injected photon must be more energetic in order for the pair-production process γ + γ BG → e + + e − to be effective. Numerically, the values of E C and E X at a given injection time t inj are estimated to be [1,2] where m e is the electron mass and T inj is the temperature of the thermal bath at t inj . The reprocessed photon spectrum in Eq. (3.1) is the spectrum which effectively contributes to the photoproduction and/or photodisintegration of light elements after the BBN epoch. In order to illustrate how the shape of this spectrum depends on the injection time t inj , we plot in Fig. 1 the function K(E, t inj ) which determines the shape of this spectrum as a function of E for several different values of the injection time t inj . Since the ultraviolet cutoff E C in the photon spectrum increases with t inj , injection at later times can initiate photoproduction and photodisintegration reactions with higher energy thresholds. The dashed vertical line, which we include for reference, represents the lowest threshold energy associated with any such reaction which can have a significant impact on the primordial abundance of any light nucleus which is tightly constrained by observation. As discussed in Sect. III B, this reaction turns out to be the deuterium-photodisintegration reaction D + γ → n + p, which has a threshold energy of roughly 2.2 MeV. Thus, the portion of the photon spectrum which lies within the gray shaded region in Fig. 1 has no effect on the abundance of any relevant nucleus. Since E C lies below this threshold for t inj ≈ 10 4 s, electromagnetic injection between the end of BBN and this timescale has essentially no effect on the abundances of light nuclei.
One additional complication that we must take into account in assessing the effect of injection on the abundances of light nuclei is that the reprocessed photon spectrum established by Class-I processes immediately after injection at t inj is subsequently degraded by Class-III processes, which slowly act to bring this reprocessed spectrum into thermal equilibrium with the background photons in the radiation bath. The timescale δt th (E, t inj ) on which these processes act on a photon of energy E is roughly As discussed in the text, this function determines the shape of the reprocessed photon spectrum for an instantaneous injection of photons of energy E at time tinj. The dashed vertical line at E ≈ 2.2 MeV indicates the lowest threshold energy associated with any photoproduction or photodisintegration reaction which can significantly alter the abundance of any light nucleus whose primordial abundance is tightly constrained by observation.
cross-section for the relevant scattering processes, which include Compton scattering and e + e − pair production off nuclei. The cross-sections for all relevant individual contributing processes can be found in Ref. [1]. Note that while at low energies σ th (E) is well approximated by the Thomson cross-section, this approximation breaks down at higher energies as other processes become relevant.
Given these observations, the spectrum of the resulting non-thermal population of photons at time t not only represents the sum of all contributions from injection at all times t inj < t but also reflects the subsequent degradation of these contributions by the Class-III processes which serve to thermalize this population of photons with the radiation bath. This overall non-thermal photon spectrum takes the form where is the Green's function which solves the differential equation The population of non-thermal photons described by Eq. (3.7) serves as the source for the initial photoproduction and photodisintegration reactions ultimately responsible for the modification of light-element abundances after BBN. We shall therefore henceforth refer to photons in this population as "primary" photons.
In what follows, we will find it useful to employ what we shall call the "uniform-decay approximation." Specifically, we shall approximate the full exponential decay of each dark-sector species χ i as if the entire population of such particles throughout the universe were to decay precisely at the same time τ i . As we shall see, this will prove critical in allowing us to formulate our ultimate analytic approximations. We shall nevertheless find that the results of our approximations are generally in excellent agreement with the results of a full numerical analysis.
Within the uniform-decay approximation, the contribution to dρ/dt(t inj ) from the decay of a single unstable particle species χ takes the form of a Dirac δ-function: where ρ χ (t inj ) is the energy density of χ at time t inj and where χ is the fraction of the energy density released by χ decays which is transferred to photons. It therefore follows that in this approximation, the primary-photon spectrum in Eq. (3.7) reduces to

B. Light-Element Production/Destruction
Generally speaking, the overall rate of change of the number density n a of a nuclear species N a due to the injection of electromagnetic energy at late times is governed by a Boltzmann equation of the form where H is the Hubble parameter and where C (p) a and C (s) a are the collision terms associated with two different classes of scattering processes which contribute to this overall rate of change. We shall describe these individual collision terms in detail below. Since n a is affected by Hubble expansion, it is more convenient to work with the corresponding comoving number density Y a ≡ n a /n B , where n B denotes the total number density of baryons. The Boltzmann equation for Y a then takes the form (3.13) The collision term C (p) a represents the collective contribution from Class-II processes directly involving the population of primary photons described by Eq. (3.11). The principal processes which contribute to C (p) a are photoproduction processes of the form N b +γ → N a +N c and photodisintegration processes of the form N a + γ → N b + N c , where N b and N c are other nuclei in the thermal plasma. The collision term associated with these processes takes the form where the indices b and c run over the nuclei present in the plasma, where σ a (E) respectively denote the cross-sections for the corresponding photoproduction and photodisintegration processes discussed above, and where E (ac) b and E (bc) a are the respective energy thresholds for these processes. Expressions for these cross-sections and values for the corresponding energy thresholds can be found, e.g., in Ref. [2].
By contrast, C (s) a represents the collective contribution from additional, secondary processes which involve not the primary photons themselves but rather a non-thermal population of energetic nuclei produced by interactions involving those primary photons. In principle, these secondary processes include both reactions that produce nuclei of species N a and reactions which destroy them. In practice, however, because the non-thermal population of any given species N b generated by processes involving primary photons is comparatively small, the effect of secondary processes on the populations of most nuclear species is likewise small. As we shall discuss in more detail in Sect. III C, the only exception is 6 Li, which is not produced in any significant amount during the BBN epoch but which can potentially be produced by secondary processes initiated by photon injection at subsequent times. Since these processes involve the production rather than the destruction of 6 Li, we focus on the effect of secondary processes on nuclei which appear in the final state rather than the initial state in what follows.
The energetic nuclei which participate in secondary processes are the products of the same kinds of reactions which lead to the collision term in Eq. (3.14). Thus, the kinetic-energy spectrum d n b (E b , t)/dE b of the nonthermal population of a nuclear species N b produced in this manner is in large part determined by the energy spectrum of the primary photons. In calculating this spectrum, one must in principle account for the fact that a photon of energy E γ can give rise to a range of possible E b values due to the range of possible scattering angles between the three-momentum vectors of the incoming photon and the excited nucleus in the center-of-mass frame. However, it can be shown [27] that a reasonable approximation for the collision term C (s) a for 6 Li is nevertheless obtained by taking E b to be a one-to-one function of E γ of the form [28] where E (bd) c is the energy threshold for the primary process N c + γ → N b + N d . In this approximation, is the energy-loss rate for N b due to Coulomb scattering with particles in the thermal background plasma. The exponential factor β b (E γ , t) accounts for the collective effect of additional processes which act to reduce the number of nuclei of species N b . The lower limit of integration in Eq. (3.16) is given by is the inverse of the function defined in Eq. (3.15). In other words, E −1 (E b ) is the photon energy which corresponds to a kinetic energy E b for the excited nucleus.
In principle, the processes which contribute to β a (E γ , t) include both decay processes (in the case in which N b is unstable) and photodisintegration processes of the form N b +γ → N c +N d involving a primary photon. In practice, the photodisintegration rate due to these processes is much slower that the energy-loss rate due to Coulomb scattering for any species of interest. Moreover, as we shall see in Sect. III C, the only nuclear species whose non-thermal population has a significant effect on the 6 Li abundance are tritium (T) and the helium isotope 3 He. Because these two species are mirror nuclei, the secondary processes in which they participate affect the 6 Li abundance in the same way and have almost identical cross-sections and energy thresholds. Thus, in terms of their effect on the production of 6 Li, the populations of T and 3 He may effectively be treated together as if they were the population of a single nuclear species. Although tritium is unstable and decays via beta decay to 3 He with a lifetime of τ T ≈ 5.6 × 10 8 s, these decays have no impact on the combined population of T and 3 He. We may therefore safely approximate β b (E γ , t) ≈ 0 for this combined population of excited nuclei in what follows.
The most relevant processes through which an energetic nucleus N b in this non-thermal population can alter the abundance of another nuclear species N b are scattering processes of the form N b + N f → N a + X, in which an energetic nucleus N b from the non-thermal population generated by primary processes scatters with a background nucleus N f , resulting in the production of a nucleus of species N a and some other particle X (which could be either an additional nucleus or a photon). In principle, processes of the form N b + N a → N f + X can also act to reduce the abundance of N a . However, as discussed above, this reduction has a negligible impact on Y a for any nuclear species N a which already has a sizable comoving number density at the end of BBN. Thus, we focus here on production rather than destruction when assessing the impact of secondary processes on the primordial abundances of light nuclei.
With this simplification, the collision term associated with secondary production processes takes the form where v(E b ) is the (non-relativistic) relative velocity of nuclei N b and N f in the background frame, where is the crosssection for the scattering process N b +N f → N a +X with corresponding threshold energy E (aX) bf , and where E(E C ) is the cutoff in d n b (E b , t)/dE b produced by primary processes.

C. Constraints on Primordial Light-Element Abundances
The nuclear species whose primordial abundances are the most tightly constrained by observation -and which are therefore relevant for constraining the late decays of unstable particles -are D, 4 He, 6 Li, and 7 Li. The abundance of 3 He during the present cosmological epoch has also been constrained by observation [29,30]. However, uncertainties in the contribution to this abundance from stellar sources make it difficult to translate the results of these measurements into bounds on the primordial 3 He abundance [31,32]. The effect of these uncertainties can be mitigated in part if we consider the ratio (D + 3 He)/H rather than 3 He/H, as the former is expected to be largely unaffected by stellar processing [33,34]. In this paper, we focus our attention on D, 4 He, 6 Li, and 7 Li, as the relationship between the measured abundances of these nuclei and their corresponding primordial abundances is more transparent.
The observational constraints on the primordial abundances of these nuclei can be summarized as follows. Bounds on the primordial 4 He abundance are typically phrased in terms of the primordial helium mass fraction Y p ≡ (ρ4 He /ρ B ) p , where ρ4 He is the mass density of 4 He, where ρ B is the total mass density of baryonic matter, and where the subscript p signifies that it is only the primordial contribution to ρ4 He which is used in calculating Y p , with subsequent modifications to this quantity due to stellar synthesis, etc., ignored. The 2σ limits on Y p are [35] 0.2369 < Y p < 0.2529 . (3.18) The observational 2σ limits on the 7 Li abundance are [36] 1.0 × 10 −10 < 7 Li H p < 2.2 × 10 −10 , (3.19) where the symbols 7 Li and H denote the primordial number densities of the corresponding nuclear species.
Constraining the primordial abundance of D is complicated by a mild tension which currently exists between the observational results for D/H derived from measurements of the line spectra of low-metallicity gas clouds [37] and the results obtained from numerical analysis of the Boltzmann equations for BBN [38] with input from Planck data [39], which predict a slightly lower value for this ratio. We account for these tensions by choosing our central value and lower limit on D/H in accord with the central value and 2σ lower limit from numerical calculations, while at the same time adopting the 2σ observational upper limit as our own upper limit on this ratio. Thus, we take our bounds on the D abundance to be While this upper bound is identical to the corresponding constraint quoted in Ref. [2], this is a numerical accident resulting from a higher estimate of the 6 Li/ 7 Li ratio (due to the recent detection of additional 6 Li in low-metallicity stars) and a reduction in the upper bound on the 7 Li/H ratio.
Having assessed the observational constraints on D, 4 He, 6 Li, and 7 Li, we now turn to consider the effect that the late-time injection of electromagnetic radiation has on the abundance of each of these nuclear species relative to its initial abundance at the conclusion of the BBN epoch. Of all these species, 4 He is by far the most abundant. For this reason, reactions involving 4 He nuclei in the initial state play an outsize role in the production of other nuclear species. Moreover, since the abundances of all other such species in the thermal bath are far smaller than that of 4 He, reactions involving these other nuclei in the initial state have a negligible impact on the 4 He abundance. Photodisintegration processes initiated directly by primary photons are therefore the only processes which have an appreciable effect on the primordial abundance of 4 He. A number of individual such processes contribute to the overall photodisintegration rate of 4 He, all of which have threshold energies E thresh 20 MeV.
The primordial abundance of 7 Li, like that of 4 He, evolves in response to photon injection primarily as a result of photodisintegration processes initiated directly by primary photons. At early times, when the energy ceiling E C in Eq. (3.5) for the spectrum of these photons is relatively low, the process 7 Li + γ → 4 He + T, which has a threshold energy of only ∼ 2.5 MeV, dominates the photodisintegration rate. By contrast, at later times, additional processes with higher threshold energies, such as 7 Li + γ → 6 Li + n and 7 Li + γ → 4 He + 2n + p, become relevant.
While the reactions which have a significant impact on the 4 He and 7 Li abundances all serve to reduce these abundances, the reactions which have an impact on the D abundance include both processes which create deuterium nuclei and processes which destroy them. At early times, photodisintegration processes initiated by primary photons -and in particular the process D + γ → n + p, which has a threshold energy of only E thresh ≈ 2.2 MeV -dominate and serve to deplete the initial D abundance. At later times, however, additional processes with higher energy thresholds turn on and serve to counteract this initial depletion. The dominant such process is the photoproduction process 4 He + γ → D + n + p, which has a threshold energy of E thresh ≈ 25 MeV. Unlike 4 He, 7 Li, and D, the nucleus 6 Li is not generated to any significant degree by BBN. However, a population of 6 Li nuclei can be generated after BBN as a result of photon injection at subsequent times. The most relevant processes are 7 Li + γ → 6 Li + n and the secondary production processes 4 He + 3 He → 6 Li + p and 4 He + T → 6 Li + n, where T denotes a tritium nucleus. All of these processes have energy thresholds E thresh ≈ 7 MeV. The abundances of 3 He and T, which serve as reactants in these secondary processes, are smaller at the end of BBN than the abundance of 4 He by factors of O(10 4 ) and O(10 6 ), respectively (for reviews, see, e.g., Ref. [41]). At the same time, the non-thermal populations of 3 He and T generated via the photodisintegration of 4 He are much larger than the nonthermal population of 4 He, which is generated via the photodisintegration of other, far less abundant nuclei. Thus, to a very good approximation, the reactions which contribute to the secondary production of 6 Li involve an excited 3 He or T nucleus and a "background" 4 He nucleus in thermal equilibrium with the radiation bath.
In Table I, we provide a list of the relevant reactions which can serve to alter the abundances of light nuclei as a consequence of photon injection at late times, along with their corresponding energy thresholds. Expressions for the cross-sections for these processes are given in Ref. [2]. While there exist additional nuclear processes beyond those listed in Table I 6 Li, and 7 Li in scenarios involving photon injection at late times, along with the corresponding energy threshold E thresh for each process. We note that the excited T and 3 He nuclei which participate in the secondary production of 6 Li are generated primarily by the processes listed above which contribute to the destruction of 4 He.
to the collision terms in Eq. (3.13), these processes do not have a significant impact on the Y a of any relevant nucleus when the injected energy density is small and can therefore be neglected. In should be noted that the population of excited T and 3 He nuclei which participate in the secondary production of 6 Li are generated primarily by the same processes which contribute to the destruction of 4 He. We note that we have not included processes which contribute to the destruction of 6 Li. The reason is that the collision terms for these processes are proportional to Y6 Li itself and thus only become important in the regime in which the rate of electromagnetic injection from unstable-particle decays is large. By contrast, for reasons that shall be discussed in greater detail below, we focus in what follows primarily on the regime in which injection is small and 6 Li-destruction processes are unimportant. However, we note that these processes can have an important effect on Y6 Li in the opposite regime, rendering the bound in Eq. (3.21) essentially unconstraining for sufficiently large injection rates [2].
Finally, we note that the rates and energy thresholds for 4 He + 3 He → 6 Li + p and 4 He + T → 6 Li + n are very similar, as are the rates and energy thresholds for the 4 He-destruction processes which produce the nonthermal populations of 3 He and T [2]. In what follows, we shall make the simplifying approximation that 3 He and T are "interchangeable" in the sense that we treat these rates -and hence also the non-thermal spectra d n a (E, t)/dE of 3 He and T -as identical. Thus, although T decays via beta decay to 3 He on a timescale τ T ≈ 10 8 s, we neglect the effect of the decay kinematics on the resulting non-thermal 3 He spectrum. As we shall see, these simplifying approximations do not significantly impact our results.

D. Towards an Analytic Approximation: Linearization and Decoupling
In order to assess whether a particular injection history is consistent with the constraints discussed in the previous section, we must evaluate the overall change δY a (t) ≡ Y a (t)−Y init a in the comoving number density Y a of a given nucleus at time t = t now , where Y init a denotes the initial value of Y a at the end of BBN. In principle, this involves solving a system of coupled differential equations, one for each nuclear species present in the thermal bath, each of the form given in Eq. (3.13).
In practice, however, we can obtain reasonably reliable estimates for the δY a without having to resort to a full numerical analysis. This is possible ultimately because observational constraints require |δY a | to be quite small for all relevant nuclei, as we saw in Sect. III C. The equations governing the evolution of the Y a are coupled due to feedback effects in which a change in the comoving number density of one nuclear species N a alters the reaction rates associated with the production of other nuclear species. However, if the change in Y a is sufficiently small for all relevant species, these feedback effects can be neglected and the evolution equations effectively decouple.
In order for the evolution equations for a particular nuclear species N a to decouple, the linearity criterion |δY b | Y b must be satisfied for any other nuclear species N b = N a which serves as a source for reactions that significantly affect the abundance of N a at all times after the conclusion of the BBN epoch. In principle, there are two ways in which this criterion could be enforced by the observational constraints and consistency conditions discussed in Sect. III C. The first is simply that the applicable bound on each Y b which serves as a source for N a is sufficiently stringent that this bound is always violated before the linearity criterion |δY b | Y b fails. The second possibility is that while the direct bound on Y b may not in and of itself require that |δY b | be small, the comoving number densities Y a and Y b are nevertheless directly related in such a way that the applicable bound on Y a is always violated before the linearity criterion |δY b | Y b fails. If one of these two conditions is satisfied for every species N b which serves as a source for N a , we may treat the evolution equation for N a as effectively decoupled from the equations which govern the evolution of all other nuclear species. We emphasize that N a itself need not satisfy the linearity criterion in order for its evolution equation to decouple in this way.
We now turn to examine whether and under what circumstances our criterion for the decoupling of the evolution equations is satisfied in practice for all relevant nuclear species. In Fig. 2 we illustrate the network of reactions which can have a significant effect on the values of Y a for these species. The nuclei which appear in the initial state of one of the primary production processes The nodes in the network represent different nuclear species: nodes represented by solid circles correspond to nuclei whose primordial comoving number densities Ya are reliably constrained by data, while nodes represented by dashed circles correspond to other nuclear species involved in these reactions. An arrow pointing from one node to another indicates that the nucleus associated with the node from which the arrow originates serves as a source for the nucleus associated with the node to which the arrow points. A solid arrow indicates that the corresponding reaction can potentially have a non-negligible impact on the abundance of the product nucleus, while a dashed arrow indicates that the effect of the corresponding reaction is always negligible. Note that 1 H nuclei -i.e., protons -are generated as a byproduct of many of the other reactions shown. However, the impact of these contributions to the 1 H abundance is negligible and the corresponding arrows have been omitted for clarity. A checked box superimposed on the arrows emerging from a node indicates that observational constraints enforce the linearity criterion |δYa| Ya for the corresponding nucleus. By contrast, an open box indicates that this criterion is not satisfied. The Boltzmann equation for a given nucleus Na effectively decouples, in the sense that feedback effects can be neglected in calculating Ya, whenever the linearity criterion is satisfied for all N b = Na which serve as a source for Na (though not necessarily for Na itself). Since no other nucleus serves as a source for 4 He or 7 Li, the Boltzmann equations for both of these species trivially decouple. The Boltzmann equation for D also decouples because 4 He, the one nucleus which serves as a significant source for D, is required to satisfy the linearity criterion. However, the Boltzmann equation for 6 Li does not decouple, as one of its source nuclei -in particular, 7 Li -does not satisfy the linearity criterion. Further details are described in the text. Table I are 4 He, which serves as a source for D and 6 Li, and 7 Li, which serves as a source for 6 Li. We note that while 3 He and T each appear in the initial state of one of the secondary production processes for 6 Li, it is the non-thermal population of each nucleus which plays a significant role in these reactions. Since the nonthermal populations of both 3 He and T are generated primarily as a byproduct of 4 He destruction, requiring that the linearity criterion be satisfied for 4 He and 7 Li is sufficient to ensure that the evolution equations for all relevant nuclear species decouple.

listed in
We begin by assessing whether the direct constraints on Y4 He and Y7 Li themselves are sufficient to enforce the linearity criterion. We take the initial values of these comoving number densities at the end of the BBN epoch to be those which correspond to the central observational values for Y p and 7 Li/H quoted in Ref. [42], namely Y p = 0.2449 and 7 Li/H = 1.6 × 10 −10 . Since neither 4 He nor 7 Li is produced at a significant rate by interactions involving other nuclear species, the evolution equation in Eq. (3.13) for each of these nuclei takes the form represents the rate at which Y b is depleted as a result of photodisintegration processes. This depletion rate varies in time, but depends neither on Y b nor on the comoving number density of any other nucleus. Since Y b decreases monotonically in this case, it follows that if this comoving quantity lies within the observationally-allowed range today, it must also lie within this range at all times since the end of the BBN epoch.
The bound on Y4 He which follows from Eq. (3.18) is sufficiently stringent that |δY4 He | Y4 He is indeed required at all times since the end of BBN for consistency with observation. Thus, our linearity criterion is always satisfied for 4 He. By contrast, the bound on Y7 Li in Eq. (3.19) is far weaker in the sense that |δY7 Li | need not necessarily be small in relation to Y7 Li itself. Moreover, while the contribution to δY6 Li from primary production is indeed directly related to δY7 Li , we find that the observational bound on Y6 Li is not always violated before the linearity criterion |δY7 Li | Y7 Li fails -even if we assume Y6 Li ≈ 0 at the end of BBN. The reason is that not every 7 Li nucleus destroyed by primary photodisintegration processes produces a 6 Li nucleus. Indeed, 7 Li + γ → 4 He + T and 7 Li + γ → 4 He + 2n + p contribute to the depletion of Y7 Li as well.
Since the energy threshold for 7 Li + γ → 4 He + T is lower than the threshold for 7 Li + γ → 6 Li + n, there will be a range of t inj within which injection contributes to the destruction of 7 Li without producing 6 Li at all. Furthermore, even at later injection times t inj 4.9 × 10 5 s, when the primary-photon spectrum from injection includes photons with energies above the threshold for 7 Li + γ → 6 Li + n, the 7 Li-photodisintegration rates associated with this process and the rates associated with 7 Li + γ → 4 He + T and 7 Li + γ → 4 He + 2n + p are comparable. Consequently, only around 30% of 7 Li nuclei destroyed by primary photodisintegration for t inj 4.9 × 10 5 s produce a 6 Li nucleus in the process. Thus, a large |δY7 Li | invariably results in a much smaller contribution to δY6 Li . Thus, 7 Li does not satisfy the linearity criterion when serving as a source for 6 Li.
That said, while our linearity criterion is not truly satisfied for 7 Li, we can nevertheless derive meaningful constraints on decaying particle ensembles from observational bounds on 6 Li by neglecting feedback effects on Y7 Li in calculating Y6 Li . Since the Boltzmann equation for Y7 Li takes the form given in Eq. (3.22), Y7 Li is always less than or equal to its initial value Y init 7 Li at the end of BBN. This in turn implies that the collision term C (p) 6 Li in the Boltzmann equation for 6 Li is always less than or equal to the value that it would have had if the linearity criterion for 7 Li had been satisfied. It therefore follows that the contribution to δY6 Li from primary production which we would obtain if we were to approximate Y7 Li by Y init 7 Li at all times subsequent to the end of BBN is always an overestimate. In this sense, then, the bound on electromagnetic injection which we would obtain by invoking this linear approximation for Y7 Li represents a conservative bound. Moreover, it turns out that because of the relationship between Y7 Li and Y6 Li , the bound on decaying ensembles from the destruction of 7 Li is always more stringent than the bound from the primary production of 6 Li. Thus, adopting the linear approximation for 7 Li in calculating C (p) 6 Li does not artificially exclude any region of parameter space for such ensembles once the combined constraints from all relevant nuclear species are taken into account.
Motivated by these considerations, in what follows we shall therefore adopt the linear approximation in which in calculating the collision terms C (p) a and C (s) a for any nuclear species N a = N b for which N b serves as a source. As we have seen, this approximation is valid for all species except for 7 Li, which serves as a source for primary 6 Li production. Moreover, adopting this approximation for 7 Li in calculating C (p) 6 Li yields a conservative bound on electromagnetic injection from decaying particle ensembles.
As discussed above, the advantage of working within the linear approximation is that the Boltzmann equations for all relevant N a effectively decouple and may be solved individually in order to yield analytic approximations for δY a . In the simplest case, in which the collision terms C where t 0 represents the time at the conclusion of the BBN epoch beyond which the initial abundance Y a (t 0 ) = Y init a generated by standard primordial nucleosynthesis remains essentially fixed in the absence of any subsequent injection. Moreover, even in cases in which C (p) a and C (s) a include both source and sink terms, we may still evaluate δY a in this way, provided that observational constraints restrict δY a to the region |δY a | Y a and therefore allow us to ignore feedback effects and approximate Y a ≈ Y init a as a constant on the right side of Eq. (3.13).
Since the Boltzmann equation for 6 Li contains no nonnegligible sink terms, and since observational constraints require that |δY a | Y a for both 4 He and D, it follows that δY a is well approximated by Eq. (3.23) for these species. Indeed, the only relevant nucleus which does not satisfy these criteria for direct integration is 7 Li. Nevertheless, since 7 Li is destroyed by a number of primary photodisintegration processes but not produced in any significant amount, the Boltzmann equation for this nucleus takes the particularly simple form specified in Eq. (3.22). This first-order differential equation can easily be solved for Y a , yielding an expression for the comoving number density at any time t ≥ t 0 : When this relation is expressed in terms of δY a rather than Y a , we find that it may be recast in the more revealing form In situations in which the linearity criterion |δY a (t)| Y init a is satisfied for the nucleus N a itself at all times t ≥ t 0 , Taylor expansion of the left side of this equation yields which is also the result obtained by direct integration of the Boltzmann equation for N a in the approximation that Y a ≈ Y init a . Comparing Eqs. (3.25) and (3.26), we see that if we were to neglect feedback and take Y a ≈ Y init a when evaluating δY a for a species for which this approximation is not particularly good, the naïve result that we would obtain for δY a would in fact correspond to the value of the quantity Y init Thus, given that a dictionary exists between the value of δY a obtained from Eq. (3.25) and the value obtained from Eq. (3.26), for simplicity in what follows we shall derive our analytic approximation for δY a using Eq. (3.26) and simply note that the appropriate substitution should be made for the case of 7 Li. That said, we also find that the constraint on Ω χ that we would derive from Eq. (3.26) in single-particle injection scenarios from the observational bound on Y7 Li differs from the constraint that we would derive from the more accurate approximation in Eq. (3.25) by only O(10%). Thus, results obtained by approximating δY7 Li by the expression in Eq. (3.26) are nevertheless fairly reliable in such scenarios -and indeed can be expected to be reasonably reliable in scenarios involving decaying ensembles as well.

E. Analytic Approximation: Contribution from Primary Processes
Having discussed how the Boltzmann equations for the relevant N a effectively decouple in the linear regime, we now proceed to derive a set of approximate analytic expressions for δY a from these decoupled equations. We begin by considering the contribution to δY a that arises from primary photoproduction or photodisintegration processes. The contribution from secondary processes, which is relevant only for 6 Li, will be discussed in Sect. III F.
Our ultimate goal is to derive an approximate analytic expression for the total contribution to δY a due the injection of photons from an entire ensemble of decaying states. However, our first step in this derivation shall be to consider the simpler case in which the injection is due to the decay of a single unstable particle species χ with a lifetime τ χ . We shall work within the uniformdecay approximation, in which the non-thermal photon spectrum takes the particularly simple form in Eq. (3.11). In this approximation, the lower limit of integration in Eq. (3.23) may be replaced by τ χ , while the upper limit can be taken to be any time well after photons at energies above the thresholds E (ac) b and E (bc) a for all relevant photoproduction and photodisintegration processes have thermalized. Thus, we may approximate the change in the comoving number density of each relevant nuclear species as In evaluating each δY a , we may also take advantage of the fact that the rates for the relevant reactions discussed in Sect. III C turn out to be such that the first (source) and second (sink) terms in Eq. (3.14) are never simultaneously large for any relevant nuclear species. Indeed, the closest thing to an exception occurs during a very small time interval within which the source and sink terms for D are both non-negligible. Thus, depending on the value of τ χ and its relationship to the timescales associated with these reactions, we may to a very good approximation treat the effect of injection from a single decaying particle as either producing or destroying N a .
With these approximations, the integral over t in Eq. (3.27) may be evaluated in closed form. In particular, when the source term in Eq. (3.14) dominates, we find that δY a is given by where we have defined .
By contrast, when the sink term dominates, we find that δY a is given by While the expressions for δY a in Eqs. (3.30) and (3.28) pertain to the case of a single unstable particle within the uniform-decay approximation, it is straightforward to generalize these results to more complicated scenarios. Indeed, within the linear approximation, the total change δY a in Y a which results from multiple instantaneous injections over an extended time interval is well approximated by the sum of the individual contributions from these injections. In the limit in which this set of discrete injections becomes a continuous spectrum, this sum becomes an integral over the injection time t inj . Thus, in the continuum limit, δY a is well approximated by where dY a /dt inj is the differential change in Y a due to an infinitesimal injection of energy in the form of photons at time t inj . The approximation in Eq. (3.31) allows us to account for the full exponential time-dependence of the electromagnetic injection due to particle decay in calculating δY a for any given nucleus. By extension, Eq. (3.31) gives us the ability to compare the results for δY a obtained both with and without invoking the uniform-decay approximation, thereby providing us insight into how reliably δY a can be computed with this approximation.
As an example, let us consider the effect of a single decaying particle χ with lifetime τ χ on the comoving number density of 4 He. Within the uniform-decay approximation, δY4 He is given by Eq. (3.30) because the sink term in Eq. (3.14) dominates. By contrast, when the full exponential nature of χ decay is taken into account, the corresponding result is where dρ χ (t inj )/dt inj denotes the rate of change in the energy density of χ per unit time t inj . Prior to t MRE , the energy density of an unstable particle with an extrapolated abundance Ω χ may be written is the critical density of the universe at present time and where t MRE is the time of matterradiation equality. The corresponding rate of change in the energy density, properly evaluated in the comoving frame and then transformed to the physical frame, is The corresponding expressions for continuum injection in cases in which the source term in Eq. (3.14) dominates are completely analogous and can be derived in a straightforward way.
In Fig. 3, we compare the results obtained for δY4 He within the uniform-decay approximation to the results obtained with the full exponential time-dependence of χ decay taken into account. In particular, within the (Ω χ , τ χ ) plane, we display contours of the corresponding 4 He mass fraction Y p obtained within the uniformdecay approximation (upper panel) and through the full exponential calculation (lower panel). For concreteness, in calculating these contours we have assumed an initial value Y p = 0.2449 for the 4 He mass fraction at the end of BBN, following Ref. [42].
Comparing the two panels of Fig. 3, we see that the results for Y p obtained within the uniform-decay approximation are indeed very similar throughout most of the parameter space shown to the results obtained through the full calculation. Nevertheless, we observe that discrepancies do arise. For example, we note that the constraints obtained within the uniform-decay approximation are slightly stronger for particles with lifetimes within the range 10 7 s τ χ 10 8 s and slightly weaker for τ χ above this range than the constraints obtained through the full calculation. These discrepancies are ultimately due to the fact that the reprocessed photon spectrum in Eq. (3.1) depends on the temperature of the radiation bath. Thus, there is a timescale t σ,max for which the reaction rate for a particular photoproduction or photodisintegration process is maximized for fixed dρ(t inj )/dt inj . Since all of the energy density ρ χ initially associated with the decaying particle is injected precisely at τ χ within the uniform-decay approximation, this approximation yields slightly more stringent constraints than those obtained through the full calculation when τ χ ∼ t σ,max . By the same token, the constraints obtained within the uniform-decay approximation are slightly less stringent than those obtained through the full calculation when τ χ differs significantly from t σ,max .
We also observe that for lifetimes τ χ 10 7 s, the uniform-decay approximation likewise yields constraints that are weaker than those obtained through the full exponential calculation. The reason for this is that the upper energy cutoff E C in the reprocessed photon spectrum is proportional to t 1/2 inj . For sufficiently early injection times, E C lies below the threshold energy E (bc) a for the photodisintegration reactions which contribute to δY4 He . Injection at such early times therefore has essentially no impact on Y p . This implies that within the uniform-decay approximation, a particle with a lifetime in this regime likewise has no effect on Y p . By contrast, within the full exponential calculation, injection occurs at a significant rate well after τ χ , leading to a non-negligible change in Y p even when the lifetime of the decaying particle is short.
In summary, the results shown in Fig. 3 indicate that other than in the regime where τ χ is sufficiently short that E C lies below the threshold energy for 4 He destruction, the result for δY4 He obtained within the uniform-decay approximation is very similar to the result obtained through the full exponential calculation for the same τ χ and Ω χ . We find this to be the case for the other relevant nuclei as well. Thus, having shown that the results for δY a obtained within the uniform-decay approximation accord well with those obtained through the full exponential calculation, at least for sufficiently long τ χ , we shall adopt this approximation in deriving our constraints on ensembles of electromagnetically-decaying particles. As we shall see, the advantage of working within the uniform-decay approximation is that within this approximation it is possible to write down a simple analytic expression for δY a . However, as we shall discuss in more detail below, we shall adopt an alternative strategy for approximating δY a within the regime in which τ χ is short and the results obtained within the uniform-decay approximation do not agree with those obtained through the full exponential calculation.
In order to write down our analytic expressions for δY a , we shall make one additional approximation: we shall treat the ratio of cross-sections S (bd) c (E) as not varying significantly as a function of energy between the threshold energy E (bc) a and E C . When this approximation holds, we may treat this ratio as a constant and pull it outside the integral over photon energies. In order to justify this approximation, we begin by noting that the cross-sections σ (bc) a (E) for primary processes typically peak at a value slightly above E (bc) a but fall precipitously with E beyond that point (see, e.g., Ref. [2]). In the vicinity of the peak, however, the variation of σ  a (E). Moreover, because K(E, τ χ ) falls rapidly with E, the dominant contribution to the energy integral in Eq. (3.30) arises from photons just above threshold even within the approximation that S (bc) a (E) is constant. Thus, to a good approximation, we may replace S (bd) c (E) by a constant on the order of its peak value and take this quantity outside the energy integral.
Within this approximation, it is now possible to analytically evaluate the integral in Eq. (3.28). The form of the result depends on the relationship between the threshold energy E (bc) a for the scattering process and the energy scales E C (t inj ) and E X (t inj ) which determine the shape of the photon spectrum at time t inj = τ χ . There are three cases of interest: the case in which E a . Moreover, the relations in Eq. (3.5) imply that each of these cases corresponds to a specific range of τ χ . In particular, the respective lifetime regimes are τ χ < t Ca , t (bc) (3.35) Within the τ χ < t (bc) Ca regime, none of the photons produced by χ decay exceed the threshold for the photodisintegration process. The contribution to δY a within the uniform-decay approximation is therefore formally zero. As discussed above, this is the regime in which the uniform-decay approximation fails to reproduce the results obtained through the full exponential calculation. Thus, in order to derive a meaningful bound on decaying particles with lifetimes in this regime, we instead model the injection of photons using the full continuum expression in Eq. (3.32) with dρ(t inj )/dt inj given by Eq. (3.34). However, in order to arrive at a simple analytic expression for δY a , we include only the contribution from injection times in the range t (bc) Ca beyond leading order in the resulting expression. With these approximations, in this regime we find where the proportionality constant A (bc) a for each contributing reaction is independent of the properties of the decaying particle. This treatment ensures that we obtain a more reliable estimate for the contribution to δY a from particles with τ χ in this regime.
Within the remaining two lifetime regimes, the contribution to δY a within the uniform-decay approximation is non-vanishing. Thus, within these regimes, we obtain our approximation for δY a by integrating Eq. (3.28), as discussed above. For t (bc) is a proportionality constant and where β ≡ E X /E C ≈ 0.27 is the τ χ -independent ratio of the energy scales in Eq. (3.5). Likewise, for t (bc) We emphasize that the proportionality constant B Strictly speaking, Eq. (3.38) does not hold for arbitrarily large τ χ , since photons produced by extremely late decays are not efficiently reprocessed by Class-I processes into the spectrum in Eq. (3.1). In order to account for this in what follows, we shall consider each term in the sum in Eq. (3.38) to be valid only for injection times t inj < t f a is a characteristic cutoff timescale associated with the reaction. Photons injected after this cutoff timescale are assumed to have no effect on Y a . For most reactions, it is appropriate to take t (bc) f a ≈ 10 12 s, as this is the timescale beyond which certain crucial Class-I processes effectively begin to shut off and the reprocessed photon spectrum is no longer reliably described by Eq. (3.1).

F. Analytic Approximation: Contribution from Secondary Processes
The approximate analytic expressions for δY a which we have derived in Sect. III E are applicable to all of the primary photoproduction or photodisintegration processes relevant for constraining electromagnetic injection from unstable-particle decays after BBN. However, since secondary production can contribute non-negligibly to the production of 6 Li, we must derive analogous expressions for δY a in the case of secondary production as well. Moreover, since secondary production is fundamentally different from primary production in terms of particle kinematics, there is no reason to expect that these expressions should have the same functional dependence on τ χ as that exhibited by the expressions in Eqs. (3.36)- (3.38). Indeed, as we shall see, they do not.
We begin this undertaking by observing that within the uniform-decay approximation, the contribution to δY a from secondary production is given by .

(3.39)
We may simplify this expression by noting that the timedependence of the energy-loss rate b(E b , t) of excited nuclei due to Coulomb scattering is primarily due to the dilution of the number density of electrons n e (t).
Since n e (t) scales with t in the same manner as n c (t), the quantity is, to a good approximation, a comoving quantity, and hence independent of t. Thus, Υ bcf (E) may be pulled outside the time integral in Eq. (3.39). Within this approximation, the contribution to δY a reduces to In order to proceed further, we must first assess the dependence of the quantities S bcf (E b ) on the respective energy scales E γ and E b . However, we find that Υ (aX) bcf (E b ) varies reasonably slowly over the relevant range of E b [2,28]. Thus, to a good approximation, this quantity may also be pulled outside the integral over E b .
By contrast, we find that S (ab) c (E γ ) cannot reliably be approximated as a constant over the range of E γ relevant for secondary production. This deserves further comment -especially because we have approximated S (ab) c (E γ ) as a constant in deriving the expressions in Eqs. (3.36)-(3.38) for primary production. As we shall now make clear, there are important differences between the kinematics of primary and secondary production which enable us to approximate S (ab) c (E γ ) as independent of E γ in the former case but not in the latter.
For primary production, as discussed in Sect. III E, the rapid decrease of K(E γ , τ χ ) with E γ suppresses the contribution to δY a from photons with E γ E (ac) b . The dominant contribution to δY a therefore comes from a narrow region of the spectrum just above threshold within which S have little collective impact on δY a . Thus, in approximating the overall contribution to δY a from primary production, it is reasonable to treat S (ac) b (E γ ) as a constant. By contrast, for secondary production -or at least for the secondary production of 6 Li, the one nuclear species in our analysis for which secondary production can have a significant impact on δY a -photons with E γ well above the threshold energy E (bd) c for any relevant primary process play a more important role. One reason for this is that the energy threshold E (aX) bc for each secondary processes which contributes meaningfully to 6 Li production corresponds to a primary-photon energy ] well above the associated primary-process threshold E (bd) c . Indeed, the kinetic-energy thresholds for 4 He + 3 He → 6 Li + p and 4 He + T → 6 Li + n given in Table I  (E γ ). Thus, it is really the energy threshold for the secondary process which sets the minimum value of E γ relevant for the secondary production of 6 Li. Since S (bd) c (E γ ) varies more rapidly with E γ at these energies than it does around its peak value, it follows that this variation cannot be neglected in determining the overall dependence of δY a on τ χ in this case.
There is, however, another reason why the variation of S (bd) c (E γ ) with E γ cannot be neglected in the case of secondary production -a deeper reason which is rooted more in fundamental differences between primary and secondary production than in the values of the particular energy thresholds associated with processes pertaining to 6 Li. The overall contribution to δY a from primary production in Eq. (3.28) involves a single integral over E γ . Thus, for primary production, the fall-off in K(E γ , τ χ ) itself with E γ is sufficient to suppress the partial contribution to δY a from photons with E γ E (bd) c . By contrast, the overall contribution to δY a from secondary production involves integration not only over E γ but also over E b . In this case, the fall-off in K(E γ , τ χ ) with E γ is not sufficient to suppress the partial contribution to δY a from photons with energies well above threshold. Thus, for any secondary production process, an accurate estimate for δY a can only be obtained when the variation of S (bd) c (E γ ) with E γ -a variation which can be quite significant at such energies -is taken into account.
We must therefore explicitly incorporate the functional dependence of S (bd) c (E γ ) on E γ into our calculation of δY a . We recall that S (bd) c (E γ ), as we have defined it in Eq. (3.29), represents the ratio of the cross-section σ (bd) c (E γ ) for the primary process N c + γ → N b + N d which produces the population of excited nuclei to the cross-section σ th (E γ ) for the Class-III processes which serve to thermalize the primary-photon spectrum. The cross-sections for the two primary processes 4 He + γ → 3 He + n and 4 He + γ → T + p relevant for secondary 6 Li production, expressed as a functions of the photon energy E γ , both take the form [2] σ (bd) c,0 ≈ 20.6 mb and σ (bd) c,0 ≈ 19.8 mb for these two processes, respectively. By contrast, σ th (E γ ) includes two individual contributions. The first contribution is due to e + e − pair-production off nuclei via Bethe-Heitler processes of the form γ + N a → X + e + + e − , where N a denotes a background nucleus and X denotes some hadronic final state. The cross-section for this process is where α is the fine-structure constant, where m e is the electron mass, and where σ T ≈ 661 mb is the Thomson cross-section. The second contribution to σ th (E γ ) is due to Compton scattering, for which the cross-section is Given the dependence of these cross-sections on E γ , we find that over the photon-energy range E (bd) c < E γ 1 GeV, the ratio S (bd) c (E γ ) for each relevant process is well approximated by a simple function of the form c,0 is a constant. In Fig. 4 we show a comparison between the exact value of S (bd) c (E γ ) for the individual process 4 He + γ → T + p and our approximation in Eq. (3.45) over this same range of E γ . We see from this figure that our approximation indeed provides a good fit to S (bd) c (E γ ) over most of this range. Moreover, while the discrepancy between these two functions becomes more pronounced as E γ ∼ O(1 GeV), we emphasize that S ] for the production of an excited T nucleus above the kineticenergy threshold for the secondary process 4 He+T → 6 Li+n. Thus, photons in the gray shaded region play no role in the secondary production of 6 Li. We note that the corresponding S (bd) c (Eγ) curves and thresholds for 4 He + γ → 3 He + n (the other primary process relevant for secondary 6 Li production) are nearly identical to those shown here. δY a , any discrepancy between Eq. (3.45) and the exact expression for S (bd) c at such energies may safely be ignored.
Armed with our approximation for S (bd) c (E γ ) in Eq. (3.45), it is now straightforward to evaluate the integral in Eq. (3.41) and obtain an approximate analytic expression for δY a . Just as it does for primary photoproduction and photodisintegration, the functional dependence of δY a on τ χ for secondary production depends on the relationship between τ χ and a pair of characteristic timescales determined by the energy thresholds for the relevant processes. As we have discussed above, it is ] rather than E (bd) c which sets the minimum E γ required of a photon in order for it to contribute to the secondary production of 6 Li. It therefore follows that the characteristic timescales for the secondary production of this nucleus are determined by . Thus, for secondary production, we define within the uniform-decay approximation, δY a is formally zero, as it is in the corresponding lifetime regime for primary production. Thus, in order to derive a meaningful constraint on unstable particles with lifetimes within this regime, we follow the same procedure as we employed in order to calculate the contribution to δY a from primary production for a particle with τ χ < t (bc) Ca . We model the injection of photons using the expression for secondary production appropriate for continuum injection, with dρ(t inj )/dt inj given by Eq. (3.34). In analogy to our treatment of the primary process in this regime, we include only the contribution from injection times in the range t Dropping terms beyond leading order in the dimensionless variable where the proportionality constant A (bcf ) a for each combination of primary and secondary processes is independent of τ χ . We note that the dependence of δY a on τ χ in this expression is exactly the same as in Eq. (3.36). By contrast, for t represents the ratio of the photon-energy threshold for primary production to the minimum photon energy needed to produce an excited nucleus of species N b with a kinetic energy above the energy threshold for the secondary process. Numerically, we find θ ≈ 0.46 for 4 He + γ → T + p followed by 4 He + T → 6 Li + n, while θ ≈ 0.52 for 4 He + γ → 3 He + n followed by 4 He + 3 He → 6 Li + p. Finally, for t (bcf ) Xa < τ χ , we find Thus, to summarize, we have derived a set of simple, analytic approximations for the change δY a in the abundance of a given nuclear species N a due to the late injection of photons by a decaying particle χ. For primary photoproduction or photodisintegration processes, we find that δY a is given by Eq. (3.36), Eq. (3.37), or Eq. (3.38), depending on the lifetime τ χ of the particle. Likewise, for secondary production, we find that δY a is given by Eq.

G. The Fruits of Linearization: Light-Element Constraints on Ensembles of Unstable Particles
We now turn to the task of extending these results to the case of an ensemble of decaying particles χ i with lifetimes τ i and extrapolated abundances Ω i . Indeed, we have seen that if the linearity criterion is satisfied both for N a itself and for all of its source nuclei N b = N a , all feedback effects on Y a can be neglected. Thus, in this regime, the overall change δY a in the abundance of a light nucleus is well approximated by the sum of the individual contributions associated with the individual χ i from each pertinent process. Indeed, it is only because we have entered a linear regime that such a direct sum is now appropriate. These, then, are the fruits of linearization.
While certain nuclei in our analysis -namely 6 Li and 7 Li -do not have the property that the linearity criterion is always satisfied for both the nucleus itself and its source nuclei, we emphasize that we are nevertheless able to derive meaningful bounds on the comoving number densities of these nuclei. As noted in Sect. III D, artificially adopting the linearity criterion for Y7 Li in the Boltzmann equation for 6 Li always yields a conservative bound on δY6 Li , regardless of the injection history. For 7 Li, the issue is that feedback effects on the photodisintegration rate due to changes in the comoving number density of 7 Li itself are not necessarily small. Thus, strictly speaking, δY7 Li is not well approximated by a direct sum of the contributions from the individual χ i . However, since this direct sum is simply an approximation of the integral in Eq. (3.25), it is equivalent to the quantity Y init 7 Li ln[1+δY7 Li (t)/Y init 7 Li ]. Thus, for all relevant N a , we find that an ensemble of decaying particles makes an overall contribution to δY a -or, in the case of 7 Li, to the quantity Y init a ln[1 + δY a (t)/Y init a ] -which can be approximated as a direct sum of the individual contributions from the individual ensemble constituents. Indeed, this yields either a reliable estimate for the true value of δY a or a reliably conservative bound on δY a .
In principle, each of these individual contributions to a given δY a can involve a large number of reactions with different energy thresholds and scattering kinematics, each with its own distinct fit parameters A (bc) Ca , etc. In practice, however, the number of reactions which contribute significantly to δY a for any one of the four relevant nuclear species is quite small, as can be seen from Table I. Moreover, it is often the case that many if not all of the reactions which have a nonnegligible impact on a given δY a have very similar energy thresholds and scattering kinematics. When this is the case, the contribution to δY a from these processes can, to a very good approximation, collectively be modeled using a single set of fit parameters. For example, the only reactions which have a significant impact on δY4 He are the primary photodisintegration processes 4 He + γ → T + p and 4 He + γ → 3 He + n, which have very similar energy thresholds and scattering kinematics. Thus, a fit involving a single set of parameters yields an accurate approximation for δY4 He . The reactions which contribute to δY7 Li differ more significantly in terms of their energy thresholds and scattering kinematics. Nevertheless, as we shall see, the collective effect of these processes is also well modeled by a single set of fit parameters. By contrast, for D there are two processes with qualitatively different energy thresholds and scattering kinematics which must be modeled using separate sets of fit parameters: primary photodisintegration via D + γ → n + p and primary photoproduction via 4 He + γ → D + n + p. Likewise, for 6 Li, two sets of fit parameters are required: one for primary photoproduction via 7 Li + γ → 6 Li + n and one for secondary production via both 4 He + 3 He → 6 Li + p and 4 He + T → 6 Li + n.
In general, then, we require at most two distinct sets of parameters in order to model the overall contribution to δY a for each relevant nucleus due to electromagnetic injection from an ensemble of decaying particles. Thus, in general, we may express this overall contribution as the sum of two terms (3.51) For each relevant nuclear species, one of these terms is associated with a primary process: primary photoproduction in the case of 6 Li and primary photodisintegration in the case of 4 He, 7 Li, and D. Thus, δY where A a , B a , t Ca , t Xa , etc., are model parameters whose assignments we shall discuss below. By contrast, the form of the second term in Eq. (3.51) differs depending on the nucleus in question. For 4 He and 7 Li, we simply have δY (2) a = 0. For D, this term is associated with primary production and thus takes exactly the same functional form as δY (1) a . In other words, we have where A * a , B * a , t * Ca , t * Xa , etc., represent an additional set of model parameters distinct from the parameters A a , B a , t Ca , t Xa , etc., in Eq. (3.52). For 6 Li, the term δY (2) a is associated with secondary production and therefore takes the form   We note that for simplicity and compactness of notation, we have implicitly taken i = 1 for each χ i in formulating the expressions appearing in Eqs. (3.52)-(3.54). However, it is straightforward to generalize these results to the case in which i differs from unity for one or more of the χ i . In particular, the corresponding expressions for δY (1) a and δY (2) a in this case may be obtained by replacing Ω i for each species with the product Ω i i . Moreover, we remind the reader that for the case of 7 Li, a more accurate estimate of δY a can be obtained by replacing δY . We now determine the values of the parameters in Eq. (3.52). We begin by noting that the parameter t Aa for each relevant nucleus N a may be assigned essentially any value between the end of the BBN epoch and the timescale at which the first reaction which contributes to δY a becomes efficient. Thus, for simplicity, we choose a universal value t Aa = 10 4 s for all N a , which corresponds to the timescale at which the first reaction in Table I namely the photodisintegration process D + γ → n + p -effectively turns on. For t f a , we likewise choose a universal value t f a = 10 12 s for all N a . As discussed in Sect. VI, the reason for this choice is that many of the Class-I processes which serve to establish the reprocessed photon spectrum in Eq. (3.1) start to become inefficient around this timescale. As a result, the photon spectrum that results from electromagnetic injection at times t inj 10 12 s differs from that in Eq. (3.1), and the analytic approximation for δY Given these results, we now determine the values of the remaining parameters in Eqs. (3.52)-(3.54) by demanding that the constraint contour we obtain from each of these equations be consistent with the contour obtained from a more complete numerical analysis of the Boltzmann system in the limiting case in which our decaying "ensemble" comprises only a single particle species χ. In this single-particle case, our analytic expression for δY a becomes a function of only two variables: the abundance Ω χ and lifetime τ χ of χ. We determine the constraint contours corresponding to the observational limits quoted in Sect. III C by surveying (Ω χ ,τ χ ) space and numerically solving the full, coupled system of Boltzmann equations for δY a at each point. The undetermined parameters in Eqs. (3.52)-(3.54) are then chosen such that our analytic expression provides a good fit to the corresponding constraint contour for each nucleus. We list the values for all parameters obtained in this manner in Table II. In deriving these constraint contours, it is necessary to translate the observational constraints in Eqs. on the corresponding subsequent change δY a in each Y a . In so doing, we must specify a set of initial values of Y a at the end of BBN. For Y4 He and Y7 Li , we take Y p = 0.2449 and 7 Li/H = 1.6 × 10 −10 , as discussed in Sect. III D. For Y D , we take the value which corresponds to the central theoretical prediction D/H = 2.413 × 10 −5 for the primordial D/H ratio derived in Ref. [37]. Finally, since 6 Li is not produced in any significant amount during BBN, we take Y6 Li = 0. The values for the observational upper and lower limits δY min a and δY max a on δY a which correspond to these choices for the initial Y a are listed in Table II alongside the values for our fit parameters.
In Fig. 5, we display contours showing the results of our numerical calculation of δY a in (Ω χ , τ χ ) space for each relevant nucleus. More specifically, we display contours of Y p for 4 He, whereas we display contours of log 10 Y a for D, 7 Li, and 6 Li. Moreover, for purposes of illustration, we display separate contours for D production and destruction -i.e., contours evaluated considering either production or destruction processes in C (p) D only, with the cross-sections for the opposite set of processes artificially set to zero. Likewise, we also display separate contours for primary and secondary 6 Li production. We emphasize, however, that in assessing our overall constraints on decaying ensembles, we consider these processes together, allowing for the possibility of a cancellation between In all panels except that pertaining to 4 He destruction, we plot contours of log 10 Ya, while for 4 He destruction we plot contours of Yp. Note that for later reference, we have plotted separate contours for D production and destruction; likewise, we show separate contours for primary and secondary 6 Li production. Regions of the (τχ, Ωχ) plane which are found by numerical analysis to be excluded are those above and/or to the right of the gray dashed line in each panel, while the solid colored curve in each panel demarcates the analogous region found to be excluded by fitting the free parameters in our analytic approximations for δYa in Eqs. (3.52)-(3.54) to our numerical results. We see that the full numerical results and our analytic approximations agree very well in all cases except that of primary 6 Li production, where our analytic approximation diverges from the numerical result at large Ωχ and thus provides an even more conservative bound on the allowed parameter space.
the two individual contributions to δY a (in the case of D) or a reinforcement (in the case of 6 Li). The gray dashed line in each panel of Fig. 5 indicates the upper bound on Ω χ obtained through a numerical calculation of δY a which follows from the observational limits in Eqs. (3.18)-(3.21). By contrast, the solid curve in each panel represents the corresponding constraint contour obtained through our analytic approximation for δY a in Eq. (3.51). We see from Fig. 5 that in all cases except for that of primary 6 Li production, the constraint contours we obtain using the linear and uniform-decay approximations are quite similar to those obtained in the more realistic case in which we account for the fact that χ has a constant decay rate and in which we account for feedback effects in the collision terms C for each relevant nucleus. Moreover, as discussed in Sect. III D, we see that for primary 6 Li production the constraint contour obtained using these approximations indeed represents a conservative bound on Ω χ . We also note that with this sole exception, the results of our analytic approximations are in good agreement with the corresponding results in Ref. [2] within the regime in which Ω χ is small and the linear approximation holds. Once again, we note that the effects of additional processes which we do not incorporate into this analysis can have important effects on the Y a in the regime in which Ω χ -and therefore the overall magnitude of electromagnetic injection -is large. These include, for example, photodisintegration processes which contribute to the depletion of a previously generated population of 6 Li as well as additional processes which contribute to the production or destruction of D [2]. However, since the regions of (Ω χ , τ χ ) space in which these effects on D and 6 Li are relevant are excluded by constraints on the primordial abundances of other nuclei, our neglecting these processes will have essentially no impact on the overall constraints we derive on ensembles of decaying particles.
In Fig. 6 we plot our analytic approximations to the constraint contours obtained for each relevant nucleus together in (Ω χ , τ χ ) space. Note that the 6 Li contour includes both primary and secondary contributions to δY6 Li . In determining the D contour, we have included contributions from both photoproduction and photodisintegration processes, the interference between which is responsible for the "funnel" region at τ χ ≈ 10 6.5 s within which the bound on Ω χ is considerably weakened. However, since the upper and lower bounds on δY D derive from different considerations, different colors are used to distinguish the part of the contour which corresponds to the upper bound (dark blue) from the part which corresponds to the lower bound (light blue). We note that in all cases, these contours are in good agreement with the corresponding results in Ref. [2] within the  7 Li destruction (green), 4 He destruction (magenta), and 6 Li production (red), the latter including both primary and secondary contributions. In addition, we also plot constraints stemming from limits on potential distortions of the CMB-photon spectrum, both µ-type (orange) and ytype (black), as will be discussed in Sect. IV. In each case, the regions of parameter space above each contour are excluded.
regime in which Ω χ is small and the linear approximation holds.

IV. SPECTRAL DISTORTIONS IN THE CMB
In addition to altering the abundances of light nuclei, electromagnetic injection from the decays of unstable particles in the early universe can also have observable effects on the spectrum of CMB photons, distorting that spectrum away from a pure blackbody profile. Observational limits on these distortions therefore also constrain such decays. Once again, while the constraints on a single decaying particle species are well known, the corresponding constraints on an ensemble of decaying particles are less well established. Our aim in this section is to derive approximate analytic expressions for the CMB distortions which can arise due to entire ensembles of decaying particles -expressions analogous to those we derived in Sect. III for the change δY a in the comoving number density in each relevant nucleus after BBN.
One of the subtleties which arises in constraining distortions in the CMB-photon spectrum from ensembles of decaying particles with a broad range of lifetimes is that both the manner in which that spectrum is distorted and the magnitude of that distortion depend on the timescale over which the injection takes place. At early times, energetic photons produced by these decays are rapidly brought into both kinetic and thermal equilibrium by the Class-III processes discussed in Sect. II. As a result, the spectrum of photons retains a blackbody shape. However, at later times, photon-numberchanging processes such as double-Compton scattering and bremsstrahlung "freeze out," in the sense that the rates for these processes drop below the expansion rate H of the universe. Once this occurs, injected photons are no longer able to attain thermal equilibrium with photons in the radiation bath. Nevertheless, they are able to attain kinetic equilibrium with these photons through elastic Compton scattering. As a result, the photon spectrum no longer retains its original blackbody shape; rather, it is distorted into a Bose-Einstein distribution characterized by a pseudo-degeneracy parameter µ.
Eventually, at an even later time t EC ≈ 9×10 9 s, elastic Compton scattering also freezes out. Thus, photons injected at times t t EC attain neither kinetic nor thermal equilibrium with photons in the radiation bath. Nevertheless, the injected photons continue to interact with electrons in the radiation bath at an appreciable rate until the time of last scattering t LS ≈ 1.19 × 10 13 s. During this window, photons in the radiation bath are up-scattered by electrons which acquire significant energy from these interactions. The ultimate result is a photon spectrum which is suppressed at low frequencies and enhanced at high frequencies. This spectral distortion is characterized by a non-zero Compton y parameter y C . Finally, after last scattering, any additional photons injected from particle decay simply contribute to the diffuse photon background.
In order to derive analytic expressions for the constraints on the parameters µ and y C in the presence of a decaying ensemble, we shall once again make use of the uniform-decay approximation as well as a linear approximation similar to the approximation we invoked in calculating the δY a in Sect. III. The validity of such an approximation follows from the fact that both µ and y C are tightly constrained by observation. The most stringent limits on these parameters, which are currently those from COBE/FIRAS data [43], are |µ| < 9.0 × 10 −5 y C < 1.2 × 10 −5 . (4.1) Moreover, proposed broad-spectrum experiments, such as PIXIE [44], could potentially increase the sensitivity to which these distortions can be probed by more than an order of magnitude. These stringent constraints on the shape of the CMB-photon spectrum imply that the contribution δρ γ to the photon energy density from injection during the relevant epoch must be quite small, in the sense that δρ γ ρ γ . Thus, to a very good approximation, we may ignore feedback effects on the overall photon number density n γ and energy density ρ γ when analyzing distortions in that spectrum. This allows us to replace these quantities by their equilibrium values.
This linear approximation significantly simplifies the analysis of the resulting CMB spectrum [5,6,45,46]. For example, in this approximation we may study the evolution of the CMB-photon spectrum using the method of Green's functions [20,47]. In this approach, we may express a generic spectral distortion which constitutes a shift ∆I ν (t now ) in the present-day intensity of photons with frequency ν relative to the expected intensity in the standard cosmology in the form where Q is the energy density injected at time t and where G(ν, t, t ) is a Green's function which relates the intensity of photons with frequency ν injected at time t to the photon intensity at that same frequency at a later time t > t. In general, G(ν, t, t now ) is well approximated by a sum of terms representing an overall temperature shift, a µ-type distortion, and y-type distortion [20]. This Green's-function formalism provides a tool for numerically estimating the overall distortion to the CMB-photon spectrum from any arbitrary injection history, provided that the overall injected energy is small. We shall draw on this formalism extensively in deriving our analytic approximations for δµ and δy C .

A. The µ Parameter
The time-evolution of the pseudo-degeneracy parameter µ is governed by an equation of the form [5,6] 3) The first term in Eq. (4.3) is a source term which arises due to the injection of energetic photons. Within the linear approximation, this source term takes the generic form where n γ and ρ γ respectively denote the overall number density and energy density of photons. We emphasize that this expression is completely general and applies even in absence of any source of photon injection, in which case the two terms in Eq. (4.4) cancel. The second term in Eq. (4.3) is a sink term associated with processes which serve to bring the population of photons in the radiation bath into thermal equilibrium and thereby erase µ. As discussed above, the dominant such processes are double-Compton scattering and bremsstrahlung, the respective thermalization rates for which are well approximated by power-law expressions of the form [5] where t MRE once again denotes the time of matterradiation equality. The coefficients in this expression are given by where h is the Hubble constant, Ω B is the present-day abundance of baryons, and T now ≈ 2.7 K is the presentday temperature of the CMB. Up to this point, our expressions governing the evolution of µ are completely general and may be applied to any injection history, provided that the total energy density injected is sufficiently small that the linear approximation is valid. In order to derive an approximate analytic expression for the overall change δµ in µ due to an ensemble of unstable particles, we proceed in essentially the same manner as we did in deriving expressions for the changes δY a in the comoving number densities of light nuclei in Sect. III. We begin by deriving an expression for δµ in the presence of a single decaying particle species χ with mass m χ and lifetime τ χ . When we take into account the full exponential nature of the decay, we find where N (γ) χ is the average number of photons produced per decay and where χ is the fraction of the energy released by χ decays which is transferred to photons. By contrast, within the uniform-decay approximation in which dρ χ /dt is given by Eq. (3.10), we have In this approximation, a simple, analytic expression for δµ may be obtained by directly integrating Eq. (4.3) from the time t e ≈ 1.69 × 10 3 s at which electron-positron annihilation freezes out to the time t EC at which elastic Compton scattering freezes out. In particular, we find that a particle species with a lifetime τ χ anywhere within this range gives rise to a µ-type distortion of size where m p is the proton mass, where Ω B is the total present-day abundance of baryons, where η ≡ n B /n γ ≈ 6.2 × 10 −10 is the baryon-to-photon ratio of the universe, and where Ω χ once again refers to the extrapolated abundance of χ. For τ χ outside this range, δµ ≈ 0. Since the bremsstrahlung term in the exponential in Eq. (4.9) is generically negligible in comparison with the double-Compton-scattering term, we may safely neglect it in what follows. Moreover, for m χ and τ χ within our regimes of interest, the second term in the prefactor in Eq. (4.9) is subleading and can likewise be neglected. Thus, we find that δµ is well approximated by  This latter quantity represents the timescale at which the damping due to double-Compton scattering becomes inefficient and contributions to δµ become effectively unsuppressed.
The dependence of δµ on τ χ has a straightforward physical interpretation. The exponential suppression factor reflects the washout of µ by double-Compton scattering at early times τ χ t µ . The additional polynomial factor reflects the fact that the equilibrium photon density grows more rapidly with time than does the matter density at earlier times. This implies that particles decaying at an earlier time yield smaller perturbations to the photon spectrum.

B. The Compton y Parameter
As discussed above, elastic Compton scattering serves to maintain kinetic equilibrium among photons in the radiation bath after double-Compton scattering and bremsstrahlung freeze out. However, elastic Compton scattering eventually freezes out as well, at time t ∼ t EC . An injection of photons which occurs after t EC but before the time t LS of last scattering contributes to the development of a non-vanishing Compton y parameter y C . Within the linear approximation, a reliable estimate for the rate of change of y C may be obtained from the relation [48] where dρ(t inj )/dt inj denotes the rate at which energy density is injected into the photon bath. We note that in contrast with the corresponding evolution equation for µ in Eq. (4.3), this equation contains no damping term. In deriving our analytic approximation for the change δy C from an ensemble of unstable particle species, we once again begin by considering the case of a single such species χ. When we take into account the full exponential nature of the decay, we obtain The corresponding result within the uniform-decay approximation is In integrating this latter expression over t in order to obtain our approximate expression for δy C , we must account for the fact that the epoch during which injection contributes to y C straddles the time t MRE of matter-radiation equality. Using the appropriate timetemperature relations 4.15) before and after the transition to matter-domination, where T MRE is the temperature at matter-radiation equality, we find (4.16) where η once again denotes the baryon-to-photon ratio of the universe.

C. The Fruits of Linearization: CMB-Distortion Constraints on Ensembles of Unstable Particles
Having derived analytic approximations for the spectral distortions δµ and δy C due to a single decaying particle species within the uniform-decay approximation, we now proceed to generalize these results to the case of an arbitrary injection history.
Within the linear approximation, the overall contribution to δµ from a decaying ensemble is simply the sum of the individual contributions from the constituent particles χ i with lifetimes in the range t e τ i t EC . Likewise, the overall contribution to y C is a sum of individual contributions from χ i with t EC τ i t LS . However, injection from constituents with lifetimes τ i ∼ t EC gives rise to intermediate-type distortions which are distinct from the purely µ-type or purely y-type distortions discussed above. While these intermediatetype distortions are not merely a superposition of µ-type and y-type distortions [46], isolating and constraining such intermediate-type distortions would require a more precise measurement of the CMB-photon spectrum than COBE/FIRAS data currently provide. We therefore find that -at least for the time being -the net effect of injection at times t ∼ t EC may indeed be modeled through such a superposition, with δµ and δy C given by modified versions of Eqs. (4.10) and (4.16). The appropriate modifications of these expressions can be inferred from the forms of the Green's functions in Ref. [20]. In particular, we model the overall contributions to δµ from the decaying ensemble by an expression of the form 1 − e −(t2µ/τi) 4αµ/3 .

(4.17)
Distortion t 0a (s) t Ba (s) t 1a (s) αa Aa Ba µ-type 3.3 × 10 6 2.4 × 10 9 1.2 × 10 10 1.0 × 10 0 5.4 × 10 −3 3.1 × 10 −1 yc-type 3.3 × 10 6 6.8 × 10 7 5.4 × 10 9 1.2 × 10 0 3.7 × 10 −5 2.7 × 10 −1  Likewise, for δy C we have In these equations we have introduced two sets of parameters A µ , B µ , α µ , t 0µ , etc., and A y , B y , α y , t 1y , etc., in order to characterize δµ and δy C , respectively. We note, however, that the timescales t 2µ ≡ t 1y are not independent model parameters but merely shorthand notation for quantities which are completely determined by t 1µ and t 1y , respectively. We also note that this parametrization accounts not only for µ-type distortions at times t t EC and y-type distortions at times t t EC , but also for subdominant contributions to δy C prior to t EC and to δµ after t EC . Finally, we note that in formulating Eqs. (4.17) and (4.18) we have once again implicitly taken i = 1 for all χ i , just as we did in Eqs. (3.52)-(3.54). The corresponding expressions for δµ and δy C may nevertheless be obtained for the case in which i differs from unity for one or more of the χ i by once again replacing Ω i for each species with the product Ω i i .
We now determine the values of the parameters in Eqs. (4.17) and (4.18) in a manner similar to that in which we determined the values of the parameters in Eqs. (3.52)-(3.54). Specifically, for each nucleus we demand that our analytic approximations yield constraint contours consistent with those obtained from numerical analysis in the case in which the ensemble comprises a single decaying particle with lifetime τ χ and abundance Ω χ . In particular, we obtain arrays of δµ and δy C values by surveying (Ω χ , τ χ ) space and computing these distortions numerically at each point using the Green'sfunction formalism of Refs. [20,47]. The undetermined parameters in Eqs. (4.17) and (4.18) are then chosen such that our analytic expressions for δµ and δy C provide good fits to the corresponding constraint contours. Our best-fit values for all of these parameters are quoted in Table III. In Fig. 7, we display contours showing the results of our numerical calculations for δµ (left panel) and δy C (right panel) within the (Ω χ , τ χ ) plane. The gray dashed line in each panel demarcates the upper boundary of the region within the (Ω χ , τ χ ) plane which is consistent with the limits on these distortions quoted in Eq. (4.1). By contrast, the solid curve in each panel represents the corresponding constraint contour obtained through our analytic approximation in Eqs. (4.17) and (4.18).
As we see, the results obtained through our analytic approximation represent a good fit to our numerical results for both δµ and δy C . The constraint contours obtained from our analytic approximations for both µ-and y-type distortions are also included in Fig. 6 in order to facilitate comparison with the constraint contours associated with the primordial abundances of light nuclei. We note that the constraint on δµ is more stringent than the constraint on y C for τ χ t EC , as expected, while the opposite is true for τ χ t EC . We also see that both of these constraints are subdominant in comparison with the constraints on the primordial abundances of light nuclei, except at late times τ χ 10 10 s.

V. MODIFICATIONS TO THE IONIZATION HISTORY OF THE UNIVERSE
As we have seen in Sect. IV, electromagnetic injection from late-decaying particles can give rise to certain distortions in the CMB-photon spectrum. However, such injections can affect the CMB in other ways as well. In particular, when such injection occurs during recombination, the resulting photons and electrons contribute to the heating and ionization of neutral hydrogen and helium as they cool, thereby expanding the surface of last scattering. This in turn leads to alterations in the pattern of CMB anisotropies, including a damping of correlations between temperature fluctuations as well as an enhancement of correlations between polarization fluctuations at low multipole moments [7]. These considerations therefore place additional constraints on decaying particle ensembles.
In a nutshell, these effects arise because the rates for many of the Class-I processes discussed in Sect. II are still sizable at times t ∼ t LS . These processes therefore continue to rapidly redistribute the energies of photons injected during the recombination epoch to lower energies, thereby allowing for efficient photoionization of neutral hydrogen and helium at a rate that exceeds the rate of cosmic expansion. However, there is only a somewhat narrow time window around t LS during which an injection of photons can have a significant impact on the surface of last scattering. At times t t LS , any additional ionization of particles in the plasma is effectively washed out. By contrast, at times t t LS , the relevant Class-I processes effectively shut off and photoionization is suppressed.
To a good approximation, then, we may regard any modification of the ionization history of the universe as being due to injection at t ∼ t LS . Thus, the crucial quantity which is constrained is the overall injection rate of energy density in the form of photons or other electromagnetically-interacting particles at t LS . Since the universe is already matter-dominated by the recombination epoch, the injection rate associated with a single decaying particle species χ with lifetime τ χ and extrapolated abundance Ω χ at t LS is given by 1) where χ once again denotes the fraction of the energy density released by χ decays which is transferred to photons. Thus, we model the constraint on such a decaying particle from its effects on the ionization history of the universe by an inequality of the form where t IH and Γ IH are taken to be free parameters. We next determine the values of Γ IH and t IH by fitting these parameters to the constraint contours obtained in Refs. [10,49] from a numerical analysis of the injection and deposition dynamics. The constraint contours obtained from subsequent numerical studies by other authors [50] are quite similar. In particular, we find Note that recent results from the final Planck data release [51] will affect the numerical values for Γ IH and t IH slightly but not dramatically. In Fig. 8, we compare the results of our analytic approximation with the numerical results of Refs. [10,49]. We observe that our approximation for the constraint contour is in good agreement with these numerical results for τ χ 10 15 s and deviates significantly from this contour only when τ χ ∼ t LS . However, as discussed in Ref. [10], the results to which we perform our fit yield an artificially conservative bound on Ω χ for lifetimes τ χ ∼ t LS , as the production of additional Lyman-α photons due to excitations from the ionizing particles was intentionally neglected. Other analyses of the ionization history [52] which include this effect yield constraint contours which more closely resemble that obtained from our analytic approximation.
The constraint in Eq. (5.2) may be generalized to an ensemble of decaying particles in a straightforward manner. Indeed, since this constraint is ultimately a bound on the overall injection rate at t ∼ t LS , the corresponding constraint on an ensemble of decaying particles χ i takes the form FIG. 8: Constraints on a decaying particle with lifetime τχ and extrapolated abundance Ωχ from modifications of the ionization history of the universe. These constraints correspond to the case in which χ = 1 and the entirety of the energy density released by χ decays is transferred to photons. The gray dashed contour represents the bound obtained in Refs. [10,49] from numerical analysis. The brown shaded region above the contour is excluded. The solid curve represents the constraint contour obtained by fitting the corresponding analytic approximation in Eq. (5.2) to these results. Dashed vertical lines indicating the time tLS of last scattering and the present age of the universe tnow have also been included for reference. We see that our approximation for the constraint contour is in good agreement with these numerical results for τχ 10 15 s and deviates significantly from this contour only when τχ ∼ tLS.
with values of Γ IH and t IH once again given in Eq. (5.3).
We emphasize that the constraint in Eq. (5.4) stems solely from the effect of electromagnetic injection from the decays of the ensemble constituents on the surface of last scattering. Since this effect arises primarily from injection at t ∼ t LS , the constraints are most relevant for ensembles in which the constituents have lifetimes roughly around this timescale. However, electromagnetic injection can affect the ionization history of the universe in other ways as well. For example, energy deposited in the intergalactic medium as a result of particle decays at timescales well after t LS serves both to elevate the gas temperature and to increase the ionization fraction of hydrogen. These effects can shift the onset of reionization to earlier times [49,53] and impact the 21-cm absorption/emission signal arising from hyperfine transitions in neutral hydrogen [54][55][56]. These considerations are relevant primarily for particles with lifetimes τ i t now -particles which would in principle contribute to the present-day dark-matter abundance. Since our primary focus in this paper is on ensembles whose constituents have lifetimes τ i t now , we do not investigate these additional considerations here. However, we note that they can play an important role in constraining ensembles of particles with longer lifetimes, such as those which arise within the context of the DDM framework.

VI. CONTRIBUTIONS TO THE DIFFUSE PHOTON BACKGROUND
At times t t LS , the universe becomes transparent to photons with a broad range of energies O(keV) E γ O(TeV) [7]. Photons within this energy range which are injected on timescales t t LS are not reprocessed by Class-I processes and do not interact with CMB photons at an appreciable rate. Thus, as discussed in Sect. IV, such photons therefore simply free-stream until the present epoch and contribute to the diffuse photon background. The total diffuse background of photons observed from the location of Earth includes both a galactic contribution from processes occurring within the Milky-Way halo (at redshifts z ≈ 0) and an extra-galactic contribution from processes which occurred within the more distant past. The decays of unstable particles with lifetimes τ i t now contribute essentially exclusively to the extra-galactic component of this total background. By contrast, particles with lifetimes τ i t nowparticles which would contribute non-negligibly to the present-day abundance of dark matter -in principle contribute to both components. Since our primary focus in this paper is on particles with lifetimes in the regime τ i t now , we focus on the extra-galactic component of the diffuse photon background. However, we also note that since the local dark-matter density within the Milky-Way halo is several orders of magnitude larger than the corresponding cosmological density, the galactic component of the diffuse photon background is typically the more important of the two for constraining particles with lifetimes τ i t now .
The overall diffuse extra-galactic contribution dΦ/dE γ to the differential flux of photons per unit energy observed by a detector on Earth from the decays of an ensemble of unstable constituent particles χ i is simply the direct sum of the individual contributions dΦ i /dE γ from these constituents. In order to evaluate these individual contributions, we begin by noting that the decays of a given species χ i , averaged over all possible decay channels, produce a characteristic injection spectrum of photons dN i /dE γ . The spectrum of photons arriving at present time at a detector on Earth is the integrated total contribution from decays occurring at different redshifts. A photon injected at redshift z with energy E γ is redshifted to energy E γ = (1+z) −1 E γ when it arrives at the location of Earth. Thus, we find that the observed flux of photons observed at a detector on Earth from the decay of an individual decaying species χ i may be written in the form where dΩ is the differential solid angle of the detector, where ρ i (z) is the energy density of χ i at redshift z, where κ(z) is the optical depth of the universe to photons emitted at redshift z, and where (z) is the line-of-sight distance away from Earth, expressed as a function of z.
We focus here primarily on photons with energies in the range O(keV) E γ O(TeV), for which the universe is effectively transparent at all times t t LS and opaque at times t t LS . Moreover, for simplicity, we shall approximate the optical depth for such photons with a function of the form where z LS ≡ z(t LS ) ≈ 1100 is the redshift at the time of last scattering. In other words, we shall assume that the universe is completely transparent to all photons within this energy range after last scattering and completely opaque beforehand. Moreover, since t LS t MRE , we approximate the line-of-sight distance (z) in a flat universe dominated by matter and dark energy as where c is the speed of light, where H now is the presentday value of the Hubble parameter, and where we have defined where Ω m and Ω Λ respectively denote the present-day abundances of massive matter and dark energy, the former including contributions from both dark and baryonic matter. The energy density ρ i (z) of a decaying particle can be expressed in terms of its extrapolated abundance Ω i as where the time t(z) which corresponds to redshift z is .
With these substitutions, Eq. (6.1) becomes Finally, in assessing our prediction for dΦ/dE γ for a decaying particle ensemble, we must account for the fact that at a realistic detector, the photon energies recorded by the instrument differ from the actual energies of the incoming photons due to the non-zero energy resolution of the detector. We account for this effect by convolving the spectrum of incoming photons with a smearing function R (E γ − E γ ), which represents the probability for the detector to register a photon energy E γ , given an actual incoming photon energy E γ . For concreteness, we consider a Gaussian smearing function of the form where is a dimensionless parameter which sets the overall, energy-dependent standard deviation σ(E γ ) = E γ of the distribution. Thus, we have our final result for the extra-galactic photon flux: with the individual fluxes dΦ i /dE γ given by Eq. (6.7) and with R (E γ − E γ ) given by Eq. (6.8).
Limits on additional contributions to the extra-galactic photon flux at energies E γ 800 keV can be derived from measurements of the total gamma-ray flux by instruments such as COMPTEL, EGRET, and Fermi-LAT. For photon energies in the range 800 keV E γ 30 MeV, COMPTEL data on the diffuse extra-galactic background spectrum [57] are well modeled by a powerlaw expression of the form MeV −1 cm −1 s −1 sr −1 . (6.10) Likewise, for energies in the range 30 MeV E γ 1.41 GeV, EGRET data [58] are well modeled by the power-law expression At even higher energies, the spectrum of the extragalactic diffuse photon background can be inferred from Fermi-LAT data. Within the energy range 100 MeV E γ 820 GeV, these data are well described by a powerlaw expression with an exponential suppression at high energies [59]: (6.12) Finally, we note that while we focus here primarily on photons with energies above 800 keV, estimates of the extra-galactic photon background at lower energies have been computed from BAT [60] and INTEGRAL [61] data and from observations of Type-I supernovae [62]. The expression in Eq. (6.9) is completely general and applicable to any ensemble of particles which decays to photons with energies within the transparency window O(keV) E γ O(TeV). To illustrate how this expression can be evaluated in practice, we consider a concrete example involving a specific set of injection spectra dN i /dE γ for the ensemble constituents. In particular, we consider a scenario in which each constituent χ i decays with a branching fraction of effectively unity to a pair of photons via the process χ i → γ + γ. Under the assumption that the χ i are "cold" (i.e., non-relativistic) by the time t LS , each of these photons has energy E γ ≈ m i /2 in the comoving frame at the moment of injection. We thus have While higher-order processes will generically also give rise to a continuum spectrum which peaks around E γ ∼ O(10 −3 m i ), we focus here on the contribution from the photons with E γ ≈ m i /2, since line-like features arising from monochromatic photon emission are far easier to detect than continuum features. Thus, for the injection spectrum in Eq. (6.13), the total contribution to the differential extra-galactic diffuse photon flux from the decays of the constituents within this ensemble is given by (6.14)

VII. APPLICATIONS: TWO EXAMPLES
In Sects. III-VI, we examined the observational limits on electromagnetic injection from particle decays after the BBN epoch and formulated a set of constraints which can be broadly applied to any generic ensemble of unstable particles. In this section, we examine the implications of these constraints by applying them to a pair of toy ensembles which exhibit a variety of scaling behaviors for the constituent masses m i , extrapolated abundances Ω i , and lifetimes τ i . For simplicity, we shall focus primarily on the "early-universe" regime in which 1 s τ i t LS for all of the χ i . Within this regime, the leading constraints on injection are those from modifications of the abundances of light nuclei and from distortions in the CMB.
We consider two different classes of mass spectra for our decaying ensemble constituents. The first consists of spectra in which the masses of these decaying particles are evenly distributed on a linear scale according to the scaling relation where k = 0, 1, ..., N − 1 and where m 0 and ∆m are taken to be free parameters which describe the mass of the lightest ensemble constituent and the subsequent mass splittings respectively. Mass spectra of this approximate linear form are realized, for example, for the KK excitations of a particle propagating in a flat extra spacetime dimension whose total length is small compared with 1/m 0 . The second class of mass spectra we consider are spectra in which the particle masses are evenly distributed on a logarithmic scale according to the scaling relation where m 0 once again denotes the mass of the lightest ensemble constituent, and where ξ is a dimensionless free parameter. A mass spectrum of this sort is expected, for example, for axion-like particles in axiverse constructions [16]. A variety of scaling relations for the decay widths Γ i of the χ i across the ensemble can likewise be realized in different scenarios. For simplicity, in this paper we focus on the case in which each χ i decays with a branching fraction of effectively unity to a pair of photons via the process χ i → γ + γ. Moreover, we shall assume that the scaling relation for Γ i takes the form where Λ is a parameter with dimensions of mass and where C i is a dimensionless coefficient -a coefficient which in principle can also have a non-trivial dependence on m i . Such a scaling relation arises naturally, for example, in situations in which each of the χ i decays via a non-renormalizable Lagrangian operator of dimension d = 5 in an effective theory for which Λ is the cutoff scale. For purposes of illustration, we focus on the case in which the coefficient C i = C χ is identical for all C i . In this case, the Γ i are specified by a single free parameter: the ratio M * ≡ Λ/ C χ . We emphasize that it is possible -and indeed in most situations even expected -that C χ 1. Thus, unlike Λ itself, the parameter M * may exceed the reduced Planck mass M P without necessarily rendering the theory inconsistent. However, we remark that for values of M * in this regime, the Weak Gravity Conjecture [63] imposes restrictions on the manner in which the decay operator in Eq. (7.3) can arise.
Finally, we consider a range of different possible scaling relations for the extrapolated abundances Ω i across the ensemble. In particular, we consider a family of scaling relations for the Ω i of the form (7.4) in which the abundance Ω 0 of the lightest ensemble constituent and the scaling exponent γ are both taken to be free parameters. In principle, γ can be either positive or negative. However, since τ i ∝ m −3 i and since the leading constraints on injection at timescales t inj t LS are typically less stringent for earlier t inj , the most interesting regime turns out to be that in which γ ≥ 0 and the initial energy density associated with each χ i is either uniform across the ensemble or decreases with increasing τ i .
In summary, then, each of our toy ensembles is characterized by a set of six parameters: {N, m 0 , ∆m, M * , Ω 0 , γ} for the mass spectrum in Eq. (7.1) and {N, m 0 , ξ, M * , Ω 0 , γ} for the mass spectrum in Eq. (7.2). Note that for any particular choice of these parameters, the lifetimes of the individual ensemble constituents span a range τ N −1 ≤ τ i ≤ τ 0 .

A. Results for Uniform Mass Splittings
We begin by examining the results for a toy ensemble with a mass spectrum given by Eq. (7.1). For such a mass spectrum, the density n m of states per unit mass is uniform across the ensemble; however, it is also useful to consider the density n τ of states per unit lifetime, which is more directly related to the injection history. In the continuum limit in which the difference between the lifetimes of sequential states χ i and χ i+1 in the ensemble is small and the lifetime τ may be approximated as a continuous variable, we find that The density of states per unit lifetime is therefore greatest within the most unstable regions of the ensemble. Moreover, we also note that the density of states per unit ln(τ ) also decreases with ln(τ ). These considerations turn out to have important implications for the bounds on such ensembles, as shall become apparent below. In Fig. 9, we show the constraints on an ensemble comprising N = 3 unstable particles in which the lightest ensemble constituent has a mass m 0 = 10 GeV. We consider this ensemble to be "low-density" in that it has a density of states per unit τ which is small throughout the range of parameters shown. (This will later be contrasted with analogous results in Fig. 10 for a "high-density" ensemble.) The contours shown in the left and center panels of Fig. 9 represent upper bounds on the total extrapolated abundance Ω tot as functions of M * for fixed ∆m. The contours shown in the left panel correspond to the choice of ∆m = 15 GeV; those in the center panel correspond to the choice of ∆m = 95 GeV. In the right panel, we show the corresponding contours as functions of ∆m for fixed M * = 10 19.5 GeV. The blue, green, and red curves shown in each panel represent the bounds on Ω tot from constraints on the primordial abundances of D, 7 Li, and 6 Li, respectively. We use a uniform color for the D contour -a contour which reflects the sum of contributions to δY D from both production and destruction processes. The dotted, dashed, and solid curves of each color correspond respectively to the choices γ = {0, 1, 2} for the scaling exponent in Eq. (7.4). The corresponding black curves shown in each panel represent the bounds arising from limits on distortions of the CMBphoton spectrum obtained by taking whichever bound (either that from δµ or that from δy C ) is more stringent at every point.
In interpreting the results shown in Fig. 9, we note that each choice of parameters specifies a particular range of lifetimes τ N −1 ≤ τ i ≤ τ 0 for the ensemble constituents. This range of lifetimes varies across the range of M * shown in the left panel, from 1.03 × 10 4 s ≤ τ i ≤ 6.58 × 10 5 s for M * = 10 16.5 GeV to 1.03 × 10 10 s ≤ τ i ≤ 6.58 × 10 11 s for M * = 10 19.5 GeV. Across this entire range of M * values, the range of τ i values for the ensemble is reasonably narrow, spanning only roughly a single order of magnitude. As a result, the bounds on Ω tot for a decaying ensemble shown in the left panel of Fig. 9 as functions of M * closely resemble the bounds on Ω χ for a single decaying particle species shown in Fig. 6 as a function of τ χ . Moreover, since the τ i for all particles in the ensemble are roughly comparable, the results are not particularly sensitive to γ, which determines how Ω tot is distributed across the ensemble. The most stringent constraint on Ω tot is the bound associated with the primordial D abundance over almost the entire range of M * shown. The only exception occurs in the regime where M * is large and the lifetimes of all of the ensemble constituents are extremely long, in which case the constraint associated with CMB distortions supersedes the D constraint. We note that for values of Ω tot which lie above the D contour, D may be either overabundant or underabundant. We also note that a "funnel" similar to that appearing in Fig. 6 is visible in the D constraint contours for all values of γ shown in the figure, and that this funnel occurs where M * is such that the lifetimes of the ensemble constituents are approximately τ i ≈ 10 6.5 s. By contrast, the results shown in the center panel of Fig. 9 correspond to the regime in which the range of lifetimes spanned by the ensemble is quite broad, extending over several orders of magnitude in τ i . In particular, the range of lifetimes varies across the range of M * shown in this panel from 1.03 × 10 4 s ≤ τ i ≤ 8.22 × 10 7 s for M * = 10 17.5 GeV to 8.22 × 10 7 s ≤ τ i ≤ 6.58 × 10 11 s for M * = 10 19.5 GeV. We note that the range of M * shown in this panel has been truncated relative to the left panel in order to ensure that the characteristic decay timescales for all ensemble constituents occur well after  Fig. 6, we see that the bounds on a decaying ensemble differ not only quantitatively but even qualitatively from the bounds on its individual constituents. the BBN epoch. For an ensemble with such a broad range of lifetimes, the effect of distributing Ω tot over multiple particle species is readily apparent. Indeed, the shapes of the constraint contours in the center panel of Fig. 9 bear little resemblance to the shapes of the contours in Fig. 6. Thus, comparing the results shown in the center panel of Fig. 9 with those in Fig. 6, we see that the bounds on a decaying ensemble indeed differ not only quantitatively but even qualitatively from the bounds on its individual constituents. For example, two separate funnels arise in the D constraint contour for γ = 2. Each of these funnels represents not merely a cancellation between the production and destruction contributions to Y D from a single particle, but rather a cancellation among the individual contributions from all constituent particle species in the ensemble. The funnel on the left arises due to a cancellation between the net negative contributions to Y D from χ 1 and χ 2 and the net positive contribution from χ 0 . By contrast, the funnel on the right arises due to a cancellation between a net negative contribution to Y D from χ 2 and positive contributions from χ 0 and χ 1 , with the latter contribution suppressed because τ 1 ≈ 10 6.5 s. The positions of these funnels, then, reflect the collective nature of the ensemble -and they have important consequences. Indeed, within certain regions of parameter space within which the D constraints are weakened by cancellations among δY D contributions from different ensemble constituents, the constraint associated with the 6 Li abundance actually represents the leading bound on Ω tot for the ensemble.
Finally, in the right panel of Fig. 9, we show how the bounds on Ω tot vary as a function of ∆m for fixed M * = 10 19.5 GeV. We choose this benchmark for M * because the CMB-distortion constraints and the leading constraints from light nuclei are comparable for this choice of M * . In interpreting the results displayed in this panel, it is important to note that τ 0 is independent of ∆m. By contrast, the τ i for all other χ i decrease with increasing ∆m. Since, broadly speaking, the impact that a decaying particle has on both CMB distortions and on the primordial abundances of light nuclei tends to decrease as τ i decreases, the bounds on decaying ensembles generally grow weaker as ∆m increases. This weakening of the constraints generally grows more pronounced as γ increases and as a larger fraction of the abundance is carried by the heavier ensemble constituents. For γ = 0, on the other hand, the fraction of Ω tot carried by the lightest ensemble constituent χ 0 is independent of ∆m. As a result, many of the constraint contours shown in the right panel of Fig. 9 become essentially flat for sufficiently large ∆m as the lifetimes of the remaining χ i become so short that the impact of their decays on the corresponding cosmological observables is negligible in comparison with those of χ 0 .
Thus far, we have focused on the regime in which the spacing between the lifetimes τ i of the ensemble constituents is large -in other words, the regime in which the density of states per unit τ is small across the ensemble. However, it is also interesting to consider the opposite regime, in which the density of states per unit τ is sufficiently large over the entire range τ N −1 ≤ τ ≤ τ 0 that the ensemble effectively acts as a continuous source of electromagnetic injection over this interval. This range of lifetimes is completely determined in our toy model by M * , m 0 , and the quantity m N −1 = (N − 1)∆m + m 0 . Thus, we may make a direct comparison between the bounds on such "high-density ensembles" and the bounds on "low-density" ensembles with the same τ i range by varying N and ∆m oppositely while holding M * , m 0 , and m N −1 fixed. Towards this end, in Fig. 10, we show the upper bounds on Ω tot for an ensemble with a higher density of states per unit lifetime. The constraint contours in the left and center panels of the figure represent the same choices of m 0 and m N −1 as in the corresponding panels of Fig. 9, but for N = 20 rather than N = 3. For these parameter choices, the ensemble is effectively in the "high-density" regime over the entire range of M * shown in each panel. In the right panel of Fig. 10, we once again show constraint contours as functions of ∆m for fixed M * .
As one might expect, the primary effect of distributing the abundance of the ensemble more evenly across the same range of lifetimes is that features in the constraint contours associated with particular χ i are in large part smoothed out. While the D funnels in the left panel are still present, this reflects the fact that the range of τ i for the ensemble is sufficiently narrow that distributing the abundance more democratically across this range has little impact on the constraint contours. A single D funnel also appears in the constraint contour for γ = 2 in the center panel of Fig. 10. The presence of this feature reflects the fact that for γ = 2, the vast majority of Ω tot is carried by the most massive -and therefore most unstable -ensemble constituents. Moreover, the density of states per unit lifetime is also much higher for these shorter-lived states than it is for the rest of the ensemble. The funnel can be understood as corresponding to the value of M * for which the lifetimes of these states, which are smaller than but roughly similar to τ N −1 , satisfy τ i ∼ τ N −1 ≈ 10 6.5 s. Once again, however, we note that the location of the funnel is shifted slightly away from this value of M * due to the collective contribution to Y D from the longer-lived constituents in the ensemble. A similar feature is also apparent in the constraint contour for γ = 2 in the right panel of Fig. 10.
We can gain further insight into the results displayed in Fig. 10 by examining the continuum limit in which ∆m → 0.
In this limit, the sums appearing in Eqs. (3.52)-(3.54) and in Eqs. (4.17) and (4.18) can be recast as integrals over the continuous lifetime variable τ -integrals which may be evaluated in a straightforward manner, yielding analytic expressions for δY a , δµ, and δy C . For example, in the case in which the τ i span the entire range from t Aa to t f a for each relevant nucleus, we find that for any γ ≥ 0 the term δY where Γ(n, z) is the incomplete gamma function, where we have defined where A a and B a are defined in Table II and where the limits of integration are where the limits of integration are given by expressions analogous to those appearing in Eq. (7.8).
Expressions for the CMB-distortion parameters δµ and δy C may be obtained in a similar fashion. In particular, for the case in which the τ i span the entire range from t e to t LS , we find that δµ takes the form x min 2µ (7.10) for all γ ≥ 0, except for the special case in which γ = 1. For this special case, the replacement should be made in the fourth line of Eq. (7.10). Likewise, we find that δy C takes the form x min 2y , (7.12) where 2 F 1 (a, b; c; z) denotes the ordinary hypergeometric function, where A a , B a , and α a are defined in Table III, and where the limits of integration are We now consider the results for a toy ensemble with a mass spectrum given by Eq. (7.2). The density of states per unit mass is inversely proportional to m in this case, and the corresponding density of states per unit lifetime is which likewise decreases with τ . However, we note that the density of states per unit ln(τ ) in this case remains constant across the ensemble.
In Fig. 11 we show the bounds on Ω tot for a "lowdensity" ensemble with a mass spectrum given by Eq. (7.2) and the same benchmark parameter choices N = 3 and m 0 = 10 GeV as in Fig. 9. In the left and center panels, the value of ξ has been chosen such that the range of τ i for the ensemble constituents is the same as in Fig. 9 at both ends of the range of M * shown. For the mass spectrum considered here, this range of τ i depends on m 0 , on M * , and on the quantity ξ N −1 . In the right panel of Fig. 11, we show how the bounds on Ω tot vary as a function of ξ for fixed M * = 10 19.5 GeV. The contours shown in Fig. 11 exhibit the same overall scaling behaviors as those in Fig. 9. However, there are salient differences which reflect the fact the the density of states per unit ln(τ ) is uniform in this case across the ensemble.
In Fig. 12, we show the corresponding results for an ensemble with m 0 = 10 GeV and a mass spectrum again given by Eq. (7.2), but with a sufficiently high n τ throughout the range of M * or ξ shown in each panel that the ensemble is effectively always within the "highdensity" regime. We find that taking N = 10 is sufficient to achieve this. The left and center panels of Fig. 12 correspond to the same ranges of M * and the same value of ξ N −1 as the left and center panels of Fig. 11. The right panel once again shows how the bounds on Ω tot vary as a function of ξ for fixed M * = 10 19.5 GeV. Once again, we see that the principal consequence of increasing n τ is that constraint contours become increasingly smooth and featureless.
Once again, as we did for the case of an ensemble with a uniform mass splitting, we may derive the corresponding expressions for δY a , δµ, and δy C in the case of an exponentially rising mass splitting in the continuum limit by direct integration of the expressions in Eqs. (3.52)-(3.54) and in Eqs. (4.17) and (4.18). For the mass spectrum in Eq. (7.2), this limit corresponds to taking ξ → 1. For the case in which the τ i span the entire range from t Aa to t f a for each relevant nucleus, we find that for γ > 0, the term δY    Fig. 11, but for "high-density" ensembles with N = 10. In the left and center panels, we plot the upper bound on Ωtot as a function of the coupling-suppression scale M * with fixed ξ = 1.2 and with fixed ξ = 1.39, respectively. In the right panel, we plot the bound as a function of ξ with fixed M * = 10 19.5 GeV. where the limits of integration are given by expressions analogous to those appearing in Eq. (7.8). Once again, as was the case with δY (1) a , the expression for δY (2) a in Eq. (7.19) requires modification for the special case in which γ = 0. In particular, the replacement specified in Eq. (7.17) should likewise be made in the last term in square brackets on the second line of Eq. (7.19). Moreover, for γ = 0, the quantity Ξ(t) is given by Eq. (7.18) rather than by Eq. (7.16).
The corresponding expressions for δµ and δy C for the case in which the τ i span the entire range from t e to t LS can likewise be computed by direct integration. In particular, we find that δµ is given by x min 2µ (7.20) for all γ ≥ 0, except for the special case in which γ = 2. For this special case, the replacement should be made in the fourth line of Eq. (7.20). Likewise, we find that δy C is given by x max 0y x min 0y + B y Ξ(t 1y ) 6x α+ 3−2γ 6 6α µ − 2γ + 3 2 F 1 1, 1 + 3 − 2γ 6α y ; 2 + 3 − 2γ 6α y ; −x αy x max 1y x min 1y + B y Ξ(t 2y ) 3x x min 2y .

(7.22)
The limits of integration in Eqs. (7.20) and (7.22) are defined as in Eq. (7.8), and Ξ(t) is defined as in Eq. (7.16) for γ > 0 and as in Eq. (7.18) for γ = 0. Once again, we note that in cases in which the τ i do not span the relevant range of timescales during which particle decays can affect one of these cosmological observables, the limits of integration in Eqs. (7.15)-(7.22) should be replaced by the values which restrict the overall range of lifetimes appropriately.
In summary, the results shown in Figs. 9-12 illustrate some of the ways in which the constraints on injection from unstable-particle decays can be modified in scenarios in which the injected energy density is distributed across an ensemble of particles with a range of lifetimes. These results also illustrate that a nontrivial interplay between the contributions from different decaying particle species within the ensemble can have unexpected and potentially dramatic effects on the upper bound on Ω tot -as when cancellations among positive and negative contributions to δY D from a broad range of particles within the ensemble result in a significant weakening of this upper bound. Thus, we see that the bounds on a decaying ensemble can exhibit collective properties and behaviors that transcend those associated with the decays of its individual constituents.

VIII. CONCLUSIONS: DISCUSSION AND SUMMARY OF RESULTS
In this paper, we have considered ensembles of unstable particle species and investigated the cosmological constraints which may be placed on such ensembles due to limits on electromagnetic injection since the conclusion of the BBN epoch. Indeed, as we have discussed, such injection has the potential to modify the primordial abundances of light nuclei established during BBN. Such injection can also give rise to spectral distortions in the CMB, alter the ionization history of the universe, and leave observable imprints in the diffuse photon background. For each of these individual considerations, we have presented an approximate analytic formulation for the corresponding constraint which may be applied to generic ensembles of particles with lifetimes spanning a broad range of timescales 10 2 s τ i t now . For ease of reference, these analytic approximations, along with the corresponding equation numbers, are compiled in Table IV. In deriving these results, we have taken advantage of certain linear and uniform-decay approximations. We have also demonstrated how these results can be applied within the context of toy scenarios in which the mass spectrum for the decaying ensemble takes one of two characteristic forms realized in certain commonly-studied

Constraint
Analytic formulation

Primordial abundances of light nuclei
Analytic approximations for δYa given in Eq. (3.51), with the corresponding limits and parameter values given in Table II.

Spectral distortions in the CMB
Analytic approximation for δµ and δy C given in Eqs. (4.17) and (4.18), with the parameter values given in Table III  Diffuse extra-galactic photon background Differential signal flux given by Eq. (6.9), with differential background fluxes given in Eqs. (6.10)-(6.12). extensions of the SM.
Several comments are in order. First, it is worth noting that while the values of the parameters in Table II were derived assuming a particular set of initial comoving number densities Y a for the relevant nuclei at the end of the BBN epoch, the corresponding parameter values for different sets of initial comoving number densities Y a may be obtained in a straightforward manner from the values in Table II. Provided that Y a and Y a are not significantly different, the characteristic timescales t Ba , t Ca , and t Xa for each process are to a good approximation unchanged, and a shift in the initial comoving number densities can be compensated by an appropriate rescaling of the relevant normalization parameters A a and B a . For example, a shift in the initial helium mass fraction at the end of BBN can be compensated by a rescaling of the fit parameters A * D and B * D associated with D production and of the parameters A4 He and B4 He associated with 4 He destruction by Y p /Y p , as well as a rescaling of the parameters A * 6 Li and B * 6 Li associated with secondary 6 Li production by (Y p /Y p ) 2 . The quadratic dependence on Y p /Y p exhibited by the last of these rescaling factors reflects the fact that the reactions through which the non-thermal populations of 3 He and T nuclei are initially produced and the reactions through which these nuclei in turn contribute to 6 Li production both involve 4 He in the initial state. Similarly, a shift in the initial 7 Li abundance can be compensated by a rescaling of the fit parameters A7 Li and B7 Li associated with 7 Li destruction and of the parameters A6 Li and B6 Li associated with primary 6 Li production by Y 7 Li /Y7 Li , while a shift in the initial D abundance can be compensated by a rescaling of the parameters A D and B D associated with D destruction by Y D /Y D .
Second, we remark that in this paper we have focused primarily on constraints related to electromagnetic injection. In situations in which a significant fraction of the energy liberated during unstable particle decays is released in the form of hadrons, the limits obtained from bounds on the primordial abundances of light nuclei may differ significantly from those we have derived here. Detailed analyses of the limits on hadronic injection in the case of a single decaying particle species have been performed by a number of authors [3,[64][65][66]. It would be interesting to generalize these results to the case of a decaying ensemble as well.
Third, as we have seen, there are two fundamentally different timescale regimes in which the corresponding physics is subject to very different leading bounds: an "early" regime 10 2 s < ∼ t inj < ∼ 10 12 s, and a "late" regime t inj > ∼ 10 12 s. This distinction is particularly important when considering the implications of our results within the context of the Dynamical Dark Matter framework [18,19]. Within this framework, the more stable particle species within the ensemble provide the dominant contribution to the dark-matter abundance at present time, while the species with shorter lifetimes have largely decayed away. However, the masses, decay widths, abundances, etc., of the constituents of a realistic DDM ensemble are governed by the same underlying set of scaling relations. Thus, in principle, one might hope to establish bounds on the ensemble as a whole within the DDM framework by constraining the properties of those ensemble constituents which decay on timescales τ t LS -i.e., "early" timescales on which the stringent limits on CMB distortions and the primordial abundances of light nuclei can be brought to bear.
In practice, however, it turns out that these considerations have little power to constrain most realistic DDM scenarios. On the one hand, both the lifetimes τ i and extrapolated abundances Ω i of the individual ensemble constituents typically decrease monotonically with i in such scenarios -either exponentially [67] or as a power law [19,68,69]. On the other hand, as indicated in Fig. 8, the constraint associated with the broadening of the surface of last scattering due to ionization from particle decays becomes increasingly stringent as the lifetime of the particle decreases down to around τ i ∼ t LS . This implies that for an ensemble of unstable particles in which the collective abundance of the cosmologically stable constituents is around Ω DM ≈ 0.26, the Ω i for those constituents with lifetimes τ i t LS are constrained to be extremely small. Thus, the constraints associated with CMB distortions and with the primordial abundances of light nuclei in the "early" regime are not particularly relevant for constraining the properties of DDM ensembles in which the abundances scale with lifetimes in this way. Nevertheless, we note that in ensembles in which τ i and Ω i do not increase monotonically with i (such as those in Ref. [69]), the constraints associated with these considerations could indeed be relevant. This is currently under study [70].