††thanks: ORCID: 0000-0003-1122-6606††thanks: ORCID:
0000-0003-2486-0681††thanks: ORCID: 0000-0001-9088-7845††thanks: ORCID:
0000-0001-6051-9216

# A Stimulating Explanation of the Extragalactic Radio Background

Andrea Caputo<EMAIL_ADDRESS>School of Physics and Astronomy, Tel-Aviv
University, Tel-Aviv 69978, Israel  Department of Particle Physics and
Astrophysics,Weizmann Institute of Science, Rehovot 7610001, Israel  Hongwan
Liu<EMAIL_ADDRESS>Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA
Department of Physics, Princeton University, Princeton, NJ 08544, USA
Siddharth Mishra-Sharma<EMAIL_ADDRESS>The NSF AI Institute for Artificial
Intelligence and Fundamental Interactions  Center for Theoretical Physics,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Department of
Physics, Harvard University, Cambridge, MA 02138, USA  Maxim Pospelov  School
of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
William I. Fine Theoretical Physics Institute, School of Physics and
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  Joshua T.
Ruderman<EMAIL_ADDRESS>Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA

###### Abstract

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.

††preprint: MIT-CTP/5441

Introduction.—The cosmic microwave background (CMB) between
$60600\text{\,}\mathrm{GHz}$ is one of the most well-studied 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.

Much less is known about the extragalactic radio background (ERB) at
frequencies $\nu\lesssim$10\text{\,}\mathrm{GHz}$$. The ARCADE 2 collaboration
combined their own observations [3] with measurements from other radio
telescopes [4, 5, 6, 7], spanning a frequency range of
$22\text{\,}\mathrm{MHz}$–$90\text{\,}\mathrm{GHz}$, and found that the ERB
can be modeled as

$T(\nu)=T_{0}+T_{\text{R}}\left(\frac{\nu}{$310\text{\,}\mathrm{MHz}$}\right)^{\beta}\,,$
(1)

where $T(\nu)$ denotes the brightness temperature. $T$ is related to the
photon spectrum via
$T(\omega)=(\pi^{2}/\omega)\mathrm{d}n_{\gamma}/\mathrm{d}\omega$, where
$\mathrm{d}n_{\gamma}/\mathrm{d}\omega$ is the number density of photons per
unit energy (taking $\hbar=c=k_{\text{B}}=1$), and $\omega\equiv 2\pi\nu$. The
power-law fit in Ref. [3] found $T_{0}=$2.725(1)\text{\,}\mathrm{K}$$,
$T_{\text{R}}=$24.1(21)\text{\,}\mathrm{K}$$ and $\beta=-2.599(36)$. For
$\nu\gtrsim$100\text{\,}\mathrm{GHz}$$, this is consistent with the CMB
temperature of $T_{0}^{\text{FIRAS}}=$2.7255(85)\text{\,}\mathrm{K}$$ measured
by FIRAS [1, 8]. At lower frequencies, however, the power law is an excellent
fit to $T_{\text{exc}}\equiv T-T_{0}$. A more recent re-analysis including
all-sky maps from the LWA1 Low Frequency Sky Survey (LLFSS) reached a similar
conclusion [9].

Explaining $T_{\text{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\lesssim\beta\lesssim-2.5$ [11, 12]. Extragalactic radio synchrotron
sources such as active galactic nuclei and star-forming galaxies can produce
$\beta\approx-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_{\text{exc}}$ [16, 18, 19, 20, 21, 22] between
$100\text{\,}\mathrm{MHz}$–$10\text{\,}\mathrm{GHz}$. Electrons that are
reaccelerated during cluster mergers have been proposed as a significant
contribution to $T_{\text{exc}}$ [23], albeit under optimistic assumptions;
this mechanism also produces a softer power law than is required.

An alternative explanation of $T_{\text{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 [27,
28, 26]. 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_{\text{exc}}$ that is fully produced at high redshifts. In particular, the
LEDA and LOFAR results constrain any cosmological radio background with
$T_{\text{exc}}\propto\omega^{-2.6}$ at $\sim 10\%$ ($13.2<z<27.4$) [29] and
$\sim 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 [40, 41, 42, 39, 43, 44]. These measurements find
that fluctuations in $T_{\text{exc}}$, $\Delta T/T_{\text{exc}}$, are $\sim
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_{\text{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].

$\mathnormal{a}$${\color[rgb]{0.29,0.33,0.89}\definecolor[named]{pgfstrokecolor}{rgb}{0.29,0.33,0.89}A^{\prime}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}$
$\displaystyle A^{\prime}$ blackbody stimulates
$A^{\prime}$$\gamma$$T\varpropto\omega^{-5/2}$$X$

Figure 1: A schematic representation of the key aspects of the proposed class
of models.

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_{\text{exc}}$. These models rely on three basic ingredients: 1) a
particle decaying into dark photons $A^{\prime}$; 2) the presence of a thermal
bath of $A^{\prime}$ which stimulates this decay; and 3) $A^{\prime}$
resonantly oscillating into radio photons. This class of models leads to
$T_{\text{exc}}\propto\omega^{-5/2}$, close to the observed power-law
dependence. Relatively low-$z$ resonant oscillations as well as
$\mathrm{d}n_{A^{\prime}}/\mathrm{d}\omega\propto\omega^{-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^{\prime}$, and $A^{\prime}$ 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^{\prime\prime}$, which then
resonantly converts twice via $A^{\prime\prime}\to A^{\prime}\to\gamma$, with
some dark states charged under a dark photon $A^{\prime\prime}$ facilitating
the $A^{\prime\prime}\to A^{\prime}$ transition. However, two resonant
conversions makes it difficult to produce a sufficiently large
$T_{\text{exc}}$ while ensuring that $T_{\text{exc}}\propto\omega^{-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_{\text{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.—A particle physics model that has the following three features can
generate an ERB with $T_{\text{exc}}\propto\omega^{-5/2}$:

  1. 1.

a cold component of DM $a$ which is stable on cosmological timescales, but
undergoes a two-body decay with lifetime $\tau_{\text{vac}}$ in vacuum (for
simplicity, we take $a$ to be all of DM);

  2. 2.

one of the daughter particles of the decay is a dark photon, $A^{\prime}$,
with an existing blackbody distribution, characterized by a temperature
$T^{\prime}(z)=T_{0}^{\prime}(1+z)$, where $T_{0}^{\prime}$ is its temperature
today; and

  3. 3.

the $A^{\prime}$ has mass $${10}^{-15}\text{\,}\mathrm{eV}$\lesssim
m_{A^{\prime}}\lesssim${10}^{-9}\text{\,}\mathrm{eV}$$ and is emitted
relativistically with energy $\omega_{A^{\prime}}$. $A^{\prime}$ mixes with
the Standard Model photon $\gamma$ with kinetic mixing parameter $\epsilon$.

These three features are summarized in Fig. 1. The existence of the thermal
population of $A^{\prime}$ enhances the decay rate of $a$ due to Bose
enhancement [52, 53], leading to a redshift-dependent effective decay lifetime
$\tau(z)$, where

$\displaystyle\tau(z)=\tau_{\text{vac}}\left[1+nf^{\text{BB}}_{A^{\prime}}(z)\right]^{-1}\,,$
(2)

with $f^{\text{BB}}_{A^{\prime}}=(e^{\omega_{A^{\prime}}/T^{\prime}}-1)^{-1}$
being the blackbody occupation number of $A^{\prime}$ with energy
$\omega_{A^{\prime}}$.

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 $\alpha=1$ or $2$ dark photons. Finally,
$\omega_{A^{\prime}}$ depends on the kinematics of the decay of $a$, with
$\omega_{A^{\prime}}=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_{\text{exc}}\propto\omega^{-5/2}$,
processes such as inverse decays of $A^{\prime}$ into $a$ cannot significantly
distort either the power-law index of the $A^{\prime}$ spectrum produced by
the decay, or the blackbody distribution of $A^{\prime}$ across all relevant
frequencies. In accordance with this fiducial model, we fix $n=2$, $\alpha=1$,
and $\omega_{A^{\prime}}=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^{\prime}$ and $\gamma$ enables resonant
conversion between the two particles whenever
$m_{A^{\prime}}^{2}=m_{\gamma}^{2}$ [54, 55], where $m_{\gamma}^{2}$ is the
effective photon plasma mass squared. This quantity is proportional to the
free electron number density, $n_{\text{e}}$. The converted photons ultimately
form the present-day $T_{\text{exc}}$. We calculate the sky-averaged
conversion probability per redshift $\mathrm{d}\langle
P_{A^{\prime}\to\gamma}\rangle/\mathrm{d}z$ taking into account
inhomogeneities using the formalism developed in Refs. [56, 55] (see also
Refs. [57, 58, 59]).

Radio background production.—The particle $a$ decays throughout cosmic
history, producing a number density of $A^{\prime}$ per unit energy given by
(see App. A for a complete derivation, which follows from Ref. [60])

$\displaystyle\frac{\mathrm{d}n_{A^{\prime}}}{\mathrm{d}\omega}(\omega,z)=\frac{\alpha\rho_{a}(z)}{m_{a}\omega\tau(z_{\star})H(z_{\star})}\Theta(z_{\star}-z)\,,$
(3)

where $1+z_{\star}\equiv\omega_{A^{\prime}}(1+z)/\omega$ is the redshift at
which a daughter $A^{\prime}$ with frequency $\omega$ at redshift $z$ was
produced, $H$ is the Hubble parameter, $\rho_{a}(z)$ is the DM energy density,
and $\Theta$ is the Heaviside step function. The total spectrum of
$A^{\prime}$ is then a sum of this decay spectrum and the blackbody spectrum
of the thermal $A^{\prime}$ 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 $\mathrm{d}\langle
P_{A^{\prime}\to\gamma}\rangle/\mathrm{d}z$, i.e.,

$\displaystyle\frac{1}{(1+z)^{3}}\frac{\mathrm{d}n_{\gamma}}{\mathrm{d}x}=\int_{z}^{\infty}\mathrm{d}z^{\prime}\,\frac{\mathrm{d}\langle
P_{A^{\prime}\to\gamma}\rangle}{\mathrm{d}z^{\prime}}\frac{1}{(1+z^{\prime})^{3}}\frac{\mathrm{d}n_{A^{\prime}}}{\mathrm{d}x}(x,z^{\prime})\,,$
(4)

where $x\equiv\omega/[T_{0}(1+z)]$. For the range of $m_{A^{\prime}}$ that is
of interest, resonant conversions occur only after recombination. In the range
$3\times 10^{-4}\lesssim x\lesssim 0.2$, which are relevant for
$T_{\text{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

$\displaystyle\frac{\mathrm{d}n_{\gamma}}{\mathrm{d}x}(x,z)=\frac{\rho_{a}(z)}{m_{a}}\frac{\alpha}{x}\underbrace{\frac{1}{\tau(z_{\star})}}_{\propto\,x^{-1}}\underbrace{\frac{1}{H(z_{\star})}}_{\propto\,x^{3/2}}\underbrace{\int_{z}^{z_{\star}}\mathrm{d}z^{\prime}\frac{\mathrm{d}\langle
P_{A^{\prime}\to\gamma}\rangle}{\mathrm{d}z^{\prime}}}_{\propto\,x^{-1}}\,.$
(5)

This gives $\mathrm{d}n_{\gamma}/\mathrm{d}x\propto x^{-3/2}$, or
$T_{\text{exc}}\propto\omega^{-5/2}$, the desired frequency dependence. Note
that $z_{\star}$ must occur during matter domination in order for
$H(z_{\star})\propto x^{-3/2}$, while $\mathrm{d}\langle
P_{A^{\prime}\to\gamma}\rangle/\mathrm{d}z\propto x^{-1}$ is derived in Refs.
[54, 55]. $\tau(z_{\star})\propto x$ follows from Eq. (2) when
$T^{\prime}\gg\omega_{A^{\prime}}$.

Figure 2: _(Left)_ Point-wise posterior for $T_{\text{exc}}$ 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 $T_{\text{eg}}$ [16] (dashed grey) for
completeness. The spectrum for a single point in parameter space is shown in
pink. Radio data, plotted as $T-T_{0}^{\text{FIRAS}}$, 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 $T_{\text{eg}}=0$ is shown by the
pink dashed line. The inset shows the same quantities scaled by a factor of
$\nu^{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}^{\prime}/T_{0}$,
$m_{A^{\prime}},\tau_{a}$, and $\epsilon$—in our fiducial model. Median and
middle-68% containment values are indicated in the subtitles, while for
$\log_{10}(m_{a}/$\text{\,}\mathrm{eV}$)$ we obtain $-3.66^{+0.31}_{-0.29}$
(not shown). The potential reaches of PIXIE [61, 56], and CMB-S4 [62] in our
parameter space are shown in blue, with projected unconstrained regions
indicated by the arrow.

Eq. (5) is the main result of this paper. For our fiducial model, we find the
following approximate parametric dependence of $T_{\text{exc}}$ at $z=0$:

$T_{\text{exc}}(\omega)\approx$10\text{\,}\mathrm{K}$\left(\frac{\nu}{$310\text{\,}\mathrm{MHz}$}\right)^{-5/2}\left(\frac{T_{0}^{\prime}}{0.2T_{0}}\right)\left(\frac{${10}^{21}\text{\,}\mathrm{s}$}{\tau_{\text{vac}}}\right)\\\
\times\left(\frac{${10}^{-4}\text{\,}\mathrm{eV}$}{\omega_{A^{\prime}}}\right)^{3/2}\left(\frac{$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$}{m_{a}}\right)\left(\frac{P_{\text{GHz}}}{10^{-5}}\right)\,,$
(6)

where $P_{\text{GHz}}$ is $\langle P_{A^{\prime}\to\gamma}\rangle$ for
$\nu=$1\text{\,}\mathrm{GHz}$$ today, which the FIRAS measurement of the CMB
energy spectrum limits to be $P_{\text{GHz}}\lesssim 10^{-2}$ [54, 56].

Fit and constraints.— We can now perform a fit of $T_{\text{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^{\prime}},\tau_{\text{vac}},T_{0}^{\prime},\epsilon\\}$ from
our new physics model, and $T_{0}$. Fitting the radio data places some
requirements on the model parameters. For $T_{\text{exc}}$ to be present
across all data points, we require $m_{a}/2\gtrsim
2\pi\times$10\text{\,}\mathrm{GHz}$$. On the other hand, photons at
$22\text{\,}\mathrm{MHz}$ must originate from decays of $a$ during matter
domination to satisfy $T_{\text{exc}}\propto\omega^{-5/2}$, leading to the
requirement that $m_{a}/2\lesssim
2\pi(1+z_{\text{eq}})\times$22\text{\,}\mathrm{MHz}$$, where $z_{\text{eq}}$
is the redshift of matter-radiation equality. Together, this means
$$8\text{\times}{10}^{-5}\text{\,}\mathrm{eV}$\lesssim
m_{a}\lesssim$6\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$$. Resonant
conversion must predominantly occur after $z_{\star}$ for
$10\text{\,}\mathrm{GHz}$ photons so that the full power law extends to at
least that frequency; the requirements on $m_{a}$ show that $z_{\star}\lesssim
6$, which favors $${10}^{-14}\text{\,}\mathrm{eV}$\lesssim
m_{A^{\prime}}\lesssim$3\text{\times}{10}^{-13}\text{\,}\mathrm{eV}$$.
Finally, $T^{\prime}$ must be sufficiently large for
$f_{A^{\prime}}^{\text{BB}}\gg 1$ and $\tau(z_{\star})\propto x$ for the
$10\text{\,}\mathrm{GHz}$ data points; this imposes
$T_{0}^{\prime}/T_{0}\gtrsim 0.06$.

Several important experimental constraints exist on the class of models under
consideration. First, $T_{0}^{\prime}\lesssim 0.4T_{0}$, in order to avoid
introducing excessive effective relativistic degrees of freedom [2]. The
kinetic mixing between $A^{\prime}$ and $\gamma$ also leads to spectral
distortions due to $\gamma\to A^{\prime}$ resonant conversions; the FIRAS
measurement of the CMB spectrum [1] leads to constraints on the order of
$\epsilon<10^{-7}$–$10^{-5}$ in the $m_{A^{\prime}}$ range of interest [56,
55].

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^{\prime}$ thermal population by the present-
day to be less than 2.2%, we find
$\tau_{\text{vac}}>${10}^{19}\text{\,}\mathrm{s}$\times\max\left[2.1,6.3\left(${10}^{-4}\text{\,}\mathrm{eV}$/\omega_{A^{\prime}}\right)\left(T_{0}^{\prime}/0.2T_{0}\right)\right]$
(see App. C for details). We also note that the decaying particle $a$ may be a
subcomponent $f_{\rm DM}$ of DM, which evades the lifetime bound altogether if
$f_{\rm DM}\lesssim 2.2\%$ [67].

21-cm power spectrum measurements from both LEDA and LOFAR are both in tension
with $T_{\text{exc}}$ being fully produced at high redshifts ($z\gtrsim 9.1$),
placing strong constraints on models that produce $T_{\text{exc}}$
primordially. In our model, however, $A^{\prime}\to\gamma$ resonant
conversions occur predominantly at $z\lesssim 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.
[56, 55] and 2) the DM lifetime limit discussed above. Priors on
$T_{0}^{\prime}$ and $T_{0}$ are also chosen to account for the
$N_{\text{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_{\text{exc}}$ from our model, we include an irreducible
contribution from unresolved extragalactic radio sources of
$T_{\text{eg}}=$0.23\text{\,}\mathrm{K}$(\nu/$\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$)^{-2.7}$
[16, 39], which is at least 3 times smaller than the measured brightness
temperature. This parametrization of $T_{\text{eg}}$ is in excellent agreement
with other independent estimates [73]. Although this power-law expression
strictly applies in the range
$100\text{\,}\mathrm{MHz}$–$10\text{\,}\mathrm{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}\simeq$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$$,
$m_{A^{\prime}}\simeq$2.5\text{\times}{10}^{-14}\text{\,}\mathrm{eV}$$,
$\tau\simeq$1.3\text{\times}{10}^{21}\text{\,}\mathrm{s}$$,
$T_{0}^{\prime}\simeq 0.22\,T_{0}$ and $\epsilon\simeq 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\times 10^{3}\lesssim\ell\lesssim 6\times 10^{4}$ in the
$48\text{\,}\mathrm{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 $\sim$$140\text{\,}\mathrm{MHz}$ for $10^{2}\lesssim\ell\lesssim 10^{4}$.
Such measurements can be challenging: observations at
$\sim$$140\text{\,}\mathrm{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_{\text{exc}}$ in our model, we
take the $48\text{\,}\mathrm{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 $\delta_{\text{e}}$
during $A^{\prime}\to\gamma$ conversions, since $m_{\gamma}^{2}\propto
n_{\text{e}}$. The decay anisotropy was found to exceed the radio anisotropy
power spectrum in the $48\text{\,}\mathrm{GHz}$ range, unless the decay that
produces these photons happens at $z_{\star}\gtrsim 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^{\prime}\to\gamma$ in two different directions in the sky. As such, it
depends on the two-point function of $\delta_{\text{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.

Figure 3: Predicted anisotropy power spectrum with normally (solid black line,
below) and log-normally (solid black line, above) distributed baryon
fluctuations. The model parameters are
$m_{A^{\prime}}=$3\text{\times}{10}^{-14}\text{\,}\mathrm{eV}$$, and
$\omega_{A^{\prime}}=$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$$. Upper
limits from VLA at $4.86\text{\,}\mathrm{GHz}$ [40] (purple triangles),
$8.4\text{\,}\mathrm{GHz}$ [41] (orange triangles) and ATCA at
$8.7\text{\,}\mathrm{GHz}$ [42] (green triangle) are shown, and have been
rescaled to $T_{\text{exc}}=$255\text{\,}\mathrm{K}$$, the expected value at
$140\text{\,}\mathrm{MHz}$. Representative data points of the power spectrum
measured by LOFAR are shown in red [44].

Fig. 3 shows the predicted anisotropy power spectrum in units of
$\text{\,}{\mathrm{K}}^{2}$, defined using the same conventions as in the CMB
power spectrum [2], for normally and log-normally distributed
$\delta_{\text{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 LOFAR data, as well as the
upper limits from VLA and ATCA at their respective frequencies $\nu$, which
have been rescaled by $($140\text{\,}\mathrm{MHz}$/\nu)^{2\beta}$, with
$\beta=-2.6$; this assumes that $\Delta T/T_{\text{exc}}$ is independent of
frequency. Our calculated anisotropy power spectrum for a fiducial choice of
$m_{A^{\prime}}=$3\text{\times}{10}^{-14}\text{\,}\mathrm{eV}$$ and
$\omega_{A^{\prime}}=$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$$ lies below
the $48\text{\,}\mathrm{GHz}$ measurements for both choices of the two-point
PDF of $\delta_{\text{e}}$, while only the log-normal two-point PDF exceeds
the low-$\ell$ measurements by LOFAR. However, significant scatter exists
between adjacent frequency bands in the data for $\ell\lesssim 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
$\epsilon$-$m_{A^{\prime}}$ plane. The thermal population of $A^{\prime}$ may
also lead to a value of $N_{\text{eff}}$ that is detectable in future CMB
experiments such as CMB-S4 [62]. The potential reach of both experiments are
shown in Fig. 2.

Acknowledgments.—We thank J. Dowell for useful correspondence, and D. Aloni,
A. Hook, J. Muñoz, N. Outmezguine, J. Pradler, M. Regis, M. Schmaltz, N.
Sehgal and A. Urbano for helpful discussions. AC is supported by the Foreign
Postdoctoral Fellowship Program of the Israel Academy of Sciences and
Humanities and also acknowledges support from the Israel Science Foundation
(Grant 1302/19), the US-Israeli BSF (Grant 2018236), the German-Israeli GIF
(Grant I-2524-303.7) and the European Research Council (ERC) under the EU
Horizon 2020 Programme (ERC-CoG-2015-Proposal n. 682676 LDMThExp). HL and JTR
are supported by NSF grant PHY-1915409. HL is also supported by the DOE under
Award Number DE-SC0007968 and the Simons Foundation. JTR is also supported by
the NSF CAREER grant PHY-1554858. MP is supported in part by U.S. Department
of Energy Grant No. DE-SC0011842. This work is supported by the National
Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI
Institute for Artificial Intelligence and Fundamental Interactions,
http://iaifi.org/). This work is partly supported by the U.S. Department of
Energy, Office of Science, Office of High Energy Physics of U.S. Department of
Energy under grant Contract Number DE-SC0012567. This work was performed in
part at the Aspen Center for Physics, which is supported by NSF grant
PHY-1607611. This work made use of the NYU IT High Performance Computing
resources, services, and staff expertise. This work was also performed using
the Princeton Research Computing resources at Princeton University which is a
consortium of groups including the Princeton Institute for Computational
Science and Engineering and the Princeton University Office of Information
Technology’s Research Computing department. This research has made use of
NASA’s Astrophysics Data System. This research made use of the ArviZ [78],
astropy [79, 80], CLASS [81], dynesty [72], HyRec [82], IPython [83], Jupyter
[84], matplotlib [85], NumPy [86], pyfftlog [87, 88, 89], seaborn [90], pandas
[91], SciPy [92], and tqdm [93] software packages.

## References

  * Fixsen  _et al._ [1996] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, The Cosmic Microwave Background spectrum from the full COBE FIRAS data set, Astrophys. J. 473, 576 (1996), arXiv:astro-ph/9605054 [astro-ph] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Fixsen  _et al._ [2011] D. J. Fixsen _et al._ , ARCADE 2 Measurement of the Extra-Galactic Sky Temperature at 3-90 GHz, Astrophys. J. 734, 5 (2011), arXiv:0901.0555 [astro-ph.CO] .
  * Roger  _et al._ [1999] R. S. Roger, C. H. Costain, T. L. Landecker, and C. M. Swerdlyk, The radio emission from the galaxy at 22 mhz, Astron. Astrophys. Suppl. Ser. 137, 7 (1999), arXiv:astro-ph/9902213 .
  * Maeda  _et al._ [1999] K. Maeda, H. Alvarez, J. Aparici, J. May, and P. Reich, A 45-MHz continuum survey of the northern hemisphere, Astronomy and Astrophysics, Supplement 140, 145 (1999).
  * Haslam  _et al._ [1981] C. G. T. Haslam, U. Klein, C. J. Salter, H. Stoffel, W. E. Wilson, M. N. Cleary, D. J. Cooke, and P. Thomasson, A 408 MHz all-sky continuum survey. I - Observations at southern declinations and for the North Polar region., Astronomy and Astrophysics 100, 209 (1981).
  * Reich and Reich [1986] P. Reich and W. Reich, A radio continuum survey of the northern sky at 1420 MHz. II, Astronomy and Astrophysics, Supplement 63, 205 (1986).
  * Fixsen [2009] D. J. Fixsen, The Temperature of the Cosmic Microwave Background, Astrophys. J. 707, 916 (2009), arXiv:0911.1955 [astro-ph.CO] .
  * Dowell and Taylor [2018] J. Dowell and G. B. Taylor, The Radio Background Below 100 MHz, Astrophys. J. Lett. 858, L9 (2018), arXiv:1804.08581 [astro-ph.CO] .
  * Singal  _et al._ [2018] J. Singal _et al._ , The Radio Synchrotron Background: Conference Summary and Report, Publ. Astron. Soc. Pac. 130, 036001 (2018), arXiv:1711.09979 [astro-ph.HE] .
  * Wright  _et al._ [1999] M. Wright, J. Dickel, B. Koralesky, and L. Rudnick, The Supernova Remnant Cassiopeia A at Millimeter Wavelengths, Astrophys. J. 518, 284 (1999).
  * Scheuer and Williams [1968] P. A. G. Scheuer and P. J. S. Williams, Radio Spectra, The Annual Review of Astronomy and Astrophysics 6, 321 (1968).
  * Platania  _et al._ [1998] P. Platania, M. Bensadoun, M. Bersanelli, G. De Amici, A. Kogut, S. Levin, D. Maino, and G. F. Smoot, A determination of the spectral index of galactic synchrotron emission in the 1-10 GHz range, Astrophys. J. 505, 473 (1998), arXiv:astro-ph/9707252 .
  * Carilli  _et al._ [1991] C. L. Carilli, R. A. Perley, J. W. Dreher, and J. P. Leahy, Multifrequency radio observations of Cygnus A - Spectral aging in powerful radio galaxies, Astrophys. J. 383, 554 (1991).
  * Protheroe and Biermann [1996] R. J. Protheroe and P. L. Biermann, A New estimate of the extragalactic radio background and implications for ultrahigh-energy gamma-ray propagation, Astropart. Phys. 6, 45 (1996), [Erratum: Astropart.Phys. 7, 181 (1997)], arXiv:astro-ph/9605119 .
  * Gervasi  _et al._ [2008] M. Gervasi, A. Tartari, M. Zannoni, G. Boella, and G. Sironi, The contribution of the Unresolved Extragalactic Radio Sources to the Brightness Temperature of the sky, Astrophys. J. 682, 223 (2008), arXiv:0803.4138 [astro-ph] .
  * Niţu  _et al._ [2021] I. C. Niţu, H. T. J. Bevins, J. D. Bray, and A. M. M. Scaife, An updated estimate of the cosmic radio background and implications for ultra-high-energy photon propagation, Astropart. Phys. 126, 102532 (2021), arXiv:2004.13596 [astro-ph.HE] .
  * Vernstrom  _et al._ [2011] T. Vernstrom, D. Scott, and J. V. Wall, Contribution to the diffuse radio background from extragalactic radio sources, Monthly Notices of the Royal Astronomical Society 415, 3641 (2011).
  * Condon  _et al._ [2012] J. J. Condon, W. D. Cotton, E. B. Fomalont, K. I. Kellermann, N. Miller, R. A. Perley, D. Scott, T. Vernstrom, and J. V. Wall, Resolving the radio source background: Deeper understanding through confusion, The Astrophysical Journal 758, 23 (2012).
  * Fornengo  _et al._ [2014] N. Fornengo, R. A. Lineros, M. Regis, and M. Taoso, The isotropic radio background revisited, JCAP 04, 008, arXiv:1402.2218 [astro-ph.CO] .
  * Vernstrom  _et al._ [2015] T. Vernstrom, R. P. Norris, D. Scott, and J. V. Wall, The deep diffuse extragalactic radio sky at 1.75 GHz, Monthly Notices of the Royal Astronomical Society 447, 2243 (2015), arXiv:1408.4160 [astro-ph.GA] .
  * Hardcastle  _et al._ [2021] M. J. Hardcastle, T. W. Shimwell, C. Tasse, P. N. Best, A. Drabent, M. J. Jarvis, I. Prandoni, H. J. A. Röttgering, J. Sabater, and D. J. Schwarz, The contribution of discrete sources to the sky temperature at 144 MHz, Astron. Astrophys. 648, A10 (2021), arXiv:2011.08294 [astro-ph.CO] .
  * Fang and Linden [2016] K. Fang and T. Linden, Cluster Mergers and the Origin of the ARCADE-2 Excess, JCAP 10, 004, arXiv:1506.05807 [astro-ph.HE] .
  * Subrahmanyan and Cowsik [2013] R. Subrahmanyan and R. Cowsik, Is there an Unaccounted for Excess in the Extragalactic Cosmic Radio Background?, Astrophys. J. 776, 42 (2013), arXiv:1305.7060 [astro-ph.CO] .
  * Orlando and Strong [2013] E. Orlando and A. Strong, Galactic synchrotron emission with cosmic ray propagation models, Mon. Not. Roy. Astron. Soc. 436, 2127 (2013), arXiv:1309.2947 [astro-ph.GA] .
  * Krause and Hardcastle [2021] M. G. H. Krause and M. J. Hardcastle, Can the Local Bubble explain the radio background?, Mon. Not. Roy. Astron. Soc. 502, 2807 (2021), arXiv:2101.05255 [astro-ph.HE] .
  * Singal  _et al._ [2010] J. Singal, L. Stawarz, A. Lawrence, and V. Petrosian, Sources of the radio background considered, Monthly Notices of the Royal Astronomical Society 409, 1172 (2010).
  * Singal  _et al._ [2015] J. Singal, A. Kogut, E. Jones, and H. Dunlap, Axial Ratio of Edge-on Spiral Galaxies as a Test for Bright Radio Halos, Astrophys. J. Lett. 799, L10 (2015), arXiv:1501.00499 [astro-ph.GA] .
  * Feng and Holder [2018] C. Feng and G. Holder, Enhanced global signal of neutral hydrogen due to excess radiation at cosmic dawn, Astrophys. J. 858, L17 (2018), arXiv:1802.07432 [astro-ph.CO] .
  * Biermann  _et al._ [2014] P. L. Biermann, B. B. Nath, L. I. Caramete, B. C. Harms, T. Stanev, and J. Becker Tjus, Cosmic backgrounds due to the formation of the first generation of supermassive black holes, Mon. Not. Roy. Astron. Soc. 441, 1147 (2014), arXiv:1403.3804 [astro-ph.CO] .
  * Ewall-Wice  _et al._ [2018] A. Ewall-Wice, T. C. Chang, J. Lazio, O. Dore, M. Seiffert, and R. A. Monsalve, Modeling the Radio Background from the First Black Holes at Cosmic Dawn: Implications for the 21 cm Absorption Amplitude, Astrophys. J. 868, 63 (2018), arXiv:1803.01815 [astro-ph.CO] .
  * Ewall-Wice  _et al._ [2020] A. Ewall-Wice, T.-C. Chang, and T. J. W. Lazio, The Radio Scream from black holes at Cosmic Dawn: a semi-analytic model for the impact of radio-loud black holes on the 21 cm global signal, Mon. Not. Roy. Astron. Soc. 492, 6086 (2020), arXiv:1903.06788 [astro-ph.GA] .
  * Mirocha and Furlanetto [2019] J. Mirocha and S. R. Furlanetto, What does the first highly-redshifted 21-cm detection tell us about early galaxies?, Mon. Not. Roy. Astron. Soc. 483, 1980 (2019), arXiv:1803.03272 [astro-ph.GA] .
  * Sharma [2018] P. Sharma, Astrophysical radio background cannot explain the EDGES 21-cm signal: constraints from cooling of non-thermal electrons, Mon. Not. Roy. Astron. Soc. 481, L6 (2018), arXiv:1804.05843 [astro-ph.HE] .
  * Bowman  _et al._ [2018] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen, and N. Mahesh, An absorption profile centred at 78 megahertz in the sky-averaged spectrum, Nature 555, 67 (2018), arXiv:1810.05912 [astro-ph.CO] .
  * Bernardi  _et al._ [2016] G. Bernardi, J. T. L. Zwart, D. Price, L. J. Greenhill, A. Mesinger, J. Dowell, T. Eftekhari, S. W. Ellingson, J. Kocz, and F. Schinzel, Bayesian constraints on the global 21-cm signal from the Cosmic Dawn, Mon. Not. Roy. Astron. Soc. 461, 2847 (2016), arXiv:1606.06006 [astro-ph.CO] .
  * Mertens  _et al._ [2020] F. G. Mertens _et al._ , Improved upper limits on the 21-cm signal power spectrum of neutral hydrogen at $\boldsymbol{z\approx 9.1}$ from LOFAR, Mon. Not. Roy. Astron. Soc. 493, 1662 (2020), arXiv:2002.07196 [astro-ph.CO] .
  * Fialkov and Barkana [2019] A. Fialkov and R. Barkana, Signature of Excess Radio Background in the 21-cm Global Signal and Power Spectrum, Mon. Not. Roy. Astron. Soc. 486, 1763 (2019), arXiv:1902.02438 [astro-ph.CO] .
  * Holder [2013] G. P. Holder, The unusual smoothness of the extragalactic unresolved radio background, The Astrophysical Journal 780, 112 (2013).
  * Fomalont  _et al._ [1988] E. B. Fomalont, K. I. Kellermann, M. C. Anderson, D. Weistrop, J. V. Wall, R. A. Windhorst, and J. A. Kristian, New Limits to Fluctuations in the Cosmic Background Radiation at 4.86 GHz Between 12 and 60 Arcsecond Resolution, The Astronomical Journal 96, 1187 (1988).
  * Partridge  _et al._ [1997] R. B. Partridge, E. A. Richards, E. B. Fomalont, K. I. Kellerman, and R. A. Windhorst, Small-Scale Cosmic Microwave Background Observations at 8.4 GHz, Astrophys. J. 483, 38 (1997).
  * Subrahmanyan  _et al._ [2000] R. Subrahmanyan, M. J. Kesteven, R. D. Ekers, M. Sinclair, and J. Silk, An Australia telescope survey for CMB anisotropies, Mon. Not. Roy. Astron. Soc. 315, 808 (2000), arXiv:astro-ph/0002467 .
  * Choudhuri  _et al._ [2020] S. Choudhuri, A. Ghosh, N. Roy, S. Bharadwaj, H. T. Intema, and S. S. Ali, All-sky angular power spectrum – I. Estimating brightness temperature fluctuations using the 150-MHz TGSS survey, Mon. Not. Roy. Astron. Soc. 494, 1936 (2020), arXiv:2003.07869 [astro-ph.CO] .
  * Offringa  _et al._ [2021] A. R. Offringa, J. Singal, S. Heston, S. Horiuchi, and D. M. Lucero, Measurement of the anisotropy power spectrum of the radio synchrotron background, Mon. Not. Roy. Astron. Soc. 509, 114 (2021), arXiv:2110.00499 [astro-ph.CO] .
  * Fornengo  _et al._ [2011] N. Fornengo, R. Lineros, M. Regis, and M. Taoso, Possibility of a Dark Matter Interpretation for the Excess in Isotropic Radio Emission Reported by ARCADE, Phys. Rev. Lett. 107, 271302 (2011), arXiv:1108.0569 [hep-ph] .
  * Hooper  _et al._ [2012] D. Hooper, A. V. Belikov, T. E. Jeltema, T. Linden, S. Profumo, and T. R. Slatyer, The Isotropic Radio Background and Annihilating Dark Matter, Phys. Rev. D 86, 103003 (2012), arXiv:1203.3547 [astro-ph.CO] .
  * Cline and Vincent [2013] J. M. Cline and A. C. Vincent, Cosmological origin of anomalous radio background, JCAP 02, 011, arXiv:1210.2717 [astro-ph.CO] .
  * Fang and Linden [2015] K. Fang and T. Linden, Anisotropy of the extragalactic radio background from dark matter annihilation, Phys. Rev. D 91, 083501 (2015), arXiv:1412.7545 [astro-ph.HE] .
  * Ajello  _et al._ [2015] M. Ajello _et al._ , The Origin of the Extragalactic Gamma-Ray Background and Implications for Dark-Matter Annihilation, Astrophys. J. Lett. 800, L27 (2015), arXiv:1501.05301 [astro-ph.HE] .
  * Lisanti  _et al._ [2016] M. Lisanti, S. Mishra-Sharma, L. Necib, and B. R. Safdi, Deciphering Contributions to the Extragalactic Gamma-Ray Background from 2 GeV to 2 TeV, Astrophys. J. 832, 117 (2016), arXiv:1606.04101 [astro-ph.HE] .
  * Vogel and Redondo [2014] H. Vogel and J. Redondo, Dark Radiation constraints on minicharged particles in models with a hidden photon, JCAP 02, 029, arXiv:1311.2600 [hep-ph] .
  * Caputo  _et al._ [2019] A. Caputo, M. Regis, M. Taoso, and S. J. Witte, Detecting the Stimulated Decay of Axions at RadioFrequencies, JCAP 1903 (03), 027, arXiv:1811.08436 [hep-ph] .
  * Bolliet  _et al._ [2021] B. Bolliet, J. Chluba, and R. Battye, Spectral distortion constraints on photon injection from low-mass decaying particles, Mon. Not. Roy. Astron. Soc. 507, 3148 (2021), arXiv:2012.07292 [astro-ph.CO] .
  * Mirizzi  _et al._ [2009] A. Mirizzi, J. Redondo, and G. Sigl, Microwave Background Constraints on Mixing of Photons with Hidden Photons, JCAP 0903, 026, arXiv:0901.0014 [hep-ph] .
  * Caputo  _et al._ [2020a] A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruderman, Modeling Dark Photon Oscillations in Our Inhomogeneous Universe, Phys. Rev. D 102, 103533 (2020a), arXiv:2004.06733 [astro-ph.CO] .
  * Caputo  _et al._ [2020b] A. Caputo, H. Liu, S. Mishra-Sharma, and J. T. Ruderman, Dark Photon Oscillations in Our Inhomogeneous Universe, Phys. Rev. Lett. 125, 221303 (2020b), arXiv:2002.05165 [astro-ph.CO] .
  * Bondarenko  _et al._ [2020] K. Bondarenko, J. Pradler, and A. Sokolenko, Constraining dark photons and their connection to 21 cm cosmology with CMB data, Phys. Lett. B 805, 135420 (2020), arXiv:2002.08942 [astro-ph.CO] .
  * Garcia  _et al._ [2020] A. A. Garcia, K. Bondarenko, S. Ploeckinger, J. Pradler, and A. Sokolenko, Effective photon mass and (dark) photon conversion in the inhomogeneous Universe, JCAP 10, 011, arXiv:2003.10465 [astro-ph.CO] .
  * Witte  _et al._ [2020] S. J. Witte, S. Rosauro-Alcaraz, S. D. McDermott, and V. Poulin, Dark photon dark matter in the presence of inhomogeneous structure, JHEP 06, 132, arXiv:2003.13698 [astro-ph.CO] .
  * Pospelov  _et al._ [2018] M. Pospelov, J. Pradler, J. T. Ruderman, and A. Urbano, Room for New Physics in the Rayleigh-Jeans Tail of the Cosmic Microwave Background, Phys. Rev. Lett. 121, 031103 (2018), arXiv:1803.07048 [hep-ph] .
  * Kogut  _et al._ [2011] A. Kogut _et al._ , The Primordial Inflation Explorer (PIXIE): A Nulling Polarimeter for Cosmic Microwave Background Observations, JCAP 1107, 025, arXiv:1105.2044 [astro-ph.CO] .
  * Abazajian  _et al._ [2019] K. Abazajian _et al._ , CMB-S4 Science Case, Reference Design, and Project Plan (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Poulin  _et al._ [2016] V. Poulin, P. D. Serpico, and J. Lesgourgues, A fresh look at linear cosmological constraints on a decaying dark matter component, JCAP 08, 036, arXiv:1606.02073 [astro-ph.CO] .
  * Nygaard  _et al._ [2021] A. Nygaard, T. Tram, and S. Hannestad, Updated constraints on decaying cold dark matter, JCAP 05, 017, arXiv:2011.01632 [astro-ph.CO] .
  * Chen  _et al._ [2021] A. Chen _et al._ (DES), Constraints on dark matter to dark radiation conversion in the late universe with DES-Y1 and external data, Phys. Rev. D 103, 123528 (2021), arXiv:2011.04606 [astro-ph.CO] .
  * Mau  _et al._ [2022] S. Mau _et al._ , Milky Way Satellite Census. IV. Constraints on Decaying Dark Matter from Observations of Milky Way Satellite Galaxies (2022), arXiv:2201.11740 [astro-ph.CO] .
  * Simon  _et al._ [2022] T. Simon, G. F. Abellán, P. Du, V. Poulin, and Y. Tsai, Constraining decaying dark matter with BOSS data and the effective field theory of large-scale structures (2022), arXiv:2203.07440 [astro-ph.CO] .
  * Alvi  _et al._ [2022] S. Alvi, T. Brinckmann, M. Gerbino, M. Lattanzi, and L. Pagano, Do you smell something decaying? Updated linear constraints on decaying dark matter scenarios (2022), arXiv:2205.05636 [astro-ph.CO] .
  * Skilling [2006] J. Skilling, Nested sampling for general Bayesian computation, Bayesian Analysis 1, 833 (2006).
  * Feroz  _et al._ [2009] F. Feroz, M. P. Hobson, and M. Bridges, MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics, Mon. Not. Roy. Astron. Soc. 398, 1601 (2009), arXiv:0809.3437 [astro-ph] .
  * Skilling [2004] J. Skilling, Nested Sampling, in  _Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering_, American Institute of Physics Conference Series, Vol. 735, edited by R. Fischer, R. Preuss, and U. V. Toussaint (2004) pp. 395–405.
  * Speagle [2020] J. S. Speagle, dynesty: a dynamic nested sampling package for estimating bayesian posteriors and evidences, Monthly Notices of the Royal Astronomical Society 493, 3132 (2020).
  * Massardi  _et al._ [2010] M. Massardi, A. Bonaldi, M. Negrello, S. Ricciardi, A. Raccanelli, and G. de Zotti, A model for the cosmological evolution of low-frequency radio sources, Monthly Notices of the Royal Astronomical Society 404, 532 (2010), arXiv:1001.1069 [astro-ph.CO] .
  * Limber [1953] D. N. Limber, The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field., Astrophys. J. 117, 134 (1953).
  * Kaiser [1992] N. Kaiser, Weak gravitational lensing of distant galaxies, Astrophys. J. 388, 272 (1992).
  * Kaiser [1998] N. Kaiser, Weak lensing and cosmology, Astrophys. J. 498, 26 (1998), arXiv:astro-ph/9610120 .
  * Fornengo and Regis [2014] N. Fornengo and M. Regis, Particle dark matter searches in the anisotropic sky, Front. Physics 2, 6 (2014), arXiv:1312.4835 [astro-ph.CO] .
  * Kumar  _et al._ [2019] R. Kumar, C. Carroll, A. Hartikainen, and O. Martin, Arviz a unified library for exploratory analysis of bayesian models in python, Journal of Open Source Software 4, 1143 (2019).
  * Price-Whelan  _et al._ [2018] A. M. Price-Whelan _et al._ , The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package, Astron. J. 156, 123 (2018), arXiv:1801.02634 .
  * Robitaille  _et al._ [2013] T. P. Robitaille _et al._ (Astropy), Astropy: A Community Python Package for Astronomy, Astron. Astrophys. 558, A33 (2013), arXiv:1307.6212 [astro-ph.IM] .
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 1107, 034, arXiv:1104.2933 [astro-ph.CO] .
  * Ali-Haimoud and Hirata [2011] Y. Ali-Haimoud and C. M. Hirata, HyRec: A fast and highly accurate primordial hydrogen and helium recombination code, Phys. Rev. D83, 043513 (2011), arXiv:1011.3758 [astro-ph.CO] .
  * Perez and Granger [2007] F. Perez and B. E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering 9, 21 (2007).
  * Kluyver  _et al._ [2016] T. Kluyver _et al._ , Jupyter notebooks - a publishing format for reproducible computational workflows, in  _ELPUB_ (2016).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing In Science & Engineering 9, 90 (2007).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Talman [1978] J. D. Talman, Numerical fourier and bessel transforms in logarithmic variables, Journal of Computational Physics 29, 35 (1978).
  * Hamilton [2000] A. J. S. Hamilton, Uncorrelated modes of the nonlinear power spectrum, Mon. Not. Roy. Astron. Soc. 312, 257 (2000), arXiv:astro-ph/9905191 .
  * Werthmüller [2020] D. Werthmüller, prisae/pyfftlog: First packaged release (2020).
  * Waskom  _et al._ [2017] M. Waskom _et al._ , mwaskom/seaborn: v0.8.1 (september 2017) (2017).
  * McKinney [2010] W. McKinney, Data structures for statistical computing in python, in  _Proceedings of the 9th Python in Science Conference_ , edited by S. van der Walt and J. Millman (2010) pp. 51 – 56.
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods https://doi.org/10.1038/s41592-019-0686-2 (2020).
  * da Costa-Luis [2019] C. O. da Costa-Luis, tqdm: A fast, extensible progress meter for python and cli, JOSS 4, 1277 (2019).
  * Grzadkowski  _et al._ [2021] B. Grzadkowski, M. Iglicki, and S. Mrówczyński, $t$-channel singularities in cosmology and particle physics (2021), arXiv:2108.01757 [hep-ph] .
  * Redondo and Postma [2009] J. Redondo and M. Postma, Massive hidden photons as lukewarm dark matter, JCAP 0902, 005, arXiv:0811.0326 [hep-ph] .
  * LoVerde and Afshordi [2008] M. LoVerde and N. Afshordi, Extended Limber Approximation, Phys. Rev. D 78, 123506 (2008), arXiv:0809.5112 [astro-ph] .
  * Kayo  _et al._ [2001] I. Kayo, A. Taruya, and Y. Suto, Probability distribution function of cosmological density fluctuations from Gaussian initial condition: comparison of one- and two-point log-normal model predictions with n-body simulations, Astrophys. J. 561, 22 (2001), arXiv:astro-ph/0105218 [astro-ph] .
  * Alsing and Handley [2021] J. Alsing and W. Handley, Nested sampling with any prior you like, arXiv e-prints , arXiv:2102.12478 (2021), arXiv:2102.12478 [astro-ph.IM] .
  * Rezende and Mohamed [2015] D. Rezende and S. Mohamed, Variational inference with normalizing flows, in  _Proceedings of the 32nd International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 37, edited by F. Bach and D. Blei (PMLR, Lille, France, 2015) pp. 1530–1538.
  * Papamakarios  _et al._ [2019] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference, arXiv preprint arXiv:1912.02762 (2019).
  * Durkan  _et al._ [2019] C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, Neural spline flows, in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Isi  _et al._ [2022] M. Isi, W. M. Farr, and K. Chatziioannou, Comparing Bayes factors and hierarchical inference for testing general relativity with gravitational waves, arXiv e-prints (2022), arXiv:2204.10742 [gr-qc] .
  * Rybicki and Lightman [1979] G. B. Rybicki and A. P. Lightman,  _Radiative processes in astrophysics_ (1979).
  * Bell [1978] A. R. Bell, The acceleration of cosmic rays in shock fronts. II., Mon. Not. Roy. Astron. Soc. 182, 443 (1978).
  * Bell [1978] A. R. Bell, The acceleration of cosmic rays in shock fronts - I., Monthly Notices of the Royal Astronomical Society 182, 147 (1978).
  * Fermi [1949] E. Fermi, On the origin of the cosmic radiation, Phys. Rev. 75, 1169 (1949).
  * Blandford and Ostriker [1980] R. D. Blandford and J. P. Ostriker, Supernova Shock Acceleration of Cosmic Rays in the Galaxy, Astrophys. J. 237, 793 (1980).
  * Blandford and Ostriker [1978] R. D. Blandford and J. P. Ostriker, Particle Acceleration by Astrophysical Shocks, Astrophys. J. Lett. 221, L29 (1978).
  * Malkov and Drury [2001] M. Malkov and L. O. Drury, Nonlinear theory of diffusive acceleration of particles by shock waves, Rept. Prog. Phys. 64, 429 (2001).
  * Domček  _et al._ [2021] V. Domček, J. Vink, J. V. Hernández Santisteban, T. Delaney, and P. Zhou, Mapping the spectral index of Cassiopeia A: evidence for flattening from radio to infrared, Mon. Not. Roy. Astron. Soc. 502, 1026 (2021), arXiv:2005.12677 [astro-ph.HE] .

oneΔ

Appendix

This appendix is organized as follows. In App. A we derive in detail the dark
photon spectrum as a function of frequency and redshift in the class of models
considered in the Letter. App. B gives an extended description of a particular
model which presents all the ingredients necessary to explain the
extragalactic radio background (ERB). App. C presents some of the important
experimental limits on our model, including a derivation of the cosmological
limits on the dark matter lifetime in the presence of stimulated decay, as
well as limits on the dark sector temperature due to the presence of extra
relativistic degrees of freedom. App. D presents the calculation of the
anisotropy of the ERB due to photon-dark photon resonant conversions. In App.
E we describe our data analysis procedure, before presenting in App. F our
extended results, including a complete corner plot of our parameter space, and
other systematics checks such as the use of different datasets and fits using
the ARCADE 2 data points only. Finally, in App. G we give a brief summary of
the spectrum of radiation expected from synchrotron sources.

##  Appendix A Dark photon spectrum

In this section, we provide a detailed derivation of the dark photon spectrum
as a function of frequency and redshift, in the presence of stimulated
emission.

Consider a two-body decay of a heavy particle $a$ into two daughter particles,
at least one of which is a dark photon $A^{\prime}$, emitted relativistically
with energy $\omega_{A^{\prime}}$. In the presence of a blackbody distribution
of $A^{\prime}$ with temperature $T^{\prime}$, $a$ undergoes stimulated decay,
with a decay width $\Gamma(z)$ that is related to the usual vacuum decay rate
$\Gamma_{\text{vac}}$ via

$\displaystyle\Gamma(z)=(1+nf_{A^{\prime}}^{\text{BB}})\Gamma_{\text{vac}}\,,$

where
$f_{A^{\prime}}^{\text{BB}}\equiv[\exp(\omega_{A^{\prime}}/T^{\prime})-1]^{-1}$
is the blackbody occupation number at energy $\omega_{A^{\prime}}$. $n$ is an
integer that depends on the details of the model; $n=2$ in the fiducial
particle physics model that we describe below. Note that for
$\omega_{A^{\prime}}\ll T^{\prime}$, we obtain
$\Gamma(z)\simeq(nT^{\prime}/\omega_{A^{\prime}})\Gamma_{\text{vac}}$, an
approximation that we will use frequently to obtain parametric estimates.

Given $\Gamma(z^{\prime})$, the number density of $A^{\prime}$ produced per
unit time at redshift $z^{\prime}$ is

$\displaystyle\frac{\mathrm{d}n_{A^{\prime}}}{\mathrm{d}t}=\alpha\Gamma(z^{\prime})\frac{\rho_{a}(z^{\prime})}{m_{a}}\,,$

where $\alpha=2$ if both daughter particles are $A^{\prime}$ and $\alpha=1$
otherwise, $\rho_{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 $\omega$, due to
a decay occuring at redshift $z^{\prime}$:

$\displaystyle\frac{\mathrm{d}n_{A^{\prime}}(z)}{\mathrm{d}\omega\,dz^{\prime}}=\frac{\mathrm{d}n_{A^{\prime}}(z^{\prime})}{\mathrm{d}t}\frac{\mathrm{d}t}{\mathrm{d}z^{\prime}}\frac{(1+z)^{3}}{(1+z^{\prime})^{3}}\delta_{\text{D}}\left(\omega\frac{1+z^{\prime}}{1+z}-\omega_{A^{\prime}}\right)\,,$

where the redshift factors rescale the number density of particles produced at
redshift $z^{\prime}$ down to redshift $z$, and the Dirac delta function
enforces the fact that the decay produces dark photons with energy
$\omega_{A^{\prime}}$. Integrating over $z^{\prime}$ gives

$\displaystyle\frac{\mathrm{d}n_{A^{\prime}}}{\mathrm{d}\omega}$
$\displaystyle=\int_{z}^{\infty}\mathrm{d}z^{\prime}\,\frac{\rho_{a}(z^{\prime})}{m_{a}}\frac{\alpha\Gamma(z^{\prime})}{H(z)(1+z)}\frac{(1+z)^{3}}{(1+z^{\prime})^{3}}\delta_{\text{D}}\left(\omega\frac{1+z^{\prime}}{1+z}-\omega_{A^{\prime}}\right)=\frac{\alpha\rho_{a}(z)}{m_{a}}\frac{\Gamma(z_{\star})}{\omega
H(z_{\star})}\Theta\left(z_{\star}-z\right)\,,$ (7)

where

$\displaystyle 1+z_{\star}\equiv\frac{\omega_{A^{\prime}}}{\omega}(1+z)$ (8)

is the redshift at which the decay of $a$ produced the dark photon at
frequency $\omega$, with $\omega\leq\omega_{A^{\prime}}$. We have also assumed
that $\rho_{a}\propto(1+z)^{3}$. This expression agrees with Ref. [60], if we
simply set $\Gamma(z_{\star})\to\Gamma_{\text{vac}}$. In terms of comoving
number density and $x\equiv\omega/T_{0}$ where $T_{0}$ is the CMB temperature
today, we find

$\displaystyle\frac{1}{(1+z)^{3}}\frac{\mathrm{d}n_{A^{\prime}}}{\mathrm{d}x}=\frac{\alpha\rho_{a,0}}{m_{a}}\frac{\Gamma(z_{\star})}{xH(z_{\star})}\Theta(z_{\star}-z)\,,$
(9)

which is almost constant in redshift, aside from the step function.
$\rho_{a,0}$ is the dark matter energy density today. Eq. (7) is the
expression we use in the Letter.

##  Appendix B Particle physics model

###  B.1 Model building considerations

As we explained in the Letter, the particle physics model we need to explain
the radio background requires three features that we repeat here:

  1. 1.

A cold component of DM $a$ which is stable on cosmological timescales, but
undergoes a two-body decay with lifetime $\tau_{\text{vac}}$ in vacuum. For
simplicity, we take $a$ to be all of DM;

  2. 2.

One of the daughter particles is a dark photon $A^{\prime}$, which has an
existing blackbody distribution of $A^{\prime}$ characterized by a temperature
$T^{\prime}(z)=T_{0}^{\prime}(1+z)$, where $T_{0}^{\prime}$ is its temperature
today; and

  3. 3.

$A^{\prime}$ has a mass $${10}^{-15}\text{\,}\mathrm{eV}$\lesssim
m_{A^{\prime}}\lesssim${10}^{-9}\text{\,}\mathrm{eV}$$, and is emitted
relativistically with energy $\omega_{A^{\prime}}$. $A^{\prime}$ mixes with
the Standard Model (SM) photon $\gamma$ with kinetic mixing parameter
$\epsilon$.

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^{\prime}$ power-law spectrum must maintain
$\mathrm{d}n_{A^{\prime}}/\mathrm{d}\omega\propto\omega^{-1/2}$ between
$3\times 10^{-4}\lesssim x\lesssim 0.2$, where $x\equiv\omega/T_{\text{CMB}}$,
from the redshift of production $z_{\star}(x)$ until resonant conversions
$A^{\prime}\to\gamma$ are mostly complete. Furthermore, the $A^{\prime}$
blackbody spectrum also cannot be significantly distorted; specifically,
$A^{\prime}$ particles from the bath at each value of $x$ must be described by
a blackbody spectrum at $z_{\star}(x)$. In the model of Ref. [60], both of
these requirements can be violated by inverse decays $A^{\prime}A^{\prime}\to
a$ or other $A^{\prime}$ 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_{\text{exc}}$, scattering between a low-energy
$A^{\prime}$ in the power-law spectrum and a blackbody $A^{\prime}$ is too
rapid to guarantee that the $A^{\prime}$ 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^{\prime}$ 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^{\prime}$ 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^{\prime}$ particle is a constant, the Boltzmann equation governing the
occupation number $f_{A^{\prime}}$ can be written

$\displaystyle\frac{\mathrm{d}f_{A^{\prime}}(x,t)}{\mathrm{d}t}=\frac{C[f_{A^{\prime}}]}{\omega}\,,$
(10)

where $C[f_{A^{\prime}}]$ is the collision integral, which we will define for
particular processes below. Here, we adopt the convention that
$f_{A^{\prime}}$ is related to the number density via

$\displaystyle
n_{A^{\prime}}=3\int\frac{\mathrm{d}^{3}\vec{p}}{(2\pi)^{3}}f_{A^{\prime}}\,,$

where the factor of 3 accounts for the degeneracy in the spin state of
$A^{\prime}$. To determine if any process causes a significant distortion to
the spectrum of blackbody $A^{\prime}$ or low-energy $A^{\prime}$ from the
decay of $a$, we divide Eq. (10) by $f_{A^{\prime}}$ and integrate with
respect to time or, equivalently, redshift, to obtain the change in $\log
f_{A^{\prime}}$, which gives a measure of the fraction of $A^{\prime}$ at that
value of $x$ that undergoes a scattering process. Our criterion for a small
distortion to the $A^{\prime}$ spectra is therefore

$\displaystyle|\Delta\log
f_{A^{\prime}}(x)|=\int_{z_{\min}}^{z_{\max}}\frac{\mathrm{d}z}{H(z)(1+z)}\frac{|C[f_{A^{\prime}}]|}{\omega
f_{A^{\prime}}(x,z)}\ll 1\,,$ (11)

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\text{\,}\mathrm{MHz}$–$10\text{\,}\mathrm{GHz}$, which we can cover by
considering $3\times 10^{-4}<x<0.2$.

The power-law component at fixed $x$ is produced by decays of $a$ at a
redshift $1+z_{\star}\equiv\omega_{A^{\prime}}/(xT_{0})$, where
$\omega_{A^{\prime}}$ is the energy of the emitted $A^{\prime}$ at the point
of decay (for $a$ decaying into two massless particles, this is simply
$\omega_{A^{\prime}}=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.111Technically, 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. 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_{\star}(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_{\text{min}}=0$ and
$z_{\text{max}}=z_{\star}(x)$, while for the blackbody spectrum,
$z_{\text{min}}=z_{\star}(x)$.

###  B.2 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_{\psi}$ and $A_{\text{osc}}$. The three requirements for producing
the radio excess as discussed in the main body of the paper are satisfied as
follows:

  1. 1.

the DM decays through the process $a\to A_{\psi}A_{\text{osc}}$;

  2. 2.

a blackbody distribution of $A_{\text{osc}}$ described by a temperature
$T^{\prime}$ is also present, leading to stimulated decay of $a$, and

  3. 3.

the dark photon $A_{\text{osc}}$ possesses a small kinetic mixing term with
the SM photon, and has a mass $${10}^{-15}\text{\,}\mathrm{eV}$\lesssim
m_{\text{osc}}\lesssim${10}^{-9}\text{\,}\mathrm{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

$\displaystyle\mathcal{L}\supset\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{1}{2}m_{a}^{2}a^{2}-\frac{a}{4f_{a}}F_{\psi,\mu\nu}\tilde{F}_{\text{osc}}^{\mu\nu}-\frac{\epsilon}{2}F^{\mu\nu}_{\text{osc}}F_{\mu\nu}+\mathcal{L}_{K}\,,$
(12)

where $\mathcal{L}_{K}$ contains the other kinetic terms of the dark photons,
including a mass of $m_{\text{osc}}$ for $A_{\text{osc}}$. Note that a
discrete symmetry under which $A_{\psi}\to-A_{\psi}$ and $a\to-a$ explicitly
forbids decays of $a$ to a pair of dark photons of the same species, as well
as mixing between $A_{\psi}$ 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 $\psi$ that scatters rapidly with
$A_{\psi}$, to ensure that $A_{\psi}$ always equilibrates into a thermal
distribution with temperature $T^{\prime}$, without destroying the low-energy
spectrum of $A_{\mathrm{osc}}$.222For simplicity, we assume that the entire
dark sector is described by a common temperature $T^{\prime}$. With this
modification, we overcome the main difficulty faced by the model in Ref. [60]:
there are now no low-energy $A_{\psi}$ to interact with the blackbody
distribution of $A_{\mathrm{osc}}$, which was the main source of distortion
for the blackbody spectrum.333While this completes the basic description of
our model, there are some requirements on $\psi$ to keep this model tractable
and avoid tracking the spectra of all of these particles as a function of
time. First, Compton $\psi A_{\psi}\to\psi A_{\psi}$ and double-Compton $\psi
A_{\psi}A_{\psi}\leftrightarrow\psi A_{\psi}$ scattering must be sufficiently
rapid to ensure that $A_{\psi}$ is always described by a simple, blackbody
distribution. This is easily satisfied, as long as $\psi$ has a large coupling
to $A_{\psi}$, and is sufficiently abundant. Second, the process
$aA_{\mathrm{osc}}\leftrightarrow\psi\overline{\psi}$ must not lead to
significant distortion of $A_{\mathrm{osc}}$. This is hard to determine
without tracking the full evolution of the spectrum of $A_{\mathrm{osc}}$, but
it is possible to avoid this entirely by choosing a sufficiently massive
$\psi$ to kinematically forbid $aA_{\mathrm{osc}}\to\psi\overline{\psi}$, and
making $\psi$ asymmetric. We find that $m_{\psi}\sim$30\text{\,}\mathrm{eV}$$
with a coupling to $A_{\psi}$ of $\alpha_{D}=1$, with an asymmetric number
density today of
$n_{\psi,0}=$2.5\text{\times}{10}^{-3}\text{\,}{\mathrm{cm}}^{-3}$$ can
provide sufficiently efficient scattering with $A_{\psi}$. This value of
$m_{\psi}$ is large enough to kinematically forbid
$aA_{\mathrm{osc}}\to\psi\overline{\psi}$. Introducing a similar fermion
$\chi$ with $m_{\chi}\sim m_{\psi}$ that couples to $A_{\mathrm{osc}}$ instead
guarantees that no distortion can occur when $T^{\prime}>m_{\psi},m_{\chi}$,
while allowing distortions to build up once $T^{\prime}<m_{\psi},m_{\chi}$ and
$\psi,\chi$ has frozen out. This also prevents $\psi\overline{\psi}\to
aA_{\mathrm{osc}}$ from being significant at $T^{\prime}<m_{\psi},m_{\chi}$,
since $\overline{\psi}$ annihilates away at $T^{\prime}\sim m_{\psi}$.

We are now ready to check that distortions to both the low-energy
$A_{\text{osc}}$ power-law spectrum and the $A_{\text{osc}}$ blackbody are
small. The most important process, which enters at order $1/f_{a}^{2}$, is
inverse decays, $A_{\psi}A_{\text{osc}}\to a$, which leads to the largest
distortion. The collision integral is given by

$\displaystyle C_{A_{\text{osc}}A_{\psi}\to
a}[f_{\text{osc}}]=-\frac{f_{\text{osc}}}{6}\int\frac{\mathrm{d}^{3}\vec{k}_{A_{\psi}}}{(2\pi)^{3}2\omega_{A_{\psi}}}\int\frac{\mathrm{d}^{3}\vec{p}_{a}}{(2\pi)^{3}2E_{a}}(2\pi)^{4}\delta_{\text{D}}^{4}(k_{\text{osc}}+k_{A_{\psi}}-p_{a})|\mathcal{M}|^{2}_{A_{\text{osc}}A_{\psi}\to
a}f_{A_{\psi}}(\vec{k}_{A_{\psi}})\,.$

We use $k$, $\vec{k}$ and $\omega$ for incoming four-momentum, three-momentum,
and energy; likewise, we have $p$, $\vec{p}$ and $E$ for outgoing states. The
subscript ‘osc’ represents quantities associated with $A_{\text{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

$\displaystyle|\mathcal{M}|^{2}_{A_{\text{osc}}A_{\psi}\to
a}=\frac{m_{a}^{4}}{2f_{a}^{2}}\,,$

while the occupation number of $A_{\psi}$ is given by the blackbody occupation
number, $f_{A_{\psi}}=[\exp(\omega_{A_{\psi}}/T^{\prime})-1]^{-1}$. These
simple expressions allow us to perform the collision integral relatively
easily, to obtain

$\displaystyle C_{A_{\text{osc}}A_{\psi}\to
a}[f_{\text{osc}}]=-\frac{f_{\text{osc}}m_{a}^{4}T^{\prime}}{96\pi\omega_{\text{osc}}f_{a}^{2}}\log\left[1-\exp\left(-x\frac{(1+z_{\star})^{2}}{(1+z)^{2}}\right)\right]\,,$
(13)

where $\omega_{\text{osc}}=xT_{\text{CMB}}$. With this expression, we can now
assess the total distortion to the blackbody and power-law distributions of
$A_{\text{osc}}$, and check that for typical model parameters, these
distortions are small, based on the criterion set out in Eq. (11). For the
blackbody distribution, we want to perform the integral in Eq. (11) starting
from $z_{\min}=z_{\star}$, allowing us to approximate the collision term as

$\displaystyle C_{A_{\text{osc}}A_{\psi}\to
a}[f_{\text{osc}}]=-\frac{f_{\text{osc}}m_{a}^{4}T^{\prime}}{96\pi\omega_{\text{osc}}f_{a}^{2}}\log\left(x\frac{(1+z_{\star})^{2}}{(1+z)^{2}}\right)\,.$
(14)

Integrating this gives

$\displaystyle|\Delta\log f_{\text{osc}}|$
$\displaystyle\sim\frac{2^{5/2}m_{a}^{3/2}\gamma x^{1/2}T_{0}^{3/2}}{96\pi
H_{0}\sqrt{\Omega_{m}}f_{a}^{2}}\sim 1.8\times
10^{-4}\bigg{(}\frac{\gamma}{0.2}\bigg{)}\left(\frac{$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$}{m_{a}}\right)^{3/2}\left(\frac{${10}^{21}\text{\,}\mathrm{s}$}{\tau_{\text{vac}}}\right)\bigg{(}\frac{x}{0.2}\bigg{)}^{1/2}\,,$
(15)

where $\gamma\equiv T_{0}^{\prime}/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. (11) from $z_{\min}=0$ and
$z_{\max}=z_{\star}$. The integral can be approximately done in two parts: i)
$0<1+z<(1+z_{\star})/\sqrt{x}$, where we can expand the expression in Eq. (13)
using $\log(1-\alpha)\simeq-\alpha$, and ii)
$(1+z_{\star})/\sqrt{x}<1+z<1+z_{\star}$, where we can use the approximation
in Eq. (14). The second part of the integral dominates, but the total
contribution to the distortion is given numerically by

$\displaystyle|\Delta\log f_{\text{osc}}|\sim
0.21\bigg{(}\frac{\gamma}{0.2}\bigg{)}\left(\frac{$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$}{m_{a}}\right)^{3/2}\left(\frac{${10}^{21}\text{\,}\mathrm{s}$}{\tau_{\text{vac}}}\right)\left(\frac{3\times
10^{-4}}{x}\right)^{3/4}\,.$ (16)

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\times 10^{-4}$
or $22\text{\,}\mathrm{MHz}$ is $\sim$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 $\tau_{\text{vac}}$, which decreases the
size of the distortion. We therefore conclude that inverse decays do not
significantly distort either component of the $A_{\text{osc}}$
spectrum.444Although elastic scattering processes should be subleading at
order $1/f_{a}^{4}$, they can, in fact, be very rapid for two reasons. First,
processes like $aA_{\text{osc}}\to aA_{\text{osc}}$ receive a large Bose
enhancement from the existing population of $A_{\text{osc}}$, which
significantly increases the scattering rate. Secondly, the matrix element of
$aA_{\text{osc}}\to aA_{\text{osc}}$ contains a divergence. Such divergences
are commonplace even in Standard Model processes, such as $e^{-}Z\to e^{-}Z$
[94]. In the context of cosmological fluids, these divergences are regulated
by thermal corrections to the self-energy of $A_{\psi}$, which impart an
imaginary contribution to the mass of $A_{\psi}$ [94]. We have checked that
even including Bose enhancement and regulating this divergence with loop
contributions from a thermal distribution of $\psi$ [95], all elastic
scattering processes are indeed subdominant to the leading inverse decay
process discussed here.

Our proposed model therefore satisfies the requirement that all distortions to
the $A^{\prime}_{\text{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.

##  Appendix C Experimental constraints

###  C.1 Dark matter lifetime constraint with stimulated decay

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_{\text{dcdm}}$ decays between
recombination and today, _Planck_ 2018 and baryon acoustic oscillation (BAO)
data limit $f_{\text{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_{\text{lim}}=0.0216$ for the proportion of dark matter that decays away by
the present day. For a stimulated emission lifetime given by
$\tau(z)=\tau_{\text{vec}}(1+2f_{A^{\prime}}^{\text{BB}})^{-1}$, which is the
expression we obtain with our fiducial model, we require555This corresponds to
an approximation of the effect of stimulation, we leave a more complete
treatment for future work..

$\displaystyle\int_{0}^{t}\frac{\mathrm{d}t}{\tau_{\text{vac}}}\left(1+\frac{2}{e^{\omega_{A^{\prime}}/T^{\prime}}-1}\right)<f_{\text{lim}}\,.$

Numerically, we find

$\displaystyle\int_{0}^{t}\frac{\mathrm{d}t}{\tau_{\text{vac}}}\left(1+\frac{2}{e^{\omega_{A^{\prime}}/T^{\prime}}-1}\right)\approx\frac{1}{\tau_{\text{vac}}}\times\begin{dcases}$1.5\text{\times}{10}^{18}\text{\,}\mathrm{s}$\left(\frac{2T^{\prime}_{0}}{\omega_{A^{\prime}}}\right)\,,&\frac{\omega_{A^{\prime}}}{2T^{\prime}_{0}}<3.2\,,\\\
H_{0}^{-1}\,,&\frac{\omega_{A^{\prime}}}{2T^{\prime}_{0}}\geq
3.2\,,\end{dcases}$

which leads to a limit on the vacuum decay lifetime of

$\displaystyle\tau_{\text{vac}}>${10}^{19}\text{\,}\mathrm{s}$\times\max\left[2.1\,,\,\,6.3\left(\frac{${10}^{-4}\text{\,}\mathrm{eV}$}{\omega_{A^{\prime}}}\right)\bigg{(}\frac{\gamma}{0.2}\bigg{)}\right]\,.$
(17)

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].

###  C.2 Constraints on relativistic degrees of freedom

We complete our discussion of the particle physics model by discussing its
effect on $N_{\text{eff}}$, and relevant constraints. $A_{\text{osc}}$,
$A_{\psi}$ contribute to $N_{\text{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_{\text{osc}}$ and $A_{\psi}$ before recombination, the
total energy density of these particles divided by the energy density of a
neutrino is

$\displaystyle\Delta N_{\text{eff}}\approx\frac{2\times 3+2\times
4\times\frac{7}{8}}{2\times\frac{7}{8}}\frac{\gamma^{4}}{\left(\frac{4}{11}\right)^{4/3}}=0.04\left(\frac{\gamma}{0.2}\right)^{4}\,.$
(18)

This is to be compared with the _Planck_ measurement of
$N_{\text{eff}}=2.99\pm 0.17$ [2], which limits $\Delta N_{\text{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 $\Delta N_{\text{eff}}\approx
0.01(\gamma/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 $\hat{n}$ at a fixed frequency is directly proportional to the total
probability of conversion $P(\hat{n})$ of dark photons into photons. We begin
by defining the fluctuation of the conversion probability in some particular
direction in the sky $\hat{n}$, defined explicitly by $\delta P(\hat{n})\equiv
P(\hat{n})-\langle P\rangle$, where $\langle P\rangle$ is the sky-averaged
conversion probability. In any particular direction $\hat{n}$, we have [56,
55]

$\displaystyle
P(\hat{n})=\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\int_{0}^{z_{\star}}\mathrm{d}z\,\frac{\delta_{\text{D}}(m_{\gamma}^{2}(\vec{r},z)-m_{A^{\prime}}^{2})}{H(z)(1+z)^{2}}\,,$
(19)

where $m_{\gamma}^{2}(\vec{r},z)$ is the effective plasma mass of photons
along the particular line-of-sight, and $\omega_{0}$ is the present-day,
observed angular frequency of the photons. $\vec{r}\equiv\hat{n}\chi(z)$,
where $\chi(z)$ is the comoving distance traveled by light between $z$ and the
present day. $\omega_{0}(1+z_{\star})$ is the energy of the daughter particle
from the DM decay; photons that have energy $\omega_{0}$ today were produced
by decays at $z_{\star}$, 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_{\gamma}^{2};z)$ of
$m_{\gamma}^{2}$, to obtain [56, 55]

$\displaystyle\langle
P\rangle=\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\int_{0}^{z_{\star}}\mathrm{d}z\int\mathrm{d}m_{\gamma}^{2}\,f_{1}(m_{\gamma}^{2};z)\frac{\delta_{\text{D}}(m_{\gamma}^{2}-m_{A^{\prime}}^{2})}{H(z)(1+z)^{2}}=\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\int_{0}^{z_{\star}}\mathrm{d}z\,\frac{f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)}{H(z)(1+z)^{2}}\,,$
(20)

so that

$\displaystyle\delta
P(\hat{n})=\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\int_{0}^{z_{\star}}\frac{\mathrm{d}z}{H(z)(1+z)^{2}}\left[\delta_{\text{D}}(m_{\gamma}^{2}(\vec{r},z)-m_{A^{\prime}}^{2})-f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)\right]\,.$
(21)

We now decompose the observed conversion probability over the whole sky into
spherical harmonics $Y_{\ell m}(\hat{n})$, with coefficients $a_{\ell m}$
given by

$\displaystyle a_{\ell m}=\frac{1}{\langle
P\rangle}\int\mathrm{d}\hat{n}\,\delta P(\hat{n})Y_{\ell m}^{*}(\hat{n})\,.$

The anisotropy power spectrum $C_{\ell}$ for the conversion probability,
defined as $C_{\ell}\equiv\langle a_{\ell m}^{*}a_{\ell m}\rangle$, can then
be computed as

$\displaystyle C_{\ell}=\frac{1}{\langle
P\rangle^{2}}\int\mathrm{d}\hat{n}\int\mathrm{d}\hat{n}^{\prime}\langle\delta
P(\hat{n})\delta P(\hat{n}^{\prime})\rangle Y_{\ell m}^{*}(\hat{n})Y_{\ell
m}(\hat{n}^{\prime})\,,$

where $\langle\cdots\rangle$ in the integral should be interpreted as an all-
sky average. Note that defined in the following manner, $C_{\ell}$ is
dimensionless, since it describes fluctuations in conversion probability.

To make further progress, we define the quantity

$\displaystyle
Q(\vec{r},z;\vec{r}^{\prime},z^{\prime})\equiv\left<\left[\delta_{\text{D}}(m_{\gamma}^{2}(\vec{r},z)-m_{A^{\prime}}^{2})-f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)\right]\left[\delta_{\text{D}}(m_{\gamma}^{2}(\vec{r}^{\prime},z^{\prime})-m_{A^{\prime}}^{2})-f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z^{\prime})\right]\right>\,,$
(22)

a two-point correlation function along two different lines-of-sight, described
by comoving coordinates $\vec{r}$ and $\vec{r}^{\prime}$, with photons passing
through each point at redshifts $z$ and $z^{\prime}$ respectively. In fact,
homogeneity and isotropy guarantee that this function does not depend on
$\vec{r}$ and $\vec{r}^{\prime}$ separately, but only on
$|\vec{r}-\vec{r}^{\prime}|$. Inserting this into the expression above for
$C_{\ell}$, we find

$\displaystyle C_{\ell}=\frac{1}{\langle
P\rangle^{2}}\left[\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\right]^{2}\int_{0}^{z_{\star}}\mathrm{d}z\,W(z)\int_{0}^{z_{\star}}\mathrm{d}z^{\prime}\,W(z^{\prime})\int\mathrm{d}\hat{n}\int\mathrm{d}\hat{n}^{\prime}\,Q(|\vec{r}-\vec{r}^{\prime}|,z,z^{\prime})Y_{\ell
m}^{*}(\hat{n})Y_{\ell m}(\hat{n}^{\prime})\,,$ (23)

where $W(z)\equiv[H(z)(1+z)^{2}]^{-1}$. We can now write $Q$ in terms of its
Fourier transform over $\vec{r}-\vec{r}^{\prime}$, $\tilde{Q}$, giving:

$\displaystyle Q(|\vec{r}-\vec{r}^{\prime}|,z,z^{\prime})$
$\displaystyle=\int\frac{\mathrm{d}^{3}\vec{k}}{(2\pi)^{3}}e^{i\vec{k}\cdot(\vec{r}-\vec{r}^{\prime})}\tilde{Q}(k,z,z^{\prime})$
$\displaystyle=(4\pi)^{2}\sum_{p,s=0}^{\infty}\sum_{q=-p}^{p}\sum_{t=-s}^{s}\int\frac{\mathrm{d}^{3}\vec{k}}{(2\pi)^{3}}\tilde{Q}(\vec{k},z,z^{\prime})i^{p-s}j_{p}(kr)j_{s}(kr^{\prime})Y_{pq}^{*}(\hat{k})Y_{pq}(\hat{n})Y_{st}(\hat{k})Y_{st}^{*}(\hat{n}^{\prime})$
$\displaystyle=\frac{2}{\pi}\sum_{p=0}^{\infty}\sum_{q=-p}^{p}\int\mathrm{d}k\,k^{2}\tilde{Q}(k,z,z^{\prime})j_{p}(kr)j_{p}(kr^{\prime})Y_{pq}(\hat{n})Y_{pq}^{*}(\hat{n}^{\prime})\,.$

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. (23) and integrating over $\hat{n}$ and
$\hat{n}^{\prime}$, once again exploiting the orthogonality of spherical
harmonics, gives

$\displaystyle C_{\ell}=\frac{1}{\langle
P\rangle^{2}}\left[\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\right]^{2}\frac{2}{\pi}\int_{0}^{z_{\star}}\mathrm{d}z\,W(z)\int_{0}^{z_{\star}}\mathrm{d}z^{\prime}\,W(z^{\prime})\int\mathrm{d}k\,k^{2}\tilde{Q}(k,z,z^{\prime})j_{\ell}(kr)j_{\ell}(kr^{\prime})\,.$

Finally, we can simplify this integral further assuming the Limber
approximation [74, 75, 76], which is a high-$\ell$ expansion that works
particularly well for the multipoles in which we are interested [96] and
allows us to approximate
$k^{2}j_{\ell}(kr)j_{\ell}(kr^{\prime})\approx(\pi/2)\delta_{\text{D}}(k-\ell/r)\delta_{\text{D}}(r-r^{\prime})/r^{2}$.
This finally gives

$\displaystyle C_{\ell}=\frac{1}{\langle
P\rangle^{2}}\left[\frac{\pi\epsilon^{2}m_{A^{\prime}}^{4}}{\omega_{0}}\right]^{2}\int_{0}^{z_{\star}}\frac{\mathrm{d}z}{r^{2}(z)}\,W^{2}(z)H(z)\tilde{Q}(k=\ell/r,z,z)\,.$
(24)

This compact expression makes the calculation of the conversion anisotropy
numerically tractable; in particular, the correlation function $Q$ defined in
Eq. (22) need only be evaluated at points along two different lines-of-sight
which have equal redshifts.

We now need an expression for $\tilde{Q}(k=\ell/r,z,z)$. First, we define
$f_{2}(\rho,m_{\gamma}^{2}(\vec{r}),m_{\gamma}^{2}(\vec{r}^{\prime});z)$ as
the two-point PDF for $m_{\gamma}^{2}$ at two different points, $\vec{r}$ and
$\vec{r}^{\prime}$, at redshift $z$, with
$\rho\equiv|\vec{r}-\vec{r}^{\prime}|$. Note that homogeneity and isotropy
guarantee that $f_{2}$ only depends on $\rho$. With this, the averaging in Eq.
(22) is performed by integrating the quantity with respect to
$m_{\gamma}^{2}(\vec{r})$ and $m_{\gamma}^{2}(\vec{r}^{\prime})$, weighted by
$f_{2}$, which gives

$\displaystyle
Q(|\vec{r}-\vec{r}^{\prime}|,z,z)=f_{2}(|\vec{r}-\vec{r}^{\prime}|,m_{\gamma}^{2}(\vec{r})=m_{A^{\prime}}^{2},m_{\gamma}^{2}(\vec{r}^{\prime})=m_{A^{\prime}}^{2};z)-[f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)]^{2}\,,$
(25)

with

$\displaystyle\tilde{Q}(k,z,z)=4\pi\int\mathrm{d}\rho\,\rho^{2}j_{0}(k\rho)Q(\rho,z,z)$

Given both the one- and two-point PDF of $m_{\gamma}^{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_{\gamma}^{2}$ are given by the one- and two-point PDFs of baryonic
fluctuations $\delta_{\text{b}}$ divided by one and two powers of
$\overline{m_{\gamma}^{2}}(z)$, the mean value of $m_{\gamma}^{2}$ at redshift
$z$, respectively. This is a good assumption for $z\lesssim 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
$\xi_{\text{b}}(|\vec{r}-\vec{r}^{\prime}|;z)$, which is related to the power
spectrum of baryons $P_{\text{b}}$ via

$\displaystyle\xi_{\text{b}}(\rho;z)=\int\frac{\mathrm{d}^{3}\vec{k}}{(2\pi)^{3}}j_{0}(k\rho)P_{\text{b}}(k;z)\,,$

with $\sigma_{\text{b}}^{2}(z)=\xi_{\text{b}}(0;z)$ being the variance of
baryonic fluctuations. For normally distributed fluctuations,

$\displaystyle f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)$
$\displaystyle=\frac{1}{\overline{m_{\gamma}^{2}}(z)\sqrt{2\pi\sigma_{\text{b}}^{2}(z)}}\exp\left[-\frac{g(m_{A^{\prime}}^{2};z)^{2}}{2\sigma_{\text{b}}^{2}(z)}\right]\,,$
$\displaystyle
f_{2}(\rho,m_{\gamma}^{2}(\vec{r})=m_{A^{\prime}}^{2},m_{\gamma}^{2}(\vec{r}^{\prime})=m_{A^{\prime}}^{2};z)$
$\displaystyle=\frac{1}{2\pi\overline{m_{\gamma}^{2}}^{2}(z)\sqrt{\sigma_{\text{b}}^{4}(z)-\xi_{\text{b}}^{2}(\rho;z)}}\exp\left(-\frac{g(m_{A^{\prime}}^{2};z)^{2}}{\sigma_{\text{b}}^{2}(z)+\xi_{\text{b}}(\rho;z)}\right)\,,$
(26)

where $g(m_{A^{\prime}}^{2};z)\equiv
m_{A^{\prime}}^{2}/\overline{m_{\gamma}^{2}}(z)-1$, while for log-normally
distributed fluctuations [97, 55],

$\displaystyle f_{1}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)$
$\displaystyle=\frac{1}{\overline{m_{\gamma}^{2}}(z)\sqrt{2\pi\Sigma^{2}(z)}}\exp\left[-\frac{L^{2}(m_{A^{\prime}}^{2};z)}{2\Sigma^{2}(z)}\right]\frac{1}{1+g(m_{A^{\prime}}^{2};z)}\,,$
$\displaystyle
f_{2}(\rho,m_{\gamma}^{2}(\vec{r})=m_{A^{\prime}}^{2},m_{\gamma}^{2}(\vec{r}^{\prime})=m_{A^{\prime}}^{2};z)$
$\displaystyle=\frac{1}{2\pi\overline{m_{\gamma}^{2}}^{2}(z)\sqrt{\Sigma^{4}(z)-X^{2}(\rho;z)}}\exp\left(-\frac{L^{2}(m_{A^{\prime}}^{2};z)}{\Sigma^{2}(z)+X(\rho;z)}\right)\frac{1}{[1+g(m_{A^{\prime}}^{2};z)]^{2}}\,,$
(27)

where $\Sigma^{2}(z)=\log[1+\sigma_{\text{b}}^{2}(z)]$,
$X(\rho;z)=\log[1+\xi_{\text{b}}(\rho;z)]$, and
$L(m_{A^{\prime}}^{2};z)=\log[1+g(m_{A^{\prime}}^{2};z)]+\Sigma^{2}(z)/2$.

In the limit that $\rho\to\infty$, we obtain $\xi_{\text{b}}\to 0$ and
$f_{2}\to f_{1}^{2}$ 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. (25) and (24) give us the following general expression for $C_{\ell}$,
which is the main result of this section:

$\displaystyle C_{\ell}=\frac{1}{\langle
P\rangle^{2}}\int_{0}^{z_{\star}}\frac{\mathrm{d}z}{r^{2}(z)}H(z)\left[\frac{\mathrm{d}\langle
P\rangle}{\mathrm{d}z}\right]^{2}\int\mathrm{d}\rho\,4\pi\rho^{2}j_{0}(\ell\rho/r)\left[\frac{f_{2}(\rho,m_{\gamma}^{2}(\vec{r})=m_{A^{\prime}}^{2},m_{\gamma}^{2}(\vec{r}^{\prime})=m_{A^{\prime}}^{2};z)}{f_{1}^{2}(m_{\gamma}^{2}=m_{A^{\prime}}^{2};z)}-1\right]\,,$
(28)

an expression that we can evaluate numerically for either choice of baryon
fluctuation PDFs. Note that the integrand over $\rho$ is finite as $\rho\to
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
$\ell(\ell+1)C_{\ell}\langle T\rangle^{2}$, where $\langle T\rangle$ is the
sky-averaged brightness temperature.

For small angular scales ($\ell\gtrsim 3000$), upper limits on the anisotropy
of the ERB have been obtained by the Very Large Array (VLA) at
$4.86\text{\,}\mathrm{GHz}$ [40] and $8.4\text{\,}\mathrm{GHz}$ [41], as well
as the Australia Telescope Compact Array (ATCA) at $8.7\text{\,}\mathrm{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 $\Delta T/T_{\text{exc}}\lesssim 10^{-2}$ for the ERB [39]. More recently,
LOFAR and TGSS observations have been used to measure the anisotropy power
spectrum at $\sim$$140\text{\,}\mathrm{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.

Figure 4: Predicted anisotropy power spectrum with normally (solid black line,
below) and log-normally (solid black line, above) distributed baryon
fluctuations. The model parameters are
$m_{A^{\prime}}=$3\text{\times}{10}^{-14}\text{\,}\mathrm{eV}$$ and
$\omega_{A^{\prime}}=$2\text{\times}{10}^{-4}\text{\,}\mathrm{eV}$$. Upper
limits from VLA at $4.86\text{\,}\mathrm{GHz}$ [40] (purple triangles),
$8.4\text{\,}\mathrm{GHz}$ [41] (orange triangles) and ATCA at
$8.7\text{\,}\mathrm{GHz}$ [42] (green triangle) are shown, and have been
rescaled to $T_{\text{exc}}=$255\text{\,}\mathrm{K}$$, the expected value at
$140\text{\,}\mathrm{MHz}$. Representative data points of the power spectrum
measured by LOFAR are shown in red [44].

Fig. 4 shows the predicted anisotropy power spectrum,
$[\ell(\ell+1)/(2\pi)]C_{\ell}\langle T\rangle^{2}$, as a function of $\ell$,
obtained by integrating Eq. (28) numerically for both normally (lower line)
and log-normally (upper line) distributed baryon fluctuations.666The
theoretical prediction from our model only has a weak dependence on the
observational frequency, coming from the integration limit $z_{\star}$; we
have only included the largest fluctuations here. The range between the two
lines, shaded in blue, gives some indication of the uncertainty of our
prediction. Representative values of
$m_{A^{\prime}}=$3\text{\times}{10}^{-14}\text{\,}\mathrm{eV}$$ and
$\omega_{A^{\prime}}=$2\text{\times}{10}^{-4}\text{\,}\mathrm{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 $49\text{\,}\mathrm{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\text{\,}\mathrm{MHz}$ [44]. We take $\langle
T\rangle=$255\text{\,}\mathrm{K}$$, the expected radio temperature at
$140\text{\,}\mathrm{MHz}$ based on the power-law fit in Eq. (1) of the
Letter. The anisotropy power reported by the high frequency measurements at
frequency $\nu_{\text{high}}$ has been rescaled by
$($140\text{\,}\mathrm{MHz}$/\nu_{\text{high}})^{2\beta}$, where $\beta=2.6$;
this assumes that the relative size of fluctuations $\delta T/T$ is
independent of frequency, and can be rescaled by $\langle T\rangle^{2}$ at
each frequency. Note that the observed upper limits in the
$49\text{\,}\mathrm{GHz}$ range are not necessarily in tension with the
$\ell\gtrsim 10^{3}$ results at $140\text{\,}\mathrm{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}\lesssim\ell\lesssim 6\times 10^{4}$, finding a power of
$1{10}^{2}\text{\,}{\mathrm{K}}^{2}$ across this range; this represents a
temperature fluctuation of $0.004\lesssim\Delta T/T_{\text{exc}}\lesssim
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 $\ell\lesssim 800$ at
$140\text{\,}\mathrm{MHz}$. However, as mentioned in the Letter, lower-$\ell$
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 $49\text{\,}\mathrm{GHz}$ radio
data may be of interest in future work.

##  Appendix E Data analysis and statistical methodology

###  E.1 Likelihood and posterior sampling

We assume the likelihood of the observed brightness temperatures
$T_{\mathrm{obs},i}(\nu_{i})$ given a model for the emission $T(\nu;\theta)$
to be Gaussian, where $i$ indexes individual data points and
$\theta=\\{m_{a},m_{A^{\prime}},\tau_{\mathrm{vac}},\epsilon,T_{0}^{\prime},T_{0}\\}$
are the parameters of interest describing the fiducial new-physics model.
Assuming uncorrelated measurements, the total likelihood for a dataset is
given by
$p\left(T_{\mathrm{obs}}\mid\theta\right)=\prod_{i}\mathcal{N}\left(T_{\mathrm{obs},i}(\nu_{i})\mid
T(\nu_{i};\theta),\sigma_{T,i}\right)$, where $\sigma_{T,i}$ are the
corresponding standard deviation uncertainties describing the noise model.

Given a prior $p(\theta)$ as defined in Sec. E.2 below, we infer an
approximation for the posterior distribution $p\left(\theta\mid
T_{\mathrm{obs}}\right)=p\left(T_{\mathrm{obs}}\mid\theta\right)p(\theta)/\mathcal{Z}$,
where $\mathcal{Z}$ is the marginal evidence
$\mathcal{Z}=\int\mathrm{d}\theta\,p(T_{\mathrm{obs}}\mid\theta)\,p(\theta)$.
This is done through Monte Carlo nested sampling implemented in dynesty [72].
$1000$ live points are used to model the posterior and sampling is performed
until the estimated expected contribution to the log-evidence is less than
$\Delta\log\mathcal{Z}=0.05$, with an otherwise default configuration of the
static nested sampler.

###  E.2 Constrained prior definition

A baseline prior is assumed such that $\log_{10}m_{a}\sim\mathcal{U}(-7,-1)$,
$\log_{10}m_{A^{\prime}}\sim\mathcal{U}(-17,-10)$,
$T_{0}^{\prime}/T_{0}\sim\mathcal{U}(0,0.4)$,
$\log_{10}\tau_{\text{vac}}\sim\mathcal{U}(16,25)$,
$\log_{10}\epsilon\sim\mathcal{U}(-11,-4)$, and
$T_{0}\sim\mathcal{N}(2.7255,0.00086)$, with the prior on $T_{0}$ motivated by
an analysis of COBE/FIRAS CMB data in the $2.27$–$21.33$ cm-1 frequency range
[1, 8], and the prior on $T_{0}^{\prime}/T_{0}$ motivated by the constraints
derived in Sec. C.2. Furthermore, for our fiducial analysis we wish to
incorporate external constraints on the considered parameter space as
mentioned in the main body of this _Letter_. Practical implementations of
nested sampling require that the prior distribution be provided as a
transformation from the unit hypercube to the target prior density. Although
such a transformation is readily available for commonly used prior
specifications _e.g._ , uniform- or Gaussian-distributed, transformations from
an arbitrary prior such as that respecting the constraints laid out in the
main _Letter_ through combinations of multiple parameters are not always
readily available or analytically tractable.

In order to construct an implicit prior for our purposes we follow a procedure
similar to that outlined in Ref. [98] and use normalizing flows [99, 100],
which use a series of bijections parameterized by neural networks with
tractable Jacobians in order to define transformations from a simple base
distribution to complex, expressive target distributions. Specifically, we
generate a large number $10^{6}$ of samples from the baseline prior, remove
the samples incompatible with the external constraints, and use a normalizing
flow consisting of 8 neural spline flows [101] to learn a transformation from
a standard Gaussian $\mathcal{N}(0,\mathbbm{1})$ into the constrained target
prior density. Samples $u$ from the unit hypercube proposal can then be
transformed to those corresponding to a standard Gaussian through an affine
transformation $u\rightarrow 0+\mathbbm{1}\cdot F^{-1}_{\mathcal{N}}(u)$,
where $F^{-1}_{\mathcal{N}}$ is the inverse cumulative distribution function
of the Gaussian distribution. These samples can then be further transformed to
those on the constrained prior density by passing them through the learned
normalizing flow.

###  E.3 Prior-predictive check

Generally in Bayesian analyses and in particular when non-trivial priors are
utilized such as in this work, it is useful to conduct prior-predictive checks
in order to instill confidence that the prior does not bias the analysis
towards “favorable” regions of the parameter space. In that case, formally
low-likelihood parameter points could still correspond to high posterior
density regions.

Figure 5 shows the 68 and 95% contours of the prior induced on the excess
temperature above the modeled extragalactic background, defined as
$T_{\text{eg}}=$0.23\text{\,}\mathrm{K}$(\nu/$\text{\,}\mathrm{GHz}$)^{-2.7}$
[16, 39]. It can be seen that the induced prior covers a large range of
possible excess temperatures and does not bias the current analysis in favour
of the observed data points (shown in black and red).

Figure 5: 68% (dark green) and 95% (light green) containment of the induced
prior on the excess temperature above the modeled extragalactic background,
$T_{\text{eg}}=$0.23\text{\,}\mathrm{K}$(\nu/$\text{\,}\mathrm{GHz}$)^{-2.7}$
[16]. Measured data points from Ref. [3] are shown for comparison. The induced
prior covers a large range of possible excess temperatures and does not bias
the analysis in favour of the observed data points.

##  Appendix F Extended results

###  F.1 Parameter posteriors

Figure 6 shows the individual and pair-wise marginal posteriors on all the
modeled parameters obtained for the baseline analysis.

Figure 6: Joint and individual marginal posterior obtained in the fiducial
model, same format as Fig. 2.

###  F.2 Systematic variations on the analysis

In this section we consider various systematic variations of our baseline
analysis. In particular we consider:

  1. 1.

results obtained using the data from Ref. [9] instead of those from Ref. [3].
Ref. [9] includes recent data points from the LWA1 Low Frequency Sky Survey
(LLFSS), which spans a frequency range of $4080\text{\,}\mathrm{MHz}$. These
new data points are consistent with the power-law fit in Ref. [3]. Ref. [9]
also independently reanalyzed the data from ARCADE 2 and the other experiments
used in Ref. [3], finding broad agreement, but with small changes to the
central values and more significant changes to the error bars, despite using a
similar data analysis method. Ref. [9] therefore provides a useful systematics
check on how the ERB spectrum is extracted;

  2. 2.

results obtained by fitting only the ARCADE 2 data $\gtrsim 3$ GHz and
ignoring the lower-frequency radio measurements. In this case, we neglect
stimulated emission, since this is not necessary to produce a good fit,
simplifying the model and reducing to that first introduced in Ref. [60].
ARCADE 2 is the only experiment in our list that is designed primarily with
absolute zero-level calibration in mind, in order to produce an accurate
determination of the sky-brightness [10], which warrants a separate
investigation; and

  3. 3.

results obtained without the modeled contribution from unresolved
extragalactic radio sources [16, 39]. This should concretely demonstrate that
our model does not rely on unresolved extragalactic radio sources for a good
fit, which can be anticipated from the fact that these sources are at least 3
times less bright than the measured $T_{\text{exc}}$.

Figure 7 shows the excess temperature posterior for all these systematic
variations, while Tab. 1 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. 1 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”).

###  F.3 Power-law fit and model comparison

As a point of comparison, we also show results using a variant of the
commonly-used power-law ansatz (equivalent to Eq. 1) for the excess
temperature over $T_{0}$,

$T(\nu)=T_{0}\left[1+A_{\mathrm{GHz}}\left(\frac{\nu}{1\,\mathrm{GHz}}\right)^{\beta}\right]$
(29)

which is parameterized through the CMB black-body temperature $T_{0}$, excess
temperature at 1 GHz $A_{\mathrm{GHz}}$, and spectral index $\beta$. The
results of these fits for the Fixsen _et al_ and Dowell & Taylor datasets are
shown in the left and right columns of Fig. 8, respectively. As previously
reported in the literature, these parameterizations provide a formally good
fit to the low-frequency radio data [38, 9]. Posterior summaries for these
fits are given in Tab. 2 in the same format as Tab. 1.

We show a quantitative comparison between the power law and new physics models
as follows. For a given model $\mathcal{M}$ (in our case corresponding to
either the new-physics model or the power-law ansatz) the model evidence, also
known as the marginal likelihood, is defined for data $x$ and parameters
$\theta$ as
$\mathcal{Z}_{\mathcal{M}}\equiv\int\mathrm{d}\theta\,p_{\mathcal{M}}(x\mid\theta)p_{\mathcal{M}}(\theta)$
where $p_{\mathcal{M}}(x\mid\theta)$ and $p_{\mathcal{M}}(\theta)$ are the
model likelihood and prior, respectively. The model evidence is a measure of
compatibility between the model and data, and implicitly penalizes excessive
model complexity—this is relevant in our case, since the new physics and
power-law models have different number of parameters (3 and 6 respectively).
The ratio of evidences between two models $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$ is known as the Bayes factor, and can be formally used in a
model comparison setting to quantify how much better of a fit
$\mathcal{M}_{1}$ is compared to $\mathcal{M}_{2}$.

We show the log-Bayes factors in favor of the power-law model in the last
column of Tab. 2. The power-law model tends to provide a formally better fit
using both the Fixsen _et al_ and Dowell & Taylor datasets, corresponding to
log-Bayes factors in favor of the power-law model of
$\Delta\log\mathcal{Z}=1.89$ and $\Delta\log\mathcal{Z}=0.76$ respectively.
This reflects the fact that a power law with $\beta\approx-2.6$ is a slightly
better fit to the data points, rather than $\beta\approx-5/2$ as predicted by
our model. We note however that given significant differences between models
considered, a principled model comparison is difficult. Specifically,
$\mathcal{O}(1)$ Bayes factors should not be taken at face value, since they
can depend sensitively on chosen prior specifications and ranges; see Ref.
[102] for a detailed discussion of this point. Finally, we note that even
though synchrotron emission from unknown sources can plausibly produce a power
law close to $\beta\approx-2.6$ (see App. G), this is highly dependent on the
electron distribution within these individual sources as a function of source
brightness and counts; it is therefore not immediately obvious that any future
model ascribing $T_{\text{exc}}$ to currently unknown synchrotron sources
would provide a better fit than our model does.

Figure 7: Systematic variations on the analysis, same format as Fig. 2.

Figure 8: Excess temperature for the power-law fit, same format as Fig. 2 Config. |  Summary |  $\log_{10}m_{a}$ |  $\log_{10}m_{A^{\prime}}$ |  $T_{0}^{\prime}/T_{0}$ |  $\log_{10}\tau_{a}$ |  $\log_{10}\epsilon$ |  $T_{0}$
---|---|---|---|---|---|---|---
- |  - |  $\log_{10}[\mathrm{eV}]$ |  $\log_{10}[\mathrm{eV}]$ |  - |  $\log_{10}[\mathrm{s}]$ |  - |  [K]
Fixsen (Baseline) |  Marginals |  $-3.66^{+0.10}_{-0.26}$ |  $-13.57^{+0.32}_{-0.68}$ |  $0.22^{+0.05}_{-0.12}$ |  $21.15^{+0.47}_{-1.10}$ |  $-6.98^{+0.25}_{-0.60}$ |  $2.7258^{+0.0003}_{-0.0007}$
|  HPDI |  $[-3.97,-3.37]$ |  $[-14.52,-13.04]$ |  $[0.08,0.32]$ |  $[19.58,21.98]$ |  $[-7.65,-6.30]$ |  $[2.7251,2.7266]$
No EG model |  Marginals |  $-3.96^{+0.11}_{-0.26}$ |  $-13.22^{+0.28}_{-0.75}$ |  $0.22^{+0.05}_{-0.11}$ |  $21.48^{+0.58}_{-1.27}$ |  $-7.23^{+0.27}_{-0.71}$ |  $2.7259^{+0.0003}_{-0.0007}$
|  HPDI |  $[-4.27,-3.75]$ |  $[-13.85,-12.47]$ |  $[0.09,0.33]$ |  $[19.96,22.71]$ |  $[-7.87,-6.43]$ |  $[2.7252,2.7267]$
Dowell |  Marginals |  $-4.12^{+0.15}_{-0.33}$ |  $-13.28^{+0.30}_{-0.74}$ |  $0.22^{+0.05}_{-0.11}$ |  $21.61^{+0.60}_{-1.26}$ |  $-7.25^{+0.27}_{-0.78}$ |  $2.7256^{+0.0003}_{-0.0007}$
|  HPDI |  $[-4.54,-3.76]$ |  $[-14.01,-12.50]$ |  $[0.15,0.40]$ |  $[20.10,22.98]$ |  $[-8.06,-6.50]$ |  $[2.7249,2.7264]$
Fixsen, no stim. ARCADE only |  Marginals |  $-3.38^{+0.13}_{-0.30}$ |  $-13.59^{+0.32}_{-0.64}$ |  - |  $20.33^{+0.34}_{-0.65}$ |  $-6.61^{+0.20}_{-0.45}$ |  $2.7256^{+0.0002}_{-0.0007}$
|  HPDI |  $[-3.78,-3.10]$ |  $[-14.54,-13.06]$ |  - |  $[19.36,20.80]$ |  $[-7.15,-6.08]$ |  $[2.7249,2.7264]$
Table 1: Posterior summaries for all the systematic variations as compared to our baseline analysis. Config. |  Summary |  $A_{\mathrm{GHz}}$ |  $\beta$ |  $T_{0}$ |  $\Delta\log\mathcal{Z}$
---|---|---|---|---|---
- |  - |  - |  - |  [K] |  [nat]
Fixsen Power law |  Marginals |  $0.43^{+0.01}_{-0.04}$ |  $-2.58^{+0.01}_{-0.04}$ |  $2.7258^{+0.0003}_{-0.0008}$ |  1.89
|  HPDI |  $[0.40,0.47]$ |  $[-2.62,-2.55]$ |  $[2.7250,2.7266]$ |
Dowell Power law |  Marginals |  $0.54^{+0.03}_{-0.06}$ |  $-2.58^{+0.02}_{-0.05}$ |  $2.7255^{+0.0003}_{-0.0008}$ |  0.76
|  HPDI |  $[0.47,0.61]$ |  $[-2.63,-2.53]$ |  $[2.7247,2.7264]$ |
Table 2: Posterior summaries for the power-law fits. The last column
additionally shows the log-Bayes factor $\Delta\log\mathcal{Z}$ in favor of
the power-law fit, with positive values corresponding to a preference for the
power-law model.

##  Appendix G A recap of synchrotron radiation

As extensively discussed in the Letter, many astrophysical solutions to the
ERB advocate for the presence of synchrotron radiation from unresolved
astrophysical sources. In this appendix, we briefly review synchrotron
emission with a particular focus on its spectral features. The basic physics
behind synchrotron radiation is straightforward: a charged particle
accelerated in a magnetic field will radiate. Nonrelativistic particles emit
cyclotron radiation, with a frequency given by the cyclotron frequency
$\omega_{\rm cyc}=qB/m$, where $q$ is the particle charge, $B$ the magnetic
field in which it moves and $m$ its mass. Relativistic particles emit what we
call synchrotron radiation: the frequency spectrum is far richer, covering a
range well away from the cyclotron frequency.

Taking into account strong relativistic effects (such as forward beaming), one
finds that the characteristic frequency of the synchrotron radiation of a
single charged particle is $\omega_{\rm syn}\sim\gamma^{2}\omega_{\rm cyc}$
[103]. The spectrum of the radiation is broad: $\omega<\omega_{\rm syn}$ the
spectrum falls off like a power law $\propto\omega^{1/3}$, while for
$\omega>\omega_{\rm syn}$, the spectrum falls off exponentially.

This is, however, very different from what we observe in astrophysical
systems. The reason is that we need to consider what happens to the total
synchrotron radiation emitted by a population of charged particles, with a
given energy distribution. Let us consider an ensemble of charged particles
with number density distribution per unit energy following a power law
$\mathrm{d}N/\mathrm{d}E\propto E^{p}$. The power per unit energy emitted by
this ensemble of charges reads

$\frac{\mathrm{d}P}{\mathrm{d}E}\propto\frac{\mathrm{d}N}{\mathrm{d}E}P_{E}\propto
E^{p}E^{2}B^{2},$

where $P_{E}$ is the power emitted by a single particle, given by the Larmor
formula $P_{E}\propto B^{2}\gamma^{2}\propto B^{2}E^{2}$ (we drop here all the
irrelevant constants given that we are only interested in the final spectral
shape of the synchrotron radiation power). Let us now assume that the single-
charge synchrotron spectrum is peaked enough to assume that all the power is
emitted at the synchrotron frequency, $\omega_{\rm syn}$. This implies a one-
to-one relation between the observed frequency and the energy of the single
particle $\nu\simeq\omega_{\rm syn}/(2\pi)\propto E^{2}B$. It then follows

$\frac{\mathrm{d}P}{\mathrm{d}\nu}=\frac{\mathrm{d}P}{\mathrm{d}E}\frac{\mathrm{d}E}{\mathrm{d}\nu}\propto
E^{p+2}B\times\nu^{-1/2}B^{-1/2}\propto
B\Big{(}\frac{\nu}{B}\Big{)}^{\frac{1+p}{2}},$

and therefore the spectral index reads $\alpha=(1+p)/2$. In terms of the
spectral index of brightness temperature, $\beta$, this corresponds to
$\beta=(p-3)/2$, since the brightness temperature is related to the intensity
via $T_{\nu}\propto I_{\nu}/\nu^{2}$. Therefore, the spectral index for
synchrotron radiation is intrinsically related to the energy distribution of
the charged particles which generate the radiation. In particular, a power law
is expected if the energy distribution is also a power law as assumed above.

One well-known example for which a power-law distribution is expected comes
from particle acceleration by astrophysical shocks [104, 105, 106, 107, 108].
For example, if one considers shock waves which propagate in the interstellar
medium outside a supernova remnant, it can be shown that at linear level
$p=-2$ [104, 105], which implies $\alpha=-0.5$ and $\beta=-2.5$. The spectral
index predicted by the simplest scenario of shock acceleration coincides with
the spectral index expected in our model, and is therefore slightly different
than the power-law best-fit obtained in Ref. [3]. Nevertheless, steeper
spectra are possible taking into account non-linear effects in the theory of
diffusive shock acceleration [109]. As a matter of fact, these effects are
crucial to match observations of the supernova remnant Cassiopeia A, which
exhibits a steeper spectral index, $\alpha\simeq-0.77$ or $\beta\simeq-2.77$
[110].