A Stimulating Explanation of the Extragalactic Radio Background

Despite an intense theoretical and experimental effort over the past decade, observations of the extragalactic radio background at multiple frequencies below 10 GHz are not understood in terms of known radio sources, and may represent a sign of new physics. In this Letter we identify a new class of dark sector models with feebly interacting particles, where dark photons oscillate into ordinary photons that contribute to the radio background. Our scenario can explain both the magnitude and the spectral index of the radio background, while being consistent with other cosmological and astrophysical constraints. These models predict new relativistic degrees of freedom and spectral distortions of the cosmic microwave background, which could be detected in the next generation of experiments.

Introduction.-Thecosmic microwave background (CMB) between 60-600 GHz is one of the most wellstudied electromagnetic signals in science.In this frequency range, the CMB monopole dominates over astrophysical backgrounds, and is consistent with a blackbody distribution, with distortions limited to less than 1 part in 10 4 [1].Spatial fluctuations in the brightness of the CMB are on the level of 1 part in 10 5 , and the power spectrum of these fluctuations forms a key pillar of modern cosmology [2].These facts constitute strong evidence for the primordial origins of the CMB.
Explaining T exc with known sources of radio emission has thus far been surprisingly difficult [10].Almost all attempts to do so have relied on synchrotron emission, since synchrotron sources typically produce −3 β −2.5 [11,12].Extragalactic radio synchrotron sources such as active galactic nuclei and star-forming galaxies can produce β ≈ −2.7 [13][14][15][16][17], and the source counts distribution is thought to be relatively well-understood.However, numerous studies have estimated the emission from extragalactic radio sources to be 3-10 times smaller than T exc [16,[18][19][20][21][22] between 100 MHz-10 GHz.Electrons that are reaccelerated during cluster mergers have been proposed as a significant contribution to T exc [23], albeit under optimistic assumptions; this mechanism also produces a softer power law than is required.
An alternative explanation of T exc is contamination from unmodeled Galactic sources of radio emission, which needs to be subtracted from radio data to obtain the extragalactic component [24][25][26].However, proposed additional sources of Galactic radio emission are strongly constrained [26][27][28].More sophisticated modeling using cosmic ray propagation models disfavors the possibility of significant Galactic contamination [20].
A more exotic class of solutions involves radio emission from hypothetical early structures [29], including black holes [30][31][32] and star-forming galaxies [33] at high redshifts.These solutions typically require large, persistent magnetic fields [32], whose origin and impact on inverse Compton cooling are debated [34].
The ERB has also been studied in conjunction with 21-cm experiments, including EDGES [35], LEDA [36] and LOFAR [37].These measurements are in tension with a T exc that is fully produced at high redshifts.In particular, the LEDA and LOFAR results constrain any cosmological radio background with T exc ∝ ω −2.6 at ∼ 10% (13.2 < z < 27.4) [29] and ∼ 46% (z = 9.1) [38] of the present-day radio excess respectively.
Another aspect of the ERB that makes astrophysical explanations difficult is the spatial smoothness of the emission [39], deduced from measurements of the anisotropy of the radio sky [39][40][41][42][43][44].These measurements find that fluctuations in T exc , ∆T /T exc , are ∼ 10 −2 across a range of angular scales and frequencies; this is smoother than expected if radio emission is correlated with the present-day dark matter distribution [39].
New physics explanations proposed thus far have primarily focused on synchrotron emission from DM annihilation and decay [45][46][47][48].However, these models also run into similar issues: they may result in nonsmooth emission; underproduce T exc , unless the magnetic fields responsible for synchrotron production have unusual or unlikely properties [47,48]; require a new, large population of faint sources [21], or require a large portion of the isotropic, extragalactic gamma-ray background to come from DM annihilation [46], which is disfavored [49,50].
In this Letter, we show that a simple class of experimentally viable new-physics models can explain the amplitude, power-law dependence and smoothness of T exc .These models rely on three basic ingredients: 1) a particle decaying into dark photons A ; 2) the presence of a thermal bath of A which stimulates this decay; and 3) A resonantly oscillating into radio photons.This class of models leads to T exc ∝ ω −5/2 , close to the observed power-law dependence.Relatively low-z resonant oscillations as well as dn A /dω ∝ ω −1/2 prior to oscillations are both crucial elements of such models; we have not been able to identify viable alternatives.One possibility is the decay of DM particles to final states charged under A , and A appearing as "dark bremsstrahlung".While having a spectrum similarly enhanced in the infrared, the DM has to be light and millicharged, with an effective electromagnetic charge exceeding current experimental bounds [51].Another alternative is DM decay into a dark photon A , which then resonantly converts twice via A → A → γ, with some dark states charged under a dark photon A facilitating the A → A transition.However, two resonant conversions makes it difficult to produce a sufficiently large T exc while ensuring that T exc ∝ ω −5/2 over nearly three decades of frequency.
The remainder of this Letter is organized as follows.We begin by introducing the class of models and deriving the expected T exc from them.Next, we discuss our fit of a particular example model to the radio data, and relevant experimental constraints.We then move on to discuss the anisotropy of the ERB produced by the model.We conclude with future prospects for confirmation of this model.We adopt Planck 2018 cosmology throughout [2].Model.-Aparticle physics model that has the following three features can generate an ERB with T exc ∝ ω −5/2 : 1. a cold component of DM a which is stable on cosmological timescales, but undergoes a two-body decay with lifetime τ vac in vacuum (for simplicity, we take a to be all of DM); 2. one of the daughter particles of the decay is a dark photon, A , with an existing blackbody distribution, characterized by a temperature T (z) = T 0 (1 + z), where T 0 is its temperature today; and 3. the A has mass 10 −15 eV m A 10 −9 eV and is emitted relativistically with energy ω A .A mixes with the Standard Model photon γ with kinetic mixing parameter .
These three features are summarized in Fig. 1.The existence of the thermal population of A enhances the decay rate of a due to Bose enhancement [52,53], leading to a redshift-dependent effective decay lifetime τ (z), where with f BB A = (e ω A /T − 1) −1 being the blackbody occupation number of A with energy ω A .
There are several parameters that depend on the specifics of the model.The value of n depends on the occupation number of the other daughter particle.In addition, a can decay into α = 1 or 2 dark photons.Finally, ω A depends on the kinematics of the decay of a, with ω A = m a /2 if it decays into a pair of relativistic final states.In App.B, we describe a concrete particle physics model that realizes all three features.To guarantee T exc ∝ ω −5/2 , processes such as inverse decays of A into a cannot significantly distort either the power-law index of the A spectrum produced by the decay, or the blackbody distribution of A across all relevant frequencies.In accordance with this fiducial model, we fix n = 2, α = 1, and ω A = m a /2 for all results that are specific to it, although we emphasize that this is only one example, out of potentially many, that we have studied in detail.
The small kinetic mixing between A and γ enables resonant conversion between the two particles whenever m 2 A = m 2 γ [54,55], where m 2 γ is the effective photon plasma mass squared.This quantity is proportional to the free electron number density, n e .The converted photons ultimately form the present-day T exc .We calculate the sky-averaged conversion probability per redshift d P A →γ /dz taking into account inhomogeneities using the formalism developed in Refs.[55,56] (see also Refs.[57][58][59]).Radio background production.-Theparticle a decays throughout cosmic history, producing a number density of A per unit energy given by (see App.A for a complete derivation, which follows from Ref. [60]) where 1 + z ≡ ω A (1 + z)/ω is the redshift at which a daughter A with frequency ω at redshift z was produced, H is the Hubble parameter, ρ a (z) is the DM energy density, and Θ is the Heaviside step function.The total spectrum of A is then a sum of this decay spectrum and the blackbody spectrum of the thermal A distribution, which is subdominant in the energy range of interest.
In comoving coordinates, the spectrum of photons produced is obtained by integrating the spectrum of dark photons times d P A →γ /dz, i.e., where x ≡ ω/[T 0 (1 + z)].For the range of m A that is of interest, resonant conversions occur only after recombination.In the range 3 × 10 −4 x 0.2, which are relevant for T exc , these photons only evolve via redshifting, with the baryonic fluid being essentially transparent to them [53].
Substituting Eq. (3) into Eq.( 4), we find This gives dn γ /dx ∝ x −3/2 , or T exc ∝ ω −5/2 , the desired frequency dependence.Note that z must occur during matter domination in order for H(z ) ∝ x −3/2 , while d P A →γ /dz ∝ x −1 is derived in Refs.[54,55].τ (z ) ∝ x follows from Eq. (2) when T ω A .Eq. ( 5) is the main result of this paper.For our fiducial model, we find the following approximate parametric dependence of T exc at z = 0: where P GHz is P A →γ for ν = 1 GHz today, which the FIRAS measurement of the CMB energy spectrum limits to be P GHz 10 −2 [54,56].Fit and constraints.-Wecan now perform a fit of T exc predicted by our model to the measured data from Ref. [3].There are six parameters in our fit: five parameters {m a , m A , τ vac , T 0 , } from our new physics model, and T 0 .Fitting the radio data places some requirements on the model parameters.For T exc to be present across all data points, we require m a /2 2π × 10 GHz.On the other hand, photons at 22 MHz must originate from decays of a during matter domination to satisfy T exc ∝ ω −5/2 , leading to the requirement that m a /2 2π(1 + z eq ) × 22 MHz, where z eq is the redshift of matter-radiation equality.Together, this means 8 × 10 −5 eV m a 6 × 10 −4 eV.Resonant conversion must predominantly occur after z for 10 GHz photons so that the full power law extends to at least that frequency; the requirements on m a show that z 6, which favors 10 −14 eV m A 3 × 10 −13 eV.Finally, T must be sufficiently large for f BB A 1 and τ (z ) ∝ x for the 10 GHz data points; this imposes T 0 /T 0 0.06.
Several important experimental constraints exist on the class of models under consideration.First, T 0 0.4T 0 , in order to avoid introducing excessive effective relativistic degrees of freedom [2].The kinetic mixing between A and γ also leads to spectral distortions due to γ → A resonant conversions; the FIRAS measurement of the CMB spectrum [1] leads to constraints on the order of < 10 −7 -10 −5 in the m A range of interest [55,56].
There are also limits on the DM decay lifetime obtained from the CMB power spectrum, large scale structure and the Milky Way satellite population [63][64][65][66][67][68]; however, these results assume a constant decay width.Conservatively requiring the total energy density of a that has decayed away with stimulated decay from the A thermal population by the presentday to be less than 2.2%, we find τ vac > 10 19 s × max 2.1, 6.3 10 −4 eV/ω A (T 0 /0.2T 0 ) (see App. C for details).We also note that the decaying particle a may be a subcomponent f DM of DM, which evades the lifetime bound altogether if f DM 2.2% [67].
21-cm power spectrum measurements from both LEDA and LOFAR are both in tension with T exc being fully produced at high redshifts (z 9.1), placing strong constraints on models that produce T exc primordially.In our model, however, A → γ resonant conversions occur predominantly at z 6.
We explore the posterior on the model parameters using nested sampling [69][70][71] implemented in dynesty [72].Our priors, described in detail in App.E, are constructed such that they have zero probability density in parameter regions incompatible with 1) the FIRAS spectral distortion limits of Refs.[55,56] and 2) the DM lifetime limit discussed above.Priors on T 0 and T 0 are also chosen to account for the N eff and FIRAS constraints on these parameters, respectively.The posterior on the excess temperature is shown on the left in Fig. 2. In computing T exc from our model, we include an irreducible contribution from unresolved extragalactic radio sources of T eg = 0.23 K(ν/GHz) −2.7 [16,39], which is at least 3 times smaller than the measured brightness temperature.This parametrization of T eg is in excellent agreement with other independent estimates [73].Although this power-law expression strictly applies in the range 100 MHz-10 GHz, we estimate the contribution outside of this frequency range by extrapolation.The marginal and joint posterior distributions are shown on the right in Fig. 2. The marginal posterior medians correspond to m a 2 × 10 −4 eV, m A 2.5 × 10 −14 eV, τ 1.3 × 10 21 s, T 0 0.22 T 0 and 10 −7 .We find an excellent fit to the data over a wide range of allowed model parameters.In App.F, we show the extended corner plot for posterior distributions of all parameters of interest, along with other systematic variations; these do not qualitatively change our main result.
Smoothness.-Upper limits on the anisotropy of the ERB have been reported for 4 × 10 3 6 × 10 4 in the 4-8 GHz range by VLA [40,41] and ATCA [42], while actual measurements of the power spectrum have been made by LOFAR [44] and TGSS [43] at ∼140 MHz for 10 2 10 4 .Such measurements can be challenging: observations at ∼140 MHz disagree with each other by a factor of 3, and face issues such as galactic foreground contamination and calibration errors [44].To assess the smoothness of T exc in our model, we take the 4-8 GHz upper limits and results from LOFAR-which reports a lower power-as approximate upper limits, noting that astrophysical sources can contribute more power [44], and that these observations are expected to improve.
There are two possible contributions to the ERB anisotropy produced by our model that are essentially independent: 1) decay anisotropy, due to DM density correlations from the point at which a decays, and 2) conversion anisotropy, due to correlations in free electron density fluctuations δ e during A → γ conversions, since m 2 γ ∝ n e .The decay anisotropy was found to exceed the radio anisotropy power spectrum in the 4-8 GHz range, unless the decay that produces these photons happens at z 5 [39,47].This is easily satisfied in the range of parameters providing a good fit.
We compute the conversion anisotropy power spectrum by first writing down the two-point correlation function  of the conversion probability of A → γ in two different directions in the sky.As such, it depends on the twopoint function of δ e , which we model as either a Gaussian or a log-normal random field.The anisotropy power spectrum can then be obtained by performing a spherical harmonics decomposition, and using the Limber approximation [74][75][76] following the method outlined in Ref. [77].Further details of our calculation can be found in App.D. Fig. 3 shows the predicted anisotropy power spectrum in units of K 2 , defined using the same conventions as in the CMB power spectrum [2], for normally and lognormally distributed δ e , with the result between these two choices shaded in blue.This can be taken as an estimate for the uncertainty associated with these distributions.This power should be compared to the LO-FAR data, as well as the upper limits from VLA and ATCA at their respective frequencies ν, which have been rescaled by (140 MHz/ν) 2β , with β = −2.6;this assumes that ∆T /T exc is independent of frequency.Our calculated anisotropy power spectrum for a fiducial choice of m A = 3 × 10 −14 eV and ω A = 2 × 10 −4 eV lies below the 4-8 GHz measurements for both choices of the two-point PDF of δ e , while only the log-normal twopoint PDF exceeds the low-measurements by LOFAR.However, significant scatter exists between adjacent frequency bands in the data for 10 3 [44].Producing a sufficiently smooth radio background is thus highly plausible in our model.
Conclusion.-In this Letter, we have provided a new physics explanation for the ERB, which is brighter than expected from known astrophysical sources.Our model consists of a DM candidate which decays into dark photons in the presence of a thermal bath of the latter; the dark photons then resonantly convert to ordinary photons, producing the ERB.The anisotropy of the signal is expected to be relatively smooth, consistent with measurements of the anisotropy of the radio background [39,44].Future experiments may corroborate the predictions of our model.PIXIE [61] may be sensitive to spectral distortions expected from our model, and is expected to almost fully cover the 95% region of our posterior distribution in the -m A plane.The thermal population of A may also lead to a value of N eff that is detectable in future CMB experiments such as CMB-S4 [62].The potential reach of both experiments are shown in Fig. 2.
where α = 2 if both daughter particles are A and α = 1 otherwise, ρ a (z) is the energy density of a at redshift z, and m a is the mass of a.With this, we can find the contribution to the dark photon number density at redshift z per unit frequency ω, due to a decay occuring at redshift z : where the redshift factors rescale the number density of particles produced at redshift z down to redshift z, and the Dirac delta function enforces the fact that the decay produces dark photons with energy ω A .Integrating over z gives where is the redshift at which the decay of a produced the dark photon at frequency ω, with ω ≤ ω A .We have also assumed that ρ a ∝ (1 + z) 3 .This expression agrees with Ref. [60], if we simply set Γ(z ) → Γ vac .In terms of comoving number density and x ≡ ω/T 0 where T 0 is the CMB temperature today, we find which is almost constant in redshift, aside from the step function.ρ a,0 is the dark matter energy density today.Eq. (A1) is the expression we use in the Letter.10 −9 eV, and is emitted relativistically with energy ω A .A mixes with the Standard Model (SM) photon γ with kinetic mixing parameter .
These conditions appear simple to meet; for example, the model proposed in Ref. [60] where the DM is a dark axion decaying into two identical dark photons which oscillates into SM photons, could be supplemented by an additional thermal population of dark photons.However, there are several challenges to building a successful model.In order to produce the right power law in the ERB, the A power-law spectrum must maintain dn A /dω ∝ ω −1/2 between 3 × 10 −4 x 0.2, where x ≡ ω/T CMB , from the redshift of production z (x) until resonant conversions A → γ are mostly complete.Furthermore, the A blackbody spectrum also cannot be significantly distorted; specifically, A particles from the bath at each value of x must be described by a blackbody spectrum at z (x).In the model of Ref. [60], both of these requirements can be violated by inverse decays A A → a or other A scattering processes.Moreover, it is difficult to maintain thermal equilibrium in the blackbody spectrum through scattering with another particle in the dark sector bath, since such scattering processes likely distort the power-law spectrum significantly.For typical parameters required to produce the full T exc , scattering between a low-energy A in the power-law spectrum and a blackbody A is too rapid to guarantee that the A blackbody spectrum remains undistorted at all relevant times.While it is possible that even with significant distortion to the blackbody spectrum, a reasonable fit to the ERB can still result, checking this possibility would require us to integrate the Boltzmann equation over a wide range of A frequencies.In this paper, we avoid this computational challenge by building a slightly different fiducial model without significant distortion.
Before discussing the details of our fiducial model, we will first discuss how to check the size of any A spectral distortion due to a single process.We begin by writing down the Boltzmann equation with the relevant process contributing to the collision term.Neglecting for simplicity the existence of entropy dumps in the Standard Model so that x for any A particle is a constant, the Boltzmann equation governing the occupation number f A can be written where C[f A ] is the collision integral, which we will define for particular processes below.Here, we adopt the convention that f A is related to the number density via where the factor of 3 accounts for the degeneracy in the spin state of A .To determine if any process causes a significant distortion to the spectrum of blackbody A or low-energy A from the decay of a, we divide Eq. (B1) by f A and integrate with respect to time or, equivalently, redshift, to obtain the change in log f A , which gives a measure of the fraction of A at that value of x that undergoes a scattering process.Our criterion for a small distortion to the A spectra is therefore where z min and z max are the lowest and highest redshifts respectively for which the collision term is important.Since we are only interested in producing the radio excess over a finite frequency range, we only need to check that distortions are small over a limited range of x.The radio frequency data points with an excess temperature over the CMB temperature span the frequency range 22 MHz-10 GHz, which we can cover by considering 3 × 10 −4 < x < 0.2.
The power-law component at fixed x is produced by decays of a at a redshift 1 + z ≡ ω A /(xT 0 ), where ω A is the energy of the emitted A at the point of decay (for a decaying into two massless particles, this is simply ω A = m a /2), and T 0 is the blackbody CMB temperature today.After it is produced, it must stay in this power law without undergoing significant distortions until the present day. 1 For the blackbody component to effectively produce the stimulated emission that we need, on the other hand, the blackbody spectrum must not be significantly distorted prior to redshift z (x); distortions after this redshift are unimportant.This means that for the power-law spectrum, for each value of x, we need to consider z min = 0 and z max = z (x), while for the blackbody spectrum, z min = z (x).

Fiducial model
We will now describe our fiducial particle physics model, which has all the three properties required to produce the radio background laid out in the previous section.We will show that both the low-energy and blackbody distribution of the dark photons do not undergo any significant distortion.
Our fiducial model is a modified version of the model proposed in Ref. [60].The important particles in this model are: i) an axion-like dark matter a with mass m a and decay constant f a , and ii) two dark photons, which we label A ψ and A .The three requirements for producing the radio excess as discussed in the main body of the paper are satisfied as follows: 1. the DM decays through the process a → A ψ A osc ; 2. a blackbody distribution of A osc described by a temperature T is also present, leading to stimulated decay of a, and 3. the dark photon A osc possesses a small kinetic mixing term with the SM photon, and has a mass 10 −15 eV m osc 10 −9 eV, so that it undergoes resonant conversions into the SM photon after recombination.
The terms in the dark sector Lagrangian that are relevant to us are therefore where L K contains the other kinetic terms of the dark photons, including a mass of m osc for A osc .Note that a discrete symmetry under which A ψ → −A ψ and a → −a explicitly forbids decays of a to a pair of dark photons of the same species, as well as mixing between A ψ and the SM photon, naturally leading to the Lagrangian shown above.By allowing a to decay into two different dark photons, we can now introduce a fermion ψ that scatters rapidly with A ψ , to ensure that A ψ always equilibrates into a thermal distribution with temperature T , without destroying the low-energy spectrum of A osc . 2 With this modification, we overcome the main difficulty faced by the model in Ref. [60]: there are now no low-energy A ψ to interact with the blackbody distribution of A osc , which was the main source of distortion for the blackbody spectrum. 3 We are now ready to check that distortions to both the low-energy A osc power-law spectrum and the A osc blackbody are small.The most important process, which enters at order 1/f 2 a , is inverse decays, A ψ A osc → a, which leads to the largest distortion.The collision integral is given by We use k, k and ω for incoming four-momentum, three-momentum, and energy; likewise, we have p, p and E for outgoing states.The subscript 'osc' represents quantities associated with A osc .There is no Bose enhancement in the outgoing state, since we can treat the entire population of a as having zero momentum.This also allows us to neglect the contribution from the backward reaction to the collision integral.The squared matrix element for this process (summed over initial and final states) is Technically, the dark photons can become significantly distorted after all resonant oscillations into photons have been completed, but for simplicity we set this more stringent requirement. 2For simplicity, we assume that the entire dark sector is described by a common temperature T . 3While this completes the basic description of our model, there are some requirements on ψ to keep this model tractable and avoid tracking the spectra of all of these particles as a function of time.First, Compton ψA ψ → ψA ψ and double-Compton ψA ψ A ψ ↔ ψA ψ scattering must be sufficiently rapid to ensure that A ψ is always described by a simple, blackbody distribution.This is easily satisfied, as long as ψ has a large coupling to A ψ , and is sufficiently abundant.Second, the process aAosc ↔ ψψ must not lead to significant distortion of Aosc.This is hard to determine without tracking the full evolution of the spectrum of Aosc, but it is possible to avoid this entirely by choosing a sufficiently massive ψ to kinematically forbid aAosc → ψψ, and making ψ asymmetric.We find that m ψ ∼ 30 eV with a coupling to A ψ of α D = 1, with an asymmetric number density today of n ψ,0 = 2.5 × 10 −3 cm −3 can provide sufficiently efficient scattering with A ψ .This value of m ψ is large enough to kinematically forbid aAosc → ψψ.Introducing a similar fermion χ with mχ ∼ m ψ that couples to Aosc instead guarantees that no distortion can occur when T > m ψ , mχ, while allowing distortions to build up once T < m ψ , mχ and ψ, χ has frozen out.This also prevents ψψ → aAosc from being significant at T < m ψ , mχ, since ψ annihilates away at T ∼ m ψ .
while the occupation number of A ψ is given by the blackbody occupation number, These simple expressions allow us to perform the collision integral relatively easily, to obtain where ω osc = xT CMB .With this expression, we can now assess the total distortion to the blackbody and power-law distributions of A osc , and check that for typical model parameters, these distortions are small, based on the criterion set out in Eq. (B2).For the blackbody distribution, we want to perform the integral in Eq. (B2) starting from z min = z , allowing us to approximate the collision term as Integrating this gives where γ ≡ T 0 /T 0 , a redshift invariant quantity.This indicates that the total distortion to the blackbody distribution is small.For the low-energy spectrum, we integrate Eq. (B2) from z min = 0 and z max = z .The integral can be approximately done in two parts: √ x, where we can expand the expression in Eq. (B4) using log(1 − α) −α, .

(B7)
This distortion appears to be somewhat large for the fiducial values shown here, but the relative uncertainty on the data point at x = 3.8 × 10 −4 or 22 MHz is ∼25%, and so is enough to absorb the distortion obtained here.Moreover, the posterior distribution from our fits include larger values of m a and τ vac , which decreases the size of the distortion.We therefore conclude that inverse decays do not significantly distort either component of the A osc spectrum. 4 Our proposed model therefore satisfies the requirement that all distortions to the A osc power law and blackbody spectrum are small, and is a viable model for explaining the ERB.We stress that many of the required conditions are invoked in order to simplify the analysis and to allow an unambiguous prediction of the ERB in this model.Some of the conditions on the distortions, for example, may be relaxed under a complete analysis, which would include numerically solving the Boltzmann equations for each mode.Cosmological constraints on the decay lifetime of dark matter have been obtained for decays without stimulated emission [63][64][65][66][67][68].These results show that if a subcomponent f dcdm decays between recombination and today, Planck 2018 and baryon acoustic oscillation (BAO) data limit f dcdm < 0.0216 [67] at the 95% confidence level.This result is consistent with the uncertainty in the dark matter energy density reported by Planck 2018 [2].From this, we can deduce a limit on the decay lifetime with stimulated emission, by adopting a limit of f lim = 0.0216 for the proportion of dark matter that decays away by the present day.For a stimulated emission lifetime given by τ (z) = τ vec (1+2f BB A ) −1 , which is the expression we obtain with our fiducial model, we require 5 .
4 Although elastic scattering processes should be subleading at order 1/f 4 a , they can, in fact, be very rapid for two reasons.First, processes like aAosc → aAosc receive a large Bose enhancement from the existing population of Aosc, which significantly increases the scattering rate.Secondly, the matrix element of aAosc → aAosc contains a divergence.Such divergences are commonplace even in Standard Model processes, such as e − Z → e − Z [94].In the context of cosmological fluids, these divergences are regulated by thermal corrections to the self-energy of A ψ , which impart an imaginary contribution to the mass of A ψ [94].We have checked that even including Bose enhancement and regulating this divergence with loop contributions from a thermal distribution of ψ [95], all elastic scattering processes are indeed subdominant to the leading inverse decay process discussed here. 5This corresponds to an approximation of the effect of stimulation, we leave a more complete treatment for future work.
Numerically, we find which leads to a limit on the vacuum decay lifetime of This limit is roughly an order of magnitude stronger than the lifetime limits on DM decays without stimulated emission for our fiducial choice of parameters [67].

Constraints on relativistic degrees of freedom
We complete our discussion of the particle physics model by discussing its effect on N eff , and relevant constraints.A osc , A ψ contribute to N eff as relativistic degrees of freedom.In addition, other fermions that couple to these bosons can contribute as well if they are sufficiently light.Assuming the existence of two such Dirac fermions on top of A osc and A ψ before recombination, the total energy density of these particles divided by the energy density of a neutrino is This is to be compared with the Planck measurement of N eff = 2.99 ± 0.17 [2], which limits ∆N eff < 0.34 at the 95% confidence level.Including the irreducible contribution from the thermal distribution of just one dark photon, which must be present in the class of models considered in this work, gives ∆N eff ≈ 0.01(γ/0.2) 4 .

Appendix D: Anisotropy due to dark photon conversions
In this appendix, we obtain the conversion anisotropy power spectrum, which originates from variations in electron density fluctuations along two lines-of-sight separated in the sky by some angle.Our discussion follows Ref. [77] closely, but with some novel and peculiar features originating from the resonance conversion process; for clarity, we repeat many of the same calculations presented in that reference.
For dark photon conversions, the observed temperature in any direction in the sky n at a fixed frequency is directly proportional to the total probability of conversion P (n) of dark photons into photons.We begin by defining the fluctuation of the conversion probability in some particular direction in the sky n, defined explicitly by δP (n) ≡ P (n) − P , where P is the sky-averaged conversion probability.In any particular direction n, we have [55,56] where m 2 γ ( r, z) is the effective plasma mass of photons along the particular line-of-sight, and ω 0 is the present-day, observed angular frequency of the photons.r ≡ nχ(z), where χ(z) is the comoving distance traveled by light between z and the present day.ω 0 (1 + z ) is the energy of the daughter particle from the DM decay; photons that have energy ω 0 today were produced by decays at z , which is the maximum redshift we integrate up to.
We can compute the sky-averaged conversion probability by integrating over the one-point probability density function (PDF) f 1 (m 2 γ ; z) of m 2 γ , to obtain [55,56] We now decompose the observed conversion probability over the whole sky into spherical harmonics Y m (n), with coefficients a m given by The anisotropy power spectrum C for the conversion probability, defined as C ≡ a * m a m , can then be computed as where • • • in the integral should be interpreted as an all-sky average.Note that defined in the following manner, C is dimensionless, since it describes fluctuations in conversion probability.
To make further progress, we define the quantity a two-point correlation function along two different lines-of-sight, described by comoving coordinates r and r , with photons passing through each point at redshifts z and z respectively.In fact, homogeneity and isotropy guarantee that this function does not depend on r and r separately, but only on | r − r |.Inserting this into the expression above for C , we find where We can now write Q in terms of its Fourier transform over r − r , Q, giving: In the second line, we have used the Rayleigh expansion for plane waves, and j p is the spherical Bessel function of order p.In the last line, we use the orthogonality of spherical harmonics to integrate over solid angle.Substituting this expression into Eq.(D5) and integrating over n and n , once again exploiting the orthogonality of spherical harmonics, gives Finally, we can simplify this integral further assuming the Limber approximation [74][75][76], which is a high-expansion that works particularly well for the multipoles in which we are interested [96] and allows us to approximate k 2 j (kr)j (kr ) ≈ (π/2)δ D (k − /r)δ D (r − r )/r 2 .This finally gives This compact expression makes the calculation of the conversion anisotropy numerically tractable; in particular, the correlation function Q defined in Eq. (D4) need only be evaluated at points along two different lines-of-sight which have equal redshifts.We now need an expression for Q(k = /r, z, z).First, we define f 2 (ρ, m 2 γ ( r), m 2 γ ( r ); z) as the two-point PDF for m 2 γ at two different points, r and r , at redshift z, with ρ ≡ | r − r |.Note that homogeneity and isotropy guarantee that f 2 only depends on ρ.With this, the averaging in Eq. (D4) is performed by integrating the quantity with respect to m 2 γ ( r) and m 2 γ ( r ), weighted by f 2 , which gives Given both the one-and two-point PDF of m 2 γ , we can now proceed to calculate the anisotropy power spectrum.We demonstrate that the conversion anisotropy signal can satisfy existing experimental bounds by assuming that 1) the free-electron number density is exactly equal to the baryon number density, so that the one-and two-point PDFs of m 2 γ are given by the one-and two-point PDFs of baryonic fluctuations δ b divided by one and two powers of m 2 γ (z), the mean value of m 2 γ at redshift z, respectively.This is a good assumption for z 6, where resonant converions generally take place, since reionization has been completed, and all baryons are ionized; and 2) the baryonic fluctuations follow analytic one-and two-point PDFs.We consider two analytic forms for the one-and two-points PDFs: normally distributed and log-normally distributed fluctuations, to give an indication of the dependence on the PDF.In both cases, the PDFs are fully specified by the correlation function ξ b (| r − r |; z), which is related to the power spectrum of baryons P b via where g(m 2 A ; z) ≡ m 2 A /m 2 γ (z) − 1, while for log-normally distributed fluctuations [55,97], A ; z)] + Σ 2 (z)/2.In the limit that ρ → ∞, we obtain ξ b → 0 and f 2 → f 2 1 for both distributions, indicating that the two-point PDF factorizes into the product of two, independent one-point PDFs at sufficiently large separations, when correlations are unimportant.
Eqs. (D7) and (D6) give us the following general expression for C , which is the main result of this section: an expression that we can evaluate numerically for either choice of baryon fluctuation PDFs.Note that the integrand over ρ is finite as ρ → 0. Following the convention used in CMB anisotropy power spectrum analyses as well as Ref. [44], the temperature anisotropy power spectrum is defined to be ( + 1)C T 2 , where T is the sky-averaged brightness temperature.For small angular scales ( 3000), upper limits on the anisotropy of the ERB have been obtained by the Very Large Array (VLA) at 4.86 GHz [40] and 8.4 GHz [41], as well as the Australia Telescope Compact Array (ATCA) at 8.7 GHz [42].These measurements were made with the intention of looking for CMB anisotropies at the arcminute scales, and all set an approximate upper limit of ∆T /T exc 10 −2 for the ERB [39].More recently, LOFAR and TGSS observations have been used to measure the anisotropy power spectrum at ∼140 MHz [43,44].This latest result confirms the fact that a currently unknown population of dim but numerous synchrotron sources is required to explain the ERB; they also find that these sources must exhibit some nontrivial clustering to produce the right power spectrum.
Fig. 4 shows the predicted anisotropy power spectrum, [ ( + 1)/(2π)]C T 2 , as a function of , obtained by integrating Eq. (D10) numerically for both normally (lower line) and log-normally (upper line) distributed baryon fluctuations. 6The range between the two lines, shaded in blue, gives some indication of the uncertainty of our prediction.Representative values of m A = 3 × 10 −14 eV and ω A = 2 × 10 −4 eV have been chosen, but the results are qualitatively similar for other values of these parameters.This result should be compared with existing measurements of the anisotropy power, including 1) observed upper limits of the CMB anisotropy power spectrum at small scales between 4-9 GHz by the VLA [40,41] and the ATCA, subsequently reinterpreted in Ref. [39] as an upper limit on the anisotropy power of the excess radio power above the CMB and known point sources, and 2) the LOFAR anisotropy power spectrum of the radio background at 140 MHz [44].We take T = 255 K, the expected radio temperature at 140 MHz based on the power-law fit in Eq. ( 1) of the Letter.The anisotropy power reported by the high frequency measurements at frequency ν high has been rescaled by (140 MHz/ν high ) 2β , where β = 2.6; this assumes that the relative size of fluctuations δT /T is independent of frequency, and can be rescaled by T 2 at each frequency.Note that the observed upper limits in the 4-9 GHz range are not necessarily in tension with the 10 3 results at 140 MHz, since the anisotropy can in fact be frequency dependent.
The experimental results for the anisotropy power spectrum shown in Fig. 4 cover 10 2 6 × 10 4 , finding a power of 1-10 2 K 2 across this range; this represents a temperature fluctuation of 0.004 ∆T /T exc 0.04, a result that seems unusually smooth if radio emission were correlated with large scale structure at low redshifts [39].Our predicted conversion anisotropy power spectrum for the fiducial values shown here lies mostly below the measured anisotropy power, and is therefore consistent with these experimental results, except for 800 at 140 MHz.However, as mentioned in the Letter, lower-measurements of the anisotropy are difficult, and within the LOFAR dataset, significant scatter can be observed in the anisotropy power between adjacent frequency bands, which are subsequently combined to produce the final result shown here [44].We have thus demonstrated that our model can in principle produce a sufficiently smooth conversion anisotropy power spectrum; more detailed studies involving more realistic one-and two-point PDFs, as well as potentially a more careful analysis of the 4-9 GHz radio data may be of interest in future work.Figure 7 shows the excess temperature posterior for all these systematic variations, while Tab.I shows the posterior summaries for each case considered.The parameters ranges compatible with observations show minor variations from case to case, while providing a visually good fit to the data in all cases.
We summarize the inferred posteriors in Tab.I through the median and 68% containment of the individual marginal parameter posteriors (labeled "Marginal") and the 68% highest posterior density intervals, describing the shortest interval in the higher-dimensional joint parameter space where 68% of the posterior mass is contained (labeled "HDPI").

2 XFIG. 1 .
FIG.1.A schematic representation of the key aspects of the proposed class of models.

2 66 FIG. 2 .
FIG. 2. (Left) Point-wise posterior for Texc within our proposed model, showing the middle-68% and 95% regions (dark and light blue regions, respectively).We include a subdominant but irreducible contribution from unresolved extragalactic sources Teg [16] (dashed grey) for completeness.The spectrum for a single point in parameter space is shown in pink.Radio data, plotted as T − T FIRAS 0 , include measurements from ARCADE 2 are shown in red [3], with results from other telescopes compiled in the same reference shown in black.A fit to the ARCADE 2 data points assuming no stimulated emission and Teg = 0 is shown by the pink dashed line.The inset shows the same quantities scaled by a factor of ν 5/2 in order to more clearly expose the posterior in relation to the measured data.(Right) The inferred marginal and joint posterior distributions over a subset of parameters-T 0 /T0, m A , τa, and -in our fiducial model.Median and middle-68% containment values are indicated in the subtitles, while for log 10 (ma/eV) we obtain −3.66 +0.31 −0.29 (not shown).The potential reaches of PIXIE [56, 61], and CMB-S4 [62] in our parameter space are shown in blue, with projected unconstrained regions indicated by the arrow.

Appendix C: Experimental constraints 1 .
Dark matter lifetime constraint with stimulated decay

FIG. 6 .
FIG. 6. Joint and individual marginal posterior obtained the fiducial model, same format as Fig. 2.
where we can use the approximation in Eq. (B5).The second part of the integral dominates, but the total contribution to the distortion is given numerically by