# Was Entropy Conserved between BBN and Recombination?

Alexander C. Sobotka<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of North Carolina at Chapel Hill,
Phillips Hall CB3255, Chapel Hill, North Carolina 27599, USA  Adrienne L.
Erickcek<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of North Carolina at Chapel Hill,
Phillips Hall CB3255, Chapel Hill, North Carolina 27599, USA  Tristan L. Smith
<EMAIL_ADDRESS>Department of Physics and Astronomy, Swarthmore
College, Swarthmore, PA 19081, USA

(August 30, 2024)

###### Abstract

We test the assumption of entropy conservation between Big Bang
nucleosynthesis (BBN) and recombination by considering a massive particle that
decays into a mixture of photons and other relativistic species. We employ
Planck temperature and polarization anisotropies, COBE/FIRAS spectral
distortion bounds, and the observed primordial deuterium abundance to
constrain these decay scenarios. If between 56% and 71% of the decaying
particle’s energy is transferred to photons, then $N_{\text{eff}}$ at
recombination is minimally altered, and Planck data alone allows for
significant entropy injection. If photons are injected by the decay, the
addition of spectral distortion bounds restricts the decay rate of the
particle to be
$\Gamma_{Y}>$1.91\text{\times}{10}^{-6}\text{\,}{\mathrm{s}}^{-1}$$ at 95%
confidence level (CL). We find that constraints on the energy density of the
decaying particle are significantly enhanced by the inclusion of bounds on the
primordial deuterium abundance, allowing the particle to contribute at most
2.35% (95% CL) of the energy density of the universe before decaying.

##  I Introduction

An underlying assumption of the standard cosmological model is that the
comoving entropy density of relativistic species is conserved. However, many
alternative scenarios include entropy injection into the thermal bath of
relativistic particles. Extreme entropy injection is needed to completely
repopulate the relativistic bath after an early matter-dominated era [1, 2, 3,
4]. Minor entropy injections that only make small alterations to the radiation
density have been proposed to alter the relic dark matter abundance after
thermal freeze-out [5, 6], change the pre-recombination expansion history [7,
8, 9], and even relax the Hubble tension [10, 11], which is a discrepancy
between the present-day expansion rate inferred from local measurements and
that inferred from the cosmic microwave background (CMB) [12, 13, 14].

Significant entropy injections, such as those used to transition out of an
early matter-dominated era, must be complete before neutrino decoupling in
order to avoid altering the light-element abundances predicted by standard Big
Bang nucleosynthesis (BBN) and impacting the CMB anisotropies. Therefore,
significant entropy injection is generally constrained to have completed at
temperatures hotter than about 4 MeV [15, 16].

Minor entropy injections can occur during or after BBN, but are heavily
constrained by their influence on the expansion rate, baryon-to-photon ratio,
and photodisintegration of light nuclei [17, 18, 19, 20, 21, 22, 23]. Such
scenarios have been thoroughly considered in the context of decaying axion-
like particles [24, 25, 26, 27] for which observations of primordial light-
element abundances and the CMB provide stringent constraints. However, these
investigations focus on specific axion models and therefore cannot place
comprehensive constraints on general entropy injection.

In this work, we test the assumption of entropy conservation between BBN and
recombination by considering a generic massive particle that decays into a
mixture of photons and other relativistic species (e.g. dark radiation). We
explore the bounds that can be placed on the particle decay rate,
$\Gamma_{Y}$, and the contribution that the particle makes to the energy
composition of the universe before it decays. Since these constraints depend
on what relativistic species this massive particle decays into, we also
explore what bounds can be placed on $f_{\gamma}$, the fraction of the
decaying particle’s energy that is transferred to photons.

We restrict our analysis to particles with rest-energies less than the binding
energy of deuterium ($m<$2.2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$).
Photons generated by decaying particles with larger masses will
photodisintegrate light nuclei and alter the abundances of primordial elements
such as deuterium and helium. Photodisintegration has been shown to place very
stringent constraints on such entropy injections [22, 23, 28, 29, 30, 31, 27].

For a particle that is in kinetic equilibrium with the Standard Model to not
contribute to the photodisintegration of nuclei and be non-relativistic when
decaying between BBN and recombination, it must have a mass in the range of
$$1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$\lesssim
m\lesssim$2.2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. However, this tight
mass window can be relaxed if the decaying particle is part of a hidden sector
that is not thermally coupled to the Standard Model [32, 33, 34, 35, 36].
Hidden sector particles could decay into dark radiation [37, 38] or Standard
Model particles [39, 40, 41].

If the hidden sector was originally in kinetic equilibrium with the visible
sector of Standard Model particles and then decoupled while all the Standard
Model was still relativistic, the temperature of the hidden sector would be
colder than that of the visible sector. Then, a hidden sector particle with a
mass of $m<$1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ would be non-
relativistic by the end of BBN and be light enough such that any entropy
injection from the decay would not result in photodisintegration. It is also
possible that the hidden sector was never in kinetic equilibrium with the
Standard Model. We assume the hidden sector to be sufficiently cold such that
the decaying particle is non-relativistic throughout BBN.

To constrain a generic decaying particle we employ Planck measurements of the
CMB temperature and polarization anisotropies [42], the most recent
measurement of the primordial deuterium abundance [43], and the COBE/FIRAS
limit on spectral distortions [44, 45]. Entropy injection before recombination
will alter the effective number of relativistic species, $N_{\text{eff}}$,
resulting in a change in the expansion rate. Alterations to the pre-
recombination expansion history may be constrained by observations of CMB
anisotropies [46]. Specifically, altering the expansion rate via
$N_{\text{eff}}$ affects the amplitude of small scale perturbations [47, 48]
and introduces a phase shift to the baryon-photon acoustic oscillations [49,
50, 51]. Observations of the primordial deuterium abundance are an excellent
way to provide additional constraining power on nonstandard models, with
measurements now reaching 1% precision [43]. Finally, scenarios that inject
photons before recombination may be constrained by CMB spectral distortions
[52, 53, 54].

We combine the constraining power of CMB anisotropy data, deuterium
measurements, and CMB spectral distortions by employing a Markov-Chain Monte
Carlo (MCMC) analysis. In doing so, we derive bounds on the amount and type of
radiation that can be injected by a decaying hidden sector particle.

In many ways, our investigation generalizes the work of Millea  _et al._ [25]
who considered non-photodisintegrating injections from axion-like particles
that decay solely into photons between BBN and recombination. The specific
coupling of these axion-like particles to photons restricts the possible
parameter space that can be explored. We remain agnostic regarding a specific
particle model. In doing so, this work has three key distinctions from Millea
_et al._ [25]: 1) we are able to derive a broadly applicable bound on the
maximum level of entropy injection allowed between BBN and recombination, 2)
we can consider a particle that decays into photons and dark radiation, and 3)
our full analysis of CMB anisotropies allows us to investigate the effects
that entropy injection has on standard cosmological parameters, including
$H_{0}$. We determine that minor entropy injection from a hidden sector
particle is not a promising avenue for addressing the Hubble tension. A more
exotic change to the Universe’s radiation content, such as the inclusion of
strongly interacting radiation [11], is needed to significantly increase the
value of $H_{0}$ inferred from observations of the CMB.

The structure of this paper is as follows. In Section II, we describe the
model for our decaying particle. In Section III, we explore the primary
effects that the decay has on $N_{\text{eff}}$ and the primordial abundances
of helium and deuterium, and we review the expected constraints from spectral
distortions. We outline the choice of priors and likelihood functions for our
MCMC analysis in Section IV, and Section V discusses the results of these
analyses. We conclude with a summary in Section VI, and we include an Appendix
that contains a derivation of our analytical model and technical details of
our modifications to the public code known as CLASS (A), a derivation for the
post-decay $N_{\text{eff}}$ (B), and supplemental MCMC results (C).

##  II Decaying Particle Model

We consider a sub-dominant hidden sector $Y$ particle that decays into photons
and other relativistic particles (e.g. dark radiation) sometime between BBN
and recombination. The equations governing the background evolution for the
energy density of the $Y$ particle ($\rho_{Y}$), photons ($\rho_{\gamma}$),
and other ultra-relativistic species ($\rho_{ur}$), are

$\displaystyle\frac{d}{dt}\rho_{Y}+3H\rho_{Y}$
$\displaystyle=-\Gamma_{Y}\rho_{Y}\,,$ (1)
$\displaystyle\frac{d}{dt}\rho_{\gamma}+4H\rho_{\gamma}$
$\displaystyle=f_{\gamma}\Gamma_{Y}\rho_{Y}\,,$ (2)
$\displaystyle\frac{d}{dt}\rho_{ur}+4H\rho_{ur}$
$\displaystyle=(1-f_{\gamma})\Gamma_{Y}\rho_{Y}\,,$ (3)

where $\Gamma_{Y}$ is the decay rate of the $Y$ particle, $H\equiv\dot{a}/a$
is the Hubble parameter (with the over-dot referring to a proper time
derivative), and $f_{\gamma}$ dictates what fraction of the decaying
particle’s energy is transferred to photons.111More precisely, the terms on
the right-hand side of these equations are proportional to $m_{Y}n_{Y}$ [55],
but $m_{Y}n_{Y}=\rho_{Y}$ if the $Y$ particle has negligible kinetic energy.
We define the reheat temperature, $T_{\mathrm{RH}}$, by

$\Gamma_{Y}\equiv\sqrt{\frac{8\pi}{3m^{2}_{\text{pl}}}\left(\frac{\pi^{2}}{30}g_{*}T^{4}_{\mathrm{RH}}\right)}\,,$
(4)

where $m_{\text{pl}}$ is the Planck mass and $g_{*}$ is the effective degrees
of freedom of Standard Model relativistic particles. Since we are only
considering reheat temperatures after BBN
($T_{\mathrm{RH}}<$0.01\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$),
$g_{*}=3.38$.

We assume neutrinos are composed of two massless species and one massive
species with $m_{3}=0.06$ eV. We include the two massless neutrinos in
$\rho_{ur}$ and denote the energy density of the massive neutrino as
$\rho_{ncdm}$ (non-cold dark matter). Since the massive neutrino is
relativistic prior to recombination, the total radiation energy density is
$\rho_{r}=\rho_{\gamma}+\rho_{ur}+\rho_{ncdm}$ for the reheat temperatures
that we consider.

We parameterize how much energy density the $Y$ particle contributes to the
universe via the maximum of $\rho_{Y}/\rho_{r}$, denoted as
$max(\rho_{Y}/\rho_{r})$. Under the assumption of radiation domination, this
ratio can be expressed as

$\displaystyle
max\left(\frac{\rho_{Y}}{\rho_{r}}\right)=e^{\frac{1}{2}(\tilde{\Gamma}_{Y}-1)}\frac{\rho_{Y,i}}{\rho_{r,i}}\,\,\tilde{\Gamma}_{Y}^{-1/2}\,,$
(5)

where $\rho_{x,i}\equiv\rho_{x}(a_{i})$ for species $x$,
$\widetilde{\Gamma}_{Y}\equiv\Gamma_{Y}/H_{i}$, $H_{i}\equiv H(a_{i})$, and
$a_{i}$ is the initial scale factor of our numerical solution (see Appendix A
for a derivation).

We set $\rho_{\gamma,i}$ such that the energy density of photons after
$Y$-induced entropy injection corresponds to the fiducial $\Omega_{\gamma,0}$
determined by the measured CMB temperature of
$T_{0}=$2.7255\text{\,}\mathrm{K}$$. This initial condition for
$\rho_{\gamma}$ is set by the parameter $\Omega^{\prime}_{\gamma,0}$ such that
the initial photon energy density is
$\rho_{\gamma,i}=\Omega^{\prime}_{\gamma,0}\rho_{\text{crit},0}\,a_{i}^{-4}$.
The value of $\Omega^{\prime}_{\gamma,0}$ can be calculated directly from our
decay parameters $max(\rho_{Y}/\rho_{r})$ and $f_{\gamma}$ (see Appendix A).
We choose $a_{i}$ to occur after electron-positron annihilation so the initial
values for $\rho_{\gamma,i}$ and $\rho_{ur,i}$ are related by

$\rho_{ur,i}+\rho_{ncdm,i}=\frac{7}{8}(3.044)\left(\frac{4}{11}\right)^{4/3}\rho_{\gamma,i}\,,$
(6)

where we enforce the effective number of neutrinos to be 3.044 pre-decay. The
value of $\rho_{ncdm,i}$ is entirely determined by $\rho_{\gamma,i}$ (see
Appendix A), and the initial Hubble rate will be dominated by
$\rho_{r,i}=\rho_{\gamma,i}+\rho_{ur,i}+\rho_{ncdm,i}$. From this, a value for
$\rho_{Y,i}$ can be directly calculated from the parameters
$max(\rho_{Y}/\rho_{r})$ and $\Gamma_{Y}$ via equation (5). Therefore, we can
fully describe the decay of the $Y$ particle with the parameters
$max(\rho_{Y}/\rho_{r})$, $f_{\gamma}$, and $\Gamma_{Y}$. Figure 1 depicts a
solution to equations (1)-(3) for select values of $\widetilde{\Gamma}_{Y}$,
$f_{\gamma}$, and $max(\rho_{Y}/\rho_{r})$.

Figure 1: Numerical solution of equations (1)-(3) for
$\widetilde{\Gamma}_{Y}=10^{-6}$, $f_{\gamma}=0.6$, and
$max(\rho_{Y}/\rho_{r})=0.1$. The $Y$ particle, taken to be non-relativistic,
initially evolves as $\rho_{Y}\propto a^{-3}$ until the expansion rate equals
the decay rate near $a/a_{i}\approx 10^{3}$. Both $\rho_{\gamma}$ and
$\rho_{ur}$ scale as $\rho\propto a^{-4}$ up until there is significant
injection from the $Y$ decay at $a/a_{i}\approx 10^{3}$. Photons and other
ultra-relativistic species gain 60% and 40% of the decay products,
respectively. Then, $\rho_{\gamma}$ and $\rho_{ur}$ continue again as $a^{-4}$
after entropy injection. The comoving energy densities ($\rho a^{4}$) for
photons and other relativisitc particles are plotted in the bottom panel to
emphasize this evolution.

##  III Effects of decaying particle

The decaying $Y$ particle’s imprint on the CMB anisotropy spectrum is
primarily determined by two effects: the decay products change
$N_{\text{eff}}$ at recombination and they change the primordial helium
($Y_{\text{He}}$) and deuterium (D/H) abundances by altering the expansion
rate and baryon-to-photon ratio at BBN. In this section, we discuss these
effects as well as distortions in the CMB frequency spectrum induced by the
$Y$ decay. We compute the CMB temperature anisotropy power spectrum and
spectral distortions with the Cosmic Linear Anisotropy Solving System (CLASS)
[56]. For a discussion of the modifications made to CLASS during the
implementation of our model, see Appendix A.

Figure 2: Relationship between $max(\rho_{Y}/\rho_{r})$ defined by equation
(5) and the ratio of comoving radiation energy densities before and after $Y$
entropy injection (black line). This relationship is independent of the reheat
temperature of the $Y$ particle. We fit this trend with a fourth order
polynomial $g(x)=0.0274x^{4}-0.1525x^{3}+0.6458x^{2}+2.0727x+1$, where
$x=max(\rho_{Y}/\rho_{r})$ for shorthand (red dashed line).

###  III.1 Post-decay $N_{\text{eff}}$

Figure 3: Post-decay $N_{\text{eff}}$ contours determined by $f_{\gamma}$ and
$max(\rho_{Y}/\rho_{r})$ via equation (11) are shown by the dotted red lines.
The fiducial $N_{\text{eff}}$ value of $3.044$ is maintained as long as the
photon fraction is $f_{\gamma}=0.5913$ (solid red line). We include the 2018
Planck TT,TE,EE+lowE 95% confidence limits on $N_{\text{eff}}$ (orange shaded
region) to illustrate the $Y$ particle parameter space that is naively
consistent with CMB observations. Additionally, constraints on
$\eta_{\text{BBN}}$ derived from observations of the deuterium abundance [43]
are shown by the blue shaded region. The dashed purple lines show the values
of $\eta_{\text{BBN}}$ calculated via equation (12).

Since the decay of the $Y$ particle injects free-streaming radiation before
recombination, the effective number of relativistic species, $N_{\text{eff}}$,
will deviate from the fiducial value of 3.044. After electron-positron
annihilation, entropy conservation dictates that the temperature ratio between
neutrinos and photons is $T_{\nu}/T_{\gamma}=(4/11)^{1/3}$. Entropy injection
from the $Y$ decay will alter this temperature ratio by adding either photons
or other relativistic species. We can parameterize this change in
$T_{\nu}/T_{\gamma}$ as a change in $N_{\text{eff}}$ defined by

$\rho_{ur}+\rho_{ncdm}=\frac{7}{8}N_{\text{eff}}\left(\frac{4}{11}\right)^{4/3}\rho_{\gamma}\,.$
(7)

To relate the post-decay $N_{\text{eff}}$ to our model parameters, we first
define

$g\equiv\frac{\rho_{r,f}a_{f}^{4}}{\rho_{r,i}a_{i}^{4}}\,,$ (8)

where $\rho_{r,i}a_{i}^{4}$ and $\rho_{r,f}a_{f}^{4}$ are the comoving
radiation energy densities before and after the $Y$-induced entropy injection,
respectively. One would expect the change in comoving radiation energy density
as a result of the $Y$ decay to be directly dependent on $\rho_{Y}$ when
$H\sim\Gamma_{Y}$. Indeed, $g$ is a simple function of
$max(\rho_{Y}/\rho_{r})$. This can be seen in Figure 2 where we numerically
calculate $g$ for various values of $max(\rho_{Y}/\rho_{r})$ defined by
equation (5) and fit the function with a fourth order polynomial.222If all of
the $Y$ energy density was instantaneously transferred to $\rho_{r}$, then $g$
would equal $1+max(\rho_{Y}/\rho_{r})$. However, as can be seen from Figure 1,
most of the change in $\rho_{r}a^{4}$ occurs after $\rho_{Y}/\rho_{r}$ is
maximized. As a result, $g-1>max(\rho_{Y}/\rho_{r})$.

We can employ $g$ to determine the ratio of comoving energy densities before
and after the $Y$ decay for the individual decay species:

$\displaystyle g_{\gamma}$
$\displaystyle\equiv\frac{\rho_{\gamma,f}a_{f}^{4}}{\rho_{\gamma,i}a_{i}^{4}},$
(9)
$\displaystyle=1+f_{\gamma}\left(g-1\right)\left[1+\frac{7}{8}(3.044)\left(\frac{4}{11}\right)^{4/3}\right];$
$\displaystyle g_{ur}$
$\displaystyle\equiv\frac{\rho_{ur,f}a_{f}^{4}}{\rho_{ur,i}a_{i}^{4}},$ (10)
$\displaystyle=1+\frac{(1-f_{\gamma})\left(g-1\right)}{N_{ur}}\left[3.044+\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\right]\,\,.$

Here, $N_{ur}=2.0308$ is the effective number of ultra-relativistic species
before the $Y$ decay, excluding the single massive neutrino (see Appendix A).

The single massive neutrino species is still relativistic at recombination and
contributes to $N_{\text{eff}}$. As shown in Appendix B, it follows that the
post-decay $N_{\text{eff}}$ is

$N_{\text{eff}}=\frac{1}{g_{\gamma}}\left[3.044+N_{ur}(g_{ur}-1)\right]\,.$
(11)

This relationship is shown in Figure 3 where we plot contours of the post-
decay $N_{\text{eff}}$ for a range of decay scenarios described by
$f_{\gamma}$ and $max(\rho_{Y}/\rho_{r})$ (red dotted lines). Most notably,
there is a “sweet spot” fraction of $f_{\gamma}\simeq 0.59$ that maintains
$N_{\text{eff}}=3.044$ after the decay for all values of
$max(\rho_{Y}/\rho_{r})$ (shown by the solid red line). Figure 3 shows that
$f_{\gamma}<0.59$ corresponds to an increase in $N_{\text{eff}}$ since the
decay of the $Y$ particle increases the relative energy density of ultra-
relativistic species. Conversely, $f_{\gamma}>0.59$ decreases
$\rho_{ur}/\rho_{\gamma}$ and thereby results in a decrease in
$N_{\text{eff}}$.

Some combinations of $f_{\gamma}$ and $max(\rho_{Y}/\rho_{r})$ will not be
favored by Planck data. This is demonstrated by the orange contour in Figure
3, showing the 2018 Planck TT,TE,EE+lowE reported bounds on $N_{\text{eff}}$
[42]. For scenarios in which $f_{\gamma}$ deviates from the “sweet spot”
value, Planck $N_{\text{eff}}$ bounds require smaller values of
$max(\rho_{Y}/\rho_{r})$. However, it is important to note that naively
applying these Planck $N_{\text{eff}}$ bounds to our model parameters would be
incorrect, as doing so would ignore the fact that there is a degeneracy
between $N_{\text{eff}}$ and the primordial helium abundance, $Y_{\text{He}}$,
which is also affected by the $Y$ particle and its decay products.

To understand this degeneracy with $Y_{\text{He}}$, we must first understand
how changes in $N_{\text{eff}}$ affect the CMB anisotropy spectrum. Increasing
$N_{\text{eff}}$ increases the pre-recombination expansion rate and therefore
decreases the sound horizon $r_{s}\propto 1/H$. In order to maintain the
precisely measured angular size of the sound horizon, $\theta_{s}$, one must
increase $H_{0}$ in order to decrease the angular diameter distance to the CMB
such that $\theta_{s}$ is fixed. Meanwhile, the photon diffusion length scales
as $r_{D}\propto\sqrt{1/H}$ and so increasing $N_{\text{eff}}$ while
simultaneously keeping $\theta_{s}$ fixed leads to an increase in the angular
size of the diffusion length, $\theta_{D}$ [47, 48]. In summary, an increase
in $N_{\text{eff}}$ results in more Silk damping on small angular scales (and
visa versa).

Altering the helium abundance can mitigate the effects of changing
$N_{\text{eff}}$ on the dampening tail. Decreasing $Y_{\text{He}}$ results in
more free electrons at recombination, which decreases the Compton mean free
path, resulting in less photon diffusion and thus less damping of small-scale
anisotropies. This behavior is depicted in Figure 4. The dotted orange line
shows an increase in $N_{\text{eff}}$ resulting in more damping on small
scales compared to $\Lambda$CDM (for fixed $\omega_{b}$, $z_{eq}$, and
$\theta_{s}$, where $\omega_{b}\equiv\Omega_{b,0}h^{2}$ is the present-day
baryon energy density and $z_{eq}$ is the redshift of matter-radiation
equality). Simultaneously decreasing $Y_{\text{He}}$ as $N_{\text{eff}}$
increases, shown by the blue dashed line, can mitigate the damping. The phase
shift observed in Figure 4 is a lingering effect of changing $N_{\text{eff}}$;
free-streaming relativistic particles generate a unique phase shift in the
acoustic peaks [49, 50, 51].

Figure 4: Effects on the CMB temperature power spectrum from varying
$N_{\text{eff}}$ and $Y_{\text{He}}$ at fixed $\omega_{b}$, $z_{eq}$, and
$\theta_{s}$. Increasing $N_{\text{eff}}$ compared to $\Lambda$CDM results in
more damping on small scales (orange dotted line). This effect can be
mitigated by decreasing $Y_{\text{He}}$ (blue dashed line). The beat frequency
seen in these residuals is a direct consequence of the free-streaming nature
of neutrinos producing a unique phase shift in the CMB acoustic peaks.

Therefore, even though our model parameters $f_{\gamma}$ and
$max(\rho_{Y}/\rho_{r})$ conspire to produce a non-zero $\Delta
N_{\text{eff}}$, altering $Y_{\text{He}}$ in a inverse manner (i.e. raising
$Y_{\text{He}}$ as $N_{\text{eff}}$ decreases, or visa versa) can mitigate the
effects of Silk damping that a changing $N_{\text{eff}}$ causes.
Interestingly, our model already necessitates an increase in $Y_{\text{He}}$
by altering the expansion rate and baryon-to-photon ratio at BBN.

###  III.2 Change in primordial abundances

Since the injection of photons from the $Y$ decay necessitates a decrease in
$\rho_{\gamma}$ pre-decay, any $f_{\gamma}\neq 0$ decay alters the baryon-to-
photon ratio at the time of BBN, $\eta_{\text{BBN}}$. If $\eta_{0}$ is the
baryon-to-photon ratio at recombination inferred by $\eta_{0}=6.0\times
10^{-10}\left(\omega_{b}/0.022\right)$, then the increased ratio prior to the
decay of the $Y$ particle is

$\displaystyle\eta_{\text{BBN}}=g_{\gamma}^{3/4}\eta_{0}\,.$ (12)

Given that $g_{\gamma}$ is completely determined by $f_{\gamma}$ and
$max(\rho_{Y}/\rho_{r})$ via equation (9), we are able to place constraints on
our decay parameters via constraints on $\eta_{\text{BBN}}$. This relationship
is shown in Figure 3 where we calculate $\eta_{\text{BBN}}$ as a function of
$f_{\gamma}$ and $max(\rho_{Y}/\rho_{r})$ for a fixed value of
$\omega_{b}=0.02238$ (dashed purple lines). The blue contour shows the
parameter space that is allowed by the bounds $\eta_{\text{BBN}}=(6.119\pm
0.100)\times 10^{-10}$ derived from observations of the primordial deuterium
abundance [43]. Since large $f_{\gamma}$ values correspond to significant
photon injection, any increase in $max(\rho_{Y}/\rho_{r})$ at a large
$f_{\gamma}$ leads to a significant change in $\eta_{\mathrm{BBN}}$. On the
other hand, we can see that large values of $max(\rho_{Y}/\rho_{r})$ are
allowed by bounds on $\eta_{\mathrm{BBN}}$ as $f_{\gamma}\rightarrow 0$ since
fewer photons are being injected by the decay and so $g_{\gamma}$ approaches
unity (i.e. $\eta_{\text{BBN}}\rightarrow\eta_{0}$).

Much like the bounds on $N_{\text{eff}}$ discussed in Section III.1, these
bounds on $\eta_{\mathrm{BBN}}$ seen in Figure 3 are only meant to depict an
approximation for the allowed parameter space; the true bounds will be subject
to uncertainties in $\omega_{b}$, which are taken into account later.
Furthermore, the inference of $\eta_{\mathrm{BBN}}$ from deuterium assumes a
standard expansion history during BBN. Decay scenarios with large
$max(\rho_{Y}/\rho_{r})$ and short lifetimes that end soon after the end of
BBN ($\tau_{Y}\approx 10^{4}$ sec) will contribute a non-negligible energy
density during BBN.

We therefore altered the BBN code PArthENoPe v3.0 [57] to include the
contribution of $\rho_{Y}$ to the expansion rate during BBN. We modified
PArthENoPe to accept a new input parameter, $\rho_{Y}^{\mathrm{BBN}}$, which
is the value of $\rho_{Y}$ at a reference temperature of $T=0.1$ MeV, and we
assume the $Y$ particle is non-relativistic throughout BBN so that
$\rho_{Y}\propto a^{-3}$. This $\rho_{Y}^{\mathrm{BBN}}$ parameter is entirely
determined by our decay parameters $max(\rho_{Y}/\rho_{r})$ and $\Gamma_{Y}$.
Using our modified version of PArthENoPe, we created lookup tables for CLASS
that read in $\eta_{\mathrm{BBN}}$ and $\rho_{Y}^{\mathrm{BBN}}$ and produce
values for the helium and deuterium abundances. These tables were created with
$\Delta N_{\text{eff}}=0$ (our initial conditions ensure that
$N_{\text{eff}}=3.044$ at BBN) and a neutron lifetime of
$\tau=$879.4\text{\,}\mathrm{s}$$.

Since the $Y$ particle is non-relativistic during BBN and evolves as
$\rho_{Y}\propto a^{-3}$, its contribution to $H$ cannot be modeled with a
constant $\Delta N_{\text{eff}}$. Instead, one can think of the contribution
of $\rho_{Y}$ to the expansion rate during BBN as an evolving $\Delta
N_{\text{eff}}$; at any given temperature, there is a nonzero $\Delta
N_{\text{eff}}$ that matches the non-standard Hubble rate due to the inclusion
of $\rho_{Y}$ [58]. This behavior is demonstrated in Figure 5. The black solid
lines track the abundances of deuterium (D/H) and helium ($Y_{\text{He}}$) as
a function of $\eta_{\mathrm{BBN}}$ for a specific decay scenario in which the
$Y$ particle decays away right after the end of BBN ($T_{\mathrm{RH}}=0.01$
MeV). If we then match the Hubble rate at a temperature of $T=1$ MeV using a
$\Delta N_{\text{eff}}=0.0865$ (blue dashed line), we mostly recover the
correct $Y_{\text{He}}$ since the abundance of helium is primarily set by
neutron-proton freeze-out near $T=1$ MeV. However, the abundance of deuterium
is set when helium production freezes out around $T=0.07$ MeV. So, matching at
$T=1$ MeV correctly approximates $Y_{\text{He}}$, but fails to yield the
correct D/H. In a similar manner, matching at $T=0.07$ MeV with $\Delta
N_{\text{eff}}=1.6944$ (green dotted line) results in a good approximation for
the deuterium abundance, but grossly overestimates the correct abundance of
helium. Therefore, the contribution of $\rho_{Y}$ to the Hubble rate during
BBN cannot be captured by a simple constant change to $N_{\text{eff}}$.

Figure 5: Dependence of primordial abundances on $\eta_{\mathrm{BBN}}$ for a
decay scenario that decays right after the end of BBN (black solid line).
Matching the non-standard Hubble rate during BBN as a result of $\rho_{Y}$
with a constant nonzero $\Delta N_{\text{eff}}$ does not reproduce the same
abundances. Matching Hubble at $T=1$ MeV when neutron-proton freeze-out occurs
(blue dashed line) is a good approximation for $Y_{\text{He}}$, but
underestimates D/H. Matching Hubble at $T=0.07$ MeV when helium production
freezes out (green dotted line) is a good approximation for D/H, but
overestimates $Y_{\text{He}}$.

The abundances of deuterium and helium are both sensitive to
$\eta_{\mathrm{BBN}}$ and the expansion rate. Increasing the expansion rate
during BBN leads to an increase in both $Y_{\text{He}}$ and D/H, while
increasing $\eta_{\mathrm{BBN}}$ results in larger $Y_{\text{He}}$ and smaller
D/H. Any $Y$ decay scenario in which the decay products consist of photons
sees an increase in both $\eta_{\text{BBN}}$ and $H$, and thus an increase in
$Y_{\text{He}}$. On the other hand, the abundance of deuterium is
predominately sensitive to changes in $\eta_{\mathrm{BBN}}$. So, even though
$\rho_{Y}$ increases the expansion rate during BBN, we expect an overall
decrease in D/H for any $f_{\gamma}\neq 0$ decay scenario.

Finally, we note that all $f_{\gamma}\neq 0$ scenarios result in a larger
$Y_{\text{He}}$ compared to $\Lambda$CDM while changes in $N_{\text{eff}}$ are
symmetric about $f_{\gamma}\approx 0.59$ (see Figure 3). Therefore, we expect
an asymmetry in the small-scale Silk damping resulting from decays of
different $f_{\gamma}$. Decays with $f_{\gamma}>0.59$ will have a lack of
small scale damping from a decreased $N_{\text{eff}}$, but an increased
$Y_{\text{He}}$ compensates by enhancing Silk damping. On the other hand,
scenarios with $f_{\gamma}<0.59$ will have combined damping effects of both an
increased $Y_{\text{He}}$ and $N_{\text{eff}}$. This asymmetry suggests a
preference for decay scenarios with $f_{\gamma}>0.59$.

###  III.3 Spectral Distortions

Deviations from the blackbody spectrum of the CMB, known as spectral
distortions (SDs), are the result of energy injection into the photon bath of
the CMB at a late enough time such that the photons do not thermalize before
recombination [59, 60, 61]. The form of these SDs depends on the time of
energy injection; the era of $\mu$ distortions occurs roughly in the redshift
range of $5\times 10^{4}<z<2\times 10^{6}$, while $y$-type distortions result
from injection at $z<5\times 10^{4}$ [62, 53]. For $z>2\times 10^{6}$,
thermalization processes are very efficient and so any energy injection
results in a temperature shift to the observed blackbody. Since we will be
considering reheat temperatures that correspond to a range of about
$10^{5}<z<4.5\times 10^{7}$, $\mu$ distortions will have the most constraining
power.

We calculate the distortions in CLASS using a new SD module based on the work
of Lucca  _et al._ [53]. This module calculates $\mu$ and $y$ distortions by

$\displaystyle\tilde{d}=\frac{\Delta\rho_{\gamma}}{\rho_{\gamma}}\Bigg{|}_{d}$
$\displaystyle\equiv\int\frac{dQ/dz}{\rho_{\gamma}}\cdot\mathcal{J}_{d}(z)dz\,,$
(13)

where $d$ indexes the distortion type,333Lucca  _et al._ [53] uses “$a$” to
index distortion type, but we choose to use “$d$” in order to avoid confusion
with scale factor. $dQ/dz$ is the energy injection rate, and
$\mathcal{J}_{d}(z)$ is a branching ratio that dictates the contribution of an
energy injection to a specific distortion type. The energy injection rate can
be recast as a rate with respect to proper time, $t$, via

$\displaystyle\frac{dQ/dz}{\rho_{\gamma}}$
$\displaystyle=-\frac{\dot{Q}}{(1+z)H\rho_{\gamma}}\,,$

and we define the energy injection rate of our decaying model as

$\displaystyle\dot{Q}(t)=\rho_{Y}(t)f_{\gamma}\Gamma_{Y}e^{-\Gamma_{Y}t}\,.$

The distortion parameters are then determined by equation (13) and setting
$\mu=1.4\tilde{\mu}$ and $y=\tilde{y}/4$.

The COBE/FIRAS satellite measured the blackbody spectrum of the CMB and
determined upper bounds on the distortions to be $|\mu|<9\times 10^{-5}$ and
$|y|<1.5\times 10^{-5}$ (95% CL) [44, 45]. In Figure 6, we calculate $|\mu|$
for a range of $f_{\gamma}$ and $T_{\mathrm{RH}}$ values using CLASS and
depict the parameter space allowed by the COBE/FIRAS bound by the green shaded
region. Here we see that the COBE/FIRAS bound on $\mu$ allows photon injection
for reheat temperatures hotter than about $9\times 10^{-4}$ MeV, but requires
that $f_{\gamma}$ quickly approaches zero for $T_{\mathrm{RH}}\lesssim 9\times
10^{-4}$ MeV
($\Gamma_{Y}\lesssim$3.08\text{\times}{10}^{-7}\text{\,}{\mathrm{s}}^{-1}$$ ).

Since the CMB angular scales that are accessible to Planck ($\ell\lesssim
2500$) are dominated by modes that enter the horizon at a temperatures below
$2\times 10^{-5}$ MeV, we can neglect perturbations related to the $Y$
particle for cases in which $f_{\gamma}\gtrsim 0.01$ because those values will
only be allowed by SDs when $T_{\mathrm{RH}}\gtrsim 9\times 10^{-4}$ MeV.
However, scenarios with $f_{\gamma}=0$ are not restricted in reheat
temperature by SDs, and therefore require consideration of the $Y$
perturbation equations as well as corrections to the ultra-relativistic
species perturbations. We reserve the inclusion of these perturbations for a
subsequent analysis and, in this work, focus solely on $f_{\gamma}\geq 0.01$.
Nevertheless, our results can be extended to smaller values of $f_{\gamma}$,
provided that $T_{\mathrm{RH}}>9.5\times 10^{-4}$ MeV
($\Gamma_{Y}>$3.43\text{\times}{10}^{-7}\text{\,}{\mathrm{s}}^{-1}$$ ). For
these decay rates, the posteriors for all parameters are the same when
comparing $f_{\gamma}$ fixed at either $f_{\gamma}=0.01$ or $f_{\gamma}=0$
(see Figure 11 in Appendix C).

Figure 6: Example of expected $\mu$ distortion constraints for a value of
$max(\rho_{Y}/\rho_{r})=0.01$. The green shaded region shows the parameter
space allowed by the COBE/FIRAS constraint $|\mu|<9\times 10^{-5}$. Reheat
temperatures cooler than $T_{\mathrm{RH}}\approx 9\times 10^{-4}$ MeV are only
allowed by SDs if $f_{\gamma}\rightarrow 0$.

##  IV Analysis Method

We derive constraints for $\Gamma_{Y}$, $max(\rho_{Y}/\rho_{r})$, and
$f_{\gamma}$ by conducting a Markov-Chain Monte Carlo (MCMC) analysis using
MontePython-v3 [63, 64] with our modified version of CLASS. This analysis is
implemented with a Metropolis-Hastings algorithm and flat priors on the base
six cosmological parameters {$\omega_{b}$, $\omega_{cdm}$, $\theta_{s}$,
$A_{s}$, $n_{s}$, $\tau_{reio}$}. We assume neutrinos are composed of two
massless and one massive species with $m_{3}=0.06$ eV. In order to account for
the fact that the $Y$ particle decaying into photons changes
$T_{\nu}/T_{\gamma}$, we provide CLASS with a corrected temperature of the
massive neutrino (see Appendix A). We employ Planck high-$\ell$ TTTEEE,
low-$\ell$ TT, and low-$\ell$ EE likelihood functions [65] and refer to the
combination of these data as Planck.

As discussed in Section III.2, the $Y$ decay affects the abundances of
$Y_{\text{He}}$ and D/H through its influence on $\eta_{\mathrm{BBN}}$ and the
expansion rate during BBN. Current measurements of the deuterium abundance
have reached a $\sim 1\%$ precision level with $(\mathrm{D/H})=(2.527\pm
0.030)\times 10^{-5}$ [43]. However, we note that this bound only includes
measurement uncertainty (i.e. $\sigma_{\mathrm{MEAS}}=3.0\times 10^{-7}$).
Cooke  _et al._ [43] cite uncertainties associated with BBN calculations for
two different values of the $d(p,\gamma)^{3}\text{He}$ cross-section; using
the computationally inferred cross-section leads to a $2\sigma$ discrepancy
with CMB measurements of $\eta_{\mathrm{BBN}}$, while utilizing the measured
cross-section suggested by Adelberger  _et al._ [66] yields
$\eta_{\mathrm{BBN}}=(6.119\pm 0.100)\times 10^{-10}$. Millea  _et al._ [25]
fold in uncertainty associated with nuclear reaction rates to D/H by assuming
$\sigma_{\mathrm{NUCL}}=4.5\times 10^{-7}$ and adding uncertainties in
quadrature such that $\sigma_{\mathrm{DH}}=5.4\times 10^{-7}$. It is unclear
which $d(p,\gamma)^{3}\text{He}$ cross-section is used for this uncertainty,
and Millea  _et al._ [25] assume a linear relationship between D/H and $\eta$.
Therefore, we translate the $\eta_{\mathrm{BBN}}=(6.119\pm 0.100)\times
10^{-10}$ bounds reported by Cooke  _et al._ [43] into an uncertainty in D/H.
We calculate D/H for a range of $\eta_{\mathrm{BBN}}$ with our modified
version of PArthENoPe for $\Delta N_{\text{eff}}^{\mathrm{BBN}}=0$ and a
neutron lifetime of $879.4\text{\,}\sec$. We find that
$\text{D/H}\propto\eta_{\mathrm{BBN}}^{-1.65}$ for the range of
$\eta_{\mathrm{BBN}}$ relevant to this work and that this fit is relatively
insensitive to small changes in $N_{\text{eff}}^{\mathrm{BBN}}$. We use this
power law to infer a new fractional uncertainty on D/H
($\sigma_{\text{DH}}/{\text{DH}}=1.65\times\sigma_{\eta_{\mathrm{BBN}}}/{\eta_{\mathrm{BBN}}}$).
Doing so results in $(\text{D/H})=(2.527\pm 0.068)\times 10^{-5}$. We create a
Gaussian likelihood function in MontePython for D/H for which the mean and
standard deviation are set to be $\mu_{\mathrm{DH}}=2.527\times 10^{-5}$ and
$\sigma_{\mathrm{DH}}=6.83\times 10^{-7}$, respectively. We denote this
Gaussian likelihood as D/H.

Finally, we calculate the spectral distortions $\mu$ and $y$ in CLASS
according to Section III.3. Similar to our constraint on D/H, we add a
Gaussian likelihood for both $\mu$ and $y$. These Gaussian likelihood
functions were taken to have a mean of zero and a 2$\sigma$ deviation equal to
the respective upper bounds $|\mu|<9\times 10^{-5}$ and $|y|<1.5\times
10^{-5}$ determined by COBE/FIRAS. We refer to the combination of these $\mu$
and $y$ likelihoods as SD.

We used the Gelman-Rubin [67] criterion $|R-1|<0.04$ to asses convergence of
our MCMC chains.444All $\Lambda$CDM runs have $|R-1|<0.01$. This bound of
$0.04$ was selected due to the non-Gaussian nature of the posteriors of our
decay parameters. Post-processing of chains was done using GetDist [68] and
removing the first 30% of points as burn-in.

###  IV.1 Initial priors

Figure 3 demonstrates that naive bounds from $N_{\text{eff}}$ would leave
$max(\rho_{Y}/\rho_{r})$ unconstrained if $f_{\gamma}$ is such that
$N_{\text{eff}}=3.044$ after the $Y$ decay. The addition of
$\eta_{\mathrm{BBN}}$ bounds, however, suggests that $max(\rho_{Y}/\rho_{r})$
would be limited to less than 0.05. Therefore, we choose to sample
$max(\rho_{Y}/\rho_{r})=[0,0.07]$ when including the D/H likelihood.
Otherwise, when only considering Planck or Planck+SD, we sample
$max(\rho_{Y}/\rho_{r})=[0,0.8]$ in order to sufficiently explore the
asymptotic behavior around the “sweet spot” $f_{\gamma}$ shown in Figure 3. A
value of $max(\rho_{Y}/\rho_{r})=0$ corresponds to standard $\Lambda$CDM.

As discussed in Section III.3, the smallest scale accessible to Planck enters
the horizon at a temperature of about $2\times 10^{-5}$ MeV. Therefore, a full
perturbation analysis of the $Y$ particle and its decay products is required
for reheat temperatures $T_{\mathrm{RH}}\lesssim 2\times 10^{-5}$ MeV. We
retain this for a subsequent analysis and instead restrict our current work to
$T_{\mathrm{RH}}>2\times 10^{-5}$ MeV. Furthermore, CLASS begins evolving a
given perturbation mode outside the cosmological horizon with adiabatic
initial conditions. In order to avoid the $Y$ decay influencing these initial
conditions, we choose to only consider $T_{\mathrm{RH}}\geq 9.5\times 10^{-4}$
MeV so that $Y$-induced entropy injection is complete by the time CLASS begins
evolving perturbation modes relevant to Planck
($k\lesssim$0.18\text{\,}{\mathrm{M\text{pc}}}^{-1}$$ or $\ell\lesssim 2500$).
Assuming BBN has completed by $T\approx 0.01$ MeV, this restriction on
$T_{\mathrm{RH}}$ leads to a prior of $T_{\mathrm{RH}}=[9.5\times
10^{-4},10^{-2}]$ MeV ( $\Gamma_{Y}=[3.43\times 10^{-7},3.80\times
10^{-5}]\,$\text{\,}{\mathrm{s}}^{-1}$$).

Spectral distortions will confine $f_{\gamma}$ to be very small as
$T_{\mathrm{RH}}$ approaches $9.5\times 10^{-4}$ MeV. Therefore we would like
to sample small values of $f_{\gamma}$ for colder reheat temperatures, but
explore $f_{\gamma}$ of order unity when $T_{\mathrm{RH}}\gg 9.5\times
10^{-4}$ MeV. A flat prior of $f_{\gamma}=[0,1]$ fails to sufficiently sample
small values of $f_{\gamma}$ as $T_{\mathrm{RH}}\rightarrow 9.5\times 10^{-4}$
MeV, but using a flat prior on $\log_{10}f_{\gamma}$ does not properly depict
the parameter space of large $f_{\gamma}$. Furthermore, small values of
$f_{\gamma}$, such as $f_{\gamma}=10^{-5}$, would be indicative of a finely
tuned decay scenario. To avoid these scenarios and achieve sufficient
sampling, we choose to sample $f_{\gamma}=[0.01,1]$. The initial priors used
in our analysis are summarized in Table 1.

|  Planck, Planck+SD |  Planck+SD+D/H
---|---|---
$f_{\gamma}$ | [0.01, 1]
$T_{\mathrm{RH}}$ [MeV] | $[9.5\times 10^{-4},10^{-2}]$
$max(\rho_{Y}/\rho_{r})$ | [0, 0.8] | [0, 0.07]
Table 1: Summary of priors used for decay parameters with each likelihood
combination.

##  V Results

Figure 7 shows the resulting posterior distributions of the decay parameters
as well as $\omega_{b}$, $\omega_{cdm}$, $A_{s}$, $n_{s}$, $H_{0}$,
$Y_{\text{He}}$, and D/H, with 68% and 95% CL contours for Planck, Planck+SD,
and Planck+SD+D/H data sets. Table 2 summarizes the mean and $2\sigma$ errors
for each parameter. For the posteriors of all base six cosmological
parameters, see Figure 12 in Appendix C.

In Figure 7, it can be seen that the Planck data favors changes in
$\omega_{b}$, $\omega_{cdm}$, $A_{s}$, and $H_{0}$ in order to accommodate a
$Y$ decay scenario. For cases in which $f_{\gamma}$ is greater than the “sweet
spot” value of $0.5913$, $N_{\text{eff}}$ is less than $3.044$, which results
in a larger $z_{eq}$ compared to $\Lambda$CDM. The Planck data, being very
sensitive to $z_{eq}$ via the early ISW effect, favors a decrease in
$\omega_{cdm}$ in order to restore $z_{eq}$ to its fiducial value.
Furthermore, the $Y$ induced entropy injection directly alters the pre-
recombination expansion history and therefore the size of the sound horizon,
$r_{s}$. In the case of $f_{\gamma}>0.5913$, $r_{s}$ is increased compared to
$\Lambda$CDM. The decrease in $\omega_{cdm}$ that is required to maintain a
fixed $z_{eq}$ results in a larger angular diameter distance, $d_{A}$.
However, the rate of change for $r_{s}$ is greater than that of $d_{A}$ and
thus leads to an increase in the angular size of the sound horizon
$\theta_{s}$. Since $\theta_{s}$ sets the anisotropy peak locations, the
Planck data tends to reduce this increase in $\theta_{s}$ by decreasing
$H_{0}$.

Figure 7: 1D and 2D posterior distributions of decay and cosmological
parameters for different combinations of Planck high-$\ell$ TT,TE,EE,
low-$\ell$ TT, and low-$\ell$ EE (Planck) data, CMB spectral distortions (SD)
bounds, and bounds on the observed deuterium abundance (D/H). We include the
1D posteriors for $\Lambda$CDM constrained by Planck (dashed black line). The
dotted gray line traces $f_{\gamma}=0.5913$, which maintains
$N_{\text{eff}}=3.044$ at recombination. Figure 8: Power spectra for decays
with $max(\rho_{Y}/\rho_{r})=0.2$,
$T_{\mathrm{RH}}=$3.16\text{\times}{10}^{-3}\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$,
and either $f_{\gamma}=0.4$ (dotted red line) or $f_{\gamma}=0.7826$ (dashed
blue line). Top left: power spectra of two decay scenarios compared to
spectrum of best-fit $\Lambda$CDM parameters (solid black line) with no other
parameters being altered. Top right: $\omega_{cdm}$ and $h$ are simultaneously
altered in order to fix $z_{eq}$ and $\theta_{s}$, respectively. Bottom: A
change in $A_{s}$, in addition to $\omega_{cdm}$ and $h$, brings the spectra
of these decays scenarios close to the best-fit $\Lambda$CDM spectrum. The
remaining small-scale residuals are asymmetric; the $f_{\gamma}=0.7826$ case
is aided by the increase in $Y_{\text{He}}$ and therefore has a smaller
residual than the $f_{\gamma}=0.4$ case.

This process of fixing $z_{eq}$ and $\theta_{s}$ is illustrated in Figure 8.
The top left panel of Figure 8 depicts the power spectra for two decay
scenarios described by $max(\rho_{Y}/\rho_{r})=0.2$,
$T_{\mathrm{RH}}=$3.16\text{\times}{10}^{-3}\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$,
and either $f_{\gamma}=0.4$ (dotted red line) or $f_{\gamma}=0.7826$ (dashed
blue line). These values of $f_{\gamma}$ correspond to equivalent deviations
in either direction about the sweet spot $f_{\gamma}=0.5913$. Then, in the top
right panel of Figure 8, we alter $\omega_{cdm}$ in order to fix $z_{eq}$ and
change $h\equiv
H_{0}/($100\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$)$
to compensate for the change in $\theta_{s}$. Even though $\theta_{s}$ is
fixed, there is a lingering offset in peak locations. As mentioned in Section
III.1, this phase shift is caused by the addition of extra free-streaming
radiation [49, 50, 51]. Furthermore, the alteration of $\omega_{cdm}$ results
in an amplification or suppression of all modes. Therefore, Planck data
adjusts the amplitude $A_{s}$. The bottom panel of Figure 8 shows a
continuation of the top right panel, with corrections to $\omega_{cdm}$ and
$h$, but now with an additional change in $A_{s}$ such that the first peak
height is the same as that of the best-fit $\Lambda$CDM spectrum.

The lingering small-scale damping effects of $N_{\text{eff}}$ and
$Y_{\text{He}}$ discussed in Section III.1 are now seen in the bottom panel of
Figure 8. The case of $f_{\gamma}=0.7826$ ($N_{\text{eff}}<3.044$) has a lack
of small scale damping compared to $\Lambda$CDM whereas the $f_{\gamma}=0.4$
($N_{\text{eff}}>3.044$) case has too much damping on small scales. Changing
the spectral index, $n_{s}$, can compensate for this non-standard Silk damping
[48]. Non-standard damping can also be compensated by changes in $\omega_{b}$
such that the number of free electrons increases or decreases and alters the
mean free path of photons. However, changing $\omega_{b}$ has the unintended
consequence of affecting the height ratio of odd and even peaks and is
therefore a costly approach to accommodating a $Y$ decay scenario.

The corrections to $\omega_{cdm}$, $h$, $A_{s}$, $n_{s}$, and $\omega_{b}$,
described above are apparent in the Planck 1D posteriors of Figure 7. However,
note that there is an asymmetry in these corrections. The 1D posterior of
$\omega_{cdm}$, for example, is not symmetric about the $\Lambda$CDM
distribution. This asymmetry is due to the increase in $Y_{\text{He}}$
discussed in Section III.2. For $f_{\gamma}>0.5913$, the decrease in
$N_{\text{eff}}$ is partially compensated by an increase in $Y_{\text{He}}$.
Therefore, decreasing $n_{s}$ or $\omega_{b}$ can compensate for values of
$f_{\gamma}>0.59$, making these scenarios more viable than the naive Planck
bounds on $N_{\text{eff}}$ would imply.

On the other hand, if $f_{\gamma}<0.5913$, then $N_{\text{eff}}>3.044$,
leading to too much damping on small scales. Planck will try to compensate by
increasing $n_{s}$ or $\omega_{b}$. However, $Y_{\text{He}}$ is still larger
compared to $\Lambda$CDM in this $f_{\gamma}<0.5913$ regime. So, any increase
in $n_{s}$ or $\omega_{b}$ is fighting against the combined damping effects of
increasing both $N_{\text{eff}}$ and $Y_{\text{He}}$. Therefore, Planck tends
to favor $f_{\gamma}>0.59$ rather than $f_{\gamma}<0.59$.

This asymmetric effect can also be seen in the residuals of the bottom panel
of Figure 8. The case of $f_{\gamma}=0.7826$ has smaller residuals on small
scales than that of the $f_{\gamma}=0.4$ scenario because the increased
$Y_{\text{He}}$ is helping to compensate for the lack of damping caused by a
decreased $N_{\text{eff}}$. The $f_{\gamma}-max(\rho_{Y}/\rho_{r})$ plane in
Figure 7 further demonstrates the Planck preference for $f_{\gamma}>0.5913$.
When only Planck data is considered, $max(\rho_{Y}/\rho_{r})$ is unconstrained
as long as $f_{\gamma}$ is near the “sweet spot” value that maintains
$N_{\text{eff}}\simeq 3.0$ after the $Y$-induced entropy injection. However,
the Planck 1D posterior peaks at $f_{\gamma}=0.632$ rather than
$f_{\gamma}=0.5913$ (gray dashed line), which would give
$N_{\text{eff}}=3.044$.

In theory, Planck would accept even larger $max(\rho_{Y}/\rho_{r})$ at the
sweet spot $f_{\gamma}$ until reaching a limit that corresponds to the Planck
bounds on $Y_{\text{He}}$; Planck 2018 high-$\ell$ TT,TE,EE+low-$\ell$ EE
places a 2$\sigma$ upper bound of $Y_{\text{He}}=0.283$ [42]. Indeed, we
obtain similar results when using a prior of $max(\rho_{Y}/\rho_{r})=[0,1.2]$.
However, exploring the $max(\rho_{Y}/\rho_{r})$ corresponding to this
$Y_{\text{He}}$ Planck upper bound is unnecessary since the addition of our
D/H bounds greatly constrains $Y_{\text{He}}$.555 Additionally, the accuracy
of the approximation in equation (5) begins to diminish as
$max(\rho_{Y}/\rho_{r})\rightarrow 1$. See Appendix A. Since marginalizing
over $f_{\gamma}$ leaves $max(\rho_{Y}/\rho_{r})$ unconstrained for Planck and
Planck+SD, we do not report an upper bound on $max(\rho_{Y}/\rho_{r})$ for
these likelihoods in Table 2. We also note that the steep drop-off seen at
large $Y_{\text{He}}$ in the posteriors shown in Figure 12 results from the
MCMC sampling reaching the $max(\rho_{Y}/\rho_{r})=0.8$ upper bound of our
selected prior.

The inclusion of SD constraints with Planck anisotropy data (Planck+SD) has
the added power of restricting $T_{\mathrm{RH}}$ to hotter reheat
temperatures, as expected. SDs constrain the reheat temperature to be
$T_{\mathrm{RH}}>$2.24\text{\times}{10}^{-3}\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$
($\Gamma_{Y}>$1.91\text{\times}{10}^{-6}\text{\,}{\mathrm{s}}^{-1}$$) at
$95\%$ CL. Otherwise, SDs provide no additional constraining power on
cosmological parameters.

It is clear from Figure 7 that the addition of a bound on D/H greatly
constrains $Y$ decay scenarios compared to Planck+SD. First, the D/H
constraint favors $f_{\gamma}\rightarrow 0$ to minimize changes in
$\eta_{\mathrm{BBN}}$. This agrees with our assessment in Figure 3. The
combination of this preference for $f_{\gamma}\rightarrow 0$ with the
Planck+SD 1D posterior for $f_{\gamma}$ that peaks around $f_{\gamma}\approx
0.632$ results in a relatively flat $f_{\gamma}$ posterior for Planck+SD+D/H.

The constraint on D/H disfavors large values of $max(\rho_{Y}/\rho_{r})$. Non-
zero values of $f_{\gamma}$ mean new photons will be injected by the $Y$ decay
and therefore the photon density at BBN will be decreased accordingly. D/H
prefers to minimize this effect by restricting $max(\rho_{Y}/\rho_{r})$ to be
less than 0.0235 (95% CL). This limit on $max(\rho_{Y}/\rho_{r})$ propagates
to all other parameters; large values of $Y_{\text{He}}$ are no longer
permitted and any corrections to the base six parameters that Planck favored
are no longer needed because $N_{\text{eff}}$ is not deviating much from the
standard value of $3.044$. Therefore, the full Planck+SD+D/H combination is in
excellent agreement with $\Lambda$CDM values for the base six parameters
{$\omega_{b}$, $\omega_{cdm}$, $\theta_{s}$, $A_{s}$, $n_{s}$, $\tau_{reio}$}
and $H_{0}$ (see Table 2).

Figure 9: Comparison of the bounds on $f_{\gamma}$ and
$max(\rho_{Y}/\rho_{r})$ from different implementations of constraints. The
hatched region is the overlap of the $N_{\text{eff}}$ and
$\eta_{\mathrm{BBN}}$ contours in Figure 3 which held $\omega_{b}$ fixed. The
dashed orange lines show the 68% and 95% contours from a Gaussian likelihood
constructed with the 2018 Planck bounds of
$N_{\text{eff}}=2.92^{+0.36}_{-0.37}$ and a Gaussian likelihood created with
the bounds $\eta_{\mathrm{BBN}}=(6.119\pm 0.100)\times 10^{-10}$. Using the
entire Planck+D/H likelihood (filled $68\%$ and $95\%$ contours) incorporates
any uncertainty in cosmological parameters and the asymmetric preference
Planck has for large $f_{\gamma}$.

The posterior of $\Gamma_{Y}$ for Planck+SD+D/H extends to smaller decay rates
than that of Planck+SD. This is a symptom of the D/H bounds heavily
constraining $max(\rho_{Y}/\rho_{r})$. Small decay rates ($\Gamma_{Y}\lesssim
10^{-6}\,$\text{\,}{\mathrm{s}}^{-1}$$) were disfavored by Planck+SD because
$max(\rho_{Y}/\rho_{r})$ was free to vary up to a value of 0.8. However, by
restricting $max(\rho_{Y}/\rho_{r})$ to smaller values, the addition of D/H
made smaller decay rates more probable. Therefore, the combination of SD+D/H
constraints leads to a slightly broader posterior on $\Gamma_{Y}$ than that
obtained from SDs alone. Otherwise, the addition of D/H has no influence on
$\Gamma_{Y}$. While the combination of $\Gamma_{Y}$ and
$max(\rho_{Y}/\rho_{r})$ does influence how much the $Y$ particle increases
the expansion rate during BBN, the D/H abundance is far more sensitive to
changes in $\eta_{\mathrm{BBN}}$ than changes to $H$.

While stringent, the addition of bounds on D/H does not make constraints from
Planck obsolete. Figure 9 shows a comparison of the 2D posterior bounds on
$f_{\gamma}$ and $max(\rho_{Y}/\rho_{r})$ for different implementations of
constraints. The hatched region corresponds to the overlap between the
$N_{\text{eff}}$ and $\eta_{\mathrm{BBN}}$ contours seen in Figure 3. These
limits were predicted using a fixed value of $\omega_{b}$, so this hatched
region can be considered the most naive bounds on $f_{\gamma}$ and
$max(\rho_{Y}/\rho_{r})$. The dashed orange lines illustrate the 68% and 95%
contours resulting from constraining $f_{\gamma}$ and $max(\rho_{Y}/\rho_{r})$
with a Gaussian likelihood constructed with the 2018 Planck bounds of
$N_{\text{eff}}=2.92^{+0.36}_{-0.37}$ and a Gaussian likelihood created with
the bounds $\eta_{\mathrm{BBN}}=(6.119\pm 0.100)\times 10^{-10}$. These
contours were also created for a fixed value of $\omega_{b}$. Here we see that
marginalizing the bounds on $N_{\text{eff}}$ and $\eta_{\mathrm{BBN}}$ through
an MCMC approach extends the allowed $max(\rho_{Y}/\rho_{r})-f_{\gamma}$
parameter space compared to that of the naive bounds on $N_{\text{eff}}$ and
$\eta_{\mathrm{BBN}}$. The filled contours in Figure 9 show the 2D posteriors
for $max(\rho_{Y}/\rho_{r})$ and $f_{\gamma}$ when constraining with a full
Planck likelihood rather than just a Gaussian likelihood on $N_{\text{eff}}$,
as well as Gaussian likelihood with $1\sigma$ bounds $(\text{D/H})=(2.527\pm
0.068)\times 10^{-5}$. The full Planck analysis not only accounts for
uncertainties in cosmological parameters such as $\omega_{b}$ and $n_{s}$, but
also includes the interplay between $Y_{\text{He}}$ and $N_{\text{eff}}$
discussed in Section III.1. The Planck preference for $f_{\gamma}\approx
0.632$ permits values of $f_{\gamma}$ that are otherwise ruled out by the
bounds on just $N_{\text{eff}}$ and $\eta_{\mathrm{BBN}}$. Therefore, the
effects of the $Y$ decay on $N_{\text{eff}}$ and $Y_{\text{He}}$ still
manifest in the Planck+D/H posteriors even though $max(\rho_{Y}/\rho_{r})$ is
significantly restricted compared to the Planck posteriors.

Figure 10: Comparison between $H_{0}$ posteriors of Planck+SD likelihood for
$T_{\mathrm{RH}}=[9.5\times 10^{-4},10^{-2}]$ MeV and either $f_{\gamma}=0$
(dotted maroon outline) or $f_{\gamma}=0.01$ (filled purple contour) with
$68\%$ and $95\%$ CL contours. We include the 1D $H_{0}$ posterior for
$\Lambda$CDM constrained by Planck (dashed gray line) and the vertical band
shows the $1\sigma$ and $2\sigma$ bounds determined by S$H_{0}$ES
($H_{0}=$73.30\pm
1.04\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$$).

Finally, we note that even though this analysis considered the range of
$f_{\gamma}=[0.01,1]$, our results can be extended to smaller values of
$f_{\gamma}$ for the limited range of $\Gamma_{Y}=[3.43\times
10^{-7},3.80\times 10^{-5}]\,$\text{\,}{\mathrm{s}}^{-1}$$. Figure 11 in
Appendix C demonstrates that the posteriors for all parameters are indeed
identical when comparing $f_{\gamma}=0$ and $f_{\gamma}=0.01$ for
$T_{\mathrm{RH}}=[9.5\times 10^{-4},10^{-2}]$ MeV ( $\Gamma_{Y}=[3.43\times
10^{-7},3.80\times 10^{-5}]\,$\text{\,}{\mathrm{s}}^{-1}$$). We show a subset
of these Planck+SD posteriors in Figure 10 and include the $1\sigma$ and
$2\sigma$ bounds on $H_{0}$ (vertical bands) reported by the S$H_{0}$ES
collaboration [69]. When only considering Planck+SD, both the $f_{\gamma}=0$
and $f_{\gamma}=0.01$ cases favor values of $H_{0}$ larger than the best-fit
$\Lambda$CDM value of
$H_{0}=$67.3\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$$.
While the addition of D/H bounds will restrict $max(\rho_{Y}/\rho_{r})$ in the
$f_{\gamma}=0.01$ case so that the $H_{0}$ posterior is in further
disagreement with S$H_{0}$ES, $f_{\gamma}=0$ decays will not affect
$\eta_{\mathrm{BBN}}$ and thus we do not expect the addition of D/H
constraints to substantially change the Planck+SD constraints for
$f_{\gamma}=0$ seen in Figure 10. For $f_{\gamma}=0$, the only additional
constraints that D/H would place on $max(\rho_{Y}/\rho_{r})$ would come from
the increased expansion rate during BBN changing D/H. This change in $H$,
however, has a small effect on D/H and can be reduced by decreasing
$T_{\mathrm{RH}}$. Therefore, it is likely that the Planck+SD+D/H posteriors
for $f_{\gamma}=0$ decays with colder reheat temperatures would not
significantly differ from those for Planck+SD shown in Figure 10. As mentioned
in Section III.3, investigating reheat temperatures
$T_{\mathrm{RH}}\lesssim$9\text{\times}{10}^{-4}\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$
would require consideration of perturbations. It is not yet clear if the
inclusion of $Y$ perturbations and corrections to the perturbations of ultra-
relativistic species would relax the Planck+SD+D/H bounds on an $f_{\gamma}=0$
decay.

|  $\Lambda$CDM | Planck | Planck+SD | Planck+SD+D/H
---|---|---|---|---
$f_{\gamma}$ | – | $0.66^{+0.26}_{-0.21}$ | $0.66^{+0.32}_{-0.17}$ | –
$\log_{10}(\Gamma_{Y}/$\text{\,}{\mathrm{s}}^{-1}$)$ | – | – | $>-5.72$ | $>-6.15$
$\log_{10}(T_{\mathrm{RH}}/\mathrm{MeV})$ | – | – | $>-2.65$ | $>-2.86$
$max(\rho_{Y}/\rho_{r})$ | – | – | – | $<0.0235$
$H_{0}\,[$\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$]$ | $67.3^{+1.2}_{-1.2}$ | $66.2^{+2.8}_{-3.0}$ | $66.2^{+2.9}_{-3.0}$ | $67.4^{+1.7}_{-1.5}$
$Y_{\text{He}}$ | $0.24683^{+0.00013}_{-0.00012}$ | $0.2508^{+0.0050}_{-0.0041}$ | $0.2508^{+0.0055}_{-0.0042}$ | $0.24696^{+0.00026}_{-0.00021}$
$10^{5}\mathrm{D/H}$ | $2.522^{+0.056}_{-0.056}$ | $1.52^{+0.96}_{-0.93}$ | $1.58^{+0.92}_{-0.97}$ | $2.484^{+0.077}_{-0.087}$
$10^{-2}\omega_{b}$ | $2.234^{+0.030}_{-0.029}$ | $2.225^{+0.040}_{-0.042}$ | $2.226^{+0.041}_{-0.043}$ | $2.232^{+0.030}_{-0.029}$
$\omega_{cdm}$ | $0.1202^{+0.0027}_{-0.0027}$ | $0.1172^{+0.0065}_{-0.0067}$ | $0.1172^{+0.0061}_{-0.0067}$ | $0.1208^{+0.0037}_{-0.0034}$
$100\theta_{s}$ | $1.04186^{+0.00058}_{-0.00058}$ | $1.0424^{+0.0013}_{-0.0011}$ | $1.0424^{+0.0013}_{-0.0011}$ | $1.04176^{+0.00066}_{-0.00075}$
$\ln 10^{10}A_{s}$ | $3.045^{+0.033}_{-0.031}$ | $3.036^{+0.036}_{-0.037}$ | $3.037^{+0.036}_{-0.038}$ | $3.045^{+0.033}_{-0.031}$
$n_{s}$ | $0.9642^{+0.0087}_{-0.0087}$ | $0.959^{+0.015}_{-0.016}$ | $0.959^{+0.015}_{-0.016}$ | $0.965^{+0.010}_{-0.0096}$
$\tau_{reio}$ | $0.054^{+0.016}_{-0.015}$ | $0.053^{+0.015}_{-0.015}$ | $0.054^{+0.016}_{-0.015}$ | $0.054^{+0.016}_{-0.015}$
Table 2: Mean and $2\sigma$ errors on decay and cosmological parameters from
MCMC analysis with different combinations of data sets.

##  VI Summary and Conclusions

In this work, we test the robustness of the common assumption of entropy
conservation between BBN and recombination by considering the decay of a
massive hidden sector $Y$ particle. We assume the hidden sector is
sufficiently cold compared to the visible sector of Standard Model particles
such that the particle is non-relativistic during BBN. These decay scenarios
alter the effective number of relativistic species at recombination,
$N_{\text{eff}}$, as well as the baryon-to-photon ratio and expansion rate at
BBN. By employing observations of Planck temperature and polarization
anisotropies, CMB spectral distortions, and the primordial deuterium
abundance, we determine constraints on the decay rate of the $Y$ particle
($\Gamma_{Y}$), its maximum contribution to the energy density of the universe
($max(\rho_{Y}/\rho_{r})$), and the fraction of the $Y$ particle’s energy
density that is transferred to photons ($f_{\gamma}$).

If the $Y$ particle decays into a mixture of photons and other ultra-
relativistic particles such as dark radiation, then there is a sweet spot
photon fraction, $f_{\gamma}=0.5913$, that will keep $N_{\text{eff}}$ fixed at
the fiducial 3.044. However, the injection of photons from the $Y$ particle
leads to an increase in $\eta_{\mathrm{BBN}}$ due to the requirement that
$T_{0}=$2.7255\text{\,}\mathrm{K}$$ (see Appendix A). This increase in
$\eta_{\mathrm{BBN}}$ increases $Y_{\text{He}}$ which results in Planck
observations of the CMB temperature and polarization anisotropies preferring a
larger photon fraction, with the 1D posterior for $f_{\gamma}$ peaking at
0.632. For $f_{\gamma}$ between 0.555 and 0.707, CMB anisotropies permit
significant entropy injections in which the energy density of the $Y$ particle
equals at least half of the radiation density prior to its decay. If
$f_{\gamma}$ deviates from the sweet spot value, then cosmological parameters
like $\omega_{cdm}$, $H_{0}$, $A_{s}$, $\omega_{b}$, and $n_{s}$ must be
altered from their standard $\Lambda$CDM values in order to match the observed
Planck temperature power spectrum.

Upper limits on $\mu$ and $y$ spectral distortions of the CMB restrict the
decay rate of the $Y$ particle to
$\Gamma_{Y}=$1.91\text{\times}{10}^{-6}\text{\,}{\mathrm{s}}^{-1}$$ (95% CL)
such that maximum lifetime of the $Y$ particle is
$5.25\text{\times}{10}^{5}\text{\,}\sec$. The Primordial Inflation Explorer
(PIXIE), which proposes to measure the CMB blackbody spectrum with three
orders of magnitude better sensitivity than the results of COBE/FIRAS [70],
would enhance bounds on $\mu$ and $y$ distortions and therefore place even
more restrictive bounds on the lifetime of the particle.

The deuterium abundance, D/H, proves to have the most constraining power on
$Y$-induced entropy injection. We apply bounds of $\mathrm{D/H}=(2.527\pm
0.068)\times 10^{-5}$, which we translate from bounds on the baryon-to-photon
ratio reported by Cooke  _et al._ [43]. This uncertainty in D/H includes
measurement uncertainty and uncertainties in nuclear reaction rates. This
bound on deuterium limits the energy density of the $Y$ particle to be no more
than 2.35% (95% CL) the total energy density of radiation. Using the fit
described in Figure 2, this upper bound corresponds to $g<1.05$, which
translates to a 5% increase in the comoving energy density of radiation.
Limits on D/H will likely tighten as more absorption systems are observed and
as laboratory measurements of nuclear reaction rates improve, leading to even
more restrictive bounds on the allowed level of entropy injection.

Photodisintegration of light nuclei by the injection of photons can generate
even tighter constraints on entropy injection scenarios. Balázs  _et al._ [27]
recently performed an analysis in which they explored the constraints that
could be placed on axion-like particles with masses greater than 2.2 MeV
decaying into photons between BBN and recombination by employing CMB
anisotropy and spectral distortion data as well as detailed BBN constraints
via $\eta_{\mathrm{BBN}}$ and photodisintegration of nuclei. The addition of
photodisintegration restricts the maximum energy density of the axion-like
particle, $\rho_{a}$, to be $\rho_{a}/\rho_{tot}<10^{-3}$. A decaying particle
with a mass below the binding energy of deuterium, 2.2 MeV, evades these
photodisintegration limits; the stringent BBN bounds on
$max(\rho_{Y}/\rho_{r})$ derived here only consider the effects of a $Y$ decay
on $\eta_{\mathrm{BBN}}$ and the expansion rate during BBN. We find
$max(\rho_{Y}/\rho_{r})<0.0235$ (95% CL), demonstrating that a small level of
entropy injection could be possible in the regime of a sub-$2.2$ MeV decay for
which photodisintegration bounds are irrelevant.

Even though bounds on deuterium proved to be the most restrictive on a general
entropy injection scenario, we demonstrated that there exist subtleties in the
dependence of the decay on degeneracies between $N_{\text{eff}}$ and
$Y_{\text{He}}$ that require a full CMB analysis (Figure 9). Furthermore, we
investigated the effects that entropy injection has on standard cosmological
parameters, including $H_{0}$. In doing so, we find that any injection that
includes photons is heavily restricted by observations of the CMB and D/H such
that
$H_{0}=67.4^{+1.7}_{-1.5}\,$\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$$
($95\%$ CL).

We marginalized results over the range of $0.01\leq f_{\gamma}\leq 1$ and
therefore we cannot speak with absolute confidence regarding an $f_{\gamma}=0$
injection scenario. If $f_{\gamma}=0$, then spectral distortions would not
restrict the $Y$ particle lifetime (see Figure 6), allowing for reheat
temperatures at which the scales accessible to Planck enter the horizon. It
would therefore be necessary to consider perturbation dynamics for the $Y$
particle and how its decay products affect the evolution of other
perturbations. However, we are able to extend our results to $f_{\gamma}=0$
for a limited range of reheat temperatures; when considering
$T_{\mathrm{RH}}=[9.5\times 10^{-4},10^{-2}]$ MeV ( $\Gamma_{Y}=[3.43\times
10^{-7},3.80\times 10^{-5}]\,$\text{\,}{\mathrm{s}}^{-1}$$), we can neglect
any influences the $Y$ particle has on perturbations. For this limited range
of reheat temperatures, we show that the posteriors for $f_{\gamma}=0$ and
$f_{\gamma}=0.01$ are identical when constraining with Planck anisotropies and
spectral distortions (Figure 11). While the addition of bounds on D/H would
restrict $max(\rho_{Y}/\rho_{r})$ in the $f_{\gamma}=0.01$ case, decays with
$f_{\gamma}=0$ would leave $\eta_{\mathrm{BBN}}$ unaltered and only minimally
affect the deuterium abundance by increasing the Hubble rate during BBN.
Therefore, scenarios with the $Y$ particle decaying fully into dark radiation
are less restricted than the $f_{\gamma}\neq 0$ cases studied in this work.
Such decays into dark radiation increase $H_{0}$ and begin to mitigate the
Hubble tension. When considering only Planck anisotropies and spectral
distortions, $f_{\gamma}=0$ decays in this limited range of decay rates result
in
$H_{0}=68.1^{+1.9}_{-1.8}\,$\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{M\text{pc}}}^{-1}$$
($95\%$ CL), which only slightly reduces the Hubble tension. As was
demonstrated by the progression of early-dark energy models [71, 72, 73],
however, the inclusion of perturbations can lead to significant changes in
results. For long-lived $Y$ scenarios, altering perturbations could
potentially relax the bounds on the injection of dark radiation and alleviate
the Hubble tension even further. We reserve an in-depth investigation of the
$f_{\gamma}=0$ case for a future work.

We have demonstrated that the injection of new radiation between BBN and
recombination, be it photons or other ultra-relativistic species such as dark
radiation, is highly constrained: a massive particle whose decay products
include photons can make up at most $2.35\%$ of the energy density of the
universe. This stringent limit implies that cosmological parameters, including
the Hubble constant, cannot be adjusted by introducing changes in the photon
abundance between BBN and recombination. Our results uphold the standard
cosmological picture and the assumption that entropy is conserved.

## Acknowledgments

This analysis employed the Longleaf Computing Cluster owned by the University
of North Carolina at Chapel Hill. ACS and ALE are supported in part by NSF
CAREER grant PHY-1752752. ACS also acknowledges support from the NC Space
Grant Consortium and the Bahnson Fund at UNC Chapel Hill. TLS is supported by
NSF Grant No. 2009377, NASA Grant No. 80NSSC18K0728, and the Research
Corporation.

##  Appendix A Initial conditions for CLASS

Here we describe how to recast our set of descriptive parameters for the $Y$
decay model into a set of initial conditions. We start by making the
approximation that $H(a)\approx H_{i}(a/a_{i})^{-2}$ during radiation
domination (RD). This assumption, along with taking $\rho_{Y}=\rho_{Y,i}$ when
$a=a_{i}$, allows us to analytically solve equation (1):

$\displaystyle\rho_{Y}(a)=\rho_{Y,i}\left(\frac{a_{i}}{a}\right)^{3}e^{\frac{1}{2}\widetilde{\Gamma}_{Y}\left[1-(a/a_{i})^{2}\right]}\,.$
(14)

Additionally, we define the reheat scale factor by $\Gamma_{Y}\equiv
H_{i}(a_{i}/a_{\mathrm{RH}})^{2}$:
$a_{\mathrm{RH}}/a_{i}=\widetilde{\Gamma}_{Y}^{-1/2}$, where
$\widetilde{\Gamma}_{Y}\equiv\Gamma_{Y}/H_{i}$. It follows from equation (14)
that the maximum value for $\rho_{Y}/\rho_{r}$ occurs at
$a_{\mathrm{RH}}/a_{i}$. Therefore, an analytical expression for
$max(\rho_{Y}/\rho_{r})$ can be found by evaluating equation (14) at
$a_{\mathrm{RH}}$ and taking
$\rho_{r}(a_{\mathrm{RH}})=\rho_{r,i}(a_{i}/a_{\mathrm{RH}})^{4}$:

$\displaystyle
max\left(\frac{\rho_{Y}}{\rho_{r}}\right)=e^{\frac{1}{2}(\tilde{\Gamma}_{Y}-1)}\frac{\rho_{Y,i}}{\rho_{r,i}}\,\,\tilde{\Gamma}_{Y}^{-1/2}\,.$
(15)

Equation (15) defines the $max(\rho_{Y}/\rho_{r})$ decay parameter used
throughout this work. This analytical expression for $max(\rho_{Y}/\rho_{r})$
deviates from the exact numerical maximum of $\rho_{Y}/\rho_{r}$ as
$max(\rho_{Y}/\rho_{r})$ approaches unity. Comparing equation (15) to the
numerically calculated $max(\rho_{Y}/\rho_{r})$ for various decay scenarios,
we find that equation (15) is accurate within $1\%$ for
$max(\rho_{Y}/\rho_{r})<0.079$ and within $10\%$ for
$max(\rho_{Y}/\rho_{r})<0.945$. Seeing that a value of
$max(\rho_{Y}/\rho_{r})>1$ translates to a period $Y$ domination, it is
consistent that our RD approximation would begin fail as
$max(\rho_{Y}/\rho_{r})$ approaches 1.

In order to determine a value for $\rho_{Y,i}$ from equation (15), we need to
find the initial condition
$\rho_{r,i}=\rho_{\gamma,i}+\rho_{ur,i}+\rho_{ncdm,i}$. For decay scenarios
with $f_{\gamma}\neq 0$, we must decrease the photon energy density pre-decay
such that the photon injection from the $Y$ decay results in a present day
temperature of $T_{0}=$2.7255\text{\,}\mathrm{K}$$. Since we only consider
reheat temperatures after BBN
($T_{\mathrm{RH}}\lesssim$0.01\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$), we
must decrease $\rho_{\gamma}$ at the time of neutrino decoupling and therefore
we must also re-scale the initial energy density of ultra-relativistic
species, $\rho_{ur,i}$. To do so, we parameterize the initial energy densities
by $\Omega^{\prime}_{x,0}$ such that
$\rho_{x,i}=\Omega^{\prime}_{x,0}\rho_{\text{crit},0}\,a_{i}^{-4}$ for species
$x$. In equation (8), we define the ratio of comoving radiation energy density
before and after decay as

$\displaystyle g\equiv\frac{\rho_{r,f}a_{f}^{4}}{\rho_{r,i}a_{i}^{4}}\,,$

where the $i$ and $f$ subscripts denote before and after decay, respectively.
This $g$ parameter is entirely dependent on $max(\rho_{Y}/\rho_{r})$ defined
by equation (15). Let us denote the comoving energy density of species $x$ as
$\bar{\rho}_{x,j}=\rho_{x,j}a_{j}^{4}$. Then the change in comoving radiation
density is

$\displaystyle\Delta\bar{\rho}_{r}=\bar{\rho}_{r,f}-\bar{\rho}_{r,i}=\left(g-1\right)\bar{\rho}_{r,i}\,.$

From this, we can specify the change in individual species with the photon
fraction as

$\displaystyle\Delta\bar{\rho}_{\gamma}$
$\displaystyle=f_{\gamma}\left(g-1\right)\bar{\rho}_{r,i}\,;$
$\displaystyle\Delta\bar{\rho}_{ur}$
$\displaystyle=\left(1-f_{\gamma}\right)\left(g-1\right)\bar{\rho}_{r,i}\,.$

Now we can determine the $g$ factor for each species that receives energy
injection from the decay:

$\displaystyle g_{\gamma}$
$\displaystyle=\frac{\bar{\rho}_{\gamma,f}}{\bar{\rho}_{\gamma,i}}=\frac{\bar{\rho}_{\gamma,i}+\Delta\bar{\rho}_{\gamma}}{\bar{\rho}_{\gamma,i}}=1+\frac{\Delta\bar{\rho}_{\gamma}}{\bar{\rho}_{\gamma,i}}=1+f_{\gamma}\left(g-1\right)\frac{\rho_{r,i}}{\rho_{\gamma,i}}\,,$
$\displaystyle g_{ur}$
$\displaystyle=\frac{\bar{\rho}_{ur,f}}{\bar{\rho}_{ur,i}}=\frac{\bar{\rho}_{ur,i}+\Delta\bar{\rho}_{ur}}{\bar{\rho}_{ur,i}}=1+\frac{\Delta\bar{\rho}_{ur}}{\bar{\rho}_{ur,i}}=1+\left(1-f_{\gamma}\right)\left(g-1\right)\frac{\rho_{r,i}}{\rho_{ur,i}}\,.$

We assume the minimal convention in which neutrinos are composed of two
massless species ($ur$) and one massive with $m_{ncdm}=0.06$ eV. The massive
neutrino, or non-cold dark matter (ncdm), will be relativistic at early times
and therefore contribute to the initial radiation energy density. If we denote
$N_{ur}$ as the contribution to $N_{\text{eff}}$ from the ultra-relativistic
species and $N_{ncdm}$ as the contribution from the massive neutrino, then we
enforce that $N_{ur}+N_{ncdm}=3.044$ at the time of BBN. With these
definitions, we can write

$\displaystyle\frac{\rho_{r,i}}{\rho_{\gamma,i}}$
$\displaystyle=\frac{\rho_{ncdm,i}+\rho_{ur,i}+\rho_{\gamma,i}}{\rho_{\gamma,i}}=\frac{\rho_{ncdm,i}}{\rho_{\gamma,i}}+\frac{\rho_{ur,i}}{\rho_{\gamma,i}}+1=\frac{7}{8}N_{ncdm}\left(\frac{4}{11}\right)^{4/3}\\!\\!\\!+\frac{7}{8}N_{ur}\left(\frac{4}{11}\right)^{4/3}\\!\\!\\!+1=\frac{7}{8}(3.044)\left(\frac{4}{11}\right)^{4/3}\\!\\!\\!+1\,\,,$

and

$\displaystyle\frac{\rho_{r,i}}{\rho_{ur,i}}=\frac{\rho_{ncdm,i}+\rho_{ur,i}+\rho_{\gamma,i}}{\rho_{ur,i}}=\frac{\rho_{ncdm,i}}{\rho_{ur,i}}+1+\frac{\rho_{\gamma,i}}{\rho_{ur,i}}$
$\displaystyle=\frac{N_{ncdm}}{N_{ur}}+1+\frac{8}{7}\frac{1}{N_{ur}}\left(\frac{11}{4}\right)^{4/3}\,\,;$
$\displaystyle=\frac{3.044-N_{ur}}{N_{ur}}+1+\frac{8}{7}\frac{1}{N_{ur}}\left(\frac{11}{4}\right)^{4/3}\,\,;$
$\displaystyle=\frac{1}{N_{ur}}\left[3.044+\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\right]\,\,.$

It follows that the ratio of comoving energy density for each species is

$\displaystyle g_{\gamma}$
$\displaystyle=1+f_{\gamma}\left(g-1\right)\left[1+\frac{7}{8}(3.044)\left(\frac{4}{11}\right)^{4/3}\right]\,\,;$
(16) $\displaystyle g_{ur}$
$\displaystyle=1+\frac{(1-f_{\gamma})\left(g-1\right)}{N_{ur}}\left[3.044+\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\right]\,\,.$
(17)

If $\Omega_{\gamma,0}$ is the photon energy density corresponding to
$T_{0}=$2.7255\text{\,}\mathrm{K}$$, then the initial values are found by
$\rho_{x,i}=\Omega^{\prime}_{x,0}\rho_{\text{crit},0}\,a_{i}^{-4}$ and

$\displaystyle\Omega^{\prime}_{\gamma,0}$
$\displaystyle=\frac{\Omega_{\gamma,0}}{g_{\gamma}}=\frac{\Omega_{\gamma,0}}{1+f_{\gamma}\left(g-1\right)\left[1+\frac{7}{8}(3.044)\left(\frac{4}{11}\right)^{4/3}\right]};$
(18) $\displaystyle\Omega^{\prime}_{ur,0}$
$\displaystyle=\frac{7}{8}N_{ur}\left(\frac{4}{11}\right)^{4/3}\Omega^{\prime}_{\gamma,0}\,\,.$
(19)

Since $g$ is purely a function of $max(\rho_{Y}/\rho_{r})$ (Figure 2), we can
determine the initial conditions $\rho_{\gamma,i}$ and $\rho_{ur,i}$ directly
from our decay parameters $f_{Y}$ and $max(\rho_{Y}/\rho_{r})$. In the case of
$f_{\gamma}=0$ (no decays into photons), we do not need to decrease
$\rho_{\gamma}$ at early times. Indeed, plugging $f_{\gamma}=0$ into equation
(18) leads to $\Omega^{\prime}_{\gamma,0}=\Omega_{\gamma,0}$ and the initial
condition is simply $\rho_{\gamma,i}=\Omega_{\gamma,0}\rho_{crit}a_{i}^{-4}$.
The contribution to $N_{\text{eff}}$ from the ultra-relativistic neutrino
species before decay is determined by $N_{ur}=3.044-N_{ncdm}$, where
$N_{ncdm}$ is the contribution that the massive neutrino makes to
$N_{\text{eff}}$ prior to recombination.

Massive neutrino (ncdm) calculations performed by CLASS derive the number
density, energy density, and pressure of massive neutrinos based on the input
temperature $T_{ncdm,0}/T_{\gamma,0}$. By default, CLASS assumes this
temperature to be $T_{ncdm,0}=0.71611\,T_{\gamma,0}$ in order to obtain a
mass-to-density ratio of $m/\omega_{ncdm}=93.14$ eV. However, this default
value inherently assumes that both $T_{ncdm}$ and $T_{\gamma}$ scale as
$T\propto a^{-1}$ after electron-positron annihilation. While this scaling is
indeed true for $T_{ncdm}$, any $Y$ decay scenario that injects photons will
result in $T_{\gamma}$ not consistently scaling as $a^{-1}$. The photon
temperature at some scale factor, $a_{i}$, which occurs between electron-
positron annihilation and energy injection from the $Y$ decay, will depend on
the initial photon energy density determined by equation (18). Specifically,

$\displaystyle T_{\gamma,i}a_{i}$
$\displaystyle=\left(\frac{\Omega^{\prime}_{\gamma,0}}{\Omega_{\gamma,0}}\right)^{1/4}\\!\\!\\!T_{\gamma,0}a_{0},$

where $a_{0}$ is the present-day scale factor. Assuming that
$T_{ncdm,i}=0.71611\,T_{\gamma,i}$ after electron-positron annihilation and
before the $Y$ particle alters the evolution of $T_{\gamma}$,

$\displaystyle T_{ncdm,0}$
$\displaystyle=0.71611\left(\frac{\Omega^{\prime}_{\gamma,0}}{\Omega_{\gamma,0}}\right)^{1/4}\\!\\!\\!T_{\gamma,0}\,\,.$

If there is no photon injection from the decay, then
$\Omega^{\prime}_{\gamma,0}=\Omega_{\gamma,0}$ and we recover the default
assumption of CLASS.

The contribution that the massive neutrino makes to $N_{\text{eff}}$, which we
denote as $N_{ncdm}$, is

$N_{ncdm}=\left(\frac{0.71611}{(4/11)^{1/3}}\right)^{4}=1.0132\,\,.$

This value is used to determine the number of ultra-relativistic species pre-
decay, $N_{ur}$, by enforcing $3.044=N_{ur}+N_{ncdm}$ at BBN.

##  Appendix B Post-decay $N_{\text{eff}}$

We define the post-decay $N_{\text{eff}}$ to be the number of relativistic
species immediately after entropy injection from the $Y$ decay has completed.
For the reheat temperatures that we consider in this work, this means that
massive neutrinos are still relativistic and therefore contribute to
$N_{\text{eff}}$ immediately after decay. As discussed in Section A, the
contribution that the massive neutrino makes to $N_{\text{eff}}$ before the
decay is $N_{ncdm}=(11/4)^{4/3}(0.71611)^{4}=1.0132$. While the $Y$ particle
does not inject any new massive neutrinos, the ratio of
$\rho_{ncdm}/\rho_{\gamma}$ will change due to the decay creating new photons.
Therefore, the contribution that the massive neutrino makes to
$N_{\text{eff}}$ will change post-decay. This post-decay number of massive
neutrinos is

$\displaystyle
N_{ncdm}^{post}=\frac{N_{ncdm}}{g_{\gamma}}=\left(\frac{11}{4}\right)^{4/3}\frac{(0.71611)^{4}}{g_{\gamma}}\,\,,$

where $g_{\gamma}$ is defined by equation (16). Additionally, the effective
number of ultra-relativistic species will change due to the $Y$ decay. If
$N_{ur}$ is the effective number of relativistic species before decay, then we
denote the contribution of the ultra-relativistic species to the post-decay
$N_{\text{eff}}$ as $N_{ur}^{post}$. It follows that

$\displaystyle\Omega_{ur,0}$ $\displaystyle=\Omega^{\prime}_{ur}g_{ur}$
$\displaystyle\frac{7}{8}N_{ur}^{post}\left(\frac{4}{11}\right)^{4/3}\Omega_{\gamma,0}$
$\displaystyle=\frac{7}{8}N_{ur}\left(\frac{4}{11}\right)^{4/3}\Omega^{\prime}_{\gamma,0}g_{ur}$
$\displaystyle N_{ur}^{post}$
$\displaystyle=\frac{\Omega^{\prime}_{\gamma,0}}{\Omega_{\gamma,0}}g_{ur}N_{ur}$
$\displaystyle N_{ur}^{post}$
$\displaystyle=\frac{g_{ur}}{g_{\gamma}}N_{ur}\,\,,$

where $g_{ur}$ is defined by equation (17). Therefore, the total post-decay
$N_{\text{eff}}$ is

$N_{\text{eff}}^{post}=N_{ncdm}^{post}+N_{ur}^{post}=\frac{1}{g_{\gamma}}\left[\left(\frac{11}{4}\right)^{4/3}(0.71611)^{4}+g_{ur}N_{ur}\right]=\frac{1}{g_{\gamma}}\left[3.044+N_{ur}(g_{ur}-1)\right]\,\,,$

where $N_{ur}=2.0308$. Since both $g_{\gamma}$ and $g_{ur}$ are ultimately
functions of $max(\rho_{Y}/\rho_{r})$ and $f_{\gamma}$, we are able to
calculate the post-decay $N_{\text{eff}}$ directly from our decay parameters.

##  Appendix C Additional MCMC results

Figure 11: Comparison between posteriors of Planck+SD likelihood for
$T_{\mathrm{RH}}=[9.5\times 10^{-4},10^{-2}]$ MeV and either $f_{\gamma}=0$
(dashed outline) or $f_{\gamma}=0.01$ (filled) with $68\%$ and $95\%$ CL
contours. The posteriors of $f_{\gamma}=0$ vs. those of $f_{\gamma}=0.01$ are
identical in all parameters except for a slight difference in $\Gamma_{Y}$.
The $f_{\gamma}=0$ case is completely unconstrained by SDs and therefore
equally favors long and short particle lifetimes whereas the $f_{\gamma}=0.01$
case disfavors long lifetimes due to the constraints placed by SDs.  Figure
12: 1D and 2D posterior distributions of decay and full set of cosmological
parameters for different combinations of Planck high-$\ell$ TT,TE,EE,
low-$\ell$ TT, and low-$\ell$ EE (Planck) data, CMB spectral distortions (SD)
bounds, and constraints on the deuterium abundance (D/H). We include the 1D
posteriors for $\Lambda$CDM constrained by Planck (dashed black line). The
dotted gray line traces $f_{\gamma}=0.5913$, which maintains
$N_{\text{eff}}=3.044$ at recombination.

## References

  * Pantano and Riotto [1993] O. Pantano and A. Riotto, Physics Letters B 307, 49 (1993), arXiv:hep-ph/9212214 .
  * Chung  _et al._ [1998] D. J. H. Chung, E. W. Kolb, and A. Riotto, Phys. Rev. D 59, 023501 (1998), arXiv:hep-ph/9802238 .
  * Giudice  _et al._ [2001] G. F. Giudice, E. W. Kolb, and A. Riotto, Phys. Rev. D 64, 023508 (2001), arXiv:hep-ph/0005123 .
  * Allahverdi and Drees [2002] R. Allahverdi and M. Drees, Phys. Rev. Lett. 89, 091302 (2002), arXiv:hep-ph/0203118 .
  * Asaka  _et al._ [2006] T. Asaka, M. Shaposhnikov, and A. Kusenko, Physics Letters B 638, 401 (2006), arXiv:hep-ph/0602150 .
  * Patwardhan  _et al._ [2015] A. V. Patwardhan, G. M. Fuller, C. T. Kishimoto, and A. Kusenko, Phys. Rev. D 92, 103509 (2015), arXiv:1507.01977 .
  * Ichikawa  _et al._ [2007] K. Ichikawa, M. Kawasaki, K. Nakayama, M. Senami, and F. Takahashi, J. Cosmol. Astropart. Phys. 2007 (05), 008, arXiv:hep-ph/0703034 .
  * Fischler and Meyers [2011] W. Fischler and J. Meyers, Phys. Rev. D 83, 063520 (2011), arXiv:1011.3501 .
  * Hasenkamp [2012] J. Hasenkamp, Physics Letters B 707, 121 (2012), arXiv:1107.4319 .
  * Aboubrahim  _et al._ [2022] A. Aboubrahim, M. Klasen, and P. Nath, J. Cosmol. Astropart. Phys. 2022 (04), 042, arXiv:2202.04453 .
  * Aloni  _et al._ [2022] D. Aloni, A. Berlin, M. Joseph, M. Schmaltz, and N. Weiner, Phys. Rev. D 105, 123516 (2022), arXiv:2111.00014 .
  * Bernal  _et al._ [2016] J. L. Bernal, L. Verde, and A. G. Riess, J. Cosmol. Astropart. Phys. 2016 (10), 019, arXiv:1607.05617 .
  * Schöneberg  _et al._ [2021] N. Schöneberg, G. F. Abellán, A. P. Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, arXiv e-prints (2021), arXiv:2107.10291 .
  * Di Valentino  _et al._ [2021] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Classical and Quantum Gravity 38, 153001 (2021).
  * de Salas  _et al._ [2015] P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor, and O. Pisanti, Phys. Rev. D 92, 123534 (2015), arXiv:1511.00672 .
  * Hasegawa  _et al._ [2019] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram, and S. Hannestad, J. Cosmol. Astropart. Phys. 2019 (12), 012, arXiv:1908.10189 .
  * Lindley [1985] D. Lindley, The Astrophysical Journal 294, 1 (1985).
  * Scherrer and Turner [1988a] R. J. Scherrer and M. S. Turner, The Astrophysical Journal 331, 33 (1988a).
  * Scherrer and Turner [1988b] R. J. Scherrer and M. S. Turner, The Astrophysical Journal 331, 19 (1988b).
  * Frieman  _et al._ [1990] J. A. Frieman, E. W. Kolb, and M. S. Turner, Phys. Rev. D 41, 3080 (1990).
  * Protheroe  _et al._ [1995] R. J. Protheroe, T. Stanev, and V. S. Berezinsky, Nuclear Physics B - Proceedings Supplements 43, 62 (1995).
  * Hufnagel  _et al._ [2018a] M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, J. Cosmol. Astropart. Phys. 2018 (11), 032, arXiv:1808.09324 .
  * Kawasaki  _et al._ [2020] M. Kawasaki, K. Kohri, T. Moroi, K. Murai, and H. Murayama, J. Cosmol. Astropart. Phys. 2020 (12), 048, arXiv:2006.14803 .
  * Cadamuro and Redondo [2012] D. Cadamuro and J. Redondo, J. Cosmol. Astropart. Phys. 2012 (02), 032, arXiv:1110.2895 .
  * Millea  _et al._ [2015] M. Millea, L. Knox, and B. Fields, Phys. Rev. D 92, 023010 (2015), arXiv:1501.04097 .
  * Depta  _et al._ [2020] P. F. Depta, M. Hufnagel, and K. Schmidt-Hoberg, J. Cosmol. Astropart. Phys. 2020 (05), 009, arXiv:2002.08370 .
  * Balázs  _et al._ [2022] C. Balázs, S. Bloor, T. E. Gonzalo, W. Handley, S. Hoof, F. Kahlhoefer, M. Lecroq, D. J. E. Marsh, J. J. Renk, P. Scott, and P. Stöcker, arXiv e-prints (2022), arXiv:2205.13549 .
  * Poulin and Serpico [2015] V. Poulin and P. D. Serpico, Phys. Rev. D 91, 103007 (2015), arXiv:1503.04852 .
  * Kawasaki  _et al._ [2018] M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, Phys. Rev. D 97, 023502 (2018), arXiv:1709.01211 .
  * Forestell  _et al._ [2019] L. Forestell, D. E. Morrissey, and G. White, J. High Energ. Phys. 2019 (1), 74, arXiv:1809.01179 .
  * Poulin  _et al._ [2017] V. Poulin, J. Lesgourgues, and P. D. Serpico, J. Cosmol. Astropart. Phys. 2017 (03), 043, arXiv:1610.10051 .
  * Kolb  _et al._ [1985] E. W. Kolb, D. Seckel, and M. S. Turner, Nature 314, 415 (1985).
  * Berezhiani  _et al._ [1996] Z. G. Berezhiani, A. D. Dolgov, and R. N. Mohapatra, Physics Letters B 375, 26 (1996), arXiv:hep-ph/9511221 .
  * Feng  _et al._ [2008] J. L. Feng, H. Tu, and H.-B. Yu, J. Cosmol. Astropart. Phys. 2008 (10), 043, arXiv:0808.2318 .
  * Adshead  _et al._ [2016] P. Adshead, Y. Cui, and J. Shelton, J. High Energ. Phys. 2016 (6), 16, arXiv:1604.02458 .
  * Adshead  _et al._ [2019] P. Adshead, P. Ralegankar, and J. Shelton, J. High Energ. Phys. 2019 (8), 151, arXiv:1906.02755 .
  * Anchordoqui [2021] L. A. Anchordoqui, Phys. Rev. D 103, 035025 (2021), arXiv:2010.09715 .
  * Nygaard  _et al._ [2021] A. Nygaard, T. Tram, and S. Hannestad, J. Cosmol. Astropart. Phys. 2021 (05), 017, arXiv:2011.01632 .
  * Zhang [2015] Y. Zhang, J. Cosmol. Astropart. Phys. 2015 (05), 008, arXiv:1502.06983 .
  * Berlin  _et al._ [2016] A. Berlin, D. Hooper, and G. Krnjaic, Phys. Rev. D 94, 095019 (2016), arXiv:1609.02555 .
  * Dror  _et al._ [2016] J. A. Dror, E. Kuflik, and W. H. Ng, Phys. Rev. Lett. 117, 211801 (2016), arXiv:1607.03110 .
  * Planck Collaboration  _et al._ [2020a] Planck Collaboration, N. Aghanim, _et al._ , A&A 641, A5 (2020a), arXiv:1907.12875 .
  * Cooke  _et al._ [2018] R. Cooke, M. Pettini, and C. C. Steidel, ApJ 855, 102 (2018), arXiv:1710.11129 .
  * Mather  _et al._ [1994] J. C. Mather, E. S. Cheng, D. A. Cottingham, R. E. Eplee, Jr., D. J. Fixsen, T. Hewagama, R. B. Isaacman, K. A. Jensen, S. S. Meyer, P. D. Noerdlinger, S. M. Read, _et al._ , The Astrophysical Journal 420, 439 (1994).
  * Fixsen  _et al._ [1996] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, ApJ 473, 576 (1996), arXiv:astro-ph/9605054 .
  * Samsing  _et al._ [2012] J. Samsing, E. V. Linder, and T. L. Smith, Phys. Rev. D 86, 123504 (2012), arXiv:1208.4845 .
  * Hou  _et al._ [2013] Z. Hou, R. Keisler, L. Knox, M. Millea, and C. Reichardt, Phys. Rev. D 87, 083008 (2013), arXiv:1104.2333 .
  * Planck Collaboration  _et al._ [2014] Planck Collaboration, P. A. R. Ade, _et al._ , A&A 571, A16 (2014), arXiv:1303.5076 .
  * Bashinsky and Seljak [2004] S. Bashinsky and U. Seljak, Phys. Rev. D 69, 083002 (2004), arXiv:astro-ph/0310198 .
  * Follin  _et al._ [2015] B. Follin, L. Knox, M. Millea, and Z. Pan, Phys. Rev. Lett. 115, 091301 (2015), arXiv:1503.07863 .
  * Baumann  _et al._ [2016] D. Baumann, D. Green, J. Meyers, and B. Wallisch, J. Cosmol. Astropart. Phys. 2016 (01), 007, arXiv:1508.06342 .
  * Chluba and Jeong [2014] J. Chluba and D. Jeong, Monthly Notices of the Royal Astronomical Society 438, 2065 (2014), arXiv:1306.5751 .
  * Lucca  _et al._ [2020] M. Lucca, N. Schöneberg, D. C. Hooper, J. Lesgourgues, and J. Chluba, J. Cosmol. Astropart. Phys. 2020 (02), 026, arXiv:1910.04619 .
  * Bolliet  _et al._ [2021] B. Bolliet, J. Chluba, and R. Battye, Monthly Notices of the Royal Astronomical Society 507, 3148 (2021), arXiv:2012.07292 .
  * Erickcek  _et al._ [2022] A. L. Erickcek, P. Ralegankar, and J. Shelton, J. Cosmol. Astropart. Phys. 2022 (01), 017, arXiv:2106.09041 .
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, J. Cosmol. Astropart. Phys. 2011 (07), 034, arXiv:1104.2933 .
  * Gariazzo  _et al._ [2022] S. Gariazzo, P. F. de Salas, O. Pisanti, and R. Consiglio, Computer Physics Communications 271, 108205 (2022), arXiv:2103.05027 .
  * Hufnagel  _et al._ [2018b] M. Hufnagel, K. Schmidt-Hoberg, and S. Wild, J. Cosmol. Astropart. Phys. 2018 (02), 044, arXiv:1712.03972 .
  * Zeldovich and Sunyaev [1969] Y. B. Zeldovich and R. A. Sunyaev, Astrophys Space Sci 4, 301 (1969).
  * Sunyaev and Zeldovich [1970] R. A. Sunyaev and Y. B. Zeldovich, Astrophys Space Sci 7, 20 (1970).
  * Illarionov and Siuniaev [1975] A. F. Illarionov and R. A. Siuniaev, Soviet Astronomy 18, 691 (1975).
  * Chluba and Sunyaev [2012] J. Chluba and R. A. Sunyaev, Monthly Notices of the Royal Astronomical Society 419, 1294 (2012), arXiv:1109.6552 .
  * Audren  _et al._ [2013] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, J. Cosmol. Astropart. Phys. 2013 (02), 001, arXiv:1210.7183 .
  * Brinckmann and Lesgourgues [2018] T. Brinckmann and J. Lesgourgues, arXiv e-prints (2018), arXiv:1804.07261 .
  * Planck Collaboration  _et al._ [2020b] Planck Collaboration, N. Aghanim, _et al._ , A&A 641, A6 (2020b), arXiv:1807.06209 .
  * Adelberger  _et al._ [2011] E. G. Adelberger, A. García, R. G. H. Robertson, K. A. Snover, A. B. Balantekin, K. Heeger, M. J. Ramsey-Musolf, D. Bemmerer, A. Junghans, C. A. Bertulani, _et al._ , Rev. Mod. Phys. 83, 195 (2011).
  * Gelman and Rubin [1992] A. Gelman and D. B. Rubin, Statistical Science 7, 457 (1992).
  * Lewis [2019] A. Lewis, arXiv e-prints (2019), arXiv:1910.13970 .
  * Riess  _et al._ [2022] A. G. Riess, W. Yuan, L. M. Macri, D. Scolnic, D. Brout, S. Casertano, D. O. Jones, Y. Murakami, L. Breuval, T. G. Brink, _et al._ , arXiv e-prints (2022), arXiv:2112.04510 .
  * Kogut  _et al._ [2011] A. Kogut, D. J. Fixsen, D. T. Chuss, J. Dotson, E. Dwek, M. Halpern, G. F. Hinshaw, S. M. Meyer, S. H. Moseley, M. D. Seiffert, _et al._ , J. Cosmol. Astropart. Phys. 2011 (07), 025, arXiv:1105.2044 .
  * Karwal and Kamionkowski [2016] T. Karwal and M. Kamionkowski, Phys. Rev. D 94, 103523 (2016), arXiv:1608.01309 .
  * Poulin  _et al._ [2019] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 .
  * Smith  _et al._ [2020] T. L. Smith, V. Poulin, and M. A. Amin, Phys. Rev. D 101, 063523 (2020), arXiv:1908.06995 .