Cosmological Bounds on sub-GeV Dark Vector Bosons from Electromagnetic Energy Injection

New dark vector bosons that couple very feebly to regular matter can be created in the early universe and decay after the onset of big bang nucleosynthesis (BBN) or the formation of the cosmic microwave background (CMB) at recombination. The energy injected by such decays can alter the light element abundances or modify the power and frequency spectra of the CMB. In this work we study the constraints implied by these effects on a range of sub-GeV dark vectors including the kinetically mixed dark photon, and the $B-L$, $L_e-L_\mu$, $L_e-L_\tau$ , and $L_\mu-L_\tau$ dark U(1) bosons. We focus on the effects of electromagnetic energy injection, and we update previous investigations of dark photon and other dark vector decays by taking into account non-universality in the photon cascade spectrum relevant for BBN and the energy dependence of the ionization efficiency after recombination in our treatment of modifications to the CMB.


Introduction
Direct measurements of the cosmos give strong evidence for a standard ΛCDM cosmology containing the elementary particles and forces predicted by the Standard Model (SM) [1,2]. Observations of the cosmic microwave background (CMB) [3,4], baryon acoustic oscillations (BAO) [5,6,7], and supernovae distances [8,9,10] fix the parameters of the ΛCDM model very precisely. The success of big bang nucleosyntheis (BBN) at predicting the light element abundances (notwithstanding the lithium puzzles) provides further support for a ΛCDM cosmology with SM particles and interactions going back even earlier in time, up to temperatures on the order of a few MeV [11,12,13,14,15,16,17,18].
In this work we apply and extend these methods to derive cosmological bounds from BBN and the CMB on a range of sub-GeV dark vector bosons. New vector bosons are ubiquitous in extensions of the SM, including grand unification [55,56,57], superstring theory [58,59,60], and theories of dark matter [61,62,63]. By definition, dark vector bosons couple only very feebly to the SM, and can thus be much lighter than the electroweak scale without having been discovered yet [64,65,66]. For very small couplings, dark vectors are very difficult to produce in the laboratory and can develop macroscopically long lifetimes. Thermal cosmological production is also reduced in this regime and the dark vectors might never develop a full thermal abundance. Even so, dark vector bosons can still be created in the early universe through direct reheating [67,68,69,70], freeze-in processes [71,72,73,74,75], or the decays or annihilations of heavier dark states [76,77]. The later decays of the vector bosons then inject energy into the cosmological medium, possibly modifying the light element abundances and the CMB relative to SM+ΛCDM [78,79].
We investigate a range of light (but massive) U(1) vector bosons including a dark photon (DP) that connects to the SM exclusively through gauge kinetic mixing with hypercharge [80,81], as well as the anomaly-free combinations B − L [82,83], L e − L µ [84], L e − L τ [84], and L µ − L τ [85,86] with only very tiny direct couplings to the SM. Our focus is on the mass range m V ∈ [1 MeV, 1 GeV] over which these dark vectors have significant decay fractions to electromagnetic channels, and we concentrate on the impact of the electromagnetic component of the energy injected by the decays which is typically the most important effect for lifetimes greater than τ V 10 4 s.
This study expands on previous investigations of cosmological bounds on dark photons in two ways [78,79]. First, applying the technology developed in Ref. [87], we compute the full electromagnetic cascade induced by dark vector decays, including some new effects relative to Refs. [78,79]. From this we derive the resulting photon cascade spectrum that is essential for determining the impact of electromagnetic energy injection on the light element abundances. Our calculation also differs from the more common approach of using the universal photon spectrum defined in Ref. [37], based on the results of Refs. [34,88,89,90,91], to derive the cascade photon spectrum. The universal spectrum method provides a very good description of the photon cascade produced by the interactions of highly-energetic electromagnetic primaries (E 100 GeV) with the cosmological background prior to recombination, and has the major simplifying feature that it only depends on the total amount of electromagnetic energy injected rather than the detailed primary injection spectrum. However, it was shown in Refs. [87,92,93,94] that lower-energy electromagnetic injection near or below the GeV scale can lead to cascade photon spectra that deviate significantly from the predictions of the universal spectrum, and that can depend on the type and spectrum of the energy [87].
In light of this first point, the second way that we expand on previous work is by computing the detailed primary electromagnetic injection spectra from sub-GeV vector boson decays [95,96]. We study the energy spectra of electrons and photons created by the decay modes V → {e + e − , µ + µ − , π + π − , π + π − π 0 , π 0 γ}, which are expected to be the most important ones for the dark vectors of interest [97,98]. We then apply these injection spectra to computing the full electromagnetic cascade spectrum produced by a given dark vector decay for our BBN studies or use them to calculate limits from the CMB.
The outline of this paper is as follows. In Sec. 2 we compute the decay rates and fractions of the dark vector species and we find the energy spectra of the resulting photons and electrons. Next, in Sec. 3 we study the effects of electromagnetic energy injected by dark vector decays on the light elements and use current data to set limits on the pre-decay dark vector densities as a function of their lifetime and mass. We investigate the related effects of dark vector decays on the power and frequency spectra of the CMB in Sec. 4. In Sec. 5 we review the calculation of the density of dark vectors from thermal production in the early universe, and we find the constraints from BBN and the CMB on the masses and couplings of dark vectors for the corresponding thermal yields. Finally, Sec. 6 is reserved for our conclusions. Some additional details on the departure from electromagnetic universality in BBN and when it occurs are presented in Appendix A.

Dark Vector Bosons and their Decays
We begin by reviewing the dominant decay channels for the sub-GeV dark vector bosons of interest: kinetically mixed dark photon (DP), B − L, L e − L µ , L e − L τ , L µ − L τ . In all cases we assume the dark vectors decay exclusively to the SM, with no additional exotic decay channels available. With the decay fractions in hand, we turn next to computing the electron and photon spectra they produce in the rest frame of the dark vector relevant for computing cosmological bounds on them.

Vector Boson Decay Branching Fractions
The sub-GeV dark vectors under consideration interact with the SM primarily through gauge kinetic mixing for the dark photon (DP), with kinetic mixing parameter , and via direct gauge coupling to the SM for the rest, where g V 1 and we normalize to |Q L | = 1 for leptons. The direct-coupling dark vectors will also have a gauge kinetic mixing with the photon but we assume (self-consistently) that it is sub-leading relative to the direct gauge interaction. We assume further that the dark vectors obtain a mass m V in the MeV-GeV range from the Higgs [99] or Stueckelberg [100,101] mechanisms, and that no decay channels to light SM-singlet states are available. In particular, for the direct-coupling vectors we assume as in Ref. [98] that the light neutrinos are Majorana and mostly left-handed (which requires a spontaneous breaking of the underlying gauge invariance), and that any SM singlets are heavy enough to be neglected.
With these couplings, the perturbative decay widths to SM fermions f are where N c is the number of QCD colours of the fermion, α V = 2 α for the dark photon (with Q L,R = Q em ), and α V = g 2 V /4π for the direct vectors (with Q R = 0 for the light neutrinos given the assumptions stated above).
For sub-GeV vectors, we must also consider decays of the vector to spefcific hadronic states. Since these decays are strongly suppressed for the leptonic vectors, we only include them for the DP and B −L vectors here. Our approach to describing hadronic dark vector decays uses a combination of the data-driven methods of Refs. [78,97,98] and the analytic description based on vector-meson-dominance (VMD) of Ref. [102]. In the VMD picture [103,104,105,106], which also guides the methods of Refs. [78,97,98], the ρ and ω vector mesons are treated as gauge bosons of a hidden U (2) flavour symmetry. 1 Dark vector decays to hadrons in this picture occur through direct couplings -the DP to the electromagnetic current and the B−L vector to the baryonic current -as well as through an induced kinetic mixing with the ρ and ω with approximate strength [102,97] where t A = diag 1/2, ±1/2 for A = ω, ρ is the U (2) generator, g √ 3 × 4π, Q V = diag 2/3, −1/3 for DP and Q V = diag 1/3, 1/3 for B −L, and g V = −e for the DP.
The leading hadronic decay mode of the dark photon is usually V → π + π − . In the VMD picture, it arises from a combination of the direct countribution of the charged pions to the electromagnetic current and the induced kinetic mixing with the ρ. The corresponding decay width is where the leading-order form factor is [106] with g ρππ g, m ρ = 775.25 MeV, and Γ ρ (s = m ρ ) = 149 MeV. The true form factor is modified relative to Eq. (6) by the energy dependence of the ρ self-energy as well as ρ-ω mixing. In our analysis we follow Ref. [97] and extract the form factor from e + e − → π + π − (γ) measurements by the BABAR experiment [107]. Note that ρ-ω mixing also allows the B−L vector to decay to π + π − but we do not include the effect because the branching fraction for this channel is always less than a couple percent.
Mixing with the ω allows the DP and B −L vectors to decay to π + π − π 0 and π 0 γ, with the effect being strongest near the ω mass. The width for V → π 0 γ is [102] For the decay V → π + π − π 0 , the width is [102] where g 2 /4π 3.0, the relevant form factor is with m ω = 782.65 MeV and Γ ω = 8.49 MeV, and [104] where E ± are the energies of the outgoing charged pions in the decay frame, p i = (E i , p i ) are the pion 4-momenta whose components can all be expressed in terms of E ± by momentum conservation, and the integration region covers  In Fig. 1 we show the branching fractions (BR) of the dark photon (left) and B − L vector (right) over the mass range of interest. The hadronic resonance structures are evident near the ρ and ω poles, while the leptonic modes dominate at lower masses. The branching fractions of the leptonic vectors can also be obtained straightforwardly and are mostly determined by the relative charges and open decay channels. For the L µ −L τ vector, decays below the muon threshold are dominated by neutrinos with only a BR 3 × 10 −5 fraction to e + e − (via induced kinetic mixing) [98]. Hadronic decays of these vector bosons were also computed in Ref. [98], and their branchings BR had 10 −3 are too small to have a relevant cosmological effect over the ranges we study.

Electromagnetic Injection Spectra
The analysis above shows that the relevant decay channels of the dark vectors of interest are: e + e − , µ + µ − , π + π − , π + π − π 0 , π 0 γ, νν. These decays will all ultimately produce a collection of photons, electrons, and neutrinos. The resulting energy spectra of photons and electrons from these decay modes, in the lab frame where the decaying vector boson is at rest, are needed for our study of cosmological constraints to follow. Together, the full spectra per decay are where the sum runs over the decay modes i = e + e − , µ + µ − , . . ., and e refers to the sum of electrons and positrons. 2 We compute the spectra for each of the exclusive decay modes here, concentrating on the electrons and photons produced at leading non-trivial order. Our analysis is similar to that of Ref. [95] with some differences and additions. See also Ref. [96] for an alternate approach based on HERWIG 7 [108].
This channel was considered previously in Ref. [87]. The leading-order electron spectrum is trivial with two electrons, each with energy E = m V /2, In addition to electrons, photons can be produced from final-state-radiation (FSR). Following Ref. [95], we use the spectrum per decay As shown in Ref. [87], these FSR photons can have an important effect on BBN for sub-GeV injection.
At leading order, this channel ultimately produces two electrons and four neutrinos. Beyond leading order, photons can be created through FSR and radiative muon decays.
To obtain the leading-order energy spectrum of electrons in the vector boson lab frame, we focus on the dominant µ − → e −ν e ν µ decay channel. In the muon rest frame, for muon polarization P ∈ [−1, 1], the electron energy spectrum is the well-known Michel form [109,110]: where y = 2E /m µ and ϑ is the angle between the electron direction and the polarization axis in this frame. For anti-muon decays, the same spectrum applies but with the sign of P reversed. However, for V → µ + µ − the net effect of muon polarization cancels and this term can be neglected. To get the lab frame spectrum we must boost the muon rest frame with Lorentz factor γ = m V /2m µ = 1/ 1 − β 2 . Choosing the z axis to lie along the muon boost direction and assuming the electron is emitted at polar angle θ from this axis in the muon rest frame, the electron energy and momenta in the lab frame are Changing variables from (E , Ω ) to (E, Ω ) for solid angle Ω in the muon frame, we obtain the full lab frame distribution where E = E (E, Ω), the first term in the integrand is the relevant Jacobian factor, the factor of 1/4π normalizes the solid angle integral, and the factor of two counts the identical contributions from the µ − and µ + branches.
The photon spectrum from muon decays is computed following Refs. [95,111] as a sum of FSR and radiative contributions, This approach neglects interference between the two channels, but this effect is expected to be very small due to the narrow width of the muon. We use the expression of Eq. (14) with m f = m µ for the FSR part. The radiative spectrum is the sum of contributions from the µ + and µ − decays. In the muon rest frame, both are equal to [95,111,112] The radiative photon spectrum in the lab frame is then obtained by boosting as in Eq. (17).
Vector decays in this channel produce electrons primarily through the chain π − →ν µ µ − with µ − → ν µνe e − (and its conjugate). We focus on this chain exclusively since the inclusive branching BR(π + → e + ν e ) 1.2344 × 10 −4 [113] is small enough to be numerically insignificant in our analysis. Photons are also produced through FSR and intermediate radiative decays.
To obtain the leading-order electron spectrum, we focus on the π − branch with π − → ν µ µ − since an identical contribution arises for the π + branch. In the π − rest frame, the µ − is produced with Lorentz factor Aligning the z-axis with the direction of the outgoing muon, the muon obtains polarization P = +1. Defining E as the electron energy in the muon rest frame and θ as the electron direction relative to the z-axis, the energy spectrum in this frame is given by Eq. (15) (with E → E and θ → θ ). The electron energy E spectrum in the pion rest frame then follows from the boosting method described above with One more boost in an arbitrary pion directionn = (sin α cos φ, sin α sin φ, cos α) with Lorentz factor γ = m V /2m π = 1/ 1 − β 2 is needed to transform to the lab frame where the vector is at rest. The electron energy in this frame is with The lab frame electron energy distribution is then obtained by a simple change of variables using uniform distributions for cos θ ∈ [−1, 1], cos α ∈ [−1, 1], φ ∈ [0, 2π], and E given by the Michel spectrum of Eq. (15). Note that for the positron spectrum from antimuon decays, the positron polarization is P = −1 while the sign of the P term in the Michel spectrum is reversed leading to the same distribution in E and cos θ .
The photon spectrum from this vector decay channel is obtained from a combination of FSR and internal radiative decays as in Eq. (18). For the FSR part, from the V → π + π − step, we use the result of Ref. [95], The radiative contributions come from the π − → µ −ν µ and µ − → e −ν e ν µ steps of the decay chain. These were studied along with the radiative contribution from π − → e −ν e in Ref. [95] where it was found that the total contribution is nearly completely dominated by the muon decay step. The contribution to the photon spectrum then follows by applying the boosting procedure described above to the photon spectrum of Eq. (19).
This channel produces a pair of boosted photons from the π 0 decay as well as a monochromatic photon directly. The π 0 and monochromatic photon energies in the lab frame are For the photons from the π 0 decay, we have E = m π 0 /2 in the π 0 rest frame and a Lorentz factor γ = E 0 /m π 0 = 1/ 1 − β 2 . Summing these contributions and applying our previous results, the full photon spectrum is therefore where For this channel we concentrate exclusively on the direct photons and electrons that arise from the dominant π 0 → γγ and π − → µ −ν µ , µ − → e −ν e ν µ decay chains. The charged and neutral pions are created with a distribution of energies of given by Eq. (8), and d 2 Γ 3π /dE + dE − given by the same expression but with the function I(m 2 V ) replaced by the integrand of Eq. (10). The photon distribution is thus just the photon spectrum from a boosted π 0 decay found above weighted by the distribution of π 0 energies, with γ = 1/ 1 − β 2 = E 0 /m π 0 and B(E 0 ) from Eq. (27). Similarly, the electron (plus positron) spectrum is calculated using the boosted electron spectrum found previously for charged pion decay weighted by the joint distribution of E + and E − energies.

BBN Bounds on Sub-GeV Energy Injection
Energy injected into the cosmological plasma can disrupt the predictions of standard BBN by altering the ratio of neutrons to protons or by destroying/creating light elements after they are formed. In this section we study the photodissociation of light elements from electromagnetic energy injected by decays of long-lived sub-GeV dark vectors in the early universe. This is expected to be the most important effect of dark vector decays on BBN for decay lifetimes τ V > 10 4 s and masses greater than a few MeV. We apply our results to constrain the pre-decay abundances of such vectors.

Methods for Calculating the Impact on BBN
Decays of sub-GeV dark vectors produce photons, electrons, and neutrinos, both directly and through intermediate muons and pions. To compute photodissociation effects from these decays, we make use of the branching fractions and the photon and electron energy spectra computed in the previous section. Note that for the cosmological times at which photodissociation is effective, t 10 4 s, the muons and pions injected by sub-GeV vectors decay before they are slowed significantly by interactions with the cosmological background [39,114].
Electrons and photons injected into the dense cosmological plasma generally interact with the plasma before reacting with light nuclei. The scattering of hard primaries off the background plasma creates an electromagnetic (EM) cascade that remains strongly suppressed until well after the main element creation stage of BBN has completed. Since the formation of the EM cascade is fast relative to the rate of scattering with nuclei, the resulting photon cascade spectrum can be treated as the source for subsequent photodissociation reactions on nuclei.
To compute the electromagnetic cascade spectra of photons and electrons (with positrons counted as electrons here), we follow the methods of Ref. [87] based on the earlier works of Refs. [90,91]. Defining the differential number density per unit energy of photons or electrons in the cascade by the cascade spectra evolve according to where Γ a is the net damping rate for species a and energy E and S a (E) is the injection rate from all sources at this energy. The relevant damping and transfer reactions are generally fast relative to the Hubble rate and the effective photodissociation rates with light nuclei, and therefore the quasistatic limit of dN a /dt → 0 is a good approximation [37]. 3 This gives with both terms varying adiabatically as functions of time (or temperature).
The damping rate Γ a in Eq. (31) describes any reaction that transfers energy from species a at energy E to lower energies. The source terms in the electromagnetic cascade receive contributions from direct injection as well as from transfer reactions moving energy from higher up in the cascade down to energy E. Together, these take the explicit forms where R is the injection rate, dN a /dE is the primary energy spectrum per injection, E X is the maximum energy in the cascade, and K ab (E, E ) is transfer kernel describing reactions b(E ) + X BG → a(E) + X BG within the cascade. In the case of energy injection from the decays of species V , E X ≤ m V /2 and the injection rate is where τ V is the decay lifetime and n 0 V is the number density the species would have in the in the absence of decays.
Multiple reactions contribute to the damping and source terms in the cascade. For photons, the dominant contribution to damping at higher energies E ≥ E c ≡ m 2 e /22T is the photon-photon pair production (4P) reaction, γ+γ BG → e + +e − [91]. In the case of electrons, at the relevant energies damping is dominated by inverse Compton (IC) scattering e ∓ + γ BG → e ∓ + γ [115]. Together, these processes strongly suppress the electromagnetic cascade below nuclear dissociation thresholds until after the period of element formation, with E c > 2 MeV only for T 6 keV (t 7600 s). As a result, BBN bounds on electromagnetic injection from decay fall off very quickly for lifetimes τ V 10 4 s.
For very high-energy initial injection, with E X T m 2 e , the photon cascade spectrum resulting from the 4P, IC, and other reactions is found to have a universal form. Specifically, the universal photon spectrum depends only on the total amount of electromagnetic energy injected into the cosmological plasma, and not on the detailed injection spectra or whether it came in the form or photons or electrons [90,91]. However, recent studies of electromagnetic injection at lower energies have found significant deviations from universality [87,92,93,94]. These deviations typically occur when the primary injection energy falls below the 4P cutoff E c or nuclear dissociation thresholds. For the sub-GeV dark vectors of interest in this work, we find that this takes place throughout a very significant portion of the relevant parameter space. More details on deviations from universality and when it occurs are collected in Ref. [87] and Appendix A.
With the photon cascade spectrum in hand, the effect of photodissociation on the light element abundances can be described by Boltzmann equations of the form where N γ (E γ ) are the photon cascade spectra, A and the sums run over the relevant isotopes, and Y A are isotope number densities normalized to the entropy density, Reactions initiated by electrons are not included because their cascade spectrum is always strongly suppressed by IC scattering.
The nuclear species included in our analysis are hydrogen (H), deuterium (D = 2 H), tritium (T = 3 H), helium-3 ( 3 He), and helium (He = 4 He). Heavier species such as lithium have much smaller primordial abundances and their inclusion would not alter our results for the light elements we consider. For the nuclear cross sections in Eq. (35), we use the simple parametrizations collected in Ref. [37]. These cross sections all have the same general shape with a sharp rise at threshold up to a peak followed by a smooth fall off. The two most important thresholds for our study are E th 2.22 MeV for deuterium photodissociation and E th 19.81 MeV for helium, with the other relevant reaction thresholds typically in the range of E th ∼ 2.5 − 8.5 MeV.
To solve the coupled evolution equations of Eq. (35), we convert time to redshift and use initial primordial abundances expected from standard BBN predicted by PArthENoPE [116,117]: We then compare the resulting light element abundances to the following observed values, quoted with effective 1σ uncertainties (with theoretical and experimental uncertainties combined in quadrature): n3 He n H = (1.0 ± 0.5) × 10 −5 (Ref. [120]) .
The helium mass fraction Y p we use is consistent with Ref. [121] and earlier determinations but significantly lower than the result of Ref. [122]. The uncertainty on the ratio n D /n H is dominated by a theory uncertainty on the rate of photon capture on deuterium from Ref. [123]. For n3 He /n H , we use the determination of (n D + n3 He )/n H of Ref. [120] along with the value of n D /n H from Ref. [119]. The uncertainties quoted here are generous, and in the analysis below we implement exclusions at the 2σ level.
To illustrate our approach, we show in Fig. 2 the photodissociation bounds from BBN on the pre-decay yield Y V = n V /s for a decaying species V with mass m V and lifetime τ V assuming that all decays occur exclusively via one of the channels V → e + e − , µ + µ − , π + π − , π 0 π + π − , π 0 γ. Decays to e + e − produce the strongest effect throughout most of the lifetime-mass plane shown. Relative to the µ + µ − and π + π − channels, electron decays produce more energetic electrons and a greater electromagnetic injection fraction. At lower masses and larger lifetimes, in the lower right region of the plots, the injected electrons scatter off background photons in the Thomson regime where the upscattered photons receive energies well below nuclear photodissociation thresholds. The dominant sources of photons contributing to photodissociation are then FSR and radiative decays, and these also tend to be greater for e + e − than µ + µ − or π + π − channels. This obstacle does not affect the π 0 π + π − and π 0 γ channels which produce photons at the leading order and thus stronger BBN limits in the lower right region of the plots. Note that even after rescaling by the total electromagnetic fractions produced by these various modes, the resulting BBN limits are significantly different, reflecting the breaking of universality of the cascade photon spectrum, as discussed in more detail in Appendix A.

BBN Bounds on Dark Vectors
Using the dark vector decay fractions and the BBN methods described above, we can now derive BBN limits due to photodissociation on the pre-decay yield Y V = n V /s of a dark vector with mass m V and lifetime τ V . Our results are shown in Fig. 3 for the dark vector , V → π + π − (lower left), V → π + π − π 0 (lower middle), and V → π 0 γ (lower right). species DP (upper left), B −L (upper middle), L e −L µ (lower left), L e −L τ (lower middle), and L µ −L τ (lower right). The shapes of these exclusions can be understood by comparing the branching fractions shown in Fig. 1 and the bounds on the contributing decay modes in Fig. 2. For the dark photon (DP) and B−L vectors, the bounds get weaker due to hadronic contributions near the ρ and ω resonances. This effect is greater for the DP due to its strong mixing with the wide ρ resonance relative to the B −L that approximately only mixes with the narrow ω. With the exception of the DP, an increase in the bound is seen near the muon threshold at m V = 2m µ due to the larger visible mode decay branching fractions when the muon channel turns on and a smaller fraction of the decays go to neutrinos. For the L µ −L τ vector, the bounds below the muon threshold fall beneath the range covered by the figure since the decay products in this region are nearly completely dominated by neutrinos.
In addition to photodissociation, dark vector decays can also modify the outcome of standand BBN by altering the ratio of neutrons to protons and through the hadrodissociation of light elements. While these effects can be significant and further constrain the properties of dark vectors, we argue that they are also largely orthogonal to the bounds from photodissociation derived above. Lighter dark vectors with m V 20 MeV can alter the effective number of light degrees of freedom during and after neutron-proton freeze-out through their direct influence [25], from their decays following neutrino decoupling [25,124], or by equilibrating light right-handed neutrinos (which we assume to not be present) [82,125]. The couplings required for these effects to be significant correspond to dark vector lifetimes below τ V 10 4 s, or large initial densities m V Y V 10 −9 , and thus these considerations apply to earlier decays or produce weaker constraints than photodissociation. Heavier dark vectors that decay to pions can alter the neutron to proton ratio through reactions such as p+π − → n+π 0 or destroy light elements through hadrodissociation. The analyses of Refs. [78,79,114] find that these effects are strongest for lifetimes τ V ∼ 10 3 s and m V Y V 10 −11 GeV. Comparing to our results for photodissociation, these two sets of bounds largely apply independently of one another, with only a very small region of possible interference near τ V ∼ 10 4 s. Taken together, these considerations justify our earlier treatment of photodissociation bounds from BBN in isolation.

CMB Bounds on Sub-GeV Energy Injection
Energy injection in the early universe can also modify the power and freqency spectra of the CMB. In this section we describe the methods that we use to compute these effects and show the constraints they imply for dark vector decays.

Methods for the CMB Power Spectrum
Energy injected during and after recombination can ionize newly-formed atoms and broaden the surface of last scattering. These effects can modify the power spectra of CMB fluctuations relative to the standard recombination history [51]. Precision measurements of the CMB power spectra can therefore be used to constrain new sources of energy injection in this period [52,53].
The total rate of energy injection per unit volume from a decaying species of mass m V is simply m V R, where R is the decay rate per unit volume given in Eq. (34). For decays near or after recombination, there may be a delay between the initial decay and when the resulting energy is deposited within the cosmological medium [54]. We treat this as in Refs. [126,127] and write the total energy deposition rate per unit volume as where n 0 V is the number density in the absence of decays and Here, dN a /dE are the energy spectra per decay, with the sum in the numerator running over a = e, γ and the sum in the denominator over b = e, γ, ν. The T (a) (z , z, E) are transfer functions that describe the fraction of energy deposited at z for injection at energy E and redshift z ≤ z. Note that the exponential depletion due to decays has been incorporated into the deposition function f (z) as in Ref. [127].
In addition to the total efficiency of energy deposition described by f (z), similar functions f i (z) and corresponding transfer functions T (a) i can be defined for energy deposition into specific channels such as hydrogen ionization, helium ionization, heating of the cosmological medium, and photons with energies below 10.2 eV that free-stream [128]. Arrays of these transfer functions are collected in Ref. [129] and we apply them to our analysis. The effect of energy injection on the CMB power spectra is almost entirely due to the ionization of hydrogen, with much smaller contributions from helium ionization and Lyman-α photons [130]. As in Ref. [128], we define an ionization fraction χ ion by where f ion is the total efficiency of depositing energy into the ionization of hydrogen and helium.
A powerful method to connect the energy deposition function f (z) to bounds from CMB measurements was derived in Ref. [131] based on a principal component method. In this approach, orthogonal eigenfunctions with respect to the space of energy deposition histories are obtained from data (or data projections) that characterize the impact of energy injection on the CMB power spectra, with a marginalization over ΛCDM parameters included. The significance of the energy injection is then given by the orthogonal sum over the expansion coefficients of f (z) with respect to the eigenfunction basis weighted by their corresponding eigenvalues. Ordering the eigenfunctions by the sizes of their eigenvalues, only a small set of these principal components need to be included to get an excellent approximation of the full result [131] (which requires a intensive computation with tools such as CLASS [132] and CosmoMC [133] for each model of interest).
To estimate the CMB bounds from Planck [4] on dark vector decays, we make use of principal component functions and eigenvectors over the space of energy deposition histories described in Refs. [126,131,134] and collected in Ref. [129]. These were computed in Ref. [131] using the ionization fraction χ base (z) formulated in Refs. [51,135]. This was subsequently improved in Ref. [128] to account for losses to photons with energies below 10.2 eV that do not contribute to ionization. To implement this improvement, we follow the prescription of Refs. [136,137,138] by replacing f (z) in the analysis with With this prescription, we estimate the bounds from Planck and make projections for a cosmic variance limited (CVL) experiment with multipoles up to = 2500 using the principal component functions of Ref. [129]. In the case of Planck, these functions were derived based on a projection of the experimental sensitivity rather than data [129,131]. However, this projection turns out to be remarkably accurate at predicting the limits on s-wave dark matter annihilation derived from the Planck 2015 [139] and Planck 2018 [4] data sets. For decays, we also find that the projected limits for longer lifetimes τ V 10 13 s agree well with the principal components derived for dark matter decay with respect to the space of energy injection spectra based on Planck 2015 data [138].

Methods for the CMB Frequency Spectrum
Late-time energy injection can also distort the CMB frequency spectrum [48,49] from the nearly perfect blackbody that is observed [140]. For decays after the decoupling of double-Compton scattering at redshift z th 2.0 × 10 6 but before the decoupling of Compton scattering at z C 5.2 × 10 4 , the decay products equilibrate kinetically but generate an effective photon chemical potential [48,49]. Decays after Compton decoupling but before recombination at z rec 1090 produce a distortion that can be described by the Compton y parameter [48,49]. The current limits on µ and y from COBE/FIRAS are [140] µ < 9 × 10 −5 , |y| < 1.5 × 10 −5 , while the proposed PIXIE satellite is projected to have sensitivity to constrain [141] µ < 1 × 10 −8 , |y| < 2 × 10 −9 .
These limits can be used to constrain decays in the early universe.
An accurate expression for both types of distortions is [142,143,144,145] where the J i are window functions, ρ γ is the energy density of photons, and ∆ρ γ is the rate of electromagnetic energy injection per unit volume. For the window functions, we use "method C" of Ref. [145]: The energy injection profile for a decaying dark vector is ∆ρ γ = f em m V R, with R given by Eq. (34), and f em the electromagnetic energy fraction of the decay. For sub-GeV dark vectors we approximate this fraction by Note that in contrast to the BBN and CMB power spectrum bounds discussed above, this constraint on decaying particles is largely insentive to the energy spectrum of the decays.

CMB Bounds on Dark Vectors
In Fig. 4 we show the estimated bounds from deviations in the CMB power spectrum relative to Planck measurements on the pre-decay dark vector yield Y V = n V /s as a function of decay lifetime τ V and mass m V for the dark vector species DP (upper left), B −L (upper middle), L e − L µ (lower left), L e − L τ (lower middle), and L µ − L τ (lower right). We note that these bounds constrain values of m V Y V that are orders of magnitude lower than from BBN. However, the CMB constraints fall off quickly for lifetimes below τ V 10 13 s and BBN becomes more important. Let us also emphasize that the results shown for τ V 10 13 s have a very significant theoretical uncertainty that we do not attempt to estimate (see also Refs. [127,131,138]). The shape of the CMB bounds mirrors those from BBN, with direct photon and electron injection having a greater impact than muons or charged pions, and thus the emergence of muon decays and hadronic resonances with increasing dark vector mass produce distinct features in the figures.
In the left panel of Fig. 5 we show the approximately mass-independent bounds on the combination f em m V Y V from µ-and y-type distortions of the CMB frequency spectrum from COBE/FIRAS measurements as well as projected constraints from PIXIE. The µ-type distortions are dominant prior to the freezeout of double-Compton scattering, and y-type distortions become more important after that. To convert the curves in the left panel into bounds on a specific dark vector theory, it is only a matter of determining f em . In the right panel of Fig. 5 we plot f em as a function of mass m V ∈ [ MeV, GeV] for the five dark vector varieties discussed above. Relative to BBN electromagnetic effects, dark vector decays are constrained slightly more weakly by current data on the CMB frequency spectrum, but will become much more constraining with data from PIXIE. However, for τ V 10 13 s, energy injection primarily alters the CMB power spectra rather than the frequency spectrum [53,54], and the power spectra give the strongest limits.

Thermal Production by Freeze-In
A minimal population of dark vectors will be created by thermal N → 1 freeze-in reactions of the form (SM) N → V . 4 This production is described by where s is the entropy density, f V is the distribution function of the dark vector species, and C[f V ] is the standard collision term [1]. For underpopulated dark vectors with f V f eq V and neglecting Bose enhancement and Pauli blocking factors, detailed balance implies where Γ V is the total decay width in the V rest frame and the brackets refer to the thermal average of the inverse Lorentz factor.
The collision term for dark photon production was computed including thermal corrections and full Fermi-Dirac statistics in Ref. [78]. There it was found that full calculation via leptons (or perturbative quarks) is approximated very well by evaluating the collision term with a Maxwell-Boltzmann distribution and neglecting thermal corrections. Using a Maxwell-Boltzmann form for f eq V in Eq. (53), the integral can be evaluated analytically with the result [146,147] where K 1 is the modified Bessel function of the first kind.
To evaluate the thermal yield with this collision term, we follow Ref. [78] and convert time to x = m V /T and divide the production into temperatures above and below the QCD phase transition, assumed to occur at Λ QCD 157 MeV. The net yield is then with where x QCD = m V /Λ QCD andΓ V is the vector decay width into perturbative quark (and lepton) final states. For this, we consider up quarks with m u = 2.2 MeV, down quarks with m d = 4.7 MeV, and strange quarks with m s = 93 MeV [110]. With this approach, our results for the yield of the dark photon agree well with Ref. [78]. To the extent that the number of relativistic degrees of freedom is constant during production and is dominated by T < Λ QCD , the final dark vector yield is approximately [78] Y V 9 4π We use the full result of Eqs. (55,56,57) in the analysis to follow.

Cosmological Bounds on Thermal Dark Vector Parameters
By combining the thermal production yields with the BBN and CMB constraints derived previously, we can put limits on dark vectors produced thermally in the early universe. To compare the limits on various dark vector theories on an equal footing and to connect with previous work, it is convenient to define an effective coupling ef f by where α −1 137.036 is the usual fine-structure constant.
In Fig. 6  indicates the potential future sensitivity from a cosmic variance limited experiment up to = 2500, the solid yellow region shows the exclusion from COBE/FIRAS from distortions in the CMB frequency spectrum, and the dashed yellow contour shows the projected reach of PIXIE for such distortions. No significant exclusion is found for the L µ −L τ vector. Let us also emphasize that these exclusions are conservative in the sense that they only rely on the assumption of an early universe dominated by SM radiation at temperature T m V . The exclusions would be even stronger if there were other production sources for the dark vector such as particle decay or annihilation [76,77,94].
The cosmological limits we find for the dark photon are qualitatively similar to those derived in Refs. [78,79]. Compared to these works, our analysis includes several updates that affect the detailed quantitative results for the dark photon and our study extends to other types of dark vectors. The most significant update is our treatment of electromagnetic energy injection and the resulting photon cascade spectrum relevant for BBN. In Ref. [78] the universal photon spectrum was used, while in Ref. [79] the modification of the photon spectrum from e + e − injection due to the Thomson limit of inverse Compton scattering was included leading to much weaker exclusions from BBN photodissociation effects. Our result lies between these two exclusions -we confirm the suppression of the photon cascade spectrum from lower-energy e + e − injection pointed out in Ref. [79] but we also find a significant (and sometimes dominant) additional contribution from final-state photon radiation in these decays that bring our result closer to Ref. [78]. In our treatment of deviations in the CMB power spectrum, we also use the updated treatment of the ionization fraction presented in Refs. [128,130,136,138].
The parameter bounds presented in Fig. 6 are the strongest constraints on the species of (thermal) sub-GeV dark vectors studied here with 10 −12 and m V GeV. For ∼ 10 −11 and m V ∼ GeV, Refs. [78,79] also found disjoint BBN exclusions on the dark photon from neutron-proton converstions induced by charged pion and kaon injection, while for ∼ 10 −13 and m V 2 GeV they obtained BBN bounds mainly from the direct injection of neutrons from vector decays. These BBN bounds from hadronic injection are expected to carry over similarly to the B−L vector for masses above m V 1 GeV. Dark vectors can also alter the effective number of relativistic neutrino species N ef f by thermalizing light singlet neutrinos [82] or by injecting energy preferentially into either electromagnetic species or neutrinos after neutrino decoupling [21,22,23,25,124], with the corresponding bounds based on freeze-in production of a dark photon extending down to 5 × 10 −10 for m V 8 MeV. The emission of gamma-rays by the decays of dark vectors that would have been produced in supernova SN1987A [148,149] was shown to rule out couplings down to ef f ∼ 10 −12 for masses between 1-100 MeV [149], independent of the cosmological production mechanism of the dark vector. Taken together, our cosmological exclusions of dark vectors (assuming freeze-in production) are disjoint from these others.
Beyond the limits shown above, we have also estimated the gamma-ray signals of dark vector decays with lifetimes longer than the age of the universe. Translating these signals to the equivalent from the decay of a long-lived particle making up all the DM, they correspond to effective DM lifetimes greater than τ ef f 10 30 s, well beyond the current limits of τ 10 25 s for DM → e + e − in this mass region [150].

Conclusions
In this work we have investigated the cosmological bounds on a range of sub-GeV dark vectors based on the electromagnetic energy injected by their decays in the early universe. This injection can alter the SM+ΛCDM predictions for the light element abundances from BBN and the power and frequency spectra of the CMB. Our work expands on previous studies of cosmological bounds on kinetically-mixed dark photons [78,79], and extends them to four other dark vector species: B −L, L e −L µ , L e −L τ , and L µ −L τ .
The cosmological bounds we derive are based on electromagnetic energy injection. We take into account deviations from universality in electromagnetic effects on BBN from lowerenergy injection by computing explicitly the development of the electromagnetic cascade photon spectrum. In doing so, we take into consideration the detailed photon and electron energy injection spectra including radiative effects such as FSR. For the specific case of a dark photon created through thermal freeze-in, the BBN bounds we find are slightly weaker than Ref. [78] and somewhat stronger than Ref. [79]. Our results for CMB exclusions on the dark photon appear to be in agreement with Refs. [78,79], although we apply a more detailed energy injection spectrum and we update the effective ionization fraction following Ref. [136,137,138]. We also extend these methods to other dark vector species.
Relative to direct laboratory tests of dark vectors, the bounds we find from cosmology are able to constrain much smaller couplings to SM matter, with CMB bounds extending all the way to ef f ∼ 10 −18 . For the dark photon, it is challenging to obtain kinetic mixings this small [151], but they can arise when there is an approximate charge conjugation (C) invariance within the dark sector [152]. For the other direct coupling vectors, the bounds on the gauge couplings we find from cosmology are strong enough to be potentially relevant for certain superstring swampland considerations, corresponding to the set of consistent low-energy effective theories that cannot be embedded into a consistent theory of quantum gravity [153,154,155]. In particular, m V ∼ GeV and ef f = α V /α ∼ 10 −18 approach the parameter range covered by the weak gravity conjecture proposed in Ref. [154]. 1607611, for their hospitality while this work was being completed. This work is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), with DEM supported in part by Discovery Grants and LF by a CGS D scholarship. TRIUMF receives federal funding via a contribution agreement with the National Research Council of Canada.

A Appendix: Departures from EM Universality in BBN
Electromagnetic energy injection well above the weak scale has been shown to produce a universal photon cascade spectrum that to an excellent approximation depends only on the total amount of electromagnetic energy injected [37,91]. Correspondingly, the limits from BBN on very massive particle decays depend on the total electromagnetic injection (and decay lifetime) but not the detailed decay type or energy spectrum. In Refs. [87,92,93,94] it was shown that this feature breaks down for sub-GeV injections of photons or electrons, and thus it becomes necessary to track the type of electromagnetic injection (electrons or photons) and its energy spectrum.
In Fig. 7 we show the combined bounds from BBN on the mass times yield m V Y V of a massive particle V decaying to either V → γγ (left) or V → e + e − (right) as a function of the lifetime τ V and the particle injection energy E V = m V /2 from the two-body decay. This figure illustrates universality in the upper left regions of the panels with nearly identical bounds for photon and electron injection for given values of τ V and m V . However, universality is seen to break down badly in the lower right regions of the panels, with significantly different BBN bounds depending on the type of electromagnetic injection within these regions.
The breakdown of universality at lower energies can be understood in terms of the energy and temperature dependences of the photon-photon pair production (4P) and inverse Compton (IC). These two processes play a crucial role in the formation of the electromagnetic cascade, and for high-energy injection they are the dominant processes that seed the rest of the cascade and drive universality. As the injection energy falls below the weak scale, the typical interaction energy can fall below important thresholds for these reactions leading to deviations from universality.
Energetic photons injected into the cosmological plasma when the temperature is below the MeV scale interact primarily with the plasma through the 4P process. This process is very efficient relative to other processes (and the Hubble rate) provided the photon energy is greater than a critical energy E c , given by [37,91] This corresponds to the energy threshold for e + e − pair production on background photons extended by a Boltzmann tail. For photons with energies above the threshold E c , the 4P process is extremely efficient and these photons are quickly reprocessed into electrons at energies below E c . In contrast, photons with energy below E c are only reprocessed by other, less efficient process. The black line in Fig. 7 shows the boundary at which E V = E c (τ V ), with injection energies above the 4P critical energy to the upper left of this line.
Electrons injected into the cosmological plasma below MeV temperatures with energy E e T interact primarily through IC. This upscatters a background photon and reduces the electron energy. The energy spectrum of the upscattered photons depends sensitively on the combination [115] y e ≡ E e T m 2 e .
For y e 1, the scattering occurs in the Klein-Nishina regime with the outgoing photon energy E IC γ typically of the same order as the initial electron energy, E IC γ ∼ E e . This allows electrons to efficiently populate the photon cascade spectrum. In contrast, when y e 1 the scattering enters the Thomson regime in which the upscattered photon energy is only a small fraction of the intial electron energy, and typically less than [115] Thus, lower-energy electrons are not able to produce photons through IC with enough energy to photodissociate light nuclei. The white line in Fig. 7 shows the threshold at which the IC photon energy produced by an electron with energy E e = E V lies below E IC γ ≤ 2 MeV, corresponding to the dissociation threshold for deuterium (and the lowest threshold relevant to our analysis).
Comparing the BBN exclusion contours in Fig. 7 to the 4P (black) and IC (white) threshold boundaries, we see that universality holds to a very good approximation for injection energies above these lines. Below these lines, significant deviations from universality are manifest. For photon injection with E V ∼ 20-50 MeV < E c , the photon spectrum is more weakly suppressed and these photons are able to destroy 4 He efficiently. As the photon injection energy falls below E V < 20 MeV, the photons are no longer able to dissociate 4 He but can still destroy D and 3 He down to near 2 MeV. These features appear clearly in the left panel. Electron injection in the Thomson regime, corresponding to the region below the white line in the right panel, leads to much weaker limits from BBN. In this region, the dominant contribution to the photon cascade spectrum is typically hard finalstate radiation (FSR) from the primary decay with a relative fraction suppressed by a factor of α [87].