degeneracy with $h$ is allowed by BAO data, since it and the prior on
$R_{\star}$ and $x_{\mathrm{eq}}$ are only sensitive to $\alpha_{i}$ via its
effect on $a_{\star}$. Comparing these results with those including Planck
data highlights how precisely CMB observations can measure the early-time fine
structure constant. The datasets combined constrain the early-time electron
mass to deviate at only the percent level and the fine-structure constant at
the $0.2\%$ level from their values measured today.

For neither varying $\alpha_{i}$ nor $m_{e,i}$ do the posteriors in Fig. 12
using only BAO data lie precisely along $h\propto 1/a_{\star}$, as would be
expected from Section II.2.2 if $a_{m-\Lambda}$ were well constrained. The BAO
measurements we employ in reality only cover a relatively narrow range of
redshifts and are insufficient in number and precision to fully measure
$a_{m-\Lambda}$ independently of the overall normalization in Eq. 33, i.e.,
they remain slightly correlated. While eBOSS provides additional BAO
measurements from tracers beyond LRGs like quasars, emission line galaxies,
and Lyman-$\alpha$ absorption [173], cosmological constraints from these are
known to be mildly inconsistent [174, 175, 176, 177]. However, DESI DR1’s
results from various spanning redshifts $\sim 0.3$ to $2.33$ may be more
justifiably combined [156]. We compare to results with Planck combined instead
with DESI data in Fig. 13; the posterior shifts over the $h-a_{m-\Lambda}$
degeneracy to larger values of $h$ and smaller $a_{m-\Lambda}$. This shift is
driven not just by DESI’s preference for an early onset of dark-energy
domination (evident in its $\Lambda$CDM results as well [156]) but also for a
slightly larger uncalibrated amplitude $hr_{\mathrm{d}}$. However, we note
that the posterior of $hr_{\mathrm{d}}$ shifts to lower values by several
percent when excluding DESI DR1’s LRG results in redshift bins $z=0.51$ and
$0.71$; these data appear to be mild outliers with a fair amount of leverage
when interpreted within $\Lambda$CDM cosmology [156], so we caution against
interpreting this result too seriously.

Type Ia supernova datasets can also measure the scale factor of matter–dark-
energy equality $a_{m-\Lambda}$, and they do so independently of the sound
horizon and with distance measurements out to slightly higher redshifts than
LRG-derived BAO measurements. Since the degeneracy direction that preserves
$\theta_{s}$ does not hold $a_{m-\Lambda}$ fixed (i.e., has $a_{m-\Lambda}$
strongly correlated with $m_{e,i}$), these datasets should yield different
preferences for $m_{e,i}$ (and therefore $h$) if they prefer different
$a_{m-\Lambda}$.

Figure 13:  Posterior distribution over $h$ and $a_{m-\Lambda}$ (a derived
parameter) in $\Lambda$CDM (left) and varying-$m_{e}$ scenarios (right),
comparing results that use the Planck 2018 likelihoods combined with one of
the Pantheon+ (blue), SDSS+ BAO (red), and DESI DR1 BAO (gold) datasets. The
diagonal panels display kernel density estimates of the one-dimensional,
marginalized posterior relative to their peak value, and the lower panels
depict 1 and $2\sigma$ mass levels of the two-dimensional joint posterior.
Note that the left panels span a much broader range in $h$ and $a_{m-\Lambda}$
than do the right ones. The grey line in the lower left panel depicts the
expected degeneracy inferred from Eq. 27, i.e., with $\omega_{m}\propto
m_{e,i}$ and $h\propto m_{e,i}^{3.18}$.

Figure 13 shows that this is the case, to an extent: the medians of the
marginalized posteriors for $h$ vary by $3-8\%$ percent when Planck data is
combined with one of the SDSS+ BAO dataset, DESI’s DR1 BAO measurements, or
the Pantheon+ dataset. Because $h$ varies rather sensitively with $m_{e,i}$
along the posteriors, the corresponding variation in inferred values for
$m_{e,i}$ is less than $2.5\%$. Further, the posteriors in the
$h-a_{m-\Lambda}$ plane show a strong correlation as expected, affirming that
constraints from the SN datasets (which include no external calibration on the
fiducial brightness magnitude $M_{B}$) derive from their measurement of the
shape of the expansion history.

The $\Lambda$CDM posteriors in Fig. 13 are more mutually consistent across
dataset combinations in spite of being much narrower than those that vary the
electron mass, since in $\Lambda$CDM $h$ is constrained to simultaneously
satisfy Planck’s precise measurement of $\theta_{s}$. The varying-$m_{e}$
model, not so restricted, allows the posteriors on $a_{m-\Lambda}$ to be
mostly driven by late-time datasets. Indeed, the constraints on
$a_{m-\Lambda}$ in the left panels of Fig. 13 nearly match those each
individually infers for the flat $\Lambda$CDM model [96, 173, 156]. Similar
conclusions apply to the fine-structure constant, since it also affects the
scale factor of recombination, but to a lesser degree given that the CMB
constrains $\alpha_{i}$ independently of its effect on $a_{\star}$.

###  IV.2 Coupled, hyperlight scalars

We now study the effect of consistently accounting for a hyperlight scalar
hypothesized to couple to the electron or photon, realizing the fundamental-
constant variations in Section IV.1. We numerically implement the scalar
field’s dynamics as described in Section III (with further details in Appendix
A). In addition to the parameters specified in Section IV, we sample over the
scalar’s present energy density relative to that of CDM, i.e.,
${F_{\phi}}=\omega_{\phi}/\omega_{c}\sim\mathcal{U}(0.0,0.3)$. Directly
sampling over the scalar’s mass $m_{\phi}$ is challenging because observables
are only weakly sensitive to its value (within the mass range we consider);
again following Refs. [128, 140, 141, 142, 139], we therefore fix $m_{\phi}$
to a single value in parameter inference. We take a fiducial mass
$m_{\phi}=10^{-30}\leavevmode\nobreak\ \mathrm{eV}$, commenting on other
choices in Section IV.2.2. Because we only consider scalars that begin
oscillating in the matter era, ${F_{\phi}}$ uniquely determines the scalar’s
initial misalignment $\bar{\phi}_{i}$ via Eq. 49. We take purely adiabatic
initial conditions, assuming any isocurvature perturbations in $\phi$
generated during inflation are negligible.

####  IV.2.1 Planck data alone

The arguments of Section III that cosmological data should favor a hyperlight
scalar with energy contribution proportional to the increase in the electron
mass hinge upon the assumption that $\Lambda$CDM provides an effectively
optimal fit to features determined by late-time dynamics. In particular,
Section III shows that when _increasing_ $h$ to compensate for a smaller sound
horizon at larger $m_{e}$, the scalar’s contribution to the matter density can
restore the relative amount of dark energy and matter (i.e., the scale factor
at which they are equal, $a_{m-\Lambda}$). Though Planck data by itself does
not favor the parameter space in which this effect is relevant, a hyperlight
scalar still opens up additional freedom in the late-time expansion history
that could allow for larger $m_{e,i}$ without altering $\theta_{s}$.
Constraining the model with Planck 2018 likelihoods alone therefore tests the
extent to which they can tolerate hyperlight subcomponents of dark matter (or
whether any features in the data are actually better explained by hyperlight
scalar fields).

Figure 14 shows that Planck 2018 data only allows a hyperlight scalar to
marginally affect the results without a scalar (Section IV.1.1).

Figure 14:  Posterior distribution over $h$, $m_{e,i}/m_{e,0}$ (left) or
$\alpha_{i}/\alpha_{0}$ (right), and the fraction ${F_{\phi}}$ of the CDM
density contributed by a scalar field with mass $10^{-30}\leavevmode\nobreak\
\mathrm{eV}$. Depicted are results for varying constant scenarios without
(blue) and with (red) scalar fields as well as for $\Lambda$CDM with a scalar
field (grey), all using Planck 2018 likelihoods alone. The diagonal panels
display kernel density estimates of the one-dimensional, marginalized
posterior relative to their peak value, and the lower panels depict 1 and
$2\sigma$ mass levels of the two-dimensional joint posterior, again slightly
smoothed to facilitate comparison.

The posteriors over $h$ and $m_{e,i}$ peak at slightly smaller values, but the
degeneracy between $h$ and $m_{e,i}$ is largely unmodified by including a
scalar field. The posterior over ${F_{\phi}}$ in this case is also marginally
broader than that without any SM couplings, but only at lower $h$. At these
smallest values of $h$, the late ISW effect due to the onset of the dark
energy era is negligible, opening up room for a hyperlight scalar to induce
more of a ISW contribution to the CMB on large scales. The $95$th percentile
of the marginalized posteriors over ${F_{\phi}}$ thus increase from $0.0127$
when no constants vary to $0.0152$ when $m_{e,i}$ varies. (We note that our
results without varying constants are consistent with those from Ref. [139],
which obtains bounds of order $10^{-2}$ for masses $10^{-30}\lesssim
m_{\phi}/\mathrm{eV}\lesssim 10^{-28}$.) Such strong constraints demonstrate
how sensitive Planck is to the dynamics of dark matter perturbations.

For scalars with photon couplings (which are already strongly constrained by
high-multipole Planck measurements) we might expect the constraints on
${F_{\phi}}$ to correspond closely to constraints on hyperlight scalar fields
without varying constants. Figure 14 corroborates this expectation: the
marginalized posterior over ${F_{\phi}}$ is unchanged relative to
$\Lambda$CDM. The results in Fig. 14 do, on the other hand, exhibit a slight
anticorrelation between ${F_{\phi}}$ and $h$ in all cases, even for
$\Lambda$CDM with a hyperlight scalar field. This effect may in part be again
understood via the arguments of Section III.1 regarding the scalar’s effect on
the distance to last scattering: at constant $a_{\star}$, Eq. 43 shows that
fixing the size of the sound horizon $\theta_{s}$ requires $h\propto
A_{m}^{-2.09}\sim(1+{F_{\phi}})^{-2.09}$. For the varying-$m_{e}$ scenario,
the impact of this effect is negligible compared to the degeneracy between $h$
and $m_{e,i}$. As such, the CMB residuals depicted in Fig. 11 are
qualitatively unchanged by the inclusion of a scalar, and similar conclusions
may be drawn.

####  IV.2.2 Impact of low-redshift distance measurements

While Planck’s preference is for $m_{e,i}<m_{e,0}$, preempting the regime in
which the scalar’s could realize a novel degeneracy direction, the features in
the data that drive it hold only marginal constraining power. We might
therefore only expect this degeneracy to bear out when including direct
measurements of the late-time expansion history that would yield stronger
constraints on $a_{m-\Lambda}$. In particular, per Section III.1, fixing
$\theta_{s}$ requires a disproportionate increase in $\omega_{\Lambda}$
relative to that in $\omega_{c}$ and $\omega_{b}$ (imposed by the structure of
the acoustic peaks); late dark matter provides a means to restore the relative
amount of dark energy and matter at late times. Including a hyperlight scalar
field in the cosmological background thus opens up a degeneracy that can
satisfy both CMB and BAO measurements of the angular size of the sound horizon
at any value of $m_{e,i}$. On the other hand, Fig. 14 show that the CMB
disfavors a substantial contribution from a hyperlight scalar, even without
varying constants; the arguments of Section III.2 suggest this is due to its
effects on the dynamics of perturbations. We now test whether including low-
redshift distance data drives the former effect to outweigh the CMB’s
constraints on the latter and realizes the degeneracy anticipated in Section
III.1. Given that BAO data has a substantial impact on $m_{e,i}$ constraints
(Fig. 12), we also assess whether other low-redshift distance datasets
independently exhibit similar preferences.

Figure 15 indicates that low-redshift distance data do not outweigh Planck’s
constraining power on nonclustering subcomponents of dark matter: the
degeneracy of Section III.1 where $h\propto m_{e,i}$ is not realized.

Figure 15:  Posterior distribution over a subset of the parameters of the
electron-coupled scalar model (with $m_{\phi}=10^{-30}\leavevmode\nobreak\
\mathrm{eV}$) and the scale factor of matter–dark-energy equality
$a_{m-\Lambda}$. Posteriors use the likelihood for the full set of Planck 2018
data separately combined with the Pantheon+ (blue), SDSS+ BAO (red), and DESI
DR1 BAO (gold) datasets. The 1 and $2\sigma$ contours of the marginalized two-
dimensional posteriors are depicted as in Fig. 14. Grey lines depict the
degeneracy inferred from Eq. 27 for scenarios _without_ a scalar field, i.e.,
with $\omega_{m}\propto m_{e,i}$ and $h\propto m_{e,i}^{3.18}$.

The posteriors in the $h$-$m_{e,i}$ plane still lie along the same parameter
direction as for phenomenological varying constants (Figs. 12 and 13). The
joint posteriors over $m_{e,i}$ and ${F_{\phi}}$ _do_ , however, display some
degree of correlation, indicating that the relationship expected in Section
III.1 would hold more clearly if the CMB did not disfavor the effect of
hyperlight scalars on spatial perturbations. Indeed, the marginalized
posteriors over ${F_{\phi}}$ are identical among the four cases, and the
constraints are (marginally) sharpened compared to that for Planck alone (Fig.
14).

Finally, we note that our electron-coupling results are nearly insensitive to
the scalar’s mass $m_{\phi}$ in the range $10^{-31}\leavevmode\nobreak\
\mathrm{eV}$ to $10^{-29}\leavevmode\nobreak\ \mathrm{eV}$. In particular, the
marginalized posteriors over the standard $\Lambda$CDM parameters and
$m_{e,i}$ do not change to a discernible degree for any choice of dataset
combination. Constraints on ${F_{\phi}}$ are slightly more sensitive to
$m_{\phi}$, but in a manner entirely consistent with constraints on ultralight
axions of Ref. [139]. For $m_{\phi}/\mathrm{eV}=10^{-29}$, $10^{-30}$, and
$10^{-31}$, we obtain $95\%$ upper limits of $10^{2}\,{F_{\phi}}=1.01$,
$1.12$, and $1.79$, respectively, from Planck 2018 combined with the SDSS+ BAO
dataset. These results again indicate that the constraints on ${F_{\phi}}$ are
largely driven by effects that are not significantly correlated with any other
cosmological parameter.444 Ref. [139] shows that constraints from Planck 2018
data improve upon the 2015 data release due to an improved measurement of
$\tau_{\mathrm{reion}}$, which breaks a degeneracy between
$A_{s}e^{-2\tau_{\mathrm{reion}}}$ and ${F_{\phi}}$ that arises from a
hyperlight scalar’s suppression of large-scale power.  Given that the scalar’s
gravitational effect is predominantly on large scales, we have no reason to
suspect different findings for scalars coupled instead to the photon.

These results illustrate the power of the CMB to discriminate between models
of small subcomponents of dark matter, even independent of their effect on the
background cosmology. Though hyperlight scalars are not a successful
realization of late dark matter, other possibilities with different dynamics
could well be. In fact, the Standard Model provides one such candidate:
massive neutrinos, which we consider now.

####  IV.2.3 Comparison to massive neutrinos

Figure 16:  Comparison of posterior distributions in scenarios varying the
electron mass and the abundance of a late dark matter component from either
massive neutrinos (blue) or a new, hyperlight scalar with mass
$10^{-30}\leavevmode\nobreak\ \mathrm{eV}$ (red). Results for $\Lambda$CDM
with a varying neutrino mass are depicted in grey for comparison. Each
scenario takes a single massive neutrino state whose mass either is a free
parameter or, in the case that includes a scalar, is fixed to
$0.06\leavevmode\nobreak\ \mathrm{eV}$. The posterior is depicted over the
Hubble constant $h$, the early-time electron mass $m_{e,i}/m_{e,0}$, the net
abundance of late dark matter relative to baryons and CDM
$(\omega_{\nu}+\omega_{\phi})/(\omega_{c}+\omega_{b})$, and the scale factor
of matter–dark-energy equality $a_{m-\Lambda}$. All results use the likelihood
for the full set of Planck 2018 data combined with DESI DR1 BAO data. The 1
and $2\sigma$ contours of the marginalized two-dimensional posteriors are
depicted as in Fig. 14 and the diagonal panels depict the one-dimensional
marginal posteriors, each normalized to unity.

Having established that hyperlight scalars modify structure growth too
severely to make a sizeable contribution to the late-time matter abundance, we
consider another realization of late dark matter that is contained within the
SM: massive neutrinos. One could simultaneously vary the neutrino masses and
include a hyperlight scalar as a microphysical realization of varying
constants, but the results are not materially different. As discussed in
Section III.1, within $\Lambda$CDM massive neutrinos increase the matter
density after recombination, requiring a reduced dark energy density (and so a
smaller $h$) to offset the decrease in distance to last scattering (just as
for a hyperlight scalar). A BAO measurement of $hr_{\mathrm{d}}$ (see Section
II.2.2) breaks this degeneracy, yielding the strongest current constraints on
neutrino masses from cosmological data [75, 156]. DESI’s first data release,
already reaching the level precision of prior surveys, yields stronger
constraints on $\sum m_{\nu}$ that are driven by its preference for
$hr_{\mathrm{d}}$ several percent (or $\sim 1\sigma$) larger relative to
SDSS’s [156].555 Part of this shift is due to DESI LRG datapoints, which are
in up to a $3\sigma$ tension with SDSS in some redshift bins [156]; future
data releases may clarify whether this discrepancy is statistical.

Since late dark matter allows for early recombination without altering the
relative amount of dark energy and matter at late times, $h$ and $\sum
m_{\nu}$ are instead positively correlated. Figure 16 compares results in
varying-$m_{e}$ scenarios featuring late dark matter in the form of a massive
neutrino or a hyperlight scalar. We sample over the neutrino mass with a
uniform prior between $0$ and $1\leavevmode\nobreak\ \mathrm{eV}$. The tail of
the posterior density over
$(\omega_{\nu}+\omega_{\phi})/(\omega_{c}+\omega_{b})$, i.e., the increase in
the matter density at late times, demonstrates the Planck tolerates a
marginally larger abundance of late dark matter in the form of neutrinos than
in a hyperlight scalar. The posteriors for the hyperlight scalar scenario are
centered at slightly larger $h$ and $m_{e,i}$ only because this case takes a
fiducial neutrino mass $m_{\nu}=0.06\leavevmode\nobreak\ \mathrm{eV}$, which
effectively imposes a prior
$(\omega_{\nu}+\omega_{\phi})/(\omega_{c}+\omega_{b})>5\times 10^{-3}$.
Because (uncalibrated) BAO distances depend only on the background cosmology
at late times, the marginal posteriors over $a_{m-\Lambda}$ are effectively
identical. Ultimately, hyperlight scalars and massive neutrinos are comparably
constrained as late dark matter components because their dynamics are
qualitatively similar (as discussed in Section III).

The $95$th percentile of the posterior over $\sum m_{\nu}$ increases from
$0.073\leavevmode\nobreak\ \mathrm{eV}$ in $\Lambda$CDM to
$0.239\leavevmode\nobreak\ \mathrm{eV}$ when varying the early-time electron
mass, similar to that found in Ref. [23]. The latter result precisely matches
the constraint within $\Lambda$CDM from Planck 2018 data alone (i.e., without
BAO data) [75], which again establishes that the primary constraints on early
recombination with late dark matter derive from the new component’s effects on
structure and not those on the expansion history. This result also highlights
the sensitivity of cosmological inference of the neutrino masses to the
underlying cosmological model: the $\Lambda$CDM result nearly rules out the
inverted mass hierarchy, in early recombination scenarios cosmological data
provide no discriminatory power.

###  IV.3 Implications for cosmological concordance

The findings of the previous sections for varying-$m_{e}$ cosmologies—the
impact of consistently including a scalar field as a microphysical realization
of the scenario and the lack of strong concordance among datasets—have
important implications for the interpretation of contemporary tensions between
cosmological datasets. Prior literature pointed to the degeneracy of the
varying-$m_{e}$ model as a preferred means to infer a larger Hubble constant
from CMB data [20, 21, 22, 23] that could be (relatively more) consistent with
that measured via the distance ladder [24, 25, 26]. The results of Fig. 8,
however, suggest that the CMB by itself prefers negligible dark energy in
models where its abundance is not predominantly fixed by the distance to last
scattering, because such models have the additional freedom to better fit the
deficit in temperature power on large scales and the excess lensing of the
acoustic peaks (relative to $\Lambda$CDM) evident in Planck 2018 data.

On its own, Planck 2018 data therefore prefers late recombination: a lighter
early-time electron mass and an even lower value of $h$ than in $\Lambda$CDM,
regardless of whether one accounts for the gravitational effect of a scalar
field coupled to the electron. This finding, though evident in prior work [20,
22] (and to a slightly lesser extent in results using previous Planck data
releases [59, 60]), has not been emphasized in discussions of the model as a
solution to the Hubble tension, nor has its physical origin been explored.
Though the broad posteriors superficially appear to reduce the tension, in no
region of parameter space does the model _simultaneously_ improve the fit to
both CMB data and the calibrated distance-ladder datasets that prefer larger
$h$.666 Ref. [178] succinctly makes the general point that a proper Bayesian
interpretation of a model’s compatibility with multiple datasets requires
marginalizing likelihoods over parameter space—i.e., metrics beyond the
improvement to the best fit point (a point in parameter space with vanishing
posterior mass) or the number of standard deviations between inferred
parameter values (which provides little information on whether a model better
explains different datasets in the same part of parameter space).  Setting
aside questions of whether the lensing and low-$\ell$ anomalies have a
statistical or systematic origin, we conclude that a varying electron mass
alone does not meaningfully restore concordance of these datasets.

Planck’s preference for late recombination derives from features that can be
better explained by delaying the transition to dark-energy domination, but SNe
and BAO distances directly measure the shape of the expansion history (i.e.,
$a_{m-\Lambda}$) in this epoch. Figure 15 shows that even these direct probes
(each combined with Planck and without external calibration for SNe) are far
from concordant within the varying-$m_{e}$ model. Because $h$ and $m_{e,i}$
correlate so strongly with $a_{m-\Lambda}$, their marginal posteriors only
agree to the extent of the datasets’ individual preferences for
$a_{m-\Lambda}$.777 Constraints on the flat $\Lambda$CDM model are usually
quoted in terms of $\Omega_{m}$. We center our discussion on
$a_{m-\Lambda}=\sqrt[3]{\Omega_{m}/(1-\Omega_{m})}$ both because it more
directly connects to the shape of the late-time expansion history and because
a late dark matter component (like a hyperlight scalar) changes the late-time
matter abundance but not the early-time one (Section III.1).

Figure 17:  Posterior distributions over $h$, $m_{e,i}/m_{e,0}$, and
$a_{m-\Lambda}$ for the electron-coupled scalar model (with
$m_{\phi}=10^{-30}\leavevmode\nobreak\ \mathrm{eV}$), using the likelihood for
the full set of Planck 2018 data separately combined with the Pantheon 2018
(grey), Pantheon+ (blue), DES 5YR (green), SDSS+ BAO (red), and DESI DR1 BAO
(gold) datasets. Titles report the corresponding median and $\pm 1\sigma$
quantiles of each distribution. Note that results for the varying-$m_{e}$
model without a scalar field are nearly identical, with posteriors over $h$
and $m_{e,i}$ shifted marginally to lower values.

Results inferred from current SNe and BAO datasets are particularly
discrepant, as evident in Fig. 17. These divergent preferences are especially
pronounced for the very recent Union3 and DES 5YR SNe datasets, which both
prefer a transition to dark-energy domination even later than does the
Pantheon+ result, yielding correspondingly smaller $h$ and $m_{e,i}$.
Curiously, the posteriors from the 2018 Pantheon dataset in Fig. 17, though
much broader than the others, agree much better with Planck combined with
either BAO dataset.

These discrepancies in late time datasets yield highly divergent
interpretations of concordance. At one end, the naïve Gaussian tension between
Planck combined with DES ($h\approx 0.632\pm 0.018$) remains above the
$5\sigma$ level with Pantheon+ and SH0ES ($h\approx 0.735\pm 0.01$) [96],
despite the broader posteriors over $h$ afforded by early recombination
scenarios. In strong contrast, Planck combined with DESI yields $h\approx
0.704\pm 0.011$, apparently a mere $2\sigma$ discrepancy with Pantheon+ and
SH0ES. We reiterate that a substantial contribution to DESI’s preference for
larger $h$ is driven by its LRG data, which within $\Lambda$CDM are in slight
tension with SDSS’s [156]. Finally, because $h$ varies rather rapidly with the
electron mass at fixed $\theta_{s}$, the discrepancy in measurements of
$m_{e,i}$ in Fig. 17 is less severe but still impedes a quantitative
identification of a consensus constraint from cosmological datasets on
$m_{e,i}$.

Setting aside the more discrepant datasets (Fig. 17), even the improved
agreement on $h$ in early recombination scenarios between the Planck and DESI
combination and the SH0ES-calibrated distance ladder does not fully capture
the degree of the tension. These two dataset combinations are also in a
$2\sigma$ tension in their (marginalized) posteriors over $a_{m-\Lambda}$, as
seen in Fig. 17. Moreover, $h$ and $a_{m-\Lambda}$ are anticorrelated in early
recombination models as constrained by Planck, such that decreasing the
discrepancy in one of $h$ or $a_{m-\Lambda}$ only exacerbates that in the
other. Concordance must be assessed in the full two-dimensional parameter
space that characterizes the late-Universe expansion history; as shown in Fig.
18, the posteriors in the $h$-$a_{m-\Lambda}$ plane hardly overlap even in
their $3\sigma$ mass levels.

Figure 18:  Posterior distribution over $h$ and $a_{m-\Lambda}$ from Planck
2018 and DESI DR1 BAO data for the electron-coupled scalar model (gold) and
from Pantheon+ with SH0ES-calibrated Cepheids (blue). The vertical lines in
the diagonal panels depict the median and $2.5$th, $16$th, $84$th, and
$97.5$th percentiles of the marginalized posterior distributions. These one-
dimensional posteriors are kernel density estimates normalized relative to
their peak value to facilitate comparison. The median and corresponding $\pm
1\sigma$ uncertainties for each parameter are reported above the diagonal
panels. The lower panels display the the 1, 2, and $3\sigma$ contours (i.e.,
the $39.3\%$, $86.5\%$, and $98.9\%$ mass levels) of the marginalized, two-
dimensional posterior.

While a late dark matter contribution could increase $h$ without reducing
$a_{m-\Lambda}$ in principle (Section III.1), within the hyperlight scalar
field model Planck does not allow abundances large enough to do so. Even in a
scenario where this additional degeneracy is realized, the discrepant
inference of $a_{m-\Lambda}$ by BAO and SNe data (which trace the same
interval of redshift) cannot be fully reconciled—a problem for any model in
which the late-time expansion history is well described by flat $\Lambda$CDM.

For these reasons, we do not analyze the combination of Planck 2018 data with
the SH0ES-calibrated distance ladder. Posteriors for these combined datasets
would localize at parameters between their individual preferences, obfuscating
the fact that the individual likelihoods are degraded in this resulting
parameter space. Prior work has in addition admonished the usage of a SH0ES-
derived prior directly on $H_{0}$, emphasizing the importance of including the
full magnitude-redshift sample from Pantheon SNe, calibrated by SH0ES
Cepheids, to provide genuine constraints on modifications to late-time
cosmology [179, 180]. Our results demonstrate the importance of doing so more
generally, because modifications to early-time cosmology can correlate to
changes in the late-time expansion history. In early recombination models, the
disproportionate increase in the matter and dark-energy densities cause the
late-time expansion history to depart (in shape, or $a_{m-\Lambda}$) from that
directly measured by current SN datasets, underscoring the potential hazards
in interpreting analyses that combine the Cepheid-calibrated distance ladder
with CMB (and BAO) data.

As anticipated in Section II.1.6, our results bear a close resemblance to
those for the (counterfactual) cosmologies in which the present-day CMB
temperature $T_{0}$ is treated as an otherwise unconstrained parameter [38,
181]. The correlated impact on the late-time expansion history is common to
both scenarios, as is the role of CMB lensing and the ISW effect [38, 182].
The preceding considerations also apply (to a reasonable extent) to early
recombination mechanized by small-scale inhomogeneities in the baryon
distribution [78, 79, 80, 81] mentioned in Section II.1.6. These scenarios
only prefer a larger Hubble constant when combining Planck data with datasets
like SH0ES [78], indicating that (superficially) alleviating the tension comes
at the cost of degrading the Planck likelihoods.

Having assessed the joint effect of varying constants and hyperlight scalars
on the background evolution in the $h$-$a_{m-\Lambda}$ space, we now turn to
the model’s impact on matter clustering and the Hubble constant. Low-redshift
probes prefer larger $h$ but smaller $\sigma_{8}$ (or
$S_{8}\equiv\sqrt{\Omega_{m}/0.3})$, but a common feature of proposed
$\Lambda$CDM extensions is a positive correlation between, e.g., $h$ and
$\sigma_{8}$ [183, 184]. This expectation may be understood on general grounds
via the arguments of Section II: the radiation driving and early ISW effects
constrain the CDM abundance at recombination, requiring a larger matter
abundance today if recombination occurs at higher density, whether because it
occurs early or in the presence of additional new degrees of freedom. Absent
other effects, early-time modifications to cosmology that lead to inference of
larger $h$ typically also infer a larger amplitude of matter clustering.

Figure 19:  Posterior distribution over $h$ and $\sigma_{8}$ (the latter a
derived quantity) in scenarios vary $m_{e,i}/m_{e,0}$ (left) or
$\alpha_{i}/\alpha_{0}$ (right). Both cases also compare results with
$\Lambda$CDM with and without a scalar. Scenarios with a scalar field take its
mass to be $10^{-30}\leavevmode\nobreak\ \mathrm{eV}$, and all cases use
Planck 2018 likelihoods and the SDSS+ BAO datasets enumerated in the text.

On the other hand, per Section III.2 (and, e.g., Ref. [139]), a hyperlight
scalar subcomponent of dark matter in general suppresses matter clustering.
The correlation between the abundance of the scalar (${F_{\phi}}$) and the
early-time shift in the electron mass ($m_{e,i}/m_{e,0}-1$) outlined in
Section III.1 would achieve an anticorrelation between $h$ and $\sigma_{8}$ if
it had borne out in the full data analysis. Indeed, Fig. 19 evinces a
broadening of the posteriors (projected into the space of $h$ and
$\sigma_{8}$) in this direction for varying-$m_{e}$ cosmologies that include a
scalar compared to those that do not. However, Planck data itself ultimately
curtails the extent of this effect in parameter space; the results here
ultimately do not depart dramatically from those of Ref. [139], which
considered hyperlight scalars with only gravitational interactions. The degree
of the effect is comparable to the extent hyperlight scalars may reduce
$\sigma_{8}$ in standard $\Lambda$CDM, except that (as explained in Section
IV.2) ${F_{\phi}}$ and $h$ are themselves anticorrelated in this case. Coupled
hyperlight scalars therefore provide at the least an example model that can
accommodate CMB data at larger $h$ but smaller $\sigma_{8}$.

##  V Conclusions

The cosmic microwave background, long lauded as a precision probe of the early
Universe, depends sensitively on the physics of quantum electrodynamics. In
this work we consider microphysical realizations of scenarios in which the
fundamental parameters of QED—the electromagnetic fine-structure constant and
electron mass—take on different values during recombination than at the
present day. Companion work [32] discusses theories of new scalar fields whose
couplings vary fundamental constants in spacetime, emphasizing that
consistency with theoretical considerations and independent probes restricts
the viable space of models. Crucially, Ref. [32] shows that mechanizing
shifted early-time fundamental constants that only start evolving toward their
present-day values after recombination _requires_ the scalar field to have a
nonnegligible abundance. Here we investigate the cosmological impact of a new
scalar through both its direct coupling and its gravitational coupling to SM
matter.

Variations in the fine-structure constant and the electron mass, regardless of
their origin, primarily affect the release of the CMB via their well-studied
impacts on recombination, shifting the temperature and duration of last
scattering, and on Thomson scattering. These effects are not independently
constrained by data in $\Lambda$CDM cosmology, however; in Section II we
clarify the physical origin of degeneracies between the fundamental constants
and $\Lambda$CDM parameters. In particular, we enumerate how CMB anisotropies
largely probe the state of the Universe at last scattering, whenever it may
have occurred, and how the rate of recombination—namely, the redshift interval
over which visible CMB photons last scattered—impacts the generation of
polarization and the suppression of small-scale anisotropies. Sections II.1
and III.1 describe the combinations of standard parameters that directly
encode the most important physical effects that imprint in CMB anisotropies,
generalizing those of prior work [185, 77, 186]. In doing so, we uncover the
connection between modified recombination scenarios and the CMB monopole
temperature, if treated as an otherwise unmeasured parameter. In particular,
when either of the present-day CMB temperature and the early-time electron
mass is not fixed, the Hubble constant is no longer uniquely constrained to
fix the angular extent of the sound horizon precisely measured by CMB
anisotropies—the same “geometric” degeneracy characteristic to cosmologies
with nonzero spatial curvature [185].

The correlation of fundamental constants with $\Lambda$CDM parameters
inevitably affects cosmological dynamics at late times, to which the CMB
anisotropies are only weakly sensitive. Low-redshift distance measurements,
whether from type Ia supernovae or baryon acoustic oscillations, therefore
play a crucial role in breaking degeneracies by directly measuring the shape
of the late-time expansion history, even without calibrating their absolute
distances. When including low-redshift distance datasets, cosmological data
constrain electron mass variations at the percent level and fine-structure
constant variations at nearly the permille level. On the other hand, Section
III.1 identifies a novel solution that combines early recombination with “late
dark matter,” i.e., an enhancement of the matter abundance after
recombination. The SM already offers a source of late dark matter: massive
neutrinos, which become nonrelativistic well before the present but after
recombination. Moreover, Section III shows that a hyperlight scalar field,
invoked to mechanize variation of fundamental constants between recombination
and the present day, must begin oscillating and boost the matter abundance
after recombination as well. The positive correlation between $h$ and the
abundance of late dark matter (measured by the sum of the neutrino masses or
by ${F_{\phi}}$ for a scalar) along this degeneracy is especially notable in
contrast to the _anti_ correlation between the two in $\Lambda$CDM.

Though early recombination with late dark matter realizes a degeneracy in the
cosmological background, the scenarios are distinguished (and therefore
testable) through the behavior of the new component’s perturbations. Section
IV carefully studies the complementary impact of individual datasets in
constraining these scenarios and provides quantitative evidence for the
physical effects outlined in Sections II and III. Section IV.2 further
illustrates how the perturbations of a hyperlight scalar suppress the growth
of structure and source an integrated Sachs-Wolfe effect, signatures which
ultimately restrict the scalar abundance allowed by Planck data to no more
than $\mathcal{O}(1\%)$ of the dark matter’s.

A heavier electron at early times, which triggers recombination at a higher
temperature, has been touted as an extension to $\Lambda$CDM for which
cosmological data infer a Hubble constant $H_{0}$ in better agreement with
distance-ladder measurements. In Section IV.3 we critically assess this claim,
arguing that because the physical effects of Section II lead Planck data to
prefer late recombination to early recombination, a superficial reduction in
the $H_{0}$ tension comes at the cost of degrading the fit to individual
cosmological datasets. Moreover, though the degeneracy introduced by the
scalar’s contribution to the late-time matter abundance (identified in Section
III) offers the potential to simultaneously satisfy CMB and low-redshift
observations at arbitrarily large $H_{0}$, Planck data too strongly disfavors
the impact of a hyperlight scalar’s perturbations to fully realize this
effect. To the extent that Planck data _does_ allow for hyperlight scalars,
however, slightly larger electron masses and Hubble constants are allowed
compared to scenarios without a scalar. We show that similar conclusions apply
to massive neutrinos as an alternative realization of late dark matter. An
interesting question for further study is whether any motivated models of late
dark matter affect the dynamics of perturbations in a manner the CMB is more
amenable to.

Because low-redshift datasets play such a crucial role in breaking
degeneracies in varying electron mass scenarios, any discord between them is
put on full display. Sections IV.1.2 and IV.2.2 show that the diverging
preferences in the uncalibrated late-time expansion history (i.e., its shape
as encoded by the scale factor of matter–dark-energy equality $a_{m-\Lambda}$)
propagate directly to the electron mass and $h$. In particular, Section IV.3
highlights an increasing discrepancy in the cosmologies preferred by recent
BAO datasets and recent SNe datasets that, combined with Planck, vary in their
inferred $h$ between $0.704\pm 0.011$ and $0.632\pm 0.018$. The former result
derives from the recent DESI DR1 BAO measurements [156, 157, 158] and is
larger than that from prior BAO data, a shift driven in part by an increase
its inferred the dimensionless amplitude $hr_{\mathrm{d}}$. Section IV.3 also
emphasizes the importance of a holistic interpretation of the agreement
between posteriors derived from, e.g., CMB and BAO data with those from the
SH0ES-calibrated distance ladder: even with the additional freedom afforded to
early-recombination (and varying-$T_{0}$) cosmologies, no independent dataset
combination prefers a late-time expansion history that agrees with the SH0ES-
calibrated distance ladder in both calibration ($h$) and shape
($a_{m-\Lambda}$).

Future high-resolution CMB observations [48, 187] will greatly improve
measurements of temperature and polarization anisotropies deep into the
damping tail, promising substantially stronger constraints on fine-structure
constant variations. Even current and near-term datasets from ground-based
experiments [170, 171, 172] offer increased constraining potential. While no
physical effect positions high-resolution observations to better constrain
electron-mass variations, alternative CMB datasets can weigh in on the
influence of features in Planck data suspected to be systematic in origin. The
lensing anomaly in particular is absent in other datasets, including recent
reanalyses of Planck data [188, 189]. Bounds on the abundance of hyperlight
scalars (coupled to the SM or not) could too be alleviated because the same
systematic effects effectively disfavor the their impact on structure growth,
as has already been observed for bounds on the neutrino masses [189]. Future
BAO and SNe datasets will also provide substantially more precise measurements
of the late-time expansion history, better breaking degeneracies in CMB
constraints on early recombination. We will study constraints from these more
recent datasets and forecast bounds from upcoming observations in future work.

In a companion article [32], we apply the results of Section IV.2 to constrain
concrete extensions of the SM with new scalar fields coupled to the electron
and photon. Ref. [32] discusses numerous challenges in constructing models
that are both theoretically consistent and allowed by independent constraints
from other astrophysical and laboratory measurements. In viable theories, the
cosmological constraints of Section IV.2 in particular yield leading bounds on
hyperlight, quadratically coupled scalars that become matterlike in the matter
dominated era.

This work restricts its analysis to regimes that realize the phenomenological
scenarios of prior study, in which the change to the fundamental constants
remains fixed through recombination. The theoretical developments of Ref.
[32], however, motivate scenarios in which the scalar (and so the fundamental
constants) evolve before recombination under the influence of the cosmological
abundance of SM matter. Such a regime resembles attempts to address
coincidence problems in early dark energy via couplings to dark matter [190,
191] or neutrinos [192, 193] to explain why the early dark energy field
becomes dynamical around the time of recombination. Moreover, a complete
treatment of coupled scalars requires an extension to heavier masses for which
the scalar begins oscillating before recombination of its own volition. Time-
dependent fundamental constants around recombination may yield qualitatively
different CMB signatures via, for instance, the interplay of changes to the
shape of the visibility function and diffusion damping rate at earlier times.

In addition, scalars may equally well couple to other SM content like the
Higgs, quarks, and gluons, which would affect the masses of baryons and
mediate fifth forces between them. Fifth forces between electrons are present
in the scenarios we considered here, but their effect is negligible because
inhomogeneities in the the electron fluid are themselves negligible in linear
perturbation theory (see Section III). We leave a general study of the
cosmological signatures of coupled scalars, including the effect of evolution
prior to recombination and scalar-mediated forces on the CMB and large-scale
structure, to future work.

###### Acknowledgements.

We thank Nikita Blinov, David Cyncynates, Junwu Huang, Mikhail Ivanov, Hongwan
Liu, Cristina Mondino, Caio Nascimento, Maxim Pospelov, Murali Saravanan,
Sergey Sibiryakov, Neal Weiner, and Tien-Tien Yu for helpful conversations and
especially Marilena Loverde for many useful and extended discussions. M.B. is
supported by the U.S. Department of Energy Office of Science under Award
Number DE-SC0024375. M.B. and Z.J.W. are supported by the Department of
Physics and College of Arts and Science at the University of Washington. O.S.
was supported by the Department of Physics and a DARE Fellowship from the
Office of the Vice Provost for Graduate Education at Stanford University. This
work made use of the software packages emcee [194, 195, 196], corner.py [197],
NumPy [198], SciPy [199], matplotlib [200], xarray [201], ArviZ [202], SymPy
[203], and CMasher [204].

##  Appendix A Numerical implementation of cosmological scalar fields

The gravitational effects of cosmological scalar fields have been studied in
contexts including ultralight dark matter [205, 134, 206, 207, 128, 140, 141,
142, 139, 208], scalar field dark energy [209, 210], and early dark energy
[129, 130, 131, 133]. A scalar field’s potential typically introduces a new
timescale of cosmological importance, which in the simplest case of potentials
with minima is the scalar’s oscillation period. Scalars with quadratic
potentials generically oscillate with a fixed period much shorter than the age
of the Universe (unless their mass is close to the Hubble constant $H_{0}$, a
regime to which we do not restrict ourselves).

To mitigate the onerous computational cost required to accurately solve for
scalar fields’ dynamics, a variety of effective fluid treatments for scalar
fields have been utilized in prior literature [205, 134, 206, 128, 211, 130,
212, 131]. We follow Ref. [121], which systematically developed a fluid
approximation for massive scalar fields that enabled semianalytically and
empirically calibrating the effective equation of state and sound speed
describing the fluid. The scheme devised by Ref. [121] entails solving the
scalar’s equation of motion (the Klein-Gordon equation) directly for its first
few oscillations (when fluid approximations are largely insufficient) before
switching to the effective fluid approximation. The matching procedure of Ref.
[121] also improves substantially upon prior methods. We extend CLASS [51, 52]
by implementing the scheme proposed by Ref. [121] to enable an efficient and
accurate treatment of scalar fields. (The implementation is available at Ref.
[122].) For completeness, we briefly review the formalism and definitions
necessary to specify the scheme, referring to Ref. [121] for full details and
comparison with previous fluid treatments.

###  A.1 Dynamics

In a general spacetime, the action for a minimally coupled scalar field is
$S_{\phi}=\int\mathrm{d}^{4}x\,\sqrt{-g}\mathcal{L}_{\phi}$ in terms of its
Lagrangian density

$\displaystyle\mathcal{L}_{\phi}$
$\displaystyle=-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi).$ (58)

The scalar’s Euler-Lagrange equation is

$\displaystyle\nabla_{\mu}\nabla^{\mu}\phi$
$\displaystyle=\frac{\mathrm{d}V}{\mathrm{d}\phi},$ (59)

where $\nabla_{\mu}$ is the covariant derivative compatible with the metric
$g$. We parameterize a general, perturbed, conformal-time Friedmann-Lemaître-
Robertson-Walker (FLRW) metric with

$\displaystyle g_{\mu\nu}\equiv
a(\tau)^{2}\left(\eta_{\mu\nu}+h_{\mu\nu}\right)$ (60)

where $\eta_{\mu\nu}$ is the Minkowski metric with the mostly positive
signature and $h_{\mu\nu}$ a small perturbation we decompose as

$\displaystyle h_{00}$ $\displaystyle=-E$ (61a) $\displaystyle h_{i0}$
$\displaystyle=\partial_{i}F$ (61b) $\displaystyle h_{ij}$
$\displaystyle=A\delta_{ij}+\partial_{i}\partial_{j}B.$ (61c)

Because neither primordial vectors nor tensors have been detected, and because
we introduce no new physics that would source them, Eq. 61 includes only
scalar perturbations $A$, $B$, $E$, and $F$. This parameterization is general
but redundant due to gauge invariance; we discuss concrete gauge choices
below.

Expanding to linear order with
$\phi(t,\mathbf{x})=\bar{\phi}(t)+\delta\phi(t,\mathbf{x})$, the equation of
motion Eq. 59 becomes

$\displaystyle 0$
$\displaystyle=\bar{\phi}^{\prime\prime}+2\mathcal{H}\bar{\phi}^{\prime}+a^{2}\frac{\mathrm{d}V}{\mathrm{d}\phi}$
(62a)
$\displaystyle\begin{split}0&=\delta\phi^{\prime\prime}+2\mathcal{H}\delta\phi^{\prime}-\partial_{i}\partial_{i}\delta\phi+a^{2}\frac{\mathrm{d}^{2}V}{\mathrm{d}\phi^{2}}\delta\phi+\frac{\bar{\phi}^{\prime}}{2}\left(\partial_{i}\partial_{i}\left[B^{\prime}-2F\right]+3A^{\prime}-E^{\prime}\right)+a^{2}\frac{\mathrm{d}V}{\mathrm{d}\phi}E.\end{split}$
(62b)

Here primes denote derivatives with respect to conformal time $\tau$. The
conformal Newtonian gauge has $E=2\Psi$, $A=-2\Phi$, and both $F$ and $B$
zero, in which case

$\displaystyle 0$
$\displaystyle=\delta\phi^{\prime\prime}+2\mathcal{H}\delta\phi^{\prime}-\partial_{i}\partial_{i}\delta\phi+a^{2}\frac{\mathrm{d}^{2}V}{\mathrm{d}\phi^{2}}\delta\phi-\left(3\Phi^{\prime}+\Psi^{\prime}\right)\bar{\phi}^{\prime}+2a^{2}\frac{\mathrm{d}V}{\mathrm{d}\phi}\Psi.$
(63)

The synchronous gauge used in Ref. [213] sets $A=-2\eta$ and
$\partial_{i}\partial_{i}B=h+6\eta$ with $E$ and $F$ zero; Eq. 62b then reads

$\displaystyle 0$
$\displaystyle=\delta\phi^{\prime\prime}+2\mathcal{H}\delta\phi^{\prime}-\partial_{i}\partial_{i}\delta\phi+a^{2}\frac{\mathrm{d}^{2}V}{\mathrm{d}\phi^{2}}\delta\phi+\frac{1}{2}\bar{\phi}^{\prime}h^{\prime}.$
(64)

The scalar field’s contribution to the stress-energy tensor is

$\displaystyle\left(T_{\mu\nu}\right)^{\phi}=-2\frac{\partial\mathcal{L}_{\phi}}{\partial
g^{\mu\nu}}+g_{\mu\nu}\mathcal{L}_{\phi}$
$\displaystyle=\partial_{\mu}\phi\partial_{\nu}\phi+g_{\mu\nu}\left(-\frac{1}{2}\partial_{\alpha}\phi\partial^{\alpha}\phi-V(\phi)\right).$
(65)

At the background level, the effective energy density
$\bar{\rho}_{\phi}=-\left(\bar{T}_{0}^{\hphantom{0}0}\right)^{\phi}$ and
pressure $\bar{P}_{\phi}=\left(\bar{T}_{i}^{\hphantom{i}i}\right)^{\phi}/3$
are

$\displaystyle\bar{\rho}_{\phi}$
$\displaystyle=\frac{\left(\bar{\phi}^{\prime}\right)^{2}}{2a^{2}}+V(\bar{\phi})$
(66a) $\displaystyle\bar{P}_{\phi}$
$\displaystyle=\frac{\left(\bar{\phi}^{\prime}\right)^{2}}{2a^{2}}-V(\bar{\phi}).$
(66b)

We parameterize the scalar perturbations to the stress-energy tensor in terms
of density, pressure, and velocity perturbations $\delta\rho$, $\delta P$, and
$\delta u$, as well as anisotropic stress $\pi^{S}$:

$\displaystyle\delta T^{0}_{\hphantom{0}0}$ $\displaystyle=-\delta\rho$ (67a)
$\displaystyle\delta T^{0}_{\hphantom{0}i}$
$\displaystyle=\left(\bar{\rho}+\bar{P}\right)\partial_{i}\delta u$ (67b)
$\displaystyle\delta T^{i}_{\hphantom{i}j}$ $\displaystyle=\delta_{ij}\delta
P+\left(\partial_{i}\partial_{j}-\frac{1}{3}\delta_{ij}\partial_{k}\partial_{k}\right)\pi^{S}.$
(67c)

Often the velocity perturbation is written in terms of
$\theta=\partial_{i}\partial_{i}\delta u$. The corresponding contributions
from Eq. 65 (at linear order) are

$\displaystyle\delta\rho_{\phi}$
$\displaystyle=-\frac{E}{2a^{2}}\left(\bar{\phi}^{\prime}\right)^{2}+\frac{\bar{\phi}^{\prime}\delta\phi^{\prime}}{a^{2}}+\frac{\mathrm{d}V}{\mathrm{d}\phi}\delta\phi$
(68a) $\displaystyle\left(\bar{\rho}_{\phi}+\bar{P}_{\phi}\right)\delta
u_{\phi}$ $\displaystyle=-\frac{\bar{\phi}^{\prime}\delta\phi}{a^{2}}$ (68b)
$\displaystyle\delta P_{\phi}$
$\displaystyle=-\frac{E}{2a^{2}}\left(\bar{\phi}^{\prime}\right)^{2}+\frac{\bar{\phi}^{\prime}\delta\phi^{\prime}}{a^{2}}-\frac{\mathrm{d}V}{\mathrm{d}\phi}\delta\phi.$
(68c)

Anisotropic stress from a scalar field is zero to leading order in
perturbation theory. Recall again that $E=2\Psi$ in Newtonian gauge and $0$ in
synchronous gauge.

###  A.2 Effective fluid treatment

We now turn to effective fluid approximations for scalars with purely
quadratic potentials, $V(\phi)=m_{\phi}^{2}\phi^{2}/2$. Following Ref. [121],
we fix synchronous gauge for the remainder of the discussion and switch to
cosmic time $t$ defined by $\mathrm{d}t=a\mathrm{d}\tau$. With these choices,

$\displaystyle 0$
$\displaystyle=\ddot{\bar{\phi}}+3H\dot{\bar{\phi}}+m_{\phi}^{2}\bar{\phi}$
(69)
$\displaystyle\begin{split}0&=\delta\ddot{\phi}+3H\delta\dot{\phi}+\frac{k^{2}}{a^{2}}\delta\phi+m_{\phi}^{2}\delta\phi+\frac{1}{2}\dot{\bar{\phi}}\dot{h},\end{split}$
(70)

writing the perturbation equation in Fourier space. In the absence of metric
perturbations, these equations both exhibit oscillatory solutions with a
decaying envelope.

####  A.2.1 Background

Beginning at the background level Ref. [121] decomposes the solution to Eq. 69
onto its two oscillatory modes as

$\displaystyle\bar{\phi}(t)$
$\displaystyle=\bar{\varphi}_{c}(t)\cos\left(m_{\phi}t-m_{\phi}t_{\star}\right)+\bar{\varphi}_{s}(t)\sin\left(m_{\phi}t-m_{\phi}t_{\star}\right),$
(71)

where the amplitudes $\bar{\varphi}_{c}$ and $\bar{\varphi}_{s}$ evolve slowly
(i.e., on timescales order $1/H$). Equation 69 expands to

$\displaystyle\begin{split}0&=\left[\ddot{\bar{\varphi}}_{c}+2m_{\phi}\dot{\bar{\varphi}}_{s}+3H\left(\dot{\bar{\varphi}}_{c}+m_{\phi}\bar{\varphi}_{s}\right)\right]\cos\left(m_{\phi}t-m_{\phi}t_{\star}\right)\\\
&\hphantom{{}={}}+\left[\ddot{\bar{\varphi}}_{s}-2m_{\phi}\dot{\bar{\varphi}}_{c}+3H\left(\dot{\bar{\varphi}}_{s}-m_{\phi}\bar{\varphi}_{c}\right)\right]\sin\left(m_{\phi}t-m_{\phi}t_{\star}\right),\end{split}$
(72)

which may be solved by solving the individual expressions in brackets as a set
of two coupled differential equations. In brief, Ref. [121] defines an
effective energy density and pressure

$\displaystyle\bar{\rho}_{\phi,\mathrm{ef}}$
$\displaystyle\equiv\frac{1}{2}\left(\frac{1}{2}\left[\dot{\bar{\varphi}}_{c}^{2}+\dot{\bar{\varphi}}_{s}^{2}\right]+m_{\phi}\left[\dot{\bar{\varphi}}_{c}\bar{\varphi}_{s}-\dot{\bar{\varphi}}_{s}\bar{\varphi}_{c}\right]+m_{\phi}^{2}\left[\bar{\varphi}_{s}^{2}+\bar{\varphi}_{c}^{2}\right]\right)$
(73a) $\displaystyle\bar{P}_{\phi,\mathrm{ef}}$
$\displaystyle\equiv\frac{1}{2}\left(\frac{1}{2}\left[\dot{\bar{\varphi}}_{c}^{2}+\dot{\bar{\varphi}}_{s}^{2}\right]+m_{\phi}\left[\dot{\bar{\varphi}}_{c}\bar{\varphi}_{s}-\dot{\bar{\varphi}}_{s}\bar{\varphi}_{c}\right]\right)$
(73b)

that, via the Klein-Gordon equation Eq. 72, exactly satisfy the standard fluid
conservation law

$\displaystyle\dot{\bar{\rho}}_{\phi,\mathrm{ef}}$
$\displaystyle=-3H\left(\bar{\rho}_{\phi,\mathrm{ef}}+\bar{P}_{\phi,\mathrm{ef}}\right).$
(74)

In the scheme of Ref. [121], one solves the full Klein-Gordon equation Eq. 69
until some transition time $t_{\star}$ and matches onto the amplitudes
$\bar{\varphi}_{c}$ and $\bar{\varphi}_{s}$ in order to evaluate the effective
fluid energy density $\bar{\rho}_{\phi,\mathrm{ef}}$. Matching error (due to
oscillations of $\rho_{\phi}$ about its time-averaged $1/a^{3}$ redshifting)
is minimized by choosing

$\displaystyle\left.\frac{\ddot{\bar{\varphi}}_{I}}{\dot{\bar{\varphi}}_{I}}\right|_{t=t_{\star}}$
$\displaystyle=-\frac{\left\langle
H\right\rangle}{2}\left(3-\frac{2\mathrm{d}\left\langle
H\right\rangle/\mathrm{d}\tau}{a\left\langle
H\right\rangle^{2}}\right)\equiv\mathcal{M}.$ (75)

for $I=c$, $s$. Here $\left\langle H\right\rangle$ denotes the Hubble
parameter averaged over the oscillations induced by the scalar’s contribution
to FLRW expansion; we may neglect such effects in practice and simply
substitute $H$. Inserting the decomposed Klein-Gordon equation into Eq. 75 and
solving yields

$\displaystyle\begin{pmatrix}\dot{\bar{\varphi}}_{c}\\\
\dot{\bar{\varphi}}_{s}\end{pmatrix}$
$\displaystyle=-\frac{3H}{\mathcal{M}^{2}+3H\mathcal{M}+4m_{\phi}^{2}}\begin{pmatrix}2m_{\phi}&3H+\mathcal{M}\\\
-\mathcal{M}&2m_{\phi}\end{pmatrix}\begin{pmatrix}m_{\phi}\bar{\phi}\\\
\dot{\bar{\phi}}\end{pmatrix}$ (76)

and

$\displaystyle\bar{\varphi}_{c}$ $\displaystyle=\bar{\phi}$ (77a)
$\displaystyle\bar{\varphi}_{s}$
$\displaystyle=\dot{\bar{\phi}}/m_{\phi}-\dot{\bar{\varphi}}_{c}/m_{\phi}.$
(77b)

Ref. [121] finds these choices reduce the matching error between
$\bar{\rho}_{\phi}$ and $\bar{\rho}_{\phi,\mathrm{ef}}$ to
$\mathcal{O}([H_{\star}/m_{\phi}]^{3})$, where $H_{\star}=H(t_{\star})$.

After matching, one solves the dynamics for the effective fluid approximation
(efa) to the energy density $\bar{\rho}_{\phi,\mathrm{efa}}$,

$\displaystyle\dot{\bar{\rho}}_{\phi,\mathrm{efa}}$
$\displaystyle=-3H\left(1+w_{\phi,\mathrm{efa}}\right)\bar{\rho}_{\phi,\mathrm{efa}}.$
(78)

At late times, $\bar{\rho}_{\phi}$ redshifts like $a^{-3}$ to good precision,
motivating the choice $w_{\phi,\mathrm{efa}}=0$. However, optimal
computational efficiency switches at as large of $H_{\star}/m_{\phi}$ as
possible, at which point the (average) equation of state and rate of redshift
deviate appreciably from that of cold dark matter. Ref. [121] therefore
calibrated a time-dependent ansatz for
$w_{\phi,\mathrm{efa}}=3(H/m_{\phi})^{2}/2$, informed by computing Eq. 73
using numerical solutions to the Klein-Gordon equation. However, Ref. [121]
focused on scalars that begin oscillating in the radiation era
($m_{\phi}\gtrsim 10^{-28}\leavevmode\nobreak\ \mathrm{eV}$), whereas for the
scenarios relevant here, the scalars necessarily begin to oscillate in the
matter era. Using numerical solutions, Ref. [121] demonstrated that
$w_{\phi,\mathrm{ef}}\approx 3(H/m_{\phi})^{2}/2$ for $H\gg m_{\phi}$, noting
that the result could also be obtained by the analytic solution to the
homogeneous Klein-Gordon equation in the radiation era. Since we are
interested in lighter scalars (for which matching will occur in the matter era
or during the onset of dark-energy domination), we extend this result to FLRW
backgrounds with arbitrary equations of state.888 Generalization to lighter
masses is also of interest for future cosmological constraints on scalar
subcomponents of dark matter with such masses [128, 140, 141, 142, 139], given
projected biases from fluid approximations for future CMB experiments from
Ref. [212] (which considered masses greater than $10^{-27}\leavevmode\nobreak\
\mathrm{eV}$). It would also be interesting to extend the formalism of Ref.
[121] to potentials with minima steeper than quadratic in $\phi$, as relevant
to early dark energy [129, 130, 131, 132, 214]; however, the computational
impetus to do so is likely not nearly so severe because scalars in steeper
potentials oscillate with a frequency that decreases as the field amplitude
redshifts (see, e.g., Ref. [215]).

The Hubble parameter in a Universe with constant equation of state $w$ is
$H(t)=n/t$, where $n=2/3(1+w)$. In such a background, the solution to the
homogeneous Klein-Gordon equation (with constant field value at early times
$t\ll 1/m_{\phi}$) is $\phi(t)\propto J_{\nu}(m_{\phi}t)/t^{\nu}$ where
$\nu=(3n-1)/2$. Taking the asymptotic (i.e., $m_{\phi}t\gg 1$) expansion of
the energy density and pressure for this solution (at next-to-leading order,
to correctly compute the $1/t^{2}$ term) yields
$w_{\phi,\mathrm{ef}}=3n/4(m_{\phi}t)^{2}+\mathcal{O}[(m_{\phi}t)^{-4}]$. In a
multicomponent Universe (like $\Lambda$CDM), the instantaneous equation of
state may be obtained by combining the Friedmann equations to give
$w(t)=-2\dot{H}(t)/3H(t)^{2}-1$. Given that the fluid approximation is meant
to be used only when $m_{\phi}/H\gtrsim 10$, approximating expansion as a
power law at any given time should be wholly sufficient [barring exotic
cosmologies with $w(t)$ evolving on timescales much shorter than the Hubble
scale]. As such, we may promote our single-component result above to a general
Universe via

$\displaystyle w_{\phi,\mathrm{ef}}$
$\displaystyle=-\frac{3}{4}\frac{\dot{H}}{m_{\phi}^{2}}$ (79)

In practice, we find Eq. 79 works exceedingly well. Since $\dot{H}=-2H^{2}$
deep in the radiation era, it reproduces the results of Ref. [121] for
transitions occurring well before matter-radiation equality. Moreover, for
transitions occurring in the matter era, we observe even more rapid
convergence, improving with $(H_{\star}/m_{\phi})^{4}$ for $m_{\phi}\geq
10^{-28}\leavevmode\nobreak\ \mathrm{eV}$ compared to cubic convergence for
transitions in the radiation era.

####  A.2.2 Perturbations

In spite of the added complications that the equation of motion for
$\delta\phi$ depends on wave number $k$ and is sourced by metric
perturbations, the treatment for spatial perturbations proceeds similarly. We
briefly enumerate the corresponding definitions, referring to Ref. [121] for
further discussion and justification. The perturbation is again decomposed
onto oscillatory modes as

$\displaystyle\delta\phi(t)$
$\displaystyle=\delta\varphi_{c}(t)\cos\left(m_{\phi}t-m_{\phi}t_{\star}\right)+\delta\varphi_{s}(t)\sin\left(m_{\phi}t-m_{\phi}t_{\star}\right).$
(80)

The effective fluid variables defined as

$\displaystyle\begin{split}2\delta\rho_{\phi,\mathrm{ef}}&=\dot{\bar{\varphi}}_{c}\delta\dot{\varphi}_{c}+\dot{\bar{\varphi}}_{s}\delta\dot{\varphi}_{s}+m_{\phi}\left(\bar{\varphi}_{s}\delta\dot{\varphi}_{c}-\bar{\varphi}_{c}\delta\dot{\varphi}_{s}\right)+m_{\phi}\left(\dot{\bar{\varphi}}_{c}\delta\varphi_{s}-\dot{\bar{\varphi}}_{s}\delta\varphi_{c}\right)\\\
&\hphantom{{}={}}+2m_{\phi}^{2}\left[\bar{\varphi}_{c}\delta\varphi_{c}+\bar{\varphi}_{s}\delta\varphi_{s}\right]\end{split}$
(81a)
$\displaystyle\left(\bar{\rho}_{\phi,\mathrm{ef}}+\bar{P}_{\phi,\mathrm{ef}}\right)\theta_{\phi,\mathrm{ef}}$
$\displaystyle=\frac{k^{2}}{2a}\left[\left(\dot{\bar{\varphi}}_{c}+m_{\phi}\bar{\varphi}_{s}\right)\delta\varphi_{c}+\left(\dot{\bar{\varphi}}_{s}-m_{\phi}\bar{\varphi}_{c}\right)\delta\varphi_{s}\right]$
(81b) $\displaystyle\delta P_{\phi,\mathrm{ef}}$
$\displaystyle=\delta\rho_{\phi,\mathrm{ef}}-m_{\phi}^{2}\left[\bar{\varphi}_{c}\delta\varphi_{c}+\bar{\varphi}_{s}\delta\varphi_{s}\right]$
(81c)

exactly satisfy the perturbed fluid equations

$\displaystyle 0$
$\displaystyle=\delta\rho_{\phi,\mathrm{ef}}^{\prime}+3\mathcal{H}\left(\delta\rho_{\phi,\mathrm{ef}}+\delta
P_{\phi,\mathrm{ef}}\right)+\left(\bar{\rho}_{\phi,\mathrm{ef}}+\bar{P}_{\phi,\mathrm{ef}}\right)\left(\theta_{\phi,\mathrm{ef}}+\frac{h^{\prime}}{2}\right)$
(82a) $\displaystyle 0$
$\displaystyle=\partial_{\tau}\left[\left(\bar{\rho}+\bar{P}\right)\theta_{\phi,\mathrm{ef}}\right]+4\mathcal{H}\left(\bar{\rho}_{\phi,\mathrm{ef}}+\bar{P}_{\phi,\mathrm{ef}}\right)\theta_{\phi,\mathrm{ef}}-k^{2}\delta
P_{\phi,\mathrm{ef}}.$ (82b)

Matching proceeds identically to the background, i.e., using the analogous
matching condition to Eq. 75 and solving for $\delta\varphi_{c}$,
$\delta\varphi_{s}$, and their time derivatives in analogy to Eqs. 77a and 76.

After matching (and computing $\delta_{\phi}$ and $\theta_{\phi}$ at
$t_{\star}$), one solves the effective fluid equations [205]

$\displaystyle\begin{split}\delta_{\phi}^{\prime}&=-\left(1+w_{\phi,\mathrm{ef}}\right)\left(\theta_{\phi}+\frac{h^{\prime}}{2}\right)-3\left(c_{s}^{2}-w_{\phi,\mathrm{ef}}\right)\mathcal{H}\delta_{\phi}-9\mathcal{H}^{2}\left(1+w_{\phi,\mathrm{ef}}\right)\left(c_{s}^{2}-c_{a}^{2}\right)\frac{\theta_{\phi}}{k^{2}}\end{split}$
(83a) $\displaystyle\theta_{\phi}^{\prime}$
$\displaystyle=-\left(1-3c_{s}^{2}\right)\mathcal{H}\theta_{\phi}+\frac{c_{s}^{2}k^{2}}{1+w_{\phi,\mathrm{ef}}}\delta_{\phi},$
(83b)

where $c_{s}^{2}$ is the rest-frame fluid sound speed and the “adiabatic sound
speed” squared is

$\displaystyle c_{a}^{2}$
$\displaystyle=w_{\phi,\mathrm{ef}}-\frac{w_{\phi,\mathrm{ef}}^{\prime}}{3\mathcal{H}\left(1+w_{\phi,\mathrm{ef}}\right)}.$
(84)

Ref. [121] used a Wentzel-Kramers-Brillouin (WKB) approximation to compute an
effective sound speed $c_{s}^{2}=(\sqrt{1+\kappa^{2}}-1)/\kappa$ (using the
shorthand $\kappa=k/am_{\phi}$) in the limit that expansion is negligible.
Ref. [121] also calibrated an $\mathcal{O}[(H/m_{\phi})^{2}]$ correction to
the squared sound speed, finding (empirically from numerical solutions) the
coefficient to be $5/4$, different from the factor $3/2$ found for the
equation of state that one might naively expect. We observe in addition that a
coefficient $9/8$ works well in the matter era. To reproduce both regimes, we
take the ansatz

$\displaystyle c_{s}^{2}$
$\displaystyle=\frac{\sqrt{1+(k/am_{\phi})^{2}}-1}{k/am_{\phi}}-\frac{1}{4}\frac{\dot{H}}{m_{\phi}^{2}}+\frac{3}{4}\left(\frac{H}{m_{\phi}}\right)^{2},$
(85)

which we empirically observe works well across a wide range of masses.999
Curiously, the correction to the sound speed does coincide with that to
equation of state [Eq. 79] in the matter era.  We observe roughly an order-of-
magnitude improvement in precision compared to using the sound speed from Ref.
[121] (at the same $H_{\star}/m_{\phi}$) for scalars that begin oscillating in
the matter era. We defer a more thorough investigation of convergence and
possible analytic justification for Eq. 85 to future work.

###  A.3 Initial conditions

We set initial conditions in the radiation era, when the scale factor evolves
as $a(\tau)=\sqrt{\Omega_{r}}H_{0}\tau$. The homogeneous equation of motion
for the scalar has a slow-roll solution (i.e., where
$\mathrm{d}V/\mathrm{d}\phi$ is approximately constant) of

$\displaystyle\bar{\phi}^{\prime}(\tau)$
$\displaystyle=-\frac{\Omega_{r}H_{0}^{2}}{5}\frac{\mathrm{d}V}{\mathrm{d}\phi}\tau^{3}.$
(86)

In the superhorizon limit, $h=C(k\tau)^{2}$. While $\bar{\phi}$ is slowly
rolling in the radiation era, the superhorizon limit of the equation of motion
for $\delta\phi$ (i.e., at leading order in $k\tau\ll 1$) is solved by

$\displaystyle\delta\phi$ $\displaystyle\approx
C\frac{\Omega_{r}H_{0}^{2}}{210k^{4}}\frac{\mathrm{d}V}{\mathrm{d}\phi}(k\tau)^{6}.$
(87)

Plugging the slow-roll solution Eq. 86 and noting that

$\displaystyle\bar{\rho}_{\phi}+\bar{P}_{\phi}$
$\displaystyle=\frac{\left(\bar{\phi}^{\prime}\right)^{2}}{a^{2}},$ (88)

the corresponding initial conditions for the fluid variables are

$\displaystyle\frac{\delta_{\phi}}{C}$
$\displaystyle\equiv\frac{\delta\rho_{\phi}}{C\bar{\rho}_{\phi}}=-\frac{\Omega_{r}H_{0}^{2}}{1050k^{4}V(\bar{\phi})}\left(\frac{\mathrm{d}V}{\mathrm{d}\phi}\right)^{2}(k\tau)^{6}$
(89) $\displaystyle\frac{\theta_{\phi}}{C}$
$\displaystyle\equiv\frac{-k^{2}\delta
u_{\phi}}{C}=\frac{k^{2}\delta\phi}{C\bar{\phi}^{\prime}}=-\frac{k(k\tau)^{3}}{42},$
(90)

consistent with those in Ref. [131].

###  A.4 Implementation details

Finally, we mention a few nontrivial aspects of the implementation of the
above formalism in CLASS, particularly related to switching solving schemes.
To begin with, switching between solving the Klein-Gordon equation and
effective fluid equations changes the system of equations (and degrees of
freedom) that are solved for at some intermediate time, based on a condition
that itself depends on the instantaneous value of a dynamical quantity
(namely, $H/m_{\phi}$). Aside from the fact that the exact time this
transition occurs cannot be known in advanced (i.e., not without integrating
the equations in question), performing such transition while using adaptive
ordinary differential equation (ODE) integrations is tricky. In particular,
the transition condition may be reached in the middle of an integration
interval; changing the system of equations mid-step introduces discontinuities
that violate the smoothness conditions required by ODE solvers. One must
therefore execute the transition between integration steps. (ODE solvers
evaluate the system of equations at multiple subintervals of each step.)
However, when using “black box” integrators (i.e., where one passes a function
that evaluates the differential equations to a solver routine that performs
the entire integration from an initial condition until some endpoint, as
implemented in CLASS) it is nontrivial to determine whether a function
evaluation is the first evaluation of an integration step or some intermediate
evaluation. When integrating the background cosmology, we therefore implement
the switching logic in a callback method that CLASS’s integrator executes
between integration step (i.e., the background_sources method, in which
various derived quantities are tabulated). Still, we find that CLASS’s stiff
(implicit) ODE solver still does not work with the above procedure, requiring
that we use its explicit Runge-Kutta method instead. This could be due to the
resulting change in the Jacobian, which cannot be accounted for because the
user cannot decide when the Jacobian is recomputed (and it is not recomputed
every integration step). For the solution of perturbations, CLASS does have a
built-in system for scheduling and executing approximation scheme transitions
(used extensively for other sectors, e.g., tight-coupling approximations for
the photon-baryon plasma). We therefore simply add another such approximation
scheme to switch on the effective fluid treatment for scalar perturbations.

Several issues also arise due to CLASS’s method for tabulating and
interpolating the background cosmology. A standard feature of ODE solvers is
so-called dense output, which combines function evaluations used for
integration to additionally produce interpolating functions over solver steps.
The interpolants are designed to be accurate to some high order in the step
size (comparable to the accuracy of the integration itself). CLASS instead
simply tabulates the results at fixed intervals in $\ln a$ and constructs a
cubic Hermite interpolant. This interpolation is therefore sensitive to the
transition between the Klein-Gordon and effective fluid solutions, which can
induce spurious errors in the interpolant that have a nonnegligible effect on
later stages in CLASS’s execution (discussed below). However, it is not
tractable to simply specify that the background be tabulated with a smaller
$\ln a$ interval: CLASS restricts integration step sizes to be no larger than
this interval, drastically increasing the runtime of the background solution.
(Standard ODE solver libraries instead take as large of steps as the adaptive
method permits and evaluates the dense output interpolants.) Decreasing the
tabulation interval by just two orders of magnitude over the default
[$\Delta\ln a=\ln(a_{0}/a_{i})/N\approx 8\times 10^{-4}$, where
$a_{i}/a_{0}=10^{-14}$ and $N=40,000$ by default] makes the background runtime
comparable to that of the perturbations solution. Because intervals in
physical time increase as $\Delta t=\Delta\ln a/H$, this tabulation becomes
insufficient to accurately capture the scalar’s oscillatory behavior when
$m_{\phi}\Delta t=\Delta\ln a\cdot m_{\phi}/H\gtrsim 1/10$, which for the
default parameters occurs when $m_{\phi}/H\sim 100$. Because of the drastic
effect it has on subsequent calculations, we increase the default number of
subintervals to $N=10^{5}$, but even this is insufficient to efficiently and
rigorously test the accuracy of the fluid approximation by comparing choices
of $m_{\phi}/H_{\star}$ as large as, say, $10^{4}$. This choice is, however,
sufficient for a switching time $m_{\phi}/H_{\star}=10$ as we use in all
results. The relative difference in, e.g., CMB spectra compared to solutions
with larger $m_{\phi}/H_{\star}$ is better than the $10^{-3}$ level.

The errors from the inadequate tabulation of the background solution naturally
propagate to the evolution of scalar perturbations and can qualitatively alter
the dynamics if unmitigated. They also can trigger CLASS to crash when
determining the sampling points in time at which to evaluate the CMB source
functions (in the perturbations_timesampling_for_sources routine). To
determine the time sampling required to capture the late ISW effect, CLASS
estimates the metric growth rate, which factors in
$2a^{\prime\prime}/a-(a^{\prime}/a)^{2}=-a^{2}\bar{P}/M_{\mathrm{pl}}^{2}$. In
the matter era, the _only_ contribution to $\bar{P}$ is that from the scalar,
which oscillates between $\bar{\rho}$ and $-\bar{\rho}$. This metric is
therefore strongly sensitive to any interpolation errors, which over the
course of parameter sampling often triggered CLASS to crash due to invalid
time sampling determined by this routine. To mitigate these issues, we
implemented a smooth transition between the Klein-Gordon and effective fluid
solutions. Namely, over some number of oscillations (specified by a parameter)
after $H$ drops below the specified threshold $H_{\star}$, we solve both the
Klein-Gordon and effective fluid equations. In this intermediate interval,
when computing the scalar’s energy density and pressure, we average the
results of the two methods with weights determined by a $\mathrm{tanh}$-based
window function in time. In effect, this procedure (exponentially) smoothly
damps the effect of oscillations in the true (Klein-Gordon) solution before
fully switching to solving the effective fluid equations alone. Since the
perturbation matching requires knowing $\bar{\varphi}_{c}$,
$\bar{\varphi}_{s}$, their time derivatives, and
$\bar{\rho}_{\phi,\mathrm{ef}}$ as well, we record these at the midpoint of
this transition interval for use when matching perturbations. Note that we use
still $\bar{\rho}_{\phi,\mathrm{ef}}$ and $\bar{P}_{\phi,\mathrm{ef}}$ when
matching to effective fluid perturbations, _not_ the weighted energy density
and pressure, which are only used in the Friedmann equations.

Finally, a minor complication arises from the effective equation of state
$w_{\phi,\mathrm{ef}}$ [Eq. 79] depending on $\dot{H}$. The Friedmann
equations set $\dot{H}=-(\bar{\rho}+\bar{P})/2M_{\mathrm{pl}}^{2}$, but one
does not know the pressure contribution from the scalar
$\bar{P}_{\phi,\mathrm{ef}}$ without knowing the current value of
$w_{\phi,\mathrm{ef}}$. Since we only consider scenarios where the scalar
contributes a subfraction of the total matter density, we simply iteratively
compute $\dot{H}$ (beginning by neglecting the scalar field’s pressure
contribution) and $w_{\phi,\mathrm{ef}}$ until the values converge to a
relative precision of $10^{-12}$.

##  Appendix B Parameter inference

###  B.1 Likelihoods

We use the implementation of the standard 2018 Planck likelihoods [75, 151]
(PR3) in clik [216] with supporting data obtained from the Planck Legacy
Archive [217]. For the high-$\ell$ likelihoods (covering $30\leq\ell<2500$ in
temperature and $30\leq\ell<2000$ for $E$-mode polarization and temperature-
polarization cross correlation) we use the Plik_lite (plik_lite_v22_TTTEEE)
variant that is marginalized over the parameters of the foreground models. For
the low-$\ell$ likelihoods ($2\leq\ell\leq 30$) we use Commander
(commander_dx12_v3_2_29) for temperature and SimAll
(simall_100x143_offlike5_EE_Aplanck_B) for $E$-mode polarization. Finally, we
use Planck’s lensing autopower spectrum likelihood
(smicadx12_Dec5_ftl_mv2_ndclpp_p_teb_consext8) over multipoles $8\leq L\leq
400$.

We use BAO measurements from a variety of surveys. The Six-degree Field Galaxy
Survey (6dFGS) [152] measured $\theta_{\mathrm{BAO}}=0.327\pm 0.0142$ at an
effective sample resdhift $z=0.106$. Ref [152] used fitting functions from
Ref. [63] to compute $r_{\mathrm{d}}$, which for their fiducial cosmology with
$h=0.7$, $\omega_{b}=0.02227$, and $\Omega_{m}=0.27$ gives
$\mathrm{d}=154.06\leavevmode\nobreak\ \mathrm{Mpc}$, compared to
$149.94\leavevmode\nobreak\ \mathrm{Mpc}$ from CLASS. The 6dFGS measurement of
$\theta_{\mathrm{BAO}}=0.336\pm 0.015$ (both the mean and standard deviation)
is thus rescaled by $149.94/154.06$. The SDSS Main Galaxy Sample (MGS) DR7
provides a non-Gaussian likelihood for $1/\theta_{\mathrm{BAO}}$ evaluated at
$z=0.15$ with tabulated values (specified relative to a fiducial value
$D_{V}/r_{\mathrm{d}}=638.95/148.69$) available in supplementary material to
Ref. [153]. The BAO measurements from the Baryon Oscillation Spectroscopic
Survey (BOSS) DR12 galaxies (at $z=0.38$ and $0.51$) are a a multivariate
Gaussian likelihood for $\theta_{\mathrm{BAO},\perp}$ and
$\theta_{\mathrm{BAO},\parallel}$ at each redshift [218]. Likewise, luminous
red galaxies (LRGs) from the Extended Baryon Oscillation Spectroscopic Survey
(eBOSS) DR16 [154, 155] provide a multivariate Gaussian likelihood for
$1/\theta_{\mathrm{BAO},\perp}$ and $1/\theta_{\mathrm{BAO},\parallel}$ at
$z=0.70$. The latter are available on the SDSS-IV SVN Software Repository
[219].

In addition to the former “standard” set of BAO measurements, we also
implement Gaussian likelihoods for the BAO measurements from DESI DR1 [156,
157, 158] as tabulated in Table 1 of Ref. [156] (which has not yet been
published at the time of writing). Ref. [156] notes that each each redshift
bin are effectively uncorrelated and the results from individual tracers are
reasonably concordant, justifying their combination. Notably, rounding errors
in Table 1 of Ref. [156] have a nonnegligible effect on the likelihoods. In
particular, the effective redshift for the bright galaxy sample measurement of
$\theta_{\mathrm{BAO}}$ is reported as $0.30$, rounded from $0.295$. The
impact of the rounding error is notable, shifting the full DESI DR1 dataset’s
preferred $hr_{\mathrm{d}}$ upward by about half a percent and $a_{m-\Lambda}$
downward by nearly a percent. To our knowledge, at the time of writing the
actual data are only available at Ref. [220].

Finally, we use the Pantheon [95], Pantheon+ [96, 97], and DES 5YR [159]
measurements of the apparent magnitude [Eq. 34] of type Ia supernovae. We also
use a recent joint analysis of supernovae from numerous datasets, Union3
[160]. Per Ref. [96], SN Ia at redshifts less than $10^{-2}$, which are more
sensitive to bias due to peculiar velocities and other effects, are excluded
for Pantheon+. The Pantheon+ dataset also includes Cepheid host distances from
SH0ES [26] that can optionally be used as calibrators. Likelihoods are
Gaussian in the distance moduli $\mu=m-M_{B}$, with measured values and
covariance provided online [221, 222, 223, 224].

###  B.2 Sampling methods

To generate posterior samples, we use Markov Chain Monte Carlo (MCMC) methods
from the Python package emcee [194, 195, 196]. We utilize the ensemble move
proposal based on kernel density estimation [225, 226], which (for the
posteriors presented here) we find yields MCMC chains with shortest
autocorrelation times of emcee’s methods. We sample with 80 walkers, chosen so
that each step, which emcee performs in two batches, runs (unparallelized)
likelihood evaluations on each of the 40 cores on the CPUs we use.
Autocorrelation times (for sample parameters) span from $\sim 7$ (for the
best-constrained posteriors) to $\sim 30$ steps (for, e.g., the varying-
constant scenarios with a scalar, since ${F_{\phi}}$ is poorly constrained, or
for varying-$m_{e}$ when using only Planck data, for which the posteriors are
broad). We sample for 5,000 steps or for 10,000 steps for those cases with
more challenging posteriors, discard the first ten autocorrelation times’
worth of steps, and thin by one autocorrelation time. These choices ensure
that reported posteriors include at least $\sim 20,000-30,000$ independent
samples, enough to draw robust posterior mass contours even in the presence of
extremely strong correlations (e.g., Fig. 8) or parameters whose constraints
only place upper limits (e.g., ${F_{\phi}}$). These sample sizes are entirely
generous for one-dimensional marginalized posteriors and summary statistics
thereof, which are all consistent among smaller subsamples of the full
posterior.

Finally, we note that we checked the robustness of our results with CLASS’s
default precision settings by comparing likelihood evaluations with increased
precision over 1,000 samples of the posteriors for each scenario we consider.
In all cases, for $95\%$ of samples the logarithm of the high-$\ell$ Planck
likelihood differs by $\lesssim 0.15$ when increasing precision, and only
$0.3$ at most; the others deviate negligibly. Namely, the distribution of
errors is no worse than that for $\Lambda$CDM at the default precision
settings.

## References

  * Hannestad [1999] S. Hannestad, Possible constraints on the time variation of the fine structure constant from cosmic microwave background data, Phys. Rev. D 60, 023515 (1999), arXiv:astro-ph/9810102 .
  * Kaplinghat  _et al._ [1999] M. Kaplinghat, R. J. Scherrer, and M. S. Turner, Constraining variations in the fine structure constant with the cosmic microwave background, Phys. Rev. D 60, 023516 (1999), arXiv:astro-ph/9810133 .
  * Avelino  _et al._ [2000] P. P. Avelino, C. J. A. P. Martins, G. Rocha, and P. T. P. Viana, Looking for a varying alpha in the cosmic microwave background, Phys. Rev. D 62, 123508 (2000), arXiv:astro-ph/0008446 .
  * Battye  _et al._ [2001] R. A. Battye, R. Crittenden, and J. Weller, Cosmic concordance and the fine structure constant, Phys. Rev. D 63, 043505 (2001), arXiv:astro-ph/0008265 .
  * Avelino  _et al._ [2001] P. P. Avelino, S. Esposito, G. Mangano, C. J. A. P. Martins, A. Melchiorri, G. Miele, O. Pisanti, G. Rocha, and P. T. P. Viana, Early universe constraints on a time varying fine structure constant, Phys. Rev. D 64, 103505 (2001), arXiv:astro-ph/0102144 .
  * Landau  _et al._ [2001] S. J. Landau, D. D. Harari, and M. Zaldarriaga, Constraining nonstandard recombination: A worked example, Phys. Rev. D 63, 083505 (2001).
  * Martins  _et al._ [2004] C. J. A. P. Martins, A. Melchiorri, G. Rocha, R. Trotta, P. P. Avelino, and P. T. P. Viana, Wmap constraints on varying alpha and the promise of reionization, Phys. Lett. B 585, 29 (2004), arXiv:astro-ph/0302295 .
  * Rocha  _et al._ [2004] G. Rocha, R. Trotta, C. J. A. P. Martins, A. Melchiorri, P. P. Avelino, R. Bean, and P. T. P. Viana, Measuring alpha in the early universe: cmb polarization, reionization and the fisher matrix analysis, Mon. Not. Roy. Astron. Soc. 352, 20 (2004), arXiv:astro-ph/0309211 .
  * Stefanescu [2007] P. Stefanescu, Constraints on time variation of fine structure constant from WMAP-3yr data, New Astron. 12, 635 (2007), arXiv:0707.0190 [astro-ph] .
  * Nakashima  _et al._ [2008] M. Nakashima, R. Nagata, and J. Yokoyama, Constraints on the time variation of the fine structure constant by the 5-year WMAP data, Prog. Theor. Phys. 120, 1207 (2008), arXiv:0810.1098 [astro-ph] .
  * Menegoni  _et al._ [2009] E. Menegoni, S. Galli, J. G. Bartlett, C. J. A. P. Martins, and A. Melchiorri, New Constraints on variations of the fine structure constant from CMB anisotropies, Phys. Rev. D 80, 087302 (2009), arXiv:0909.3584 [astro-ph.CO] .
  * Menegoni  _et al._ [2012] E. Menegoni, M. Archidiacono, E. Calabrese, S. Galli, C. J. A. P. Martins, and A. Melchiorri, The Fine Structure Constant and the CMB Damping Scale, Phys. Rev. D 85, 107301 (2012), arXiv:1202.1476 [astro-ph.CO] .
  * Kujat and Scherrer [2000] J. Kujat and R. J. Scherrer, The Effect of time variation in the Higgs vacuum expectation value on the cosmic microwave background, Phys. Rev. D 62, 023510 (2000), arXiv:astro-ph/9912174 .
  * Ichikawa  _et al._ [2006] K. Ichikawa, T. Kanzaki, and M. Kawasaki, CMB constraints on the simultaneous variation of the fine structure constant and electron mass, Phys. Rev. D 74, 023515 (2006), arXiv:astro-ph/0602577 .
  * Landau  _et al._ [2008] S. J. Landau, M. E. Mosquera, C. G. Scoccola, and H. Vucetich, Early Universe Constraints on Time Variation of Fundamental Constants, Phys. Rev. D 78, 083527 (2008), arXiv:0809.2033 [astro-ph] .
  * Scoccola  _et al._ [2008] C. G. Scoccola, S. J. Landau, and H. Vucetich, WMAP 5-year constraints on time variation of $\alpha$ and $m_{e}$ in a detailed recombination scenario, Phys. Lett. B 669, 212 (2008), arXiv:0809.5028 [astro-ph] .
  * Nakashima  _et al._ [2010] M. Nakashima, K. Ichikawa, R. Nagata, and J. Yokoyama, Constraining the time variation of the coupling constants from cosmic microwave background: effect of \Lambda_QCD, JCAP 01, 030, arXiv:0910.0742 [astro-ph.CO] .
  * Landau and Scoccola [2010] S. J. Landau and C. G. Scoccola, Constraints on variation in $\alpha$ and $m_{e}$ from WMAP 7-year data, Astron. Astrophys. 517, A62 (2010), arXiv:1002.1603 [astro-ph.CO] .
  * Scoccola  _et al._ [2013] C. G. Scoccola _et al._ , The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: constraints on the time variation of fundamental constants from the large-scale two-point correlation function, Mon. Not. Roy. Astron. Soc. 434, 1792 (2013), arXiv:1209.1394 [astro-ph.CO] .
  * Hart and Chluba [2020] L. Hart and J. Chluba, Updated fundamental constant constraints from Planck 2018 data and possible relations to the Hubble tension, Mon. Not. Roy. Astron. Soc. 493, 3255 (2020), arXiv:1912.03986 [astro-ph.CO] .
  * Sekiguchi and Takahashi [2021] T. Sekiguchi and T. Takahashi, Early recombination as a solution to the $H_{0}$ tension, Phys. Rev. D 103, 083507 (2021), arXiv:2007.03381 [astro-ph.CO] .
  * Schöneberg  _et al._ [2022] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, The H0 Olympics: A fair ranking of proposed models, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Khalife  _et al._ [2024] A. R. Khalife, M. B. Zanjani, S. Galli, S. Günther, J. Lesgourgues, and K. Benabed, Review of Hubble tension solutions with new SH0ES and SPT-3G data, JCAP 04, 059, arXiv:2312.09814 [astro-ph.CO] .
  * Riess  _et al._ [2016] A. G. Riess _et al._ , A 2.4% Determination of the Local Value of the Hubble Constant, Astrophys. J. 826, 56 (2016), arXiv:1604.01424 [astro-ph.CO] .
  * Riess  _et al._ [2019] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM, Astrophys. J. 876, 85 (2019), arXiv:1903.07603 [astro-ph.CO] .
  * Riess  _et al._ [2022] A. G. Riess _et al._ , A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s-1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SH0ES Team, Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Fung  _et al._ [2023] L. W. H. Fung, L. Li, T. Liu, H. N. Luu, Y.-C. Qiu, and S. H. H. Tye, Hubble constant in the axi-Higgs universe, Phys. Rev. Res. 5, L022059 (2023), arXiv:2105.01631 [astro-ph.CO] .
  * Fung  _et al._ [2021] L. W. H. Fung, L. Li, T. Liu, H. N. Luu, Y.-C. Qiu, and S. H. H. Tye, Axi-Higgs cosmology, JCAP 08, 057, arXiv:2102.11257 [hep-ph] .
  * Luu [2023] H. N. Luu, Axion-Higgs cosmology: Cosmic microwave background and cosmological tensions, Phys. Rev. D 107, 023513 (2023), arXiv:2111.01347 [astro-ph.CO] .
  * Solomon  _et al._ [2022] R. Solomon, G. Agarwal, and D. Stojkovic, Environment dependent electron mass and the Hubble constant tension, Phys. Rev. D 105, 103536 (2022), arXiv:2201.03127 [hep-ph] .
  * Hoshiya and Toda [2023] K. Hoshiya and Y. Toda, Electron mass variation from dark sector interactions and compatibility with cosmological observations, Phys. Rev. D 107, 043505 (2023), arXiv:2202.07714 [astro-ph.CO] .
  * Baryakhtar, Masha and Simon, Olivier and Weiner, Zachary J. [2024] Baryakhtar, Masha and Simon, Olivier and Weiner, Zachary J., Searching for coupled, hyperlight scalars across cosmic history (2024), to appear.
  * Tohfa  _et al._ [2023] H. Tohfa, J. Crump, E. Baker, L. Hart, D. Grin, M. Brosius, and J. Chluba, A cosmic microwave background search for fine-structure constant evolution, arXiv:2307.06768 [astro-ph.CO] (2023).
  * Schöneberg  _et al._ [2023] N. Schöneberg, L. Vacher, J. D. F. Dias, M. M. C. D. Carvalho, and C. J. A. P. Martins, News from the Swampland — constraining string theory with astrophysics and cosmology, JCAP 10, 039, arXiv:2307.15060 [astro-ph.CO] .
  * Vacher  _et al._ [2023] L. Vacher, N. Schöneberg, J. D. F. Dias, C. J. A. P. Martins, and F. Pimenta, Runaway dilaton models: Improved constraints from the full cosmological evolution, Phys. Rev. D 107, 104002 (2023), arXiv:2301.13500 [astro-ph.CO] .
  * Fixsen  _et al._ [1996] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, The Cosmic Microwave Background spectrum from the full COBE FIRAS data set, Astrophys. J. 473, 576 (1996), arXiv:astro-ph/9605054 .
  * Fixsen [2009] D. J. Fixsen, The Temperature of the Cosmic Microwave Background, Astrophys. J. 707, 916 (2009), arXiv:0911.1955 [astro-ph.CO] .
  * Ivanov  _et al._ [2020a] M. M. Ivanov, Y. Ali-Haïmoud, and J. Lesgourgues, H0 tension or T0 tension?, Phys. Rev. D 102, 063515 (2020a), arXiv:2005.10656 [astro-ph.CO] .
  * Cyr-Racine  _et al._ [2022] F.-Y. Cyr-Racine, F. Ge, and L. Knox, Symmetry of Cosmological Observables, a Mirror World Dark Sector, and the Hubble Constant, Phys. Rev. Lett. 128, 201301 (2022), arXiv:2107.13000 [astro-ph.CO] .
  * Ge  _et al._ [2023] F. Ge, F.-Y. Cyr-Racine, and L. Knox, Scaling transformations and the origins of light relics constraints from cosmic microwave background observations, Phys. Rev. D 107, 023517 (2023), arXiv:2210.16335 [astro-ph.CO] .
  * Greene and Cyr-Racine [2023] K. L. Greene and F.-Y. Cyr-Racine, Thomson scattering: one rate to rule them all, JCAP 10, 065, arXiv:2306.06165 [astro-ph.CO] .
  * Akita and Yamaguchi [2020] K. Akita and M. Yamaguchi, A precision calculation of relic neutrino decoupling, JCAP 08, 012, arXiv:2005.07047 [hep-ph] .
  * Froustey  _et al._ [2020] J. Froustey, C. Pitrou, and M. C. Volpe, Neutrino decoupling including flavour oscillations and primordial nucleosynthesis, JCAP 12, 015, arXiv:2008.01074 [hep-ph] .
  * Bennett  _et al._ [2021] J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, Towards a precision calculation of $N_{\rm eff}$ in the Standard Model II: Neutrino decoupling in the presence of flavour oscillations and finite-temperature QED, JCAP 04, 073, arXiv:2012.02726 [hep-ph] .
  * Cielo  _et al._ [2023] M. Cielo, M. Escudero, G. Mangano, and O. Pisanti, Neff in the Standard Model at NLO is 3.043, Phys. Rev. D 108, L121301 (2023), arXiv:2306.05460 [hep-ph] .
  * Birrell  _et al._ [2014] J. Birrell, C.-T. Yang, and J. Rafelski, Relic Neutrino Freeze-out: Dependence on Natural Constants, Nucl. Phys. B 890, 481 (2014), arXiv:1406.1759 [nucl-th] .
  * Grohs and Fuller [2017] E. Grohs and G. M. Fuller, Insights into neutrino decoupling gleaned from considerations of the role of electron mass, Nucl. Phys. B 923, 222 (2017), arXiv:1706.03391 [astro-ph.CO] .
  * Abazajian  _et al._ [2016] K. N. Abazajian _et al._ (CMB-S4), CMB-S4 Science Book, First Edition, arXiv:1610.02743 [astro-ph.CO] (2016).
  * Abazajian  _et al._ [2019] K. Abazajian _et al._ , CMB-S4 Science Case, Reference Design, and Project Plan, arXiv:1907.04473 [astro-ph.IM] (2019).
  * Dvorkin  _et al._ [2022] C. Dvorkin _et al._ , The Physics of Light Relics, in  _Snowmass 2021_ (2022) arXiv:2203.07943 [hep-ph] .
  * Blas  _et al._ [2011] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 07, 034, arXiv:1104.2933 [astro-ph.CO] .
  * Lesgourgues [2011] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview, arXiv:1104.2932 [astro-ph.IM] (2011).
  * Ali-Haimoud and Hirata [2011] Y. Ali-Haimoud and C. M. Hirata, HyRec: A fast and highly accurate primordial hydrogen and helium recombination code, Phys. Rev. D 83, 043513 (2011), arXiv:1011.3758 [astro-ph.CO] .
  * Lee and Ali-Haïmoud [2020] N. Lee and Y. Ali-Haïmoud, HYREC-2: a highly accurate sub-millisecond recombination code, Phys. Rev. D 102, 083517 (2020), arXiv:2007.14114 [astro-ph.CO] .
  * Peebles [1968] P. J. E. Peebles, Recombination of the Primeval Plasma, Astrophys. J. 153, 1 (1968).
  * Peebles and Yu [1970] P. J. E. Peebles and J. T. Yu, Primeval adiabatic perturbation in an expanding universe, Astrophys. J. 162, 815 (1970).
  * Seager  _et al._ [1999] S. Seager, D. D. Sasselov, and D. Scott, A new calculation of the recombination epoch, Astrophys. J. Lett. 523, L1 (1999), arXiv:astro-ph/9909275 .
  * Ali-Haimoud and Hirata [2010] Y. Ali-Haimoud and C. M. Hirata, Ultrafast effective multi-level atom method for primordial hydrogen recombination, Phys. Rev. D 82, 063521 (2010), arXiv:1006.1355 [astro-ph.CO] .
  * Ade  _et al._ [2015] P. A. R. Ade _et al._ (Planck), Planck intermediate results - XXIV. Constraints on variations in fundamental constants, Astron. Astrophys. 580, A22 (2015), arXiv:1406.7482 [astro-ph.CO] .
  * Hart and Chluba [2018] L. Hart and J. Chluba, New constraints on time-dependent variations of fundamental constants using Planck data, Mon. Not. Roy. Astron. Soc. 474, 1850 (2018), arXiv:1705.03925 [astro-ph.CO] .
  * Weinberg [2008] S. Weinberg,  _Cosmology_ (2008).
  * Baumann [2022] D. Baumann,  _Cosmology_ (Cambridge University Press, 2022).
  * Eisenstein and Hu [1998] D. J. Eisenstein and W. Hu, Baryonic features in the matter transfer function, Astrophys. J. 496, 605 (1998), arXiv:astro-ph/9709112 .
  * Hu and White [1997a] W. Hu and M. J. White, The Damping tail of CMB anisotropies, Astrophys. J. 479, 568 (1997a), arXiv:astro-ph/9609079 .
  * Hu and Sugiyama [1996] W. Hu and N. Sugiyama, Small scale cosmological perturbations: An Analytic approach, Astrophys. J. 471, 542 (1996), arXiv:astro-ph/9510117 .
  * Vagnozzi [2021] S. Vagnozzi, Consistency tests of $\Lambda$CDM from the early integrated Sachs-Wolfe effect: Implications for early-time new physics and the Hubble tension, Phys. Rev. D 104, 063524 (2021), arXiv:2105.10425 [astro-ph.CO] .
  * Adil  _et al._ [2023] A. Adil, A. Albrecht, and L. Knox, Quintessential cosmological tensions, Phys. Rev. D 107, 063521 (2023), arXiv:2207.10235 [astro-ph.CO] .
  * Hu [1995] W. T. Hu,  _Wandering in the Background: A CMB Explorer_ , Other thesis (1995), arXiv:astro-ph/9508126 .
  * Zaldarriaga and Harari [1995] M. Zaldarriaga and D. D. Harari, Analytic approach to the polarization of the cosmic microwave background in flat and open universes, Phys. Rev. D 52, 3276 (1995), arXiv:astro-ph/9504085 .
  * Hu and White [1997b] W. Hu and M. J. White, A CMB polarization primer, New Astron. 2, 323 (1997b), arXiv:astro-ph/9706147 .
  * Hadzhiyska and Spergel [2019] B. Hadzhiyska and D. N. Spergel, Measuring the Duration of Last Scattering, Phys. Rev. D 99, 043537 (2019), arXiv:1808.04083 [astro-ph.CO] .
  * Hogg [1999] D. W. Hogg, Distance measures in cosmology, arXiv:astro-ph/9905116 (1999).
  * Baes  _et al._ [2017] M. Baes, P. Camps, and D. Van De Putte, Analytical expressions and numerical evaluation of the luminosity distance in a flat cosmology, Mon. Not. Roy. Astron. Soc. 468, 927 (2017), arXiv:1702.08860 [astro-ph.CO] .
  * Knox and Millea [2020] L. Knox and M. Millea, Hubble constant hunter’s guide, Phys. Rev. D 101, 043533 (2020), arXiv:1908.03663 [astro-ph.CO] .
  * Aghanim  _et al._ [2020a] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020a), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Percival  _et al._ [2002] W. J. Percival _et al._ (2dFGRS Team), Parameter constraints for flat cosmologies from CMB and 2dFGRS power spectra, Mon. Not. Roy. Astron. Soc. 337, 1068 (2002), arXiv:astro-ph/0206256 .
  * Howlett  _et al._ [2012] C. Howlett, A. Lewis, A. Hall, and A. Challinor, CMB power spectrum parameter degeneracies in the era of precision cosmology, JCAP 04, 027, arXiv:1201.3654 [astro-ph.CO] .
  * Jedamzik and Pogosian [2020] K. Jedamzik and L. Pogosian, Relieving the Hubble tension with primordial magnetic fields, Phys. Rev. Lett. 125, 181302 (2020), arXiv:2004.09487 [astro-ph.CO] .
  * Rashkovetskyi  _et al._ [2021] M. Rashkovetskyi, J. B. Muñoz, D. J. Eisenstein, and C. Dvorkin, Small-scale clumping at recombination and the Hubble tension, Phys. Rev. D 104, 103517 (2021), arXiv:2108.02747 [astro-ph.CO] .
  * Thiele  _et al._ [2021] L. Thiele, Y. Guan, J. C. Hill, A. Kosowsky, and D. N. Spergel, Can small-scale baryon inhomogeneities resolve the Hubble tension? An investigation with ACT DR4, Phys. Rev. D 104, 063535 (2021), arXiv:2105.03003 [astro-ph.CO] .
  * Lee and Ali-Haïmoud [2021] N. Lee and Y. Ali-Haïmoud, Probing small-scale baryon and dark matter isocurvature perturbations with cosmic microwave background anisotropies, Phys. Rev. D 104, 103509 (2021), arXiv:2108.07798 [astro-ph.CO] .
  * Wyman  _et al._ [2014] M. Wyman, D. H. Rudd, R. A. Vanderveld, and W. Hu, Neutrinos Help Reconcile Planck Measurements with the Local Universe, Phys. Rev. Lett. 112, 051302 (2014), arXiv:1307.7715 [astro-ph.CO] .
  * Dvorkin  _et al._ [2014] C. Dvorkin, M. Wyman, D. H. Rudd, and W. Hu, Neutrinos help reconcile Planck measurements with both the early and local Universe, Phys. Rev. D 90, 083503 (2014), arXiv:1403.8049 [astro-ph.CO] .
  * Bernal  _et al._ [2016] J. L. Bernal, L. Verde, and A. G. Riess, The trouble with $H_{0}$, JCAP 10, 019, arXiv:1607.05617 [astro-ph.CO] .
  * Bashinsky and Seljak [2004] S. Bashinsky and U. Seljak, Neutrino perturbations in CMB anisotropy and matter clustering, Phys. Rev. D 69, 083002 (2004), arXiv:astro-ph/0310198 .
  * Hou  _et al._ [2013] 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, 083008 (2013), arXiv:1104.2333 [astro-ph.CO] .
  * Baumann  _et al._ [2016] D. Baumann, D. Green, J. Meyers, and B. Wallisch, Phases of New Physics in the CMB, JCAP 01, 007, arXiv:1508.06342 [astro-ph.CO] .
  * Pan  _et al._ [2016] Z. Pan, L. Knox, B. Mulroe, and A. Narimani, Cosmic Microwave Background Acoustic Peak Locations, Mon. Not. Roy. Astron. Soc. 459, 2513 (2016), arXiv:1603.03091 [astro-ph.CO] .
  * Sachs and Wolfe [1967] R. K. Sachs and A. M. Wolfe, Perturbations of a cosmological model and angular variations of the microwave background, Astrophys. J. 147, 73 (1967).
  * Weller and Lewis [2003] J. Weller and A. M. Lewis, Large scale cosmic microwave background anisotropies and dark energy, Mon. Not. Roy. Astron. Soc. 346, 987 (2003), arXiv:astro-ph/0307104 .
  * Lewis and Challinor [2006] A. Lewis and A. Challinor, Weak gravitational lensing of the CMB, Phys. Rept. 429, 1 (2006), arXiv:astro-ph/0601594 .
  * Thepsuriya and Lewis [2015] K. Thepsuriya and A. Lewis, Accuracy of cosmological parameters using the baryon acoustic scale, JCAP 01, 034, arXiv:1409.5066 [astro-ph.CO] .
  * Schöneberg  _et al._ [2019] N. Schöneberg, J. Lesgourgues, and D. C. Hooper, The BAO+BBN take on the Hubble tension, JCAP 10, 029, arXiv:1907.11594 [astro-ph.CO] .
  * Bernal  _et al._ [2020] J. L. Bernal, T. L. Smith, K. K. Boddy, and M. Kamionkowski, Robustness of baryon acoustic oscillation constraints for early-Universe modifications of $\Lambda$CDM cosmology, Phys. Rev. D 102, 123515 (2020), arXiv:2004.07263 [astro-ph.CO] .
  * Scolnic  _et al._ [2018] D. M. Scolnic _et al._ (Pan-STARRS1), The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample, Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Brout  _et al._ [2022] D. Brout _et al._ , The Pantheon+ Analysis: Cosmological Constraints, Astrophys. J. 938, 110 (2022), arXiv:2202.04077 [astro-ph.CO] .
  * Scolnic  _et al._ [2022] D. Scolnic _et al._ , The Pantheon+ Analysis: The Full Data Set and Light-curve Release, Astrophys. J. 938, 113 (2022), arXiv:2112.03863 [astro-ph.CO] .
  * Davis  _et al._ [2011] T. M. Davis _et al._ , The Effect of Peculiar Velocities on Supernova Cosmology, Astrophys. J. 741, 67 (2011), arXiv:1012.2912 [astro-ph.CO] .
  * Davis  _et al._ [2019] T. M. Davis, S. R. Hinton, C. Howlett, and J. Calcino, Can redshift errors bias measurements of the Hubble Constant?, Mon. Not. Roy. Astron. Soc. 490, 2948 (2019), arXiv:1907.12639 [astro-ph.CO] .
  * Carr  _et al._ [2022] A. Carr, T. M. Davis, D. Scolnic, D. Scolnic, K. Said, D. Brout, E. R. Peterson, and R. Kessler, The Pantheon+ analysis: Improving the redshifts and peculiar velocities of Type Ia supernovae used in cosmological analyses, Publ. Astron. Soc. Austral. 39, e046 (2022), arXiv:2112.01471 [astro-ph.CO] .
  * Freedman  _et al._ [2019] 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, Astrophys. J. 882, 34 (2019), arXiv:1907.05922 [astro-ph.CO] .
  * Freedman  _et al._ [2020] W. L. Freedman, B. F. Madore, T. Hoyt, I. S. Jang, R. Beaton, M. G. Lee, A. Monson, J. Neeley, and J. Rich, Calibration of the Tip of the Red Giant Branch (TRGB), arXiv:2002.01550 [astro-ph.GA] (2020).
  * Freedman [2021] W. L. Freedman, Measurements of the Hubble Constant: Tensions in Perspective, Astrophys. J. 919, 16 (2021), arXiv:2106.15656 [astro-ph.CO] .
  * Freedman and Madore [2023] W. L. Freedman and B. F. Madore, Progress in direct measurements of the Hubble constant, JCAP 11, 050, arXiv:2309.05618 [astro-ph.CO] .
  * Jain and Seljak [1997] B. Jain and U. Seljak, Cosmological model predictions for weak lensing: Linear and nonlinear regimes, Astrophys. J. 484, 560 (1997), arXiv:astro-ph/9611077 .
  * Abbott  _et al._ [2016] T. Abbott _et al._ (DES), Cosmology from cosmic shear with Dark Energy Survey Science Verification data, Phys. Rev. D 94, 022001 (2016), arXiv:1507.05552 [astro-ph.CO] .
  * Hildebrandt  _et al._ [2017] H. Hildebrandt _et al._ , KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing, Mon. Not. Roy. Astron. Soc. 465, 1454 (2017), arXiv:1606.05338 [astro-ph.CO] .
  * Bekenstein [1982] J. D. Bekenstein, Fine Structure Constant: Is It Really a Constant?, Phys. Rev. D 25, 1527 (1982).
  * Damour and Donoghue [2010a] T. Damour and J. F. Donoghue, Phenomenology of the Equivalence Principle with Light Scalars, Class. Quant. Grav. 27, 202001 (2010a), arXiv:1007.2790 [gr-qc] .
  * Damour and Donoghue [2010b] T. Damour and J. F. Donoghue, Equivalence Principle Violations and Couplings of a Light Dilaton, Phys. Rev. D 82, 084033 (2010b), arXiv:1007.2792 [gr-qc] .
  * Brans and Dicke [1961] C. Brans and R. H. Dicke, Mach’s principle and a relativistic theory of gravitation, Phys. Rev. 124, 925 (1961).
  * Scherk and Schwarz [1974] J. Scherk and J. H. Schwarz, Dual Models for Nonhadrons, Nucl. Phys. B 81, 118 (1974).
  * Green  _et al._ [2012] M. B. Green, J. H. Schwarz, and E. Witten,  _Superstring Theory Vol. 2: 25th Anniversary Edition_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2012).
  * Taylor and Veneziano [1988] T. R. Taylor and G. Veneziano, Dilaton Couplings at Large Distances, Phys. Lett. B 213, 450 (1988).
  * Cho [1990] Y. M. Cho, UNIFIED COSMOLOGY, Phys. Rev. D 41, 2462 (1990).
  * Damour and Polyakov [1994] T. Damour and A. M. Polyakov, String theory and gravity, Gen. Rel. Grav. 26, 1171 (1994), arXiv:gr-qc/9411069 .
  * Dimopoulos and Giudice [1996] S. Dimopoulos and G. F. Giudice, Macroscopic forces from supersymmetry, Phys. Lett. B 379, 105 (1996), arXiv:hep-ph/9602350 .
  * Kaplan and Wise [2000] D. B. Kaplan and M. B. Wise, Couplings of a light dilaton and violations of the equivalence principle, JHEP 08, 037, arXiv:hep-ph/0008116 .
  * Gasperini  _et al._ [2002] M. Gasperini, F. Piazza, and G. Veneziano, Quintessence as a runaway dilaton, Phys. Rev. D 65, 023508 (2002), arXiv:gr-qc/0108016 .
  * Damour  _et al._ [2002] T. Damour, F. Piazza, and G. Veneziano, Runaway dilaton and equivalence principle violations, Phys. Rev. Lett. 89, 081601 (2002), arXiv:gr-qc/0204094 .
  * Passaglia and Hu [2022] S. Passaglia and W. Hu, Accurate effective fluid approximation for ultralight axions, Phys. Rev. D 105, 123529 (2022), arXiv:2201.10238 [astro-ph.CO] .
  * Zachary J. Weiner [2023] Zachary J. Weiner, CLASS-ULS (2023), branch class-uls-3.2.1, commit hash a5d05e4.
  * Senatore  _et al._ [2009] L. Senatore, S. Tassev, and M. Zaldarriaga, Cosmological Perturbations at Second Order and Recombination Perturbed, JCAP 08, 031, arXiv:0812.3652 [astro-ph] .
  * Khatri and Wandelt [2009] R. Khatri and B. D. Wandelt, Crinkles in the last scattering surface: Non-Gaussianity from inhomogeneous recombination, Phys. Rev. D 79, 023501 (2009), arXiv:0810.4370 [astro-ph] .
  * de Salas  _et al._ [2021] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, 2020 global reassessment of the neutrino oscillation picture, JHEP 02, 071, arXiv:2006.11237 [hep-ph] .
  * Esteban  _et al._ [2020] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, The fate of hints: updated global analysis of three-flavor neutrino oscillations, JHEP 09, 178, arXiv:2007.14792 [hep-ph] .
  * Capozzi  _et al._ [2017] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Global constraints on absolute neutrino masses and their ordering, Phys. Rev. D 95, 096014 (2017), [Addendum: Phys.Rev.D 101, 116013 (2020)], arXiv:2003.08511 [hep-ph] .
  * Hlozek  _et al._ [2015] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, A search for ultralight axions using precision cosmological data, Phys. Rev. D 91, 103512 (2015), arXiv:1410.2896 [astro-ph.CO] .
  * Poulin  _et al._ [2019] V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Early Dark Energy Can Resolve The Hubble Tension, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Poulin  _et al._ [2018] V. Poulin, T. L. Smith, D. Grin, T. Karwal, and M. Kamionkowski, Cosmological implications of ultralight axionlike fields, Phys. Rev. D 98, 083525 (2018), arXiv:1806.10608 [astro-ph.CO] .
  * Smith  _et al._ [2020] 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, 063523 (2020), arXiv:1908.06995 [astro-ph.CO] .
  * Hill  _et al._ [2020] J. C. Hill, E. McDonough, M. W. Toomey, and S. Alexander, Early dark energy does not restore cosmological concordance, Phys. Rev. D 102, 043507 (2020), arXiv:2003.07355 [astro-ph.CO] .
  * Poulin  _et al._ [2023] V. Poulin, T. L. Smith, and T. Karwal, The Ups and Downs of Early Dark Energy solutions to the Hubble tension: A review of models, hints and constraints circa 2023, Phys. Dark Univ. 42, 101348 (2023), arXiv:2302.09032 [astro-ph.CO] .
  * Hu  _et al._ [2000] W. Hu, R. Barkana, and A. Gruzinov, Cold and fuzzy dark matter, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 .
  * Amendola and Barbieri [2006] L. Amendola and R. Barbieri, Dark matter from an ultra-light pseudo-Goldsone-boson, Phys. Lett. B 642, 192 (2006), arXiv:hep-ph/0509257 .
  * Marsh [2016a] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643, 1 (2016a), arXiv:1510.07633 [astro-ph.CO] .
  * Bond  _et al._ [1980] J. R. Bond, G. Efstathiou, and J. Silk, Massive Neutrinos and the Large Scale Structure of the Universe, Phys. Rev. Lett. 45, 1980 (1980).
  * Lesgourgues and Pastor [2006] J. Lesgourgues and S. Pastor, Massive neutrinos and cosmology, Phys. Rept. 429, 307 (2006), arXiv:astro-ph/0603494 .
  * Rogers  _et al._ [2023] K. K. Rogers, R. Hložek, A. Laguë, M. M. Ivanov, O. H. E. Philcox, G. Cabass, K. Akitsu, and D. J. E. Marsh, Ultra-light axions and the S 8 tension: joint constraints from the cosmic microwave background and galaxy clustering, JCAP 06, 023, arXiv:2301.08361 [astro-ph.CO] .
  * Hložek  _et al._ [2017] R. Hložek, D. J. E. Marsh, D. Grin, R. Allison, J. Dunkley, and E. Calabrese, Future CMB tests of dark matter: Ultralight axions and massive neutrinos, Phys. Rev. D 95, 123511 (2017), arXiv:1607.08208 [astro-ph.CO] .
  * Hlozek  _et al._ [2018] R. Hlozek, D. J. E. Marsh, and D. Grin, Using the Full Power of the Cosmic Microwave Background to Probe Axion Dark Matter, Mon. Not. Roy. Astron. Soc. 476, 3063 (2018), arXiv:1708.05681 [astro-ph.CO] .
  * Laguë  _et al._ [2022] A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, Constraining ultralight axions with galaxy surveys, JCAP 01 (01), 049, arXiv:2104.07802 [astro-ph.CO] .
  * Abdalla  _et al._ [2022] E. Abdalla _et al._ , Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies, JHEAp 34, 49 (2022), arXiv:2203.06142 [astro-ph.CO] .
  * Vogt  _et al._ [2023] S. M. L. Vogt, D. J. E. Marsh, and A. Laguë, Improved mixed dark matter halo model for ultralight axions, Phys. Rev. D 107, 063526 (2023), arXiv:2209.13445 [astro-ph.CO] .
  * Laguë  _et al._ [2024] A. Laguë, B. Schwabe, R. Hložek, D. J. E. Marsh, and K. K. Rogers, Cosmological simulations of mixed ultralight dark matter, Phys. Rev. D 109, 043507 (2024), arXiv:2310.20000 [astro-ph.CO] .
  * Marsh [2016b] D. J. E. Marsh, WarmAndFuzzy: the halo model beyond CDM, arXiv:1605.05973 [astro-ph.CO] (2016b).
  * Bauer  _et al._ [2020] J. B. Bauer, D. J. E. Marsh, R. Hložek, H. Padmanabhan, and A. Laguë, Intensity Mapping as a Probe of Axion Dark Matter, Mon. Not. Roy. Astron. Soc. 500, 3162 (2020), arXiv:2003.09655 [astro-ph.CO] .
  * Ooguri and Vafa [2007] H. Ooguri and C. Vafa, On the Geometry of the String Landscape and the Swampland, Nucl. Phys. B 766, 21 (2007), arXiv:hep-th/0605264 .
  * Olive and Pospelov [2002] K. A. Olive and M. Pospelov, Evolution of the fine structure constant driven by dark matter and the cosmological constant, Phys. Rev. D 65, 085044 (2002), arXiv:hep-ph/0110377 .
  * Olive and Pospelov [2008] K. A. Olive and M. Pospelov, Environmental dependence of masses and coupling constants, Phys. Rev. D 77, 043524 (2008), arXiv:0709.3825 [hep-ph] .
  * Aghanim  _et al._ [2020b] N. Aghanim _et al._ (Planck), Planck 2018 results. V. CMB power spectra and likelihoods, Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 [astro-ph.CO] .
  * Beutler  _et al._ [2011] 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, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ [2015] 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, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Bautista  _et al._ [2020] J. E. Bautista _et al._ (eBOSS), The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the luminous red galaxy sample from the anisotropic correlation function between redshifts 0.6 and 1, Mon. Not. Roy. Astron. Soc. 500, 736 (2020), arXiv:2007.08993 [astro-ph.CO] .
  * Gil-Marin  _et al._ [2020] H. Gil-Marin _et al._ (eBOSS), The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: measurement of the BAO and growth rate of structure of the luminous red galaxy sample from the anisotropic power spectrum between redshifts 0.6 and 1.0, Mon. Not. Roy. Astron. Soc. 498, 2492 (2020), arXiv:2007.08994 [astro-ph.CO] .
  * Adame  _et al._ [2024a] A. G. Adame _et al._ (DESI), DESI 2024 VI: Cosmological Constraints from the Measurements of Baryon Acoustic Oscillations, arXiv:2404.03002 [astro-ph.CO] (2024a).
  * Adame  _et al._ [2024b] A. G. Adame _et al._ (DESI), DESI 2024 IV: Baryon Acoustic Oscillations from the Lyman Alpha Forest, arXiv:2404.03001 [astro-ph.CO] (2024b).
  * Adame  _et al._ [2024c] A. G. Adame _et al._ (DESI), DESI 2024 III: Baryon Acoustic Oscillations from Galaxies and Quasars, arXiv:2404.03000 [astro-ph.CO] (2024c).
  * Abbott  _et al._ [2024] T. M. C. Abbott _et al._ (DES), The Dark Energy Survey: Cosmology Results With ~1500 New High-redshift Type Ia Supernovae Using The Full 5-year Dataset, arXiv:2401.02929 [astro-ph.CO] (2024).
  * Rubin  _et al._ [2023] D. Rubin _et al._ , Union Through UNITY: Cosmology with 2,000 SNe Using a Unified Bayesian Framework, arXiv:2311.12098 [astro-ph.CO] (2023).
  * Spergel  _et al._ [2003] D. N. Spergel _et al._ (WMAP), First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Determination of cosmological parameters, Astrophys. J. Suppl. 148, 175 (2003), arXiv:astro-ph/0302209 .
  * Efstathiou  _et al._ [2010] G. Efstathiou, Y.-Z. Ma, and D. Hanson, Large-Angle Correlations in the Cosmic Microwave Background, Mon. Not. Roy. Astron. Soc. 407, 2530 (2010), arXiv:0911.5399 [astro-ph.CO] .
  * Copi  _et al._ [2010] C. J. Copi, D. Huterer, D. J. Schwarz, and G. D. Starkman, Large angle anomalies in the CMB, Adv. Astron. 2010, 847541 (2010), arXiv:1004.5602 [astro-ph.CO] .
  * Copi  _et al._ [2015] C. J. Copi, D. Huterer, D. J. Schwarz, and G. D. Starkman, Large-scale alignments from WMAP and Planck, Mon. Not. Roy. Astron. Soc. 449, 3458 (2015), arXiv:1311.4562 [astro-ph.CO] .
  * Copi  _et al._ [2016] C. J. Copi, M. O’Dwyer, and G. D. Starkman, The ISW effect and the lack of large-angle CMB temperature correlations, Mon. Not. Roy. Astron. Soc. 463, 3305 (2016), arXiv:1605.09732 [astro-ph.CO] .
  * Akrami  _et al._ [2020] Y. Akrami _et al._ (Planck), Planck 2018 results. VII. Isotropy and Statistics of the CMB, Astron. Astrophys. 641, A7 (2020), arXiv:1906.02552 [astro-ph.CO] .
  * Calabrese  _et al._ [2008] E. Calabrese, A. Slosar, A. Melchiorri, G. F. Smoot, and O. Zahn, Cosmic Microwave Weak lensing data as a test for the dark universe, Phys. Rev. D 77, 123531 (2008), arXiv:0803.2309 [astro-ph] .
  * Aghanim  _et al._ [2017] N. Aghanim _et al._ (Planck), Planck intermediate results. LI. Features in the cosmic microwave background temperature power spectrum and shifts in cosmological parameters, Astron. Astrophys. 607, A95 (2017), arXiv:1608.02487 [astro-ph.CO] .
  * Motloch and Hu [2020] P. Motloch and W. Hu, Lensinglike tensions in the $Planck$ legacy release, Phys. Rev. D 101, 083515 (2020), arXiv:1912.06601 [astro-ph.CO] .
  * Aiola  _et al._ [2020] S. Aiola _et al._ (ACT), The Atacama Cosmology Telescope: DR4 Maps and Cosmological Parameters, JCAP 12, 047, arXiv:2007.07288 [astro-ph.CO] .
  * Dutcher  _et al._ [2021] D. Dutcher _et al._ (SPT-3G), Measurements of the E-mode polarization and temperature-E-mode correlation of the CMB from SPT-3G 2018 data, Phys. Rev. D 104, 022003 (2021), arXiv:2101.01684 [astro-ph.CO] .
  * Balkenhol  _et al._ [2023] L. Balkenhol _et al._ (SPT-3G), Measurement of the CMB temperature power spectrum and constraints on cosmology from the SPT-3G 2018 TT, TE, and EE dataset, Phys. Rev. D 108, 023510 (2023), arXiv:2212.05642 [astro-ph.CO] .
  * Alam  _et al._ [2021] S. Alam _et al._ (eBOSS), Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory, Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 [astro-ph.CO] .
  * Aubourg  _et al._ [2015] E. Aubourg _et al._ (BOSS), Cosmological implications of baryon acoustic oscillation measurements, Phys. Rev. D 92, 123516 (2015), arXiv:1411.1074 [astro-ph.CO] .
  * Delubac  _et al._ [2015] T. Delubac _et al._ (BOSS), Baryon acoustic oscillations in the Ly$\alpha$ forest of BOSS DR11 quasars, Astron. Astrophys. 574, A59 (2015), arXiv:1404.1801 [astro-ph.CO] .
  * Addison  _et al._ [2018] G. E. Addison, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw, and J. L. Weiland, Elucidating $\Lambda$CDM: Impact of Baryon Acoustic Oscillation Measurements on the Hubble Constant Discrepancy, Astrophys. J. 853, 119 (2018), arXiv:1707.06547 [astro-ph.CO] .
  * Cuceu  _et al._ [2019] A. Cuceu, J. Farr, P. Lemos, and A. Font-Ribera, Baryon Acoustic Oscillations and the Hubble Constant: Past, Present and Future, JCAP 10, 044, arXiv:1906.11628 [astro-ph.CO] .
  * Cortês and Liddle [2024] M. Cortês and A. R. Liddle, On dataset tensions and signatures of new cosmological physics, Mon. Not. Roy. Astron. Soc. 531, L52 (2024), arXiv:2309.03286 [astro-ph.CO] .
  * Efstathiou [2021] G. Efstathiou, To H0 or not to H0?, Mon. Not. Roy. Astron. Soc. 505, 3866 (2021), arXiv:2103.08723 [astro-ph.CO] .
  * Camarena and Marra [2021] D. Camarena and V. Marra, On the use of the local prior on the absolute magnitude of Type Ia supernovae in cosmological inference, Mon. Not. Roy. Astron. Soc. 504, 5164 (2021), arXiv:2101.08641 [astro-ph.CO] .
  * Wen  _et al._ [2021] Y. Wen, D. Scott, R. Sullivan, and J. P. Zibin, Role of $T_{0}$ in CMB anisotropy measurements, Phys. Rev. D 104, 043516 (2021), arXiv:2011.09616 [astro-ph.CO] .
  * Carron  _et al._ [2022] J. Carron, A. Lewis, and G. Fabbian, Planck integrated Sachs-Wolfe-lensing likelihood and the CMB temperature, Phys. Rev. D 106, 103507 (2022), arXiv:2209.07395 [astro-ph.CO] .
  * Ade  _et al._ [2016] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XXIV. Cosmology from Sunyaev-Zeldovich cluster counts, Astron. Astrophys. 594, A24 (2016), arXiv:1502.01597 [astro-ph.CO] .
  * de Haan  _et al._ [2016] T. de Haan _et al._ (SPT), Cosmological Constraints from Galaxy Clusters in the 2500 square-degree SPT-SZ Survey, Astrophys. J. 832, 95 (2016), arXiv:1603.06522 [astro-ph.CO] .
  * Efstathiou and Bond [1999] G. Efstathiou and J. R. Bond, Cosmic confusion: Degeneracies among cosmological parameters derived from measurements of microwave background anisotropies, Mon. Not. Roy. Astron. Soc. 304, 75 (1999), arXiv:astro-ph/9807103 .
  * Kosowsky  _et al._ [2002] A. Kosowsky, M. Milosavljevic, and R. Jimenez, Efficient cosmological parameter estimation from microwave background anisotropies, Phys. Rev. D 66, 063007 (2002), arXiv:astro-ph/0206014 .
  * Ade  _et al._ [2019] P. Ade _et al._ (Simons Observatory), The Simons Observatory: Science goals and forecasts, JCAP 02, 056, arXiv:1808.07445 [astro-ph.CO] .
  * Rosenberg  _et al._ [2022] E. Rosenberg, S. Gratton, and G. Efstathiou, CMB power spectra and cosmological parameters from Planck PR4 with CamSpec, Mon. Not. Roy. Astron. Soc. 517, 4620 (2022), arXiv:2205.10869 [astro-ph.CO] .
  * Tristram  _et al._ [2024] M. Tristram _et al._ , Cosmological parameters derived from the final Planck data release (PR4), Astron. Astrophys. 682, A37 (2024), arXiv:2309.10034 [astro-ph.CO] .
  * Karwal  _et al._ [2022] T. Karwal, M. Raveri, B. Jain, J. Khoury, and M. Trodden, Chameleon early dark energy and the Hubble tension, Phys. Rev. D 105, 063535 (2022), arXiv:2106.13290 [astro-ph.CO] .
  * Lin  _et al._ [2023] M.-X. Lin, E. McDonough, J. C. Hill, and W. Hu, Dark matter trigger for early dark energy coincidence, Phys. Rev. D 107, 103523 (2023), arXiv:2212.08098 [astro-ph.CO] .
  * Sakstein and Trodden [2020] J. Sakstein and M. Trodden, Early Dark Energy from Massive Neutrinos as a Natural Resolution of the Hubble Tension, Phys. Rev. Lett. 124, 161301 (2020), arXiv:1911.11760 [astro-ph.CO] .
  * Carrillo González  _et al._ [2021] M. Carrillo González, Q. Liang, J. Sakstein, and M. Trodden, Neutrino-Assisted Early Dark Energy: Theory and Cosmology, JCAP 04, 063, arXiv:2011.09895 [astro-ph.CO] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publ. Astron. Soc. Pac. 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Hogg and Foreman-Mackey [2018] D. W. Hogg and D. Foreman-Mackey, Data analysis recipes: Using Markov Chain Monte Carlo, Astrophys. J. Suppl. 236, 11 (2018), arXiv:1710.06068 [astro-ph.IM] .
  * Foreman-Mackey  _et al._ [2019] D. Foreman-Mackey, W. Farr, M. Sinha, A. Archibald, D. Hogg, J. Sanders, J. Zuntz, P. Williams, A. Nelson, M. de Val-Borro, T. Erhardt, I. Pashchenko, and O. Pla, Emcee v3: A Python ensemble sampling toolkit for affine-invariant MCMC, Journal of Open Source Software 4, 1864 (2019).
  * Foreman-Mackey [2016] D. Foreman-Mackey, corner.py: Scatterplot matrices in python, Journal of Open Source Software 1, 24 (2016).
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Array programming with NumPy, Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python, Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9, 90 (2007).
  * Hoyer and Hamman [2017] S. Hoyer and J. Hamman, xarray: N-D labeled arrays and datasets in Python, Journal of Open Research Software 5, 10.5334/jors.148 (2017).
  * Kumar  _et al._ [2019] R. Kumar, C. Carroll, A. Hartikainen, and O. Martin, Arviz a unified library for exploratory analysis of bayesian models in python, Journal of Open Source Software 4, 1143 (2019).
  * Meurer  _et al._ [2017] A. Meurer _et al._ , SymPy: symbolic computing in Python, PeerJ Comput. Sci. 3, e103 (2017).
  * van der Velden [2020] E. van der Velden, CMasher: Scientific colormaps for making accessible, informative and ’cmashing’ plots, The Journal of Open Source Software 5, 2004 (2020), arXiv:2003.01069 [eess.IV] .
  * Hu [1998] W. Hu, Structure formation with generalized dark matter, Astrophys. J. 506, 485 (1998), arXiv:astro-ph/9801234 .
  * Hwang and Noh [2009] J.-c. Hwang and H. Noh, Axion as a Cold Dark Matter candidate, Phys. Lett. B 680, 1 (2009), arXiv:0902.4738 [astro-ph.CO] .
  * Marsh and Ferreira [2010] D. J. E. Marsh and P. G. Ferreira, Ultra-Light Scalar Fields and the Growth of Structure in the Universe, Phys. Rev. D 82, 103528 (2010), arXiv:1009.3501 [hep-ph] .
  * Winch  _et al._ [2023] H. Winch, R. Hlozek, D. J. E. Marsh, D. Grin, and K. Rogers, Extreme Axions Unveiled: a Novel Fluid Approach for Cosmological Modeling, arXiv:2311.02052 [astro-ph.CO] (2023).
  * Bean and Dore [2004] R. Bean and O. Dore, Probing dark energy perturbations: The Dark energy equation of state and speed of sound as measured by WMAP, Phys. Rev. D 69, 083503 (2004), arXiv:astro-ph/0307100 .
  * Ballesteros and Lesgourgues [2010] G. Ballesteros and J. Lesgourgues, Dark energy with non-adiabatic sound speed: initial conditions and detectability, JCAP 10, 014, arXiv:1004.5509 [astro-ph.CO] .
  * Ureña López and Gonzalez-Morales [2016] L. A. Ureña López and A. X. Gonzalez-Morales, Towards accurate cosmological predictions for rapidly oscillating scalar fields as dark matter, JCAP 07, 048, arXiv:1511.08195 [astro-ph.CO] .
  * Cookmeyer  _et al._ [2020] T. Cookmeyer, J. Cookmeyer, D. Grin, and T. L. Smith, How sound are our ultralight axion approximations?, Phys. Rev. D 101, 023501 (2020), arXiv:1909.11094 [astro-ph.CO] .
  * Ma and Bertschinger [1995] C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophys. J. 455, 7 (1995), arXiv:astro-ph/9506072 .
  * Ivanov  _et al._ [2020b] M. M. Ivanov, E. McDonough, J. C. Hill, M. Simonović, M. W. Toomey, S. Alexander, and M. Zaldarriaga, Constraining Early Dark Energy with Large-Scale Structure, Phys. Rev. D 102, 103502 (2020b), arXiv:2006.11235 [astro-ph.CO] .
  * Agrawal  _et al._ [2023] P. Agrawal, F.-Y. Cyr-Racine, D. Pinner, and L. Randall, Rock ‘n’ roll solutions to the Hubble tension, Phys. Dark Univ. 42, 101347 (2023), arXiv:1904.01016 [astro-ph.CO] .
  * Planck Collaboration [2018a] Planck Collaboration, clik (2018a), commit hash 8eaebd5.
  * Planck Collaboration [2018b] Planck Collaboration, Planck Legacy Archive (2018b), baseline release 3.00.
  * Alam  _et al._ [2017] S. Alam _et al._ (BOSS), 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, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * SDSS collaboration [2020] SDSS collaboration, SDSS-IV SVN Software Repository, Data Release 16 (2020), accessed 2023.
  * DESI collaboration [2024] DESI collaboration, DESI DR1 dataset (2024), commit hash af68b10.
  * Pantheon collaboration [2018] Pantheon collaboration, Pantheon data release repository (2018), commit hash 7eb29dc.
  * Pantheon and SH0ES collaborations [2022] Pantheon and SH0ES collaborations, Pantheon+SH0ES DataRelease repository (2022), commit hash 1083557.
  * DES collaboration [2024] DES collaboration, DES-SN5YR repository (2024), commit hash f254460.
  * Rubin, D. et. al. [2024] Rubin, D. et. al., Union3 dataset (2024), commit hash 261e356.
  * Farr  _et al._ [2014] B. Farr, V. Kalogera, and E. Luijten, A more efficient approach to parallel-tempered Markov-chain Monte Carlo for the highly structured posteriors of gravitational-wave signals, Phys. Rev. D 90, 024014 (2014), arXiv:1309.7709 [astro-ph.IM] .
  * Farr and Farr [2015] B. Farr and W. M. Farr, kombine: a kernel-density-based, embarrassingly parallel ensemble sampler (2015), in prep.