# A Step in Understanding the Hubble Tension

Daniel Aloni  Physics Department, Boston University, Boston, MA 02215, USA
Asher Berlin  Center for Cosmology and Particle Physics, Department of
Physics, New York University, New York, NY 10003, USA  Melissa Joseph  Physics
Department, Boston University, Boston, MA 02215, USA  Martin Schmaltz  Physics
Department, Boston University, Boston, MA 02215, USA  Neal Weiner  Center for
Cosmology and Particle Physics, Department of Physics, New York University,
New York, NY 10003, USA

###### Abstract

As cosmological data have improved, tensions have arisen. One such tension is
the difference between the locally measured Hubble constant $H_{0}$ and the
value inferred from the cosmic microwave background (CMB). Interacting
radiation has been suggested as a solution, but studies show that conventional
models are precluded by high-$\ell$ CMB polarization data. It seems at least
plausible that a solution may be provided by related models that distinguish
between high- and low-$\ell$ multipoles. When interactions of strongly-coupled
radiation are mediated by a force-carrier that becomes non-relativistic, the
dark radiation undergoes a “step” in which its relative energy density
increases as the mediator deposits its entropy into the lighter species. If
this transition occurs while CMB-observable modes are inside the horizon,
high- and low-$\ell$ peaks are impacted differently, corresponding to modes
that enter the horizon before or after the step. These dynamics are naturally
packaged into the simplest supersymmetric theory, the Wess-Zumino model, with
the mass of the scalar mediator near the eV-scale. We investigate the
cosmological signatures of such “Wess-Zumino Dark Radiation” (WZDR) and find
that it provides an improved fit to the CMB alone, favoring larger values of
$H_{0}$. If supernovae measurements from the SH0ES collaboration are also
included in the analysis, the inferred value of $H_{0}$ is yet larger, but the
preference for dark radiation and the location of the transition is left
nearly unchanged. Utilizing a standardized set of measures, we compare to
other models and find that WZDR is among the most successful at addressing the
$H_{0}$ tension and the best of those with a Lagrangian formulation.

###### pacs:

95.35.+d

##  I Introduction

The first results from WMAP harkened the arrival of the era of precision
cosmology. The many telescopes - terrestrial and orbital - since then have
allowed us to test the detailed nature of the evolution of the Universe. While
this has given us incredibly precise measurements of the parameters of
$\Lambda{\rm CDM}$, it has also allowed us to study whether small deviations
from standard cosmology might appear.

As the data have come in, some tensions have begun to emerge. The most notable
of these is the disagreement between the value of $H_{0}$ inferred from Planck
observations of the CMB Planck:2018vyg  and the value extracted from late-
universe supernovae data by the SH0ES collaboration Riess:2020fzl . Indeed, it
has been argued that the early-universe inferred value shows a clear
distinction from a broad set of late-universe measurements Riess:2020fzl ;
Freedman:2019jwv ; Freedman:2021ahq ; Yuan:2019npk ; Soltis:2020gpl ;
Wong:2019kwg ; Pesce:2020xfe ; LIGOScientific:2019zcs .111There are a variety
of measurements of $H_{0}$, including alternative calibrations of supernovae.
To directly compare to Ref. Schoneberg:2021qvd , we focus on the $H_{0}$ value
inferred from SH0ES Riess:2020fzl . At the same time, there is a well-known,
smaller tension with $\Lambda{\rm CDM}$ within the Planck CMB data itself;
compared to the low-$\ell$ multipoles, the high-$\ell$ multipoles in the TT
power spectrum prefer smaller values of $H_{0}$ (and larger values of the
matter density $\Omega_{m}h^{2}$) Addison:2015wyg ; Planck:2018vyg .

There is no clear solution to these tensions. It has been argued that the
discrepancy in $H_{0}$ is due to a deviation in the sound horizon of
$\Lambda{\rm CDM}$ Bernal:2016gxb ; Aylor:2018drw . As a consequence, many
proposed solutions involve the addition of new components to the cosmological
energy density before matter-radiation equality, which has the effect of
decreasing the sound horizon at the time of recombination, thus increasing the
Hubble rate inferred from the CMB (see, e.g., Refs. Brust:2017nmv ;
Blinov:2020hmc ; RoyChoudhury:2020dmd ; Brinckmann:2020bcn ; Kreisch:2019yzn ;
Escudero:2019gvw ; EscuderoAbenza:2020egd ; Escudero:2021rfi ; Karwal:2016vyq
; Poulin:2018cxd ; Lin:2019qug ; Smith:2019ihp ; Bansal:2021dfh ; Cyr-
Racine:2021alc  for a representative set). Additional radiation is one such
possibility Bernal:2016gxb ; Blinov:2020hmc ; RoyChoudhury:2020dmd ;
Brinckmann:2020bcn ; Kreisch:2019yzn , although at the cost of suppressing
power at high multipoles of the CMB due to enhanced Silk damping Hou:2011ec .
The situation is worsened if this radiation is free-streaming, since in this
case the acoustic peaks in the CMB angular power spectrum are additionally
suppressed in their amplitude as well as shifted in their location to smaller
multipoles by the so-called “drag effect” Bashinsky:2003tk ; Baumann:2015rya ;
Follin:2015hya . The drag effect is significantly suppressed if the radiation
is instead strongly-coupled to itself, which allows for a more sizable shift
to the sound horizon in models of interacting radiation. However, high-$\ell$
polarization measurements of the CMB disfavor models involving enough
additional radiation to substantially increase $H_{0}$ Bernal:2016gxb ;
Blinov:2020hmc ; RoyChoudhury:2020dmd ; Brinckmann:2020bcn , leaving an open
question of whether such basic frameworks can help explain the data.

Motivated by the low- versus high-$\ell$ determinations of $\Lambda{\rm CDM}$
parameters and the incompatibility of substantial amounts of interacting
radiation with high-$\ell$ polarization data, we examine in this paper an
$\ell$-dependent solution to the Hubble tension. Specifically, we seek to
understand whether cosmological and astrophysical data favor strongly-coupled
radiation with a mass-threshold near the eV-scale, which increases the
relative density in such radiation at late times. Our main result is that
including such “stepped” dark radiation does significantly improve the
combined fit to CMB, BAO, and SH0ES data while not degrading the fit to the
CMB and BAO. We find that when fit to a dataset including CMB and BAO (but no
direct late-universe measurements of $H_{0}$), our model predicts
$H_{0}=[68.1,71.6]\ \text{km}/\text{s}/\text{Mpc}$ at 90% confidence (in good
agreement with Ref. Freedman:2019jwv ) and reduces the $4.5\sigma$ tension
between $\Lambda{\rm CDM}$ and SH0ES to $2.3\sigma$. To further quantify the
ability of the model to resolve the tension, we apply the “$H_{0}$ Olympics”
measures defined in Ref. Schoneberg:2021qvd  and find that it is the only
concrete particle physics model (with a Lagrangian) that passes the
established rubric in all three Olympic criteria (thus qualifying for a “gold
medal”).

In Sec. II, we describe the overall setup and outline the basic physics of a
stepped dark sector, finding that it naturally fits into the simplest
supersymmetric model, the Wess-Zumino model. In Sec. III, we describe the
effects of a mass-threshold on the CMB. In Sec. IV, we perform various fits to
the data, both including and not including late-time determinations of
$H_{0}$. In Sec. V, we generalize this model to one with a larger set of
parameters that allows for a broader consideration of other related theories.
Finally, in Sec. VI we discuss how this model is favored by the data and
conclude. We also provide two appendices that include details of our
implementation of the mass-threshold and a full set of posterior densities and
best-fit cosmological parameters as determined by our MCMC analysis.

##  II A Model of Stepped Dark Radiation

Constructing a model of interacting radiation a priori is not difficult. A
Weyl fermion $\psi$ and complex scalar $\phi$ that interact through a Yukawa
coupling, $\lambda\,\psi^{2}\phi$, provide a compact and economical model of
interacting radiation. Scalars are naturally massive with masses near the
cutoff of the effective theory, and thus integrating out $\phi$ yields a four-
fermion operator, $\lambda^{2}\psi^{4}/m_{\phi}^{2}$. Requiring these
interactions to be faster than a Hubble time at recombination bounds the mass
of the scalar from above, $m_{\phi}/\lambda\lesssim\,{\rm MeV}$ (assuming that
the dark radiation temperature $T_{d}$ is comparable to the neutrino
temperature). The presence of a light scalar (which is not a Goldstone boson),
while efficient practically, is conceptually incomplete and calls for
additional physics to stabilize its mass.

Remarkably, the above picture follows from the simplest possible
supersymmetric model, the Wess-Zumino model Wess:1974tw , which has the
benefit of automatically controlling the scalar mass. The fermion and scalar
components of a single superfield $X$ with a superpotential coupling
$W=(\lambda/3)\,X^{3}$ yields the Lagrangian

${\cal
L}_{\text{WZ}}=\lambda\,\phi\,\psi^{2}+\lambda^{2}\,(\phi^{*}\phi)^{2}~{}.$
(1)

We refer to this specific model as “Wess-Zumino Dark Radiation” or WZDR.

Since supersymmetry cannot be exact in Nature, we expect an additional scalar
mass $m_{\phi}^{2}\,\phi^{*}\phi$, which may be generated from interactions
with, e.g., the Standard Model. Even without couplings to the Standard Model,
gravity will naturally generate scalar masses of order
$m_{\phi}\sim(M_{\text{SUSY}})^{2}/M_{\text{Pl}}$, where $M_{\text{SUSY}}$ is
the fundamental supersymmetry breaking scale and $M_{\text{Pl}}$ is the Planck
mass. For example, in low-scale gauge mediation Dine:1993yw ; Dine:1995ag
$M_{\text{SUSY}}\sim 100\ \text{TeV}$ and intriguingly the scalar mass is
comparable to the temperature of the photon-baryon plasma at recombination,
$m_{\phi}\sim 1\ \,{\rm eV}$.222If $X$ is sequestered from SUSY breaking,
anomaly mediation Randall:1998uk ; Giudice:1998xp  generates a scalar mass
$m_{\phi}\sim\lambda^{2}/(16\pi^{2})\,M_{\text{SUSY}}^{2}/M_{\text{Pl}}$.
Current and future measurements of the CMB are and will be sensitive to new
relativistic species during the epoch when the temperature of the baryon-
photon plasma was $T_{\gamma}\sim(0.5-50)\ \text{eV}$. This happy coincidence,
that supersymmetric models of interacting radiation often predict a mass-
threshold that occurs within a detectable epoch, makes them distinguishable
from conventional models of strongly-coupled radiation.

Figure 1: A representative example of the redshift dependence of the effective
number of additional neutrinos $N$ (dashed black, left-axis), the equation of
state $w$ (solid red, right-axis), and the speed of sound $c_{s}^{2}$ (dotted
red, right-axis) for a strongly-coupled WZDR fluid.

This then presents a physical picture as follows: some process (examples of
which we will discuss below) produces an amount of interacting radiation
consisting of $\psi$ and $\phi$ particles with an energy density equivalent to
$N_{\text{UV}}$ additional neutrino species after Big Bang nucleosynthesis
(BBN). When the temperature of this sector drops below $m_{\phi}$, the scalars
deposit their entropy into the lighter $\psi$ species as decays and
annihilations between $\phi$ and $\psi$ maintain chemical equilibrium. Due to
conservation of comoving entropy, the energy density of the remaining
fermions, as quantified by the effective number of additional neutrino species
$N_{\text{IR}}$, is greater than $N_{\text{UV}}$

$\frac{N_{\text{IR}}}{N_{\text{UV}}}=\left(\frac{g_{*}^{\phi}+g_{*}^{\psi}}{g_{*}^{\psi}}\right)^{1/3}=\left(\frac{15}{7}\right)^{1/3}\simeq
1.29~{},$ (2)

where $g_{*}^{\phi,\psi}$ is the effective number of relativistic degrees of
freedom in $\phi$, $\psi$, respectively. We refer to this general picture in
which the radiation density increases as a “stepped” model, because the
effective amount of interacting radiation relative to neutrinos at the same
epoch, $N(z)$, increases from $N_{\text{UV}}$ at large redshift $z$ to
$N_{\text{IR}}$ at low $z$ as the fluid transitions through the mass-
threshold. This model is parametrized by two quantities: $N_{\text{IR}}$, the
amount of interacting radiation at late times, and $z_{t}$, the redshift at
which the transition occurs (see Appendix A for the precise definition of
$z_{t}$). The amount of interacting radiation at early times $N_{\text{UV}}$
is determined by Eq. (2).

This transition from $N_{\text{UV}}\to N_{\text{IR}}$ is not instantaneous,
but rather occurs over a period of time. As the temperature $T_{d}$ drops
below the mass of $\phi$, its energy density becomes exponentially suppressed.
Entropy conservation dictates how $T_{d}$ and the energy density of the dark
sector evolve as a function of the scale factor. We note that simply mapping
between two relativistic fluids before and after the transition is not enough
to understand the detailed evolution of its perturbations. During the
transition, the fluid is a combination of both massive and massless particles;
consequently, its equation of state $w$ lies somewhere between that of a
massive ($w=0$) and massless ($w=1/3$) particle. Likewise, the deviation of
the sound speed $c_{s}$ away from that of a relativistic fluid value
($c_{s}^{2}=1/3$) during the transition leaves a non-negligible imprint at the
level of cosmological perturbations. Since the transition occurs while $\phi$
and $\psi$ are chemically coupled, we use entropy conservation to numerically
determine the dark sector temperature as a function of redshift, $T_{d}(z)$,
which can then be used as an input to calculate the redshift evolution of $N$,
$w$, and $c_{s}^{2}$, as shown in Fig. 1. Further details are provided in
Appendix A.

##  III The effects of a Step on the CMB

Before delving into a global analysis of the impact of WZDR on the data, let
us briefly summarize the imprints of a stepped fluid on the CMB (for a
discussion of stepped fluids with earlier transitions and their effects on
BBN, see Ref. Berlin:2019pbq ). Perhaps the simplest comparison to make is
between a stepped and unstepped interacting fluid in which late-universe
parameters are matched. This is illustrated in Fig. 2, which shows the
relative change in the CMB angular power spectrum $\Delta C_{\ell}/C_{\ell}$
(top-panel) and the change in the CMB peak position $\Delta\ell_{\text{peak}}$
(bottom-panel) between a model in which the energy density of interacting
fluid transitions as $N_{\text{UV}}\to N_{\text{IR}}$ at a redshift $z_{t}$,
and a reference model in which the energy density of interacting fluid is
always $N_{\text{IR}}$. In this comparison, we match all other cosmological
input parameters, such as $\Omega_{\Lambda}$, $\Omega_{m}$, and $Y_{p}$. We
refer to the reference model as “self-interacting dark radiation” or SIDR. In
the top-panel of Fig. 2, the most striking feature in $\Delta
C_{\ell}/C_{\ell}$ is a qualitative transition from rapid oscillations at high
$\ell$ to almost no oscillations at low $\ell$. We also note that this
transition occurs at higher $\ell$ for models with earlier steps (larger
$z_{t}$). These rapid oscillations are indicative of a shift in the locations
of the high-$\ell$ CMB peaks. This is made explicit in the bottom-panel of
Fig. 2 which shows the shift of the CMB peaks in the model with a step
relative to the peak locations in the SIDR reference model. Note that there is
no sizeable shift in the position of the low-$\ell$ multipoles, whereas the
relative shift to the position of the high-$\ell$ modes increases
approximately linearly with $\ell$.

Figure 2: (Top) The relative change to the EE spectrum $\Delta
C_{\ell}^{EE}/C_{\ell}^{EE}$ between the WZDR model and a reference SIDR model
for a variety of step locations $z_{t}$. (Bottom) The shift to the position of
the peaks $\Delta\ell_{\text{peak}}$ in the EE power spectrum as a function of
$\ell$ for the same choices of $z_{t}$. In each case, cosmological energy
densities, such as $\Omega_{\Lambda}$, $\Omega_{m}$, and $Y_{p}$, as well as
$N_{\text{IR}}=0.6$, are held fixed between the two models. The effects are
qualitatively similar for the TT spectrum.

The qualitatively different behavior of high vs. low $\ell$ is due to the
different evolution of modes which cross the horizon before and after the
step. As we will show below, the locations of the acoustic peaks are sensitive
to the expansion rate at the time when the corresponding modes enter the
horizon. Low-$\ell$ modes which enter after the step experience the same
Hubble rate as the corresponding modes in the reference SIDR model. In
contrast, high-$\ell$ modes which enter before the step experience a smaller
expansion rate at horizon crossing, which shifts the corresponding peaks
towards smaller $\ell$.

It is straightforward to analytically derive the shift of the peak locations
for transitions which happen well before matter-radiation equality. Readers
less interested in this derivation are encouraged to skip the remainder of
this section. The main results are Eqs. (6) and (7) for the phase shift.
Assuming transitions well before equality, we can ignore the contributions of
the matter density to both the Hubble parameter and the gravitational
potentials during horizon crossing for all modes of interest. For modes that
enter either well before or well after the transition, $N(z)$ is a constant at
horizon crossing, in which case one can show that the photon number
overdensity per coordinate volume $d_{\gamma}$ is constant before horizon
crossing (in conformal Newtonian gauge) Bashinsky:2003tk ; Baumann:2015rya .
Moreover, the Einstein equation involving the photon-baryon velocity,
supplemented with the continuity equation (see, e.g., Eqs. (23b) and (30) of
Ref. Ma:1995ey ), yields the initial condition for the conformal-time
derivative of $d_{\gamma}$. The superhorizon initial conditions for the photon
perturbation are therefore

$d_{\gamma}=-3\zeta~{}~{},~{}~{}\dot{d}_{\gamma}=-\frac{1}{2}\,k^{2}\,\mathcal{H}^{-1}\,\Phi~{},$
(3)

where $\zeta$ is the primordial curvature perturbation, $\mathcal{H}$ is the
conformal Hubble parameter, and $\Phi$ is the gravitational potential in the
notation of Ref. Baumann:2015rya . For low-$k$ modes which enter the horizon
well after the step, the initial conditions in Eq. (3) are therefore identical
to the initial conditions in the reference SIDR model. Furthermore, since the
WZDR and SIDR models employ identical cosmological parameters, the post-step
evolution of the perturbations is also identical. Therefore, there is no shift
to the peak locations for such low-$\ell$ modes when comparing these two
models.

On the other hand, large-$k$ modes enter the horizon well before the step,
i.e., at a time when the amount of fluid is less than in the reference model.
Therefore, these modes experience a slower expansion rate before and during
horizon crossing than for the same modes in the SIDR model. After horizon
crossing, the gravitational potential decays quickly, effectively suppressing
any direct impact that the stepped fluid has on the evolution of the photon
perturbations. However, the slower expansion rate encountered by such photon
perturbations means that they undergo a larger phase evolution until the step
is completed. The resulting phase shift is proportional to the wave number $k$
as well as the difference in conformal time $\Delta\tau$, as defined below in
Eq. (6).

To derive $\Delta\tau$ we begin by observing that the perturbation equations
are usually written as differential equations with derivatives with respect to
conformal time $\tau$ and coefficients which depend on $\mathcal{H}$. To solve
the equations one usually begins by substituting $\mathcal{H}=\tau^{-1}$ and
integrating. This relation follows from integrating the Friedmann equations,
which during radiation domination and for constant $N(z)$ yield

$\displaystyle\tau(z)=\mathcal{H}^{-1}(z)+C~{},$ (4)

where $C$ is an integration constant. This integration constant corresponds to
the origin of the conformal time and is arbitrary, a reflection of the time
translation invariance of the perturbation and Friedmann equations. Making the
customary choice $C=0$ then implies $\lim_{z\to\infty}\tau=0$ and
$\mathcal{H}=\tau^{-1}$. We will follow this convention for the reference SIDR
model. In order to make the comparison simple, we choose $C$ for the WZDR
model such that at late times both models have the same conformal time $\tau$
at fixed $z$. Using this boundary condition, we find

$\displaystyle\mathcal{H}^{-1}(z)$ $\displaystyle=\begin{cases}\tau(z)&z\ll
z_{t}\\\ \tau(z)+\Delta\tau&z\gg z_{t}~{},\end{cases}$ (5)

in the WZDR model, such that $\Delta\tau$ is the calculable constant

$\displaystyle\Delta\tau$
$\displaystyle=\int_{z_{*}}^{\infty}dz\left(\frac{1}{H_{\rm
WZDR}(z)}-\frac{1}{H_{\rm SIDR}(z)}\right)~{},$ (6)

where $H_{\text{WZDR, SIDR}}$ is the Hubble parameter in the WZDR or SIDR
model, respectively. Note that $H_{\text{WZDR}}=H_{\text{SIDR}}$ after the
step, so that the integrand vanishes for redshifts between recombination at
$z_{*}$ and the step. Matching onto the superhorizon solutions of Eq. (3), the
well-known solution for the photon perturbation inside the horizon takes the
form

$\displaystyle d_{\gamma}(z)$
$\displaystyle\propto\begin{cases}\cos[c_{\gamma}\,k\,\tau(z)]&(k\ll k_{t})\\\
\cos[c_{\gamma}\,k\,(\tau(z)+\Delta\tau)]&(k\gg k_{t})~{},\end{cases}$ (7)

where $c_{\gamma}\simeq 1/\sqrt{3}$ is the speed of sound of the photon fluid
and $k_{t}$ is defined as the wave number of a mode that enters the horizon
during the transition. This solution neglects the drag effect contribution
from the gravitational potentials Bashinsky:2003tk ; Baumann:2015rya , which
modifies both the phase and amplitude of $d_{\gamma}$ and gives rise to
subleading effects. From Eq. (7), we see that the photon perturbations for
small-$k$ modes ($k\ll k_{t}$) are now manifestly identical to those of the
reference SIDR model, whereas the high-$k$ modes have a phase shift at
recombination approximately given by $c_{\gamma}k\Delta\tau$. We note that
this phase shift can be removed by adjusting $H_{0}$ to keep the angular size
of the sound horizon $\theta_{s}$ fixed between the two models.333We have
verified numerically in CLASS that the change to $\theta_{s}$ from the step
satisfies $\Delta\theta_{s}/\theta_{s}\simeq\Delta\tau/\tau^{*}$. However,
this amounts to introducing a corresponding phase shift at low-$\ell$. Indeed,
the $\ell$-dependent signature of a stepped fluid is distinct from a
modification to the angular scale of the sound horizon at recombination.

The picture described above holds qualitatively even when matter effects are
non-negligible at the time of the transition. For small-$k$ modes, which enter
the horizon well after the step, the conformal time difference between horizon
entry and recombination is the same as in the SIDR model, and thus the
small-$k$ modes of both the WZDR and the SIDR models possess the same phase
evolution. On the contrary, large-$k$ modes that enter the horizon well before
the step experience a different expansion rate, and thus have more time to
evolve between horizon entry until the time of recombination.

In the next section, we directly confront the data and find that the WZDR
model is preferred compared to $\Lambda{\rm CDM}$ and previously studied
models of SIDR. In fact, a fit to the CMB and BAO prefers a step location near
$z_{t}\simeq 2\times 10^{4}$, the redshift at which multipoles of $\ell\sim
10^{3}$ reenter the horizon.

Model | $\Delta\chi^{2}$ | $N_{\text{eff,IR}}$ |  $H_{0}$ $(\text{km}/\text{s}/\text{Mpc})$
---|---|---|---
$\Lambda$CDM | 0.0 | $3.04$ | $67.6~{}\,[67.0,68.3]$
$\text{$\Lambda$CDM}+N_{\rm eff}$ | 0.0 | $3.04~{}\,[3.04,3.30]$ | $67.7~{}\,[67.2,69.4]$
SIDR | $-0.1$ | $3.07~{}\,[3.12,3.55]$ | $67.9~{}\,[67.4,70.4]$
WZDR | $-1.2$ | $3.27~{}\,[3.14,3.69]$ | $69.1~{}\,[68.1,71.6]$
Table 1: A summary of a fit to the SH0ES-independent dataset $\cal{D}$ for
various models listed in the first column. The second column gives the
$\Delta\chi^{2}$ of the best-fit compared to $\Lambda{\rm CDM}$. The third and
fourth columns display the preferred values of the effective number of
neutrinos in free-streaming and strongly-coupled radiation at late times
$N_{\text{eff,IR}}$ (i.e., including the Standard Model neutrino contribution)
and the Hubble parameter $H_{0}$, respectively. The first entry of each cell
for $N_{\text{eff,IR}}$ and $H_{0}$ gives the best-fit value. The brackets
that follow denote the $90\%$ C.L. posterior range, which is defined to be the
narrowest interval containing $90\%$ of the integrated posterior density.
These are computed directly from the posterior densities and not using
Gaussian fits to the posteriors.

##  IV Cosmological Imprints of WZDR

To understand the implications of this scenario, we perform a fit of various
models to a range of cosmological data. We use a modified version of standard
CLASS v2.9 CLASS  (implementing the fluid as described in Appendix A) combined
with a MCMC sampler (MontePython v3.4 Brinckmann:MP ) to study the
cosmological constraints on models of interacting radiation. For the WZDR
model, we adopt flat priors on the cosmological parameters
$\\{\omega_{b},\omega_{\mathrm{cdm}},\theta_{s},n_{s},A_{s},\tau_{\mathrm{reio}},N_{\text{IR}},\log_{10}z_{t}\\}$
with the following ranges $\log_{10}z_{t}\ \in[4.0,4.6]$ and
$N_{\text{IR}}\geq 0$.444We restrict the range of $z_{t}$ in this way in order
to avoid very late and very early transitions, as these points correspond to
conventional SIDR models. We assume a Standard Model neutrino sector that
undergoes a standard cosmological history, consisting of one massive neutrino
species with a mass of $0.06\ \text{eV}$ and two massless neutrino species.

As in Ref. Schoneberg:2021qvd , we consider the following datasets in our
default analysis. The dataset $\cal{D}$ includes the Planck 2018 dataset
Planck:2018vyg , including TT,TE, and EE in both low ($\ell<30$) and high
($\ell\geq 30$) multipole ranges, as well as CMB lensing Planck:2018vyg  and
the full set of nuisance parameters. It also includes late-universe
constraints, such as BAO measurements from BOSS DR12 ($z=0.38,0.51,0.61$) BOSS
and small $z$ measurements from 6dF ($z=0.106$) 6dF  and MGS ($z=0.15$) MGS
catalogs, as well as PANTHEON supernovae data Scolnic:2017caz . In addition to
our baseline dataset $\cal D$, we also consider a dataset $\cal D+$ that
additionally includes measurements from SH0ES Riess:2020fzl  via a prior on
the intrinsic magnitude of supernovae $M_{b}$. We summarize our results for a
range of models fitting to $\cal D$ and $\cal D+$ in Appendix B.

When we consider the impact of WZDR on cosmological observables, the immediate
question is: relative to what? Because the presence of new interacting
radiation changes integrated quantities (such as the sound horizon at
recombination), the change in any one parameter naturally changes others,
making direct comparisons between different models complicated. Two natural
points of comparison are to models with additional strongly-coupled or free-
streaming radiation without a step. The former is the conventional SIDR model
discussed above, and the latter is a minimal extension of $\Lambda{\rm CDM}$,
in which the effective number of free-streaming neutrino species is allowed to
take values greater than the $\Lambda{\rm CDM}$ contribution from the Standard
Model neutrino sector, i.e., $N_{\rm eff}>3.044$ (denoted as
$\Lambda\text{CDM}+N_{\rm eff}$) Froustey:2020mcq ; Bennett:2019ewm . Using
Table 1, we proceed to make these comparisons, which shows results of a fit to
the SH0ES-independent dataset $\cal{D}$. For convenience, we define the
variable $N_{\text{eff,IR}}$ as the late-time value of the effective number of
neutrinos in free-streaming and strongly-coupled radiation (including the
Standard Model neutrino contribution).

Model | $\Delta\chi^{2}$ | $N_{\text{eff,IR}}$ |  $H_{0}$ $(\text{km}/\text{s}/\text{Mpc})$
---|---|---|---
$\Lambda$CDM | 0.0 | $3.04$ | $68.2~{}\,[67.5,68.9]$
$\text{$\Lambda$CDM}+N_{\rm eff}$ | $-5.7$ | $3.37~{}\,[3.20,3.63]$ | $70.0~{}\,[68.9,71.6]$
SIDR | $-10.6$ | $3.51~{}\,[3.26,3.73]$ | $71.0~{}\,[69.4,72.5]$
WZDR | $-15.1$ | $3.63~{}\,[3.43,3.94]$ | $71.4~{}\,[70.2,73.2]$
Table 2: As in Table 1, but after fitting to the dataset $\cal{D}+$ that
includes SH0ES.

Figure 3: Comparison of the marginalized 1D and 2D posterior distributions for
the Hubble parameter $H_{0}$ and the late-time value of the effective number
of neutrinos in radiation $N_{\text{eff,IR}}$ (including the Standard Model
neutrino contribution) for the $\Lambda\text{CDM}+N_{\rm eff}$, SIDR, and WZDR
models when fitting to the dataset $\cal{D}$ (not including SH0ES) in the left
set of panels or $\cal{D}+$ (including SH0ES) in the right set of panels.

We note at the outset that the best-fit points for the dataset $\cal{D}$ in
conventional models involving additional free-streaming or strongly-
interacting radiation are either identical or very close to $\Lambda{\rm CDM}$
and do not improve the $\chi^{2}$ by more than $\mathcal{O}{(0.1)}$. This
changes with the addition of the WZDR-step. Considering only data from
$\cal{D}$, we see a slight improvement going from SIDR to WZDR, reducing the
best-fit $\chi^{2}$ by $1.2$, about as would be expected with the addition of
the one new parameter $z_{t}$, which has a posterior mean-value of
$z_{t}=1.9_{-0.6}^{+0.8}\times 10^{4}$. Thus, although the quantitative pull
is small, there is a nonzero preference in the data, independent of SH0ES, for
a step in an additional interacting fluid component.555Many previous analyses
have assumed that the new species of radiation was also present during BBN,
increasing the predicted abundance of primordial helium $Y_{p}$. Since that is
a much earlier era, we instead assume that the radiation is populated well
after BBN, corresponding to temperatures below $\sim 100\ \text{keV}$
Berlin:2019pbq . Thus, we will make comparisons between SIDR and WZDR with
this same assumption, so as not to penalize models of SIDR compared to WZDR.
In Table 4 of Appendix B, we also show results for SIDR models in which the
energy density in interacting radiation is also present during BBN.

Remarkably, the preferred value for $H_{0}$ in a fit to $\cal{D}$ is shifted
to larger values in the WZDR model, with a best-fit of $H_{0}=69.1\
\text{km}/\text{s}/\text{Mpc}$, as compared to $H_{0}=67.9\
\text{km}/\text{s}/\text{Mpc}$ for SIDR, $H_{0}=67.7\
\text{km}/\text{s}/\text{Mpc}$ for $\text{$\Lambda$CDM}+N_{\rm eff}$, and
$H_{0}=67.6\ \text{km}/\text{s}/\text{Mpc}$ for $\Lambda{\rm CDM}$. A similar
trend is evident in the 90% C.L. ranges listed in Table 1. This immediately
raises the question: can WZDR help address the existing tension between the
$\Lambda{\rm CDM}$-inferred value of $H_{0}$ and the late-universe measurement
of $H_{0}=(73.2\pm 1.3)\ \text{km}/\text{s}/\text{Mpc}$ by the SH0ES
collaboration?

The short answer to this question is: yes, the presence of the WZDR-step does
allow for a significant reduction to this tension. Earlier analyses have shown
that additional free-streaming radiation can naturally allow for a larger
value of $H_{0}$ when fitting to the dataset $\cal{D}+$ that includes SH0ES,
but only at the cost of significantly worsening the fit to the SH0ES-
independent dataset $\cal{D}$. Making this radiation interacting (as in the
SIDR model) somewhat ameliorates the issue, but this is still constrained by
the high-$\ell$ multipoles of the CMB polarization power spectrum
Bernal:2016gxb ; Blinov:2020hmc ; Brinckmann:2020bcn . In the WZDR model, this
is compensated by the the $\ell$-dependent modifications to the CMB, allowing
for additional levels of interacting radiation.

We provide additional results of our dedicated MCMC analysis in Fig. 3, which
shows the posteriors for $H_{0}$ and the late-time value of the effective
number of neutrinos in free-streaming and strongly-coupled radiation for the
$\text{$\Lambda$CDM}+N_{\rm eff}$, SIDR, and WZDR models (the full set of
posteriors for each of these scenarios is provided in Appendix B). As for
conventional early-universe solutions to the Hubble tension, additional
radiation is correlated with larger values of $H_{0}$, corresponding to an
approximately fixed angular size of the sound horizon at recombination
Hou:2011ec . Most notable in Fig. 3, however, is the fact that WZDR predicts
$H_{0}$ and $N_{\text{eff,IR}}$ posteriors that extend out to much larger
values, in the case that SH0ES is either included (right panel) or not
included (left panel) in the analysis. We show the resulting best-fit values
and posterior ranges for the full $\cal{D}+$ dataset in Table 2.

But, of course, simply predicting a larger value of $H_{0}$ is not a solution
to the tension if it simply provides an overall bad fit to the data. Recently,
Ref. Schoneberg:2021qvd  established a rubric for comparing models that could
address this tension with three basic measures: GT (Gaussian tension),
$Q_{\text{DMAP}}$ (difference of the maximum a posteriori), and
$\Delta\text{AIC}$ (Akaike information criterium). The values of these
measures are provided in Table 3 for $\Lambda{\rm CDM}$ and the three
benchmark radiation models. We now briefly summarize each of these in turn
(see Refs. Schoneberg:2021qvd ; Raveri:2018wln  for further details).

Model | GT | $Q_{\text{DMAP}}$ | $\Delta\text{AIC}$
---|---|---|---
$\Lambda$CDM | $4.5\sigma$ | $4.5\sigma$ | $0.0$
$\text{$\Lambda$CDM}+N_{\rm eff}$ | $3.7\sigma$ | $3.8\sigma$ | $-3.7$
SIDR | $3.1\sigma$ | $3.1\sigma$ | $-8.6$
WZDR | $2.3\sigma$ | $2.4\sigma$ | $-11.1$
Table 3: The GT (Gaussian tension), $Q_{\text{DMAP}}$ (difference of the
maximum a posteriori), and $\Delta\text{AIC}$ (Akaike information criterium)
measures, which together serve to address the success in resolving the Hubble
tension for various models.

GT and $Q_{\text{DMAP}}$ both address the residual level of tension between
the results of SH0ES (expressed as the measured supernova magnitudes $M_{b}$)
and a model-predicted value obtained by fitting to the SH0ES-independent
dataset $\cal{D}$. In particular, GT is defined as the difference in the
$M_{b}$ mean-value as determined from either method, in units of $\sigma$ (the
standard deviation obtained from $\cal{D}$ and SH0ES added in quadrature).
$Q_{\text{DMAP}}$, on the other hand, is defined as the square root of the
difference between the best-fit $\chi^{2}$ obtained by either fitting a model
to $\cal{D}+$ or $\cal{D}$, in units of $\sigma$. Smaller values of either
measure signify a reduced tension in the data, and for Gaussian distributed
posteriors, they are equivalent. We find that while
$\text{$\Lambda$CDM}+N_{\rm eff}$ and SIDR reduce the $4.5\sigma$ tension
between $\Lambda{\rm CDM}$ and SH0ES to the $\sim 4\sigma$, $3\sigma$ level,
respectively, WZDR significantly reduces it further still, down to $\sim
2\sigma$.

$\Delta\text{AIC}$ measures the success in fitting to the full dataset
including SH0ES, defined as the difference in the best-fit $\chi^{2}$ for
$\cal{D}+$ between a given model and $\Lambda{\rm CDM}$ with a $\chi^{2}$
penalty of $+2$ for each additional (beyond $\Lambda{\rm CDM}$) model
parameter:

$\Delta\text{AIC}=\chi^{2}-\chi_{\Lambda\rm CDM}^{2}+2\times(\rm new\
parameters)\,.$ (8)

Hence, more negative values signify increased agreement with the full dataset.
As shown in Table 3, this measure further illustrates the dramatic improvement
in a fit involving WZDR. In particular, while $\Lambda\text{CDM}+N_{\rm eff}$
and SIDR reduce the $\Delta\text{AIC}$ down to $-3.7$ and $-8.6$,
respectively, this is improved further to $-11.1$ in the WZDR model
(corresponding to a best-fit $\chi^{2}$ which is lower by $-15.1$ compared to
$\Lambda{\rm CDM}$).

With the addition of a single parameter, WZDR provides an improvement of
$\Delta\chi_{\cal{D}+}^{2}=-4.5$ when fitting to $\cal{D}+$, as compared to
SIDR (see Appendix B). Also noteworthy is the fact that when fitting to the
dataset $\cal{D}+$ including SH0ES, the fit to the SH0ES-independent dataset
$\cal{D}$ for WZDR is comparable to that of $\Lambda{\rm CDM}$. In contrast,
both $\Lambda{\rm CDM}$+$N_{\rm eff}$ and SIDR provide markedly worse fits to
$\cal{D}$ when fitting to the full $\cal{D}+$ (see Table 8).

To summarize, the WZDR model exemplifies a significant improvement in each of
the three standardized measures of Ref. Schoneberg:2021qvd , as shown in Table
3. Thus, this model provides a potential scenario to explain this discrepancy.
Moreover, as future data at high-$\ell$ arrive, such as from the Simons
Observatory or CMB-S4, a direct comparison to the CMB, even without SH0ES
data, should further clarify the degree of this success.

##  V Generalizations and Stepped Dark Sectors

Although the minimal supersymmetric model for interacting radiation (WZDR) has
two free parameters ($N_{\text{IR}}$ and $z_{t}$), one can envision
generalizing this to a scenario where the parameter
$N_{\text{IR}}/N_{\text{UV}}$ is not a fixed number, as in Eq. (2), for
instance if $\psi$ is a fundamental and $\phi$ an adjoint under a global
symmetry rather than residing in a single chiral superfield. This then makes
$N_{\text{IR}}/N_{\text{UV}}$ an adjustable parameter of the model.
Additionally, if the dark radiation is populated via late-time freeze-in from
neutrinos, the energy density in free-streaming radiation $N_{\text{fs}}$ can
be modified Berlin:2017ftj ; Berlin:2018ztp . We thus take $N_{\text{fs}}$ as
a free parameter, as well.

Thus from a physical point of view, in such models the four free model
parameters are the effective number of neutrinos in interacting species
$N_{\text{UV}}$ at early times, the number of interacting species
$N_{\text{IR}}$ at late times, the redshift at which the dark sector drops
below a mass-threshold $z_{t}$, and the amount of free-streaming radiation
$N_{\text{fs}}$. Such a parameterization is a simple extension of interacting
fluid models, motivated by the simple presence of a mass-threshold in the
sector, which is essentially required in models involving scalars. We refer to
this generalization of a stepped fluid as “General StepDR.”

We perform a fit of the General StepDR model to the datasets $\cal{D}$ and
$\cal{D}+$. The parameters for this model are:
$\\{\omega_{b},\omega_{\mathrm{cdm}},\theta_{s},n_{s},A_{s},\tau_{\mathrm{reio}},N_{\text{UV}},N_{\text{IR}},N_{\text{fs}},\log_{10}z_{t}\\}$
with the ranges $\log_{10}z_{t}\ \in[4.0,4.6]$, $N_{\text{UV}}\geq 0$, and
$N_{\text{IR}}\geq N_{\text{UV}}$. Additionally, we assume one massive
neutrino species with a mass of $0.06\ \text{eV}$. The results of this
analysis are included in Appendix B. From the full set of posterior
distributions in Appendix B, we see that the introduction of these additional
parameters widens the posterior density for $N_{\text{UV}}$ and
$N_{\text{IR}}$, while the favored region of parameter space for the General
StepDR model is consistent with the favored region in the restricted WZDR
model. We also see from Appendix B that the best-fit $\chi^{2}$ values are
reduced by $\mathcal{O}{(1)}$ when comparing the General StepDR and WZDR
models, as expected from the introduction of two additional free parameters.

##  VI Discussion

Models of interacting radiation provide a plausible and testable extension
beyond $\Lambda{\rm CDM}$. In many such scenarios there is a natural mass-
threshold present, corresponding to the scale of the force-carrier. In such
models, there is a “step” as the massive mediator entropy is transferred to
the remaining degrees of freedom. The simplest such model - a single fermion
and a single scalar - is nicely packaged into a well-known supersymmetric
theory, namely the Wess-Zumino model. Such a framework naturally allows for a
scalar mass which is light enough to be relevant during the era in which the
CMB-detectable modes are inside the horizon.

We emphasize that this framework of stepped dark radiation is a theoretically
simple one. The Standard Model has copious mass-thresholds, including the
critical late threshold of the electron mass where photons are heated in
analogous fashion. Our purpose here is to understand the implications of such
a step in the context of present-day cosmological data. There is a well-known
tension between CMB data and late-universe measurements of the Hubble
constant, and so it is important to distinguish them.

Remarkably, WZDR provides a good fit to the CMB data, alone. The best-fit
point, not including SH0ES data, prefers additional radiation
$N_{\text{IR}}=0.23$, excluding $N_{\text{IR}}=0$ at better than 90%
confidence. This contrasts conventional interacting radiation which has a
best-fit of $N_{\text{IR}}=0.03$, and free-streaming radiation which has a
best-fit of $\Delta N_{\rm eff}=0$. At the same time, the WZDR model prefers a
higher value of $H_{0}=69.1\ \text{km}/\text{s}/\text{Mpc}$ with a wide range
of allowed values at 90% confidence, $[68.1,71.6]\
\text{km}/\text{s}/\text{Mpc}$. Including SH0ES in the analysis changes the
values somewhat, but qualitatively, the picture does not shift, preferring a
best-fit value $H_{0}=71.4\ \text{km}/\text{s}/\text{Mpc}$.

Comparing WZDR to SIDR with the same late-time parameters, the change to the
CMB power spectrum is an $\ell$-dependent phase shift. Modes entering the
horizon earlier, with relatively less radiation present, experience a larger
sound horizon, thus shifting the peaks to lower $\ell$ multipoles. While our
analysis does not analyze $\chi^{2}$ values for separate $\ell$ regions (see,
e.g., Ref. Hill:2021yec , $\ell<650$), the WZDR and General StepDR models give
a good fit to all data, including high-$\ell$ data.

We should also comment on another discordance in cosmological data, the
$S_{8}$ tension, which we have so far ignored. The $\Lambda{\rm
CDM}$-preferred value of $S_{8}$, as inferred from fitting to the CMB, is
about $2\sigma$ higher than the values inferred from weak lensing surveys
Hildebrandt:2016iqg ; Hildebrandt:2018yau ; DES:2017myr ; HSC:2018mrq . This
tension motivates new physics that suppresses the matter power spectrum. We
have not included such large scale structure (LSS) data in our fits, in part
because the $S_{8}$ tension is not as severe or robust as the $H_{0}$ tension,
LSS data was not included in the fits of Ref. Schoneberg:2021qvd , and the
dark radiation models that we have investigated do not significantly alter the
prediction for $S_{8}$ compared to $\Lambda{\rm CDM}$ (see Appendix B), and so
neither improve nor worsen the fit to LSS data.

However, we remark that in the context of an interacting sector of dark
radiation, it is very natural to expect a component of the dark matter to
couple to dark radiation. Such a coupling endows the dark matter with pressure
that slows the growth of structure and therefore predicts a lower value for
$S_{8}$ Buen-Abad:2015ova ; Lesgourgues:2015wza ; Buen-Abad:2017gxg ;
Bansal:2021dfh . This pressure depends on a new parameter: the coupling
between dark matter and dark radiation. Thus, it would be straightforward to
resolve the $S_{8}$ tension in a natural extension of the stepped-fluid models
considered here. Such a model could arise by including an additional field $Y$
with a $Z_{2}$ symmetry and superpotenial term of the form $XY^{2}$. This
would naturally yield an interacting dark matter subcomponent.

In our analysis, we have assumed the standard BBN-predicted value of
primordial helium, assuming the dark radiation is not present at that time.
It’s quite easy for the radiation to appear at late times, for instance
through the decay of a relic particle, for which the allowed value of
$N_{\text{UV}}$ can be anything. If the dark radiation was in thermal contact
at early times such that $N_{\text{UV}}\sim\mathcal{O}{(0.1)}$ at BBN, it can
be boosted by a factor of a few if the dark sector has a strongly coupled
gauge group which confines with a mass gap, such as supersymmetric QCD. One
might also study this scenario assuming the radiation is present at BBN,
naturally raising the level of primordial helium, and likely worsening the fit
to a small degree.

It is also worth comparing WZDR and General StepDR to other models that invoke
qualitatively similar dynamics. For instance, Ref. Bansal:2021dfh  recently
discussed a mirror copy of the Standard Model within the context of the Hubble
tension. Although such theories involve some of the same ingredients as
discussed in our work, they also involve additional dynamics arising from
transitions from strongly-coupled to free-streaming radiation and epochs in
which the dark sector undergoes its own version of nucleosynthesis and
recombination. However, we are unable to meaningfully compare to the fits of
Ref. Bansal:2021dfh , since these were performed only with data that include
SH0ES and LSS. We also note that of all the models catalogued in Ref.
Schoneberg:2021qvd , the one involving a neutrino-coupled majoron
Escudero:2019gvw ; EscuderoAbenza:2020egd ; Escudero:2021rfi  is the most
similar to WZDR. However, WZDR is much more successful in addressing the
$H_{0}$ tension, as quantified by comparing the results of Table 3 and Ref.
Schoneberg:2021qvd .

Going forward, as data improve, there is a great hope that the CMB data alone
will give clarity on whether a scenario similar to ours provides a reasonable
description of Nature. A recent analysis by Ref. Hill:2021yec  showed that an
Early Dark Energy model could be distinguished with future ACT data at the
level of $25\sigma$ from the current best-fit $\Lambda{\rm CDM}$ point. The
WZDR model has similar sized residuals, and thus we can be hopeful that these
data might distinguish it from $\Lambda{\rm CDM}$ as well. A proper analysis
is warranted, however.

In summary, we have studied a simple, natural, theoretical model extension
beyond $\Lambda{\rm CDM}$, namely interacting radiation with a step, arising
from a straightforward mass-threshold. This minor extension provides a good
fit to the CMB data with a higher value of $H_{0}$. Including SH0ES data
strengthens the evidence for additional radiation, as well as the step,
improving the $\chi^{2}$ by 15.1 compared to $\Lambda{\rm CDM}$ and by 4.5
compared to conventional SIDR, with only a single additional parameter (the
step location). As further data appear, we may learn the Hubble tension is the
first indication of a dynamic and interacting dark sector.

## Acknowledgements

We thank Vivian Poulin, Nils Schöeneberg, and Andrea Pérez Sánchez for helpful
discussions and comparing data points from the “$H_{0}$ Olympics”
Schoneberg:2021qvd . The work of D.A., M.J., and M.S. is supported by the U.S.
Department of Energy (DOE) under Award DE-SC0015845. A.B. is supported by the
James Arthur Fellowship. N.W. is supported by NSF under award PHY-1915409, by
the BSF under grant 2018140, and by the Simons Foundation. Our MCMC runs were
performed on the Shared Computing Cluster, which is administered by Boston
University’s Research Computing Services.

## References

  * (1) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * (2) A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, “Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM,” Astrophys. J. Lett. 908 (2021) no. 1, L6, arXiv:2012.08534 [astro-ph.CO].
  * (3) W. L. Freedman et al., “The Carnegie-Chicago Hubble Program. VIII. An Independent Determination of the Hubble Constant Based on the Tip of the Red Giant Branch,” arXiv:1907.05922 [astro-ph.CO].
  * (4) W. L. Freedman, “Measurements of the Hubble Constant: Tensions in Perspective,” Astrophys. J. 919 (2021) no. 1, 16, arXiv:2106.15656 [astro-ph.CO].
  * (5) W. Yuan, A. G. Riess, L. M. Macri, S. Casertano, and D. Scolnic, “Consistent Calibration of the Tip of the Red Giant Branch in the Large Magellanic Cloud on the Hubble Space Telescope Photometric System and a Re-determination of the Hubble Constant,” Astrophys. J. 886 (2019) 61, arXiv:1908.00993 [astro-ph.GA].
  * (6) J. Soltis, S. Casertano, and A. G. Riess, “The Parallax of $\omega$ Centauri Measured from Gaia EDR3 and a Direct, Geometric Calibration of the Tip of the Red Giant Branch and the Hubble Constant,” Astrophys. J. Lett. 908 (2021) no. 1, L5, arXiv:2012.09196 [astro-ph.GA].
  * (7) K. C. Wong et al., “H0LiCOW – XIII. A 2.4 per cent measurement of H0 from lensed quasars: 5.3$\sigma$ tension between early- and late-Universe probes,” Mon. Not. Roy. Astron. Soc. 498 (2020) no. 1, 1420–1439, arXiv:1907.04869 [astro-ph.CO].
  * (8) D. W. Pesce et al., “The Megamaser Cosmology Project. XIII. Combined Hubble constant constraints,” Astrophys. J. Lett. 891 (2020) no. 1, L1, arXiv:2001.09213 [astro-ph.CO].
  * (9) LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., “A Gravitational-wave Measurement of the Hubble Constant Following the Second Observing Run of Advanced LIGO and Virgo,” Astrophys. J. 909 (2021) no. 2, 218, arXiv:1908.06060 [astro-ph.CO].
  * (10) N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, “The $H_{0}$ Olympics: A fair ranking of proposed models,” arXiv:2107.10291 [astro-ph.CO].
  * (11) G. E. Addison, Y. Huang, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw, and J. L. Weiland, “Quantifying discordance in the 2015 Planck CMB spectrum,” Astrophys. J. 818 (2016) no. 2, 132, arXiv:1511.00055 [astro-ph.CO].
  * (12) J. L. Bernal, L. Verde, and A. G. Riess, “The trouble with $H_{0}$,” JCAP 10 (2016) 019, arXiv:1607.05617 [astro-ph.CO].
  * (13) K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan, and W. L. K. Wu, “Sounds Discordant: Classical Distance Ladder & $\Lambda$CDM -based Determinations of the Cosmological Sound Horizon,” Astrophys. J. 874 (2019) no. 1, 4, arXiv:1811.00537 [astro-ph.CO].
  * (14) C. Brust, Y. Cui, and K. Sigurdson, “Cosmological Constraints on Interacting Light Particles,” JCAP 08 (2017) 020, arXiv:1703.10732 [astro-ph.CO].
  * (15) N. Blinov and G. Marques-Tavares, “Interacting radiation after Planck and its implications for the Hubble Tension,” JCAP 09 (2020) 029, arXiv:2003.08387 [astro-ph.CO].
  * (16) S. Roy Choudhury, S. Hannestad, and T. Tram, “Updated constraints on massive neutrino self-interactions from cosmology in light of the $H_{0}$ tension,” JCAP 03 (2021) 084, arXiv:2012.07519 [astro-ph.CO].
  * (17) T. Brinckmann, J. H. Chang, and M. LoVerde, “Self-interacting neutrinos, the Hubble parameter tension, and the cosmic microwave background,” Phys. Rev. D 104 (2021) no. 6, 063523, arXiv:2012.11830 [astro-ph.CO].
  * (18) C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, “Neutrino puzzle: Anomalies, interactions, and cosmological tensions,” Phys. Rev. D 101 (2020) no. 12, 123505, arXiv:1902.00534 [astro-ph.CO].
  * (19) M. Escudero and S. J. Witte, “A CMB search for the neutrino mass mechanism and its relation to the Hubble tension,” Eur. Phys. J. C 80 (2020) no. 4, 294, arXiv:1909.04044 [astro-ph.CO].
  * (20) M. Escudero Abenza and S. J. Witte, “Could the Hubble Tension be Pointing Towards the Neutrino Mass Mechanism?,” in Prospects in Neutrino Physics.  4, 2020.  arXiv:2004.01470 [hep-ph].
  * (21) M. Escudero and S. J. Witte, “The hubble tension as a hint of leptogenesis and neutrino mass generation,” Eur. Phys. J. C 81 (2021) no. 6, 515, arXiv:2103.03249 [hep-ph].
  * (22) T. Karwal and M. Kamionkowski, “Dark energy at early times, the Hubble parameter, and the string axiverse,” Phys. Rev. D 94 (2016) no. 10, 103523, arXiv:1608.01309 [astro-ph.CO].
  * (23) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, “Early Dark Energy Can Resolve The Hubble Tension,” Phys. Rev. Lett. 122 (2019) no. 22, 221301, arXiv:1811.04083 [astro-ph.CO].
  * (24) M.-X. Lin, G. Benevento, W. Hu, and M. Raveri, “Acoustic Dark Energy: Potential Conversion of the Hubble Tension,” Phys. Rev. D 100 (2019) no. 6, 063542, arXiv:1905.12618 [astro-ph.CO].
  * (25) T. L. Smith, V. Poulin, and M. A. Amin, “Oscillating scalar fields and the Hubble tension: a resolution with novel signatures,” Phys. Rev. D 101 (2020) no. 6, 063523, arXiv:1908.06995 [astro-ph.CO].
  * (26) S. Bansal, J. H. Kim, C. Kolda, M. Low, and Y. Tsai, “Mirror Twin Higgs Cosmology: Constraints and a Possible Resolution to the $H_{0}$ and $S_{8}$ Tensions,” arXiv:2110.04317 [hep-ph].
  * (27) F.-Y. Cyr-Racine, F. Ge, and L. Knox, “A Symmetry of Cosmological Observables, and a High Hubble Constant as an Indicator of a Mirror World Dark Sector,” arXiv:2107.13000 [astro-ph.CO].
  * (28) Z. Hou, R. Keisler, L. Knox, M. Millea, and C. Reichardt, “How Massless Neutrinos Affect the Cosmic Microwave Background Damping Tail,” Phys. Rev. D 87 (2013) 083008, arXiv:1104.2333 [astro-ph.CO].
  * (29) S. Bashinsky and U. Seljak, “Neutrino perturbations in CMB anisotropy and matter clustering,” Phys. Rev. D 69 (2004) 083002, arXiv:astro-ph/0310198.
  * (30) D. Baumann, D. Green, J. Meyers, and B. Wallisch, “Phases of New Physics in the CMB,” JCAP 01 (2016) 007, arXiv:1508.06342 [astro-ph.CO].
  * (31) B. Follin, L. Knox, M. Millea, and Z. Pan, “First Detection of the Acoustic Oscillation Phase Shift Expected from the Cosmic Neutrino Background,” Phys. Rev. Lett. 115 (2015) no. 9, 091301, arXiv:1503.07863 [astro-ph.CO].
  * (32) J. Wess and B. Zumino, “Supergauge Transformations in Four-Dimensions,” Nucl. Phys. B 70 (1974) 39–50.
  * (33) M. Dine and A. E. Nelson, “Dynamical supersymmetry breaking at low-energies,” Phys. Rev. D 48 (1993) 1277–1287, arXiv:hep-ph/9303230.
  * (34) M. Dine, A. E. Nelson, Y. Nir, and Y. Shirman, “New tools for low-energy dynamical supersymmetry breaking,” Phys. Rev. D 53 (1996) 2658–2669, arXiv:hep-ph/9507378.
  * (35) L. Randall and R. Sundrum, “Out of this world supersymmetry breaking,” Nucl. Phys. B 557 (1999) 79–118, arXiv:hep-th/9810155.
  * (36) G. F. Giudice, M. A. Luty, H. Murayama, and R. Rattazzi, “Gaugino mass without singlets,” JHEP 12 (1998) 027, arXiv:hep-ph/9810442.
  * (37) A. Berlin, N. Blinov, and S. W. Li, “Dark Sector Equilibration During Nucleosynthesis,” Phys. Rev. D 100 (2019) no. 1, 015038, arXiv:1904.04256 [hep-ph].
  * (38) C.-P. Ma and E. Bertschinger, “Cosmological perturbation theory in the synchronous and conformal Newtonian gauges,” Astrophys. J. 455 (1995) 7–25, arXiv:astro-ph/9506072.
  * (39) J. Lesgourgues, “The Cosmic Linear Anisotropy Solving System (CLASS) III: Comparision with CAMB for LambdaCDM,”Journal of Cosmology and Astroparticle Physics 2011 (Jul, 2011) 034–034, arXiv:1104.2934 [astro-ph.CO].
  * (40) T. Brinckmann and J. Lesgourgues, “MontePython 3: boosted MCMC sampler and other features,” Phys. Dark Univ. 24 (2019) 100260, arXiv:1804.07261 [astro-ph.CO].
  * (41) BOSS Collaboration, S. Alam et al., “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample,” Mon. Not. Roy. Astron. Soc. 470 (2017) no. 3, 2617–2652, arXiv:1607.03155 [astro-ph.CO].
  * (42) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, “The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant,” Mon. Not. Roy. Astron. Soc. 416 (2011) 3017–3032, arXiv:1106.3366 [astro-ph.CO].
  * (43) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, “The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$,” Mon. Not. Roy. Astron. Soc. 449 (2015) no. 1, 835–847, arXiv:1409.3242 [astro-ph.CO].
  * (44) D. M. Scolnic et al., “The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample,” Astrophys. J. 859 (2018) no. 2, 101, arXiv:1710.00845 [astro-ph.CO].
  * (45) J. Froustey, C. Pitrou, and M. C. Volpe, “Neutrino decoupling including flavour oscillations and primordial nucleosynthesis,” JCAP 12 (2020) 015, arXiv:2008.01074 [hep-ph].
  * (46) J. J. Bennett, G. Buldgen, M. Drewes, and Y. Y. Y. Wong, “Towards a precision calculation of the effective number of neutrinos $N_{\rm eff}$ in the Standard Model I: the QED equation of state,” JCAP 03 (2020) 003, arXiv:1911.04504 [hep-ph]. [Addendum: JCAP 03, A01 (2021)].
  * (47) M. Raveri and W. Hu, “Concordance and Discordance in Cosmology,” Phys. Rev. D 99 (2019) no. 4, 043506, arXiv:1806.04649 [astro-ph.CO].
  * (48) A. Berlin and N. Blinov, “Thermal Dark Matter Below an MeV,” Phys. Rev. Lett. 120 (2018) no. 2, 021801, arXiv:1706.07046 [hep-ph].
  * (49) A. Berlin and N. Blinov, “Thermal neutrino portal to sub-MeV dark matter,” Phys. Rev. D 99 (2019) no. 9, 095030, arXiv:1807.04282 [hep-ph].
  * (50) J. C. Hill et al., “The Atacama Cosmology Telescope: Constraints on Pre-Recombination Early Dark Energy,” arXiv:2109.04451 [astro-ph.CO].
  * (51) H. Hildebrandt et al., “KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing,” Mon. Not. Roy. Astron. Soc. 465 (2017) 1454, arXiv:1606.05338 [astro-ph.CO].
  * (52) H. Hildebrandt et al., “KiDS+VIKING-450: Cosmic shear tomography with optical and infrared data,” Astron. Astrophys. 633 (2020) A69, arXiv:1812.06076 [astro-ph.CO].
  * (53) DES Collaboration, T. M. C. Abbott et al., “Dark Energy Survey year 1 results: Cosmological constraints from galaxy clustering and weak lensing,” Phys. Rev. D 98 (2018) no. 4, 043526, arXiv:1708.01530 [astro-ph.CO].
  * (54) HSC Collaboration, C. Hikage et al., “Cosmology from cosmic shear power spectra with Subaru Hyper Suprime-Cam first-year data,” Publ. Astron. Soc. Jap. 71 (2019) no. 2, 43, arXiv:1809.09148 [astro-ph.CO].
  * (55) M. A. Buen-Abad, G. Marques-Tavares, and M. Schmaltz, “Non-Abelian dark matter and dark radiation,” Phys. Rev. D 92 (2015) no. 2, 023531, arXiv:1505.03542 [hep-ph].
  * (56) J. Lesgourgues, G. Marques-Tavares, and M. Schmaltz, “Evidence for dark matter interactions in cosmological precision data?,” JCAP 02 (2016) 037, arXiv:1507.04351 [astro-ph.CO].
  * (57) M. A. Buen-Abad, M. Schmaltz, J. Lesgourgues, and T. Brinckmann, “Interacting Dark Sector and Precision Cosmology,” JCAP 01 (2018) 008, arXiv:1708.09406 [astro-ph.CO].

##  Appendix A Background Evolution of a Stepped Fluid

We assume that the interaction rate in the dark fluid is large enough such
that local thermal equilibrium is maintained. At the background (homogeneous)
level, this implies that entropy conservation is sufficient to track the
evolution of the fluid parameters

$\displaystyle N\equiv\frac{\rho}{\rho^{1\nu}}\,,\quad
w\equiv\frac{p}{\rho}\,,\quad
c_{s}^{2}\equiv\frac{dp/dT_{d}}{d\rho/dT_{d}}\,,$ (9)

through the step, where $\rho^{1\nu}$ is the energy density of a single
neutrino in $\Lambda{\rm CDM}$. $N$, $w$, $c_{s}^{2}$ are functions of the
dark sector temperature $T_{d}$, which depends on the scale factor $a$.
Determining $N(a)$, $w(a)$ and $c_{s}^{2}(a)$ is the purpose of this Appendix.

The dark energy density $\rho(T_{d})$ and pressure $p(T_{d})$ of the
interacting dark fluid have contributions from both massive and massless
particles. For a single particle of mass $m$ with $g$ internal degrees of
freedom and temperature $T_{d}$ in equilibrium, the energy density and
pressure are

$\displaystyle\rho(T_{d})$
$\displaystyle=g\int\frac{d^{3}p}{(2\pi)^{3}}\sqrt{p^{2}+m^{2}}\frac{1}{e^{\sqrt{p^{2}+m^{2}}/T_{d}}\pm
1}$ $\displaystyle p(T_{d})$
$\displaystyle=\frac{g}{3}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{2}}{\sqrt{p^{2}+m^{2}}}\frac{1}{e^{\sqrt{p^{2}+m^{2}}/T_{d}}\pm
1}\,,$ (10)

where $\pm 1$ corresponds to fermions/bosons, respectively. It will be
convenient to rewrite this by factoring out the energy density and pressure of
a single relativistic boson degree of freedom,
$\rho_{B}(T_{d})=3p_{B}(T_{d})=T_{d}^{4}\pi^{2}/30$. This gives

$\displaystyle\rho(T_{d})$ $\displaystyle\equiv
g_{*}\,\rho_{B}(T_{d})\,\hat{\rho}(x)$ $\displaystyle p(T_{d})$
$\displaystyle\equiv g_{*}\,p_{B}(T_{d})\,\hat{p}(x)\,,$ (11)

where $x=m/T_{d}$, and we have defined the effective degrees of freedom,
$g_{*}=g$ for bosons and $g_{*}=(7/8)\,g$ for fermions. Note that the
dimensionless integrals $\hat{\rho}$ and $\hat{p}$ depend on the temperature
and mass only through their ratio $x$. $\hat{\rho}(0)=\hat{p}(0)=1$ and both
decay exponentially for large $x$. $\hat{\rho}$ and $\hat{p}$ are quite
similar for bosons and fermions and we approximate them by dropping the $\pm
1$ terms in the distribution functions. This approximation works to better
than 8% for all $x$, and to better than 4% for $x$ for which $\hat{\rho}$ and
$\hat{p}$ are not exponentially suppressed and therefore irrelevant.

The (Maxwell-Boltzmann) integrals can now be evaluated analytically

$\displaystyle\hat{\rho}(x)$
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\\!dq\,q^{2}\sqrt{q^{2}+x^{2}}\
e^{-\sqrt{q^{2}+x^{2}}}\bigg{/}\int_{0}^{\infty}\\!dq\,q^{3}\ e^{-q}$
$\displaystyle=\frac{x^{2}}{2}K_{2}(x)+\frac{x^{3}}{6}K_{1}(x)$
$\displaystyle\hat{p}(x)$
$\displaystyle=\int_{0}^{\infty}\\!\\!\\!\\!dq\,\frac{q^{4}}{\sqrt{q^{2}+x^{2}}}\
e^{-\sqrt{q^{2}+x^{2}}}\bigg{/}\int_{0}^{\infty}\\!\\!\\!\\!dq\,q^{3}\,e^{-q}$
$\displaystyle=\frac{x^{2}}{2}K_{2}(x)\,,$ (12)

in terms of modified Bessel functions of the second kind, $K_{i}$.

The energy density of a fluid which contains $g_{*}^{\text{IR}}$ massless and
$g_{*}^{\text{UV}}-g_{*}^{\text{IR}}$ massive particles can then be written
as666For simplicity, we assume a single mass scale $m$. It is straightforward
to generalize to several particles with different masses.

$\displaystyle\rho(T_{d})$ $\displaystyle=$ $\displaystyle
g_{*}^{\text{IR}}\rho_{B}(T_{d})(1+r_{g}\hat{\rho}(x))\,,$ (13)

where we defined the step size
$r_{g}\equiv(g_{*}^{\text{UV}}-g_{*}^{\text{IR}})/g_{*}^{\text{IR}}$. The
pressure is

$\displaystyle p(T_{d})$ $\displaystyle=$ $\displaystyle
g_{*}^{\text{IR}}p_{B}(T_{d})(1+r_{g}\hat{p}(x))\ .$ (14)

Eqs. (13) and (14) again include a factor of $7/8$ in the definition of
$g_{*}^{\text{IR}}$ and $g_{*}^{\text{UV}}$ for each fermion. Therefore, the
equation of state is

$\displaystyle
w(x)=\frac{1}{3}\frac{1+r_{g}\hat{p}(x)}{1+r_{g}\hat{\rho}(x)}=\frac{1}{3}-\frac{r_{g}}{3}\,\frac{\hat{\rho}(x)-\hat{p}(x)}{1+r_{g}\hat{\rho}(x)}\,,$
(15)

and the speed of sound is

$\displaystyle c_{s}^{2}(x)$
$\displaystyle=\frac{dp/dx}{d\rho/dx}=\frac{1}{3}\frac{1+r_{g}(\hat{p}(x)-\frac{x}{4}\hat{p}^{\prime}(x))}{1+r_{g}(\hat{\rho}(x)-\frac{x}{4}\hat{\rho}^{\prime}(x))}$
$\displaystyle=\frac{1}{3}-\frac{r_{g}}{36}\,\frac{x^{2}\hat{p}(x)}{1+r_{g}(\frac{3}{4}\hat{\rho}(x)+(\frac{1}{4}+\frac{x^{2}}{12})\hat{p}(x))}\,,$
(16)

where in the second line we used Bessel function recursion relations to
rewrite $\hat{\rho}^{\prime}$ and $\hat{p}^{\prime}$ in terms of $\hat{\rho}$
and $\hat{p}$. We see that both $w$ and $c_{s}^{2}$ are 1/3 minus a term
proportional to the step size $r_{g}$ which vanishes for both small and large
$x$ (see Fig. 1).

The expressions above are written in terms of $x=m/T_{d}$, but we will need
them as functions of scale factor $a$. Therefore, we determine $x(a)$ from
entropy conservation in the dark fluid

$\displaystyle S\propto
a^{3}\frac{\rho(T_{d})+p(T_{d})}{T_{d}}=a_{0}^{3}\frac{\rho(T_{d0})+p(T_{d0})}{T_{d0}}\,,$
(17)

where $a_{0}=1$ is the scale factor today and $T_{d0}$ is the dark sector
temperature today. Substituting the expressions for $\rho$ and $p$ from Eqs.
(13) and (14), assuming that the massive particles have annihilated away today
so that $\hat{\rho}(m/T_{d0})=\hat{p}(m/T_{d0})=0$, and defining the
transition scale factor $a_{t}\equiv T_{d0}/m$, we obtain

$\displaystyle\left(\frac{xa_{t}}{a}\right)^{3}=1+\frac{r_{g}}{4}(3\hat{\rho}(x)+\hat{p}(x))\,,$
(18)

which can be solved numerically for $x(a)$ for any given set of model
parameters $r_{g}$ and $a_{t}$. Note that well after the transition,
$\hat{\rho}(x)$ and $\hat{p}(x)$ are negligibly small. In this case, Eq. (18)
simplifies to $xa_{t}=a$, or equivalently $T_{d}=T_{d0}/a$. At early times,
well before the transition, $x\ll 1$ so that $\hat{\rho}(0)=\hat{p}(0)=1$ and
therefore $xa_{t}=(1+r_{g})^{1/3}a$ or $T_{d}=(1+r_{g})^{-1/3}T_{d0}/a$. Thus,
the temperature of the dark sector in the UV is smaller than the naive $1/a$
scaling by a factor of $(1+r_{g})^{-1/3}$ because of the step.

Finally, we can calculate the “effective number of neutrino species” from Eq.
(9)

$\displaystyle N(x)$ $\displaystyle=$ $\displaystyle
N_{\text{IR}}\frac{1+r_{g}\hat{\rho}(x)}{(1+r_{g}(\frac{3}{4}\hat{\rho}(x)+\frac{1}{4}\hat{p}(x)))^{4/3}}\,,$
(19)

where we used $\rho^{1\nu}=\frac{7}{4}\frac{\pi^{2}}{30}(\frac{T_{\nu
0}}{a})^{4}$, substituted Eq. (18), and identified

$\displaystyle N_{\text{IR}}$ $\displaystyle=$
$\displaystyle\frac{g_{*}^{\text{IR}}}{7/4}\left(\frac{T_{d0}}{T_{\nu
0}}\right)^{4}\,,$ (20)

by taking the IR ($x\rightarrow\infty$) limit. Taking the UV ($x\rightarrow
0$) limit gives the UV endpoint of the step

$\displaystyle N_{\text{UV}}$ $\displaystyle=$
$\displaystyle\frac{N_{\text{IR}}}{(1+r_{g})^{1/3}}=N_{\text{IR}}\left(\frac{g_{*}^{\text{UV}}}{g_{*}^{\text{IR}}}\right)^{1/3}\,.$
(21)

This equation confirms Eq. (2) and relates the parameter $r_{g}$ to
$N_{\text{IR}}/N_{\text{UV}}$ by

$\displaystyle r_{g}=\left(\frac{N_{\text{IR}}}{N_{\text{UV}}}\right)^{3}-1\
.$ (22)

In summary, the background quantities describing the fluid are $w(x)$,
$c_{s}^{2}(x)$, and $N(x)$, which are given in Eqs. (15), (A), and (19).
$x(a)$ is obtained by numerically solving (18), $r_{g}$ is given in (22), and
$\hat{\rho}(x)$ and $\hat{p}(x)$ in (A). These are the equations we implement
in CLASS, along with the standard perturbation equations for a strongly-
coupled fluid.

The parameters of the stepped dark fluid are the transition scale $a_{t}$ (or
equivalently, $z_{t}=1/a_{t}-1$), the amount of dark radiation in the IR,
$N_{\text{IR}}$, and the ratio $r_{g}$ determining the size of the step (or
equivalently, $N_{\text{IR}}/N_{\text{UV}}$). In the WZDR model, $r_{g}=8/7$.

##  Appendix B Triangle Plots and Parameter Values

In this Appendix, we provide triangle plots and tables (shown below), which
show the full set of posteriors and parameter values as determined from our
analysis, respectively. In Fig. 4, we compare the posteriors for the WZDR and
General StepDR dark fluid models when fitting to the dataset $\cal{D}+$
including SH0ES. In Fig. 5, we compare the posteriors for the WZDR model when
fitting to a dataset that does ($\cal{D}+$) or does not ($\cal{D}$) include
SH0ES. In Fig. 6, we compare the posteriors for the WZDR, General StepDR,
SIDR, $\Lambda\text{CDM}+N_{\rm eff}$, and $\Lambda{\rm CDM}$ models when
fitting to the dataset $\cal{D}+$. Table 4 shows $\chi^{2}$ values as well as
other fit parameters and measures relevant to evaluating the success of
various models in resolving the Hubble tension. Furthermore, Tables 5 and 6
provide posterior mean-value and $\pm 1\sigma$ ranges for various cosmological
parameters when fitting to either $\cal{D}$ or $\cal{D}+$, whereas Tables 7
and 8 show best-fit values.

Figure 4: A comparison of the posteriors for the minimal stepped fluid (WZDR in blue) and generalized stepped fluid (General StepDR in red) when fitting to the dataset $\cal{D}+$ that includes SH0ES. The dark and light shaded regions correspond to $68.3\%$ and $95.4\%$ C.L., respectively. It is interesting to note that while the introduction of two additional free parameters in the General StepDR model widens the posterior density for $N_{\text{UV}}$ and $N_{\text{IR}}$, the favored region of parameter space is approximately unchanged compared to the restricted WZDR model. Figure 5: A comparison of the posteriors for the minimal stepped fluid (WZDR) when fitting to either the SH0ES-independent dataset $\cal{D}$ (red) or the dataset $\cal{D}+$ that includes SH0ES (blue). The dark and light shaded regions correspond to $68.3\%$ and $95.4\%$ C.L., respectively. Comparing the two sets of posteriors, we see that a fit to data including SH0ES prefers a larger value of $H_{0}$ and a correspondingly larger energy density in dark radiation, as expected, while leaving the preferred location of the transition $z_{t}$ nearly unchanged. Since an enhanced radiation density strengthens the effect of the step on the CMB, the data is increasingly sensitive to the location of the transition, as can be noted from the narrower posterior for $z_{t}$ when fitting to $\cal{D}+$. Along with the increase in radiation density, the fit to $\cal{D}+$ also requires increases in $\omega_{\rm cdm}$, $\omega_{b}$, as well as $\theta_{s}$ and $n_{s}$. Figure 6: A comparison of the posteriors for models of a minimal stepped fluid (WZDR in green), conventional self-interacting dark radiation (SIDR in blue), additional free-streaming radiation ($\Lambda\text{CDM}+N_{\rm eff}$ in red), and $\Lambda{\rm CDM}$ (gray) when fitting to the dataset $\cal{D}+$ that includes SH0ES. The dark and light shaded regions correspond to $68.3\%$ and $95.4\%$ C.L., respectively. Most notable is the fact that WZDR predicts $H_{0}$ and $N_{\text{eff,IR}}$ posteriors that extend out to much larger values, reflecting its ability to more successfully resolve the Hubble tension. Model | $\Delta N_{param}$ | $M_{b}+1\sigma$ | GT | $Q_{\text{DMAP}}$ | min $\chi^{2}_{\cal{D}+}$ | min $\chi^{2}_{\cal{D}}$ | $\Delta\chi^{2}_{{\cal D+}}$ | $\Delta\text{AIC}$ |  |  $H_{0}^{BF}$ $\cal{D}$ |  $H_{0}^{Mean}$ $\cal{D}$ |  $H_{0}^{BF}$ ${\cal D+}$ |  $H_{0}^{Mean}$ ${\cal D+}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\Lambda{\rm CDM}$ | 0 | $-19.418+0.011$ | 4.5 | 4.5 | 3828.3 | 3808.5 | 0.00 | 0.00 |  | 67.6 | 67.6 | 68.2 | 68.2
$\Lambda\text{CDM}+N_{\rm eff}$ | 1 | $-19.396+0.017$ | 3.7 | 3.8 | 3822.6 | 3808.5 | $-5.7$ | $-3.7$ |  | 67.7 | 68.4 | 70.0 | 70.2
SIDR (standard BBN) | 1 | $-19.381+0.021$ | 3.2 | 3.3 | 3819.1 | 3808.4 | $-9.2$ | $-7.2$ |  | 67.8 | 68.9 | 70.6 | 70.8
SIDR | 1 | $-19.380+0.023$ | 3.1 | 3.1 | 3817.7 | 3808.4 | $-10.6$ | $-8.6$ |  | 67.9 | 68.9 | 71.0 | 71.0
WZDR | 2 | $-19.351+0.028$ | 2.3 | 2.4 | 3813.2 | 3807.3 | $-15.1$ | $-11.1$ |  | 69.1 | 69.9 | 71.4 | 71.7
General StepDR | 4 | $-19.382+0.041$ | 2.5 | 2.4 | 3811.7 | 3805.9 | $-16.6$ | $-8.6$ |  | 67.9 | 68.8 | 71.5 | 71.5
Table 4: A summary of the fit results for various models. The two versions of SIDR refer to different assumptions about when the extra radiation is populated. “SIDR (standard BBN)” refers to the model in which the energy density of interacting radiation is present during BBN, whereas “SIDR” assumes that this radiation was not yet present. In the last four columns, $H_{0}$ is given in units of $\text{km}/\text{s}/\text{Mpc}$. | $\Lambda{\rm CDM}$ | $\Lambda{\rm CDM}+N_{\rm eff}$ | SIDR | WZDR | General StepDR
---|---|---|---|---|---
100$\theta_{s}$ | $1.04195^{+0.00028}_{-0.00027}$ | $1.04166^{+0.00039}_{-0.00034}$ | $1.04208^{+0.00032}_{-0.00031}$ | $1.04260^{+0.00039}_{-0.00043}$ | $1.04339^{+0.00080}_{-0.00095}$
$\Omega_{b}h^{2}$ | $0.02240^{+0.00013}_{-0.00014}$ | $0.02248^{+0.00015}_{-0.00016}$ | $0.02255^{+0.00016}_{-0.00018}$ | $0.02257^{+0.00016}_{-0.00017}$ | $0.02249^{+0.00019}_{-0.00020}$
$\Omega_{\rm cdm}h^{2}$ | $0.11931^{+0.00091}_{-0.00094}$ | $0.1213^{+0.0013}_{-0.0022}$ | $0.1226^{+0.0019}_{-0.0032}$ | $0.1296^{+0.0026}_{-0.0035}$ | $0.1227^{+0.0037}_{-0.0042}$
$\ln 10^{10}A_{s}$ | $3.049^{+0.013}_{-0.016}$ | $3.054^{+0.015}_{-0.016}$ | $3.046\pm 0.015$ | $3.051^{+0.014}_{-0.015}$ | $3.038^{+0.021}_{-0.019}$
$n_{s}$ | $0.9662^{+0.0040}_{-0.0035}$ | $0.9703^{+0.0045}_{-0.0055}$ | $0.9670\pm 0.0038$ | $0.9737^{+0.0049}_{-0.0057}$ | $0.9664^{+0.0083}_{-0.0092}$
$\tau_{\rm reio}$ | $0.0571^{+0.0069}_{-0.0078}$ | $0.0570^{+0.0071}_{-0.0077}$ | $0.0578^{+0.0075}_{-0.0080}$ | $0.0578^{+0.0072}_{-0.0074}$ | $0.0573^{+0.0069}_{-0.0078}$
$N_{\text{fs}}$ | - | $3.17^{+0.03}_{-0.12}$ | $3.044$ | $3.044$ | $2.75^{+0.25}_{-0.25}$
$N_{\text{UV}}$ | - | - | - | $0.29^{+0.11}_{-0.16}$ | $0.39^{+0.15}_{-0.30}$
$N_{\text{IR}}$ | - | - | $0.19^{+0.05}_{-0.19}$ | $0.38^{+0.14}_{-0.21}$ | $0.48^{+0.18}_{-0.30}$
$\log_{10}(z_{t})$ | - | - | - | $4.27^{+0.14}_{-0.17}$ | $4.35^{+0.18}_{-0.11}$
$M_{b}$ | ${-19.418}^{+0.011}_{-0.012}$ | ${-19.396}^{+0.017}_{-0.025}$ | ${-19.380}^{+0.023}_{-0.036}$ | ${-19.351}^{+0.028}_{-0.037}$ | ${-19.382}^{+0.041}_{-0.045}$
$H_{0}$ [km/Mpc/s] | $67.6\pm 0.4$ | $68.4^{+0.6}_{-0.8}$ | $68.9^{+0.8}_{-1.2}$ | $69.9^{+0.9}_{-1.3}$ | $68.8^{+1.3}_{-1.5}$
$S_{8}$ | $0.825\pm 0.010$ | $0.828\pm 0.011$ | $0.824\pm 0.011$ | $0.830\pm 0.011$ | $0.823\pm 0.013$
Table 5: Mean and $\pm 1\sigma$ values for a fit to dataset $\cal{D}$. | $\Lambda{\rm CDM}$ | $\Lambda{\rm CDM}+N_{\rm eff}$ | SIDR | WZDR | General StepDR
---|---|---|---|---|---
100$\theta_{s}$ | $1.04207^{+0.00029}_{-0.00031}$ | $1.04120\pm 0.00042$ | $1.04233^{+0.00029}_{-0.00030}$ | $1.04309^{+0.00037}_{-0.00035}$ | $1.04324^{+0.00088}_{-0.00096}$
$\Omega_{b}h^{2}$ | $0.02252\pm 0.00014$ | $0.02274^{+0.00016}_{-0.00015}$ | $0.02282\pm 0.00016$ | $0.02272^{+0.00014}_{-0.00016}$ | $2.274^{+0.018}_{-0.017}$
$\Omega_{\rm cdm}h^{2}$ | $0.11817^{+0.00097}_{-0.00092}$ | $0.1245^{+0.0026}_{-0.0025}$ | $0.1269^{+0.0029}_{0.0028}$ | $0.1296^{+0.0029}_{-0.0031}$ | $0.1294^{+0.0032}_{-0.0033}$
$\ln 10^{10}A_{s}$ | $3.054^{+0.015}_{-0.016}$ | $3.067^{+0.016}_{-0.017}$ | $3.044^{+0.015}_{-0.017}$ | $3.052\pm 0.015$ | $3.045\pm 0.020$
$n_{s}$ | $0.9691^{+0.0037}_{-0.0039}$ | $0.9802^{+0.0057}_{-0.0055}$ | $0.9694^{+0.0037}_{-0.0036}$ | $0.9801^{+0.0047}_{-0.0049}$ | $0.9752^{+0.0082}_{-0.0086}$
$\tau_{\rm reio}$ | $0.0606^{+0.0071}_{-0.0086}$ | $0.0598^{+0.0075}_{-0.0082}$ | $0.0599^{+0.0070}_{-0.0084}$ | $0.0585^{+0.0072}_{-0.0081}$ | $0.0585^{+0.0071}_{-0.0078}$
$N_{\text{fs}}$ | - | $3.41^{+0.13}_{-0.15}$ | $3.044$ | $3.044$ | $2.93\pm 0.27$
$N_{\text{UV}}$ | - | - | - | $0.51^{+0.12}_{-0.13}$ | $0.601^{+0.24}_{-0.27}$
$N_{\text{IR}}$ | - | - | $0.47^{+0.15}_{-0.14}$ | $0.65\pm 0.16$ | $0.74^{+0.26}_{-0.29}$
$\log_{10}(z_{t})$ | - | - | - | $4.26^{+0.12}_{-0.11}$ | $4.30^{+0.09}_{-0.13}$
$M_{b}$ | ${-19.404}\pm 0.012$ | ${-19.342}\pm 0.026$ | ${-19.320}^{+0.027}_{-0.028}$ | ${-19.299}^{+0.027}_{-0.026}$ | ${-19.302}^{+0.031}_{-0.029}$
$H_{0}$ [km/Mpc/s] | $68.2\pm 0.4$ | $70.2\pm 0.9$ | $71.0^{+0.9}_{-1.0}$ | $71.7^{+1.0}_{-0.9}$ | $71.5^{+1.1}_{-1.0}$
$S_{8}$ | $0.814\pm 0.011$ | $0.827\pm 0.012$ | $0.818\pm 0.010$ | $0.829\pm 0.012$ | $0.823\pm 0.013$
Table 6: Mean and $\pm 1\sigma$ values for a fit to dataset $\cal{D}+$. | $\Lambda{\rm CDM}$ | $\Lambda{\rm CDM}+N_{\rm eff}$ | SIDR | WZDR | General StepDR
---|---|---|---|---|---
100$\theta_{s}$ | 1.04192 | 1.04193 | 1.04196 | 1.04236 | 1.04299
$\Omega_{b}h^{2}$ | 0.022406 | 0.022410 | 0.022450 | 0.022533 | 0.022416
$\Omega_{\rm cdm}h^{2}$ | 0.11930 | 0.11930 | 0.11980 | 0.12320 | 0.12037
$\ln 10^{10}A_{s}$ | 3.049 | 3.049 | 3.049 | 3.049 | 3.045
$n_{s}$ | 0.9666 | 0.9666 | 0.9674 | 0.9724 | 0.9693
$\tau_{\rm reio}$ | 0.0571 | 0.0571 | 0.0572 | 0.0565 | 0.0573
$N_{\text{fs}}$ | 3.044 | 3.044 | 3.044 | 3.044 | 2.85
$N_{\text{UV}}$ | - | - | - | 0.18 | 0.16
$N_{\text{IR}}$ | - | - | 0.03 | 0.23 | 0.23
$\log_{10}(z_{t})$ | - | - | - | 4.30 | 4.46
$M_{b}$ | ${-19.418}$ | ${-19.417}$ | ${-19.411}$ | ${-19.375}$ | ${-19.409}$
$H_{0}$ [km/Mpc/s] | 67.6 | 67.7 | 67.9 | 69.1 | 67.9
$S_{8}$ | 0.825 | 0.825 | 0.825 | 0.828 | 0.827
$\chi^{2}_{\rm tot}$ | 3808.5 | 3808.5 | 3808.4 | 3807.3 | 3805.9
Table 7: Best-fit values for a fit to dataset $\cal{D}$. | $\Lambda{\rm CDM}$ | $\Lambda{\rm CDM}+N_{\rm eff}$ | SIDR | WZDR | General StepDR
---|---|---|---|---|---
100$\theta_{s}$ | $1.04204$ | $1.04129$ | $1.04234$ | $1.04303$ | $1.04331$
$\Omega_{b}h^{2}$ | $0.022531$ | $0.022712$ | $0.022817$ | $0.022718$ | $0.022726$
$\Omega_{\rm cdm}h^{2}$ |  $0.1183$0 | $0.12372$ | $0.12681$ | $0.12880$ | $0.12952$
$\ln 10^{10}A_{s}$ | $3.052$ | $3.066$ | $3.043$ | $3.053$ | $3.045$
$n_{s}$ | $0.9692$ | $0.9798$ | $0.9693$ | $0.9801$ | $0.9777$
$\tau_{\rm reio}$ | $0.0597$ | $0.0599$ | $0.0594$ | $0.0587$ | $0.0575$
$N_{\text{fs}}$ | $3.044$ | $3.37$ | $3.044$ | $3.044$ | $2.94$
$N_{\text{UV}}$ | - | - | - | $0.46$ | $0.57$
$N_{\text{IR}}$ | - | - | $0.47$ | $0.59$ | $0.72$
$\log_{10}(z_{t})$ | - | - | - | $4.27$ | $4.33$
$M_{b}$ | ${-19.404}$ | ${-19.348}$ | ${-19.320}$ | ${-19.307}$ | ${-19.301}$
$H_{0}$ [km/Mpc/s] | $68.2$ | $70.0$ | 71.0 | 71.4 | 71.5
$S_{8}$ | $0.815$ | $0.825$ | $0.817$ | $0.830$ | $0.825$
$\chi^{2}_{\cal{D}}$ | $3809.9$ | $3814.5$ | $3813.5$ | $3810.3$ | $3809.37$
$\chi^{2}_{\rm tot}$ | $3828.3$ | $3822.6$ | $3817.7$ | $3813.2$ | $3811.7$
Table 8: Best-fit values for a fit to dataset $\cal{D}+$. The value of
$\chi^{2}_{\cal{D}}$ (in the second to last row) corresponds to the best-fit
points from a fit to dataset $\cal{D}+$ but only includes $\chi^{2}$
contributions from dataset $\cal{D}$.