# Dark Radiation Isocurvature from Cosmological Phase Transitions

Matthew R. Buckley  NHETC, Department of Physics and Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  Peizhi Du  NHETC, Department of Physics
and Astronomy, Rutgers University, Piscataway, NJ 08854, USA  Nicolas
Fernandez  NHETC, Department of Physics and Astronomy, Rutgers University,
Piscataway, NJ 08854, USA  Mitchell J. Weikert  NHETC, Department of Physics
and Astronomy, Rutgers University, Piscataway, NJ 08854, USA

###### Abstract

Cosmological first order phase transitions are typically associated with
physics beyond the Standard Model, and thus of great theoretical and
observational interest. Models of phase transitions where the energy is mostly
converted to dark radiation can be constrained through limits on the dark
radiation energy density (parameterized by $\Delta N_{\rm eff}$). However, the
current constraint ($\Delta N_{\rm eff}<0.3$) assumes the perturbations are
adiabatic. We point out that a broad class of non-thermal first order phase
transitions that start during inflation but do not complete until after
reheating leave a distinct imprint in the scalar field from bubble nucleation.
Dark radiation inherits the perturbation from the scalar field when the phase
transition completes, leading to large-scale isocurvature that would be
observable in the CMB. We perform a detailed calculation of the isocurvature
power spectrum and derive constraints on $\Delta N_{\rm eff}$ based on CMB+BAO
data. For a reheating temperature of $T_{\rm rh}$ and a nucleation temperature
$T_{*}$, the constraint is approximately $\Delta N_{\rm eff}\lesssim
10^{-5}(T_{*}/T_{\rm rh})^{-4}$, which can be much stronger than the adiabatic
result. We also point out that since perturbations of dark radiation have a
non-Gaussian origin, searches for non-Gaussianity in the CMB could place a
stringent bound on $\Delta N_{\rm eff}$ as well.

##  I Introduction

Extensions of the Standard Model of particle physics often result in first
order phase transitions (FOPT) in the early Universe. For example, the phase
transition (PT) in models of electroweak baryogenesis is predicted to be first
order Kuzmin  _et al._ (1985); Shaposhnikov (1986, 1987); Cohen  _et al._
(1993); Trodden (1999); Riotto and Trodden (1999); Dine and Kusenko (2003);
Morrissey and Ramsey-Musolf (2012). New physics that postulates new confining
gauge groups – such as in composite Higgs/Randall Sundrum models which address
the Planck-Weak Hierarchy problem Creminelli  _et al._ (2002); Randall and
Servant (2007); Nardini  _et al._ (2007); Konstandin  _et al._ (2010);
Konstandin and Servant (2011); Bunk  _et al._ (2018); Baratella  _et al._
(2019); Bruggisser  _et al._ (2018); Megías  _et al._ (2018); Agashe  _et al._
(2020); Fujikura  _et al._ (2020); Agashe  _et al._ (2021); Ares  _et al._
(2020); Agrawal and Nee (2021); Levi  _et al._ (2023); Csáki  _et al._ (2023);
Eröncel  _et al._ (2023); Mishra and Randall (2023, 2024), or in dark matter
models with dark QCD or dark $SU(N)$ gauge theories Bai and Schwaller (2014);
Kribs and Neil (2016); Cline (2022); Shelton and Zurek (2010); Dutta and Kumar
(2011); Holthausen  _et al._ (2013); Schwaller (2015); Hall  _et al._ (2023);
Bai and Korwar (2022); Bottaro  _et al._ (2021); Pasechnik  _et al._ (2023,
2023) – can also result in FOPT. A class of Early Dark Energy models Karwal
and Kamionkowski (2016); Mörtsell and Dhawan (2018); Poulin  _et al._ (2019),
introduced to relieve the tension between the early- and late-time
measurements of $H_{0}$ Aghanim  _et al._ (2020); Riess  _et al._ (2022),
generally require PTs in their dark sector Niedermann and Sloth (2021, 2020,
2022). Probing cosmological FOPTs is therefore a problem of wide applicability
to physics beyond the Standard Model.

Cosmological FOPTs create gravitational waves through bubble collisions, shock
waves, and turbulence produced during the transition Kosowsky  _et al._
(1992a, b); Kosowsky and Turner (1993); Kamionkowski  _et al._ (1994); Caprini
_et al._ (2016, 2020); Caldwell  _et al._ (2022). For phase transitions at
${\cal O}(10-100)~{}{\rm MeV}$ temperatures, these gravitational waves have
frequencies accessible to the Pulsar Timing Array (PTA) Agazie  _et al._
(2023); Afzal  _et al._ (2023); Antoniadis  _et al._ (2023); Reardon  _et al._
(2023); Xu  _et al._ (2023). Future experiments, including LISA Amaro-Seoane
_et al._ (2017); Auclair  _et al._ (2023), DECIGO Kawamura  _et al._ (2006,
2021), and BBO Harry  _et al._ (2006), will be able to probe transition
temperature around the TeV scale or above.

In addition to these direct probes of the stochastic gravitational wave
background generated by the FOPT, precision measurements of the Cosmic
Microwave Background (CMB) can constrain phase transitions that generate
massless or low-mass dark radiation (DR) Bai and Korwar (2022).111In this
context, gravitational waves are also a form of dark radiation. The energy
density of DR is typically parameterized in terms of the extra number of
effective neutrino species $\Delta N_{\rm eff}$ above the Standard Model
expectation $N_{\rm eff}=3.044$ Froustey  _et al._ (2020); Bennett  _et al._
(2021); Akita and Yamaguchi (2020). For free-streaming DR, the CMB sets the
upper bound $\Delta N_{\rm eff}<0.3$ at 95% confidence level (CL) Aghanim  _et
al._ (2020).

Importantly, this constraint on the number of light degrees of freedom assumes
the initial density perturbations of the dark radiation are adiabatic: with
over- and under-densities proportional to those of the other radiation
species. If the DR is sourced by a FOPT that undergoes violent bubble
nucleation, the initial spatial distribution of the perturbations will contain
additional isocurvature modes.222For general discussions of isocurvature, see
for example Bucher  _et al._ (2000); Wands  _et al._ (2000); Gordon  _et al._
(2000); Lyth  _et al._ (2003); Malik and Wands (2005); Wands (2008). However,
for FOPTs with energy scales far above the $\sim{\rm eV}$ scales of
recombination and decoupling, the isocurvature perturbations are typically
highly suppressed for wavenumbers large enough to leave observable imprints
within the CMB Freese and Winkler (2023); Elor  _et al._ (2023). As a result –
even though isocurvature perturbations are tightly constrained by CMB
measurements Akrami  _et al._ (2020a); Ghosh  _et al._ (2022)333There are also
constraints on dark radiation isocurvature from Big Bang nucleosynthesis
Adshead  _et al._ (2020). – these limits typically cannot be applied to the
high-scale FOPTs, and the dominant constraint comes from the standard
measurements on the number of relativistic degrees of freedom.

In this paper, we consider a broad class of high-scale FOPT models which
generate large-scale isocurvature perturbations whose imprint would be visible
within the CMB. These models have non-thermal (vacuum) phase transitions,
where the scalar potential that triggers the transition is independent of
temperature and couples to the inflaton field. As a result of this coupling,
the FOPT can begin during inflation, with a bubble nucleation rate per volume
$\Gamma_{\rm PT}$ which we assume to satisfy

$\Gamma_{\rm PT}\ll H_{\rm inf}^{4},$ (1)

during the inflationary epoch (when the Hubble parameter is $H_{\rm inf}$).
This hierarchy of scales can be chosen so that the PT does not complete by the
end of inflation and instead completes only after reheating. In the scenario
we consider in this paper, the energy density in the scalar field undergoing
the FOPT is mostly converted into dark radiation when the PT completes, with
perturbations inherited from the scalar field.

Bubbles nucleated in the scalar field during inflation will expand to the
horizon size in a Hubble time and remain frozen in comoving coordinates
afterwards. As the comoving horizon is rapidly changing during inflation, this
results in a wide range of bubble sizes, with the largest bubbles
corresponding to the earliest nucleation events. The bubble distribution will
later be imprinted in DR when the PT completes. A schematic of cosmological
history during and after inflation is shown in Figure 1. For phase transitions
that start when the inflationary horizon is of order the comoving scale at
recombination or larger, stochastic perturbations in DR can leave an
observable imprint on the CMB.444For earlier works on effects of large bubbles
from phase transitions during inflation, see Refs. Guth and Weinberg (1983);
Liddle and Wands (1991); Turner  _et al._ (1992); Copeland  _et al._ (1994);
Baccigalupi and Perrotta (2000); Barir  _et al._ (2023). As the bubbles are
randomly distributed, the DR perturbations are uncorrelated with the
fluctuations of the inflaton and behave as isocurvature. Since the CMB is
consistent with adiabatic initial conditions with zero isocurvature Akrami
_et al._ (2020a), we can place constraints on $\Delta N_{\rm eff}$ from this
class of FOPTs that are much stronger than in the adiabatic scenario.

Figure 1: A schematic plot of comoving horizon size $(aH)^{-1}$ (black line)
as a function of scale factor $a$. Three scale factors have been listed:
$a_{i}$ denotes the start of the FOPT, $a_{e}$ shows the end of inflation, and
$a_{*}$ corresponds to the nucleation temperature $T_{*}$ when phase
transition completes. The shaded orange region shows the time window when
modes that CMB can probe enter the horizon. The upper part of the plot shows
the evolution of perturbations in the $\chi$ field due to nucleation of
bubbles. $r_{i}\equiv(a_{i}H_{\rm inf})^{-1}$ is the size of the earliest and
largest bubble, nucleated at $a_{i}$. After the completion of the PT around
$a_{*}$, $\chi$ converts into DR with the same large-scale perturbations,
which will imprint DR isocurvature signals on the CMB after re-entering the
horizon. (We note that the distribution of bubbles shown in this plot are
schematic. We will present more realistic bubble distributions from
simulations in Section III.)

In this work, we present a detailed calculation of the DR isocurvature power
spectrum from non-thermal PTs that start during inflation. We implement this
model in CLASS Blas  _et al._ (2011); Lesgourgues (2011); Lesgourgues and Tram
(2011), appropriately modified to include DR isocurvature. Using numerical
simulations, we show how the CMB power spectra would change in the presence of
DR isocurvature. We then perform a Markov Chain Monte Carlo (MCMC) scan using
CMB and Baryon Acoustic Oscillation (BAO) data to constrain the DR energy
density (written in terms of $\Delta N_{\rm eff}$) and other parameters
related to the phase transition. We find that when the temperature $T_{*}$ at
which the phase transition completes (estimated from $\Gamma_{\rm
PT}=H(T_{*})^{4}$), is close to the reheating temperature $T_{\rm rh}$, the
isocurvature perturbations set a limit of $\Delta N_{\rm eff}\lesssim
10^{-5}(T_{*}/T_{\rm rh})^{-4}$. This limit weakens as $T_{*}/T_{\rm rh}$
decreases and approaches the adiabatic constraint of $\Delta N_{\rm eff}<0.3$
when $T_{*}\ll T_{\rm rh}$. In addition to the isocurvature, a FOPT which
begins during inflation will also result in non-Gaussianity in the CMB. While
we defer a full analysis of such constraints to a future work, here we
estimate the expected limits by analogy with neutrino density isocurvature
non-Gaussianity Akrami  _et al._ (2020b); Montandon  _et al._ (2021).

This paper is organized as follows. In Section II, we introduce a model of
non-thermal phase transitions and calculate the transition rate. In Section
III we describe the dynamics of the phase transition during inflation and the
resulting distribution of nucleated bubbles. In Section IV, we calculate the
two quantities that will lead to observational constraints on the model: the
DR isocurvature power spectrum and bispectrum. Next, we implement the model in
CLASS and show the effects of DR isocurvature on the CMB angular power spectra
in Section V. In Section VI we describe the datasets we use to set our limits
and the details of our MCMC parameter scans. The resulting constraints on the
DR isocurvature parameters from the CMB angular power spectra and bispectra
are shown in Section VII. We conclude in Section VIII.

##  II Models of Non-thermal Phase Transitions

In this work, we consider models containing a scalar field $\chi$ that
undergoes a FOPT that imprints significant isocurvature perturbations into the
CMB. The minimum requirements for such models are:

  * i.

The FOPT starts during inflation, at a moment when the comoving horizon
($r_{i}$) is above the minimum comoving length scale that the CMB can probe.

  * ii.

The PT rate remains small during inflation ($\gamma_{\rm PT}\equiv\Gamma_{\rm
PT}/H_{\rm inf}^{4}\ll 1$).

  * iii.

After reheating, the PT completes (at a temperature $T_{*}$, defined by
$\Gamma_{\rm PT}/H(T_{*})^{4}=1$) and the energy density in the $\chi$ field
converts into DR.

The interesting isocurvature effects we consider in this paper are present as
long as these three requirements are met.

As a specific implementation, we realize these features with a scalar
potential that is non-thermal – that is, independent of the temperature $T$ of
the Universe. This non-thermality can be simply achieved if there are no
sizeable couplings between the scalar field $\chi$ and the Standard Model
fields. Without temperature dependence in the potential, the parameters of the
model can be chosen such that the phase transition rate $\Gamma_{\rm PT}$ is
constant and satisfies $\Gamma_{\rm PT}\ll H_{\rm inf}^{4}$. The phase
transition in our study is triggered by the slow roll of the inflaton field
$\phi$, through couplings between $\phi$ and $\chi$. We further assume that
the energy density of the $\chi$ field is subdominant to that of $\phi$ during
inflation, and couplings between $\chi$ and other light fields in the dark
sector allow the eventual conversion of vacuum energy into dark radiation
after the phase transition completes. The example potential we consider is
described in Section II.1 along with the resulting phase transition rate. The
phase transition trigger during inflation is discussed in Section II.2.

Figure 2: Schematic plot of the scalar potential $V(\chi,\phi)$ in Eq. (6).
When $\phi<\phi_{i}$, $\chi$ stays in the true vacuum near $\chi_{+}$. Soon
after $\phi>\phi_{i}$, it turns into the false vacuum because $\mu(\phi)$
changes sign. Therefore, the FOPT starts and $\chi$ tunnels to true vacuum
near $\chi_{-}$.

###  II.1 Temperature-Independent Phase Transition

The temperature-independence of the scalar potential $V(\chi)$ can be achieved
by requiring that $V$ is dominated by a bare tree-level potential. For this
paper, we consider the simple example:

$\displaystyle
V(\chi)=-\frac{1}{2}m^{2}\chi^{2}+\frac{\mu}{3}\chi^{3}+\frac{\lambda}{4}\chi^{4}.$
(2)

The $Z_{2}$ symmetry of the potential is broken by $\mu$, which we assume to
be small ($\mu\ll m$) and positive. This condition implies that the two vacuum
states at $\chi_{\pm}\approx\pm m/\sqrt{\lambda}$ are nearly degenerate. That
is, the vacuum energy difference $\Delta V\sim\mu m^{3}/\lambda^{3/2}$ is much
smaller than the potential barrier $V_{\rm max}\sim m^{4}/\lambda$.

We assume the field is in the false vacuum state $\chi_{+}$ early in the
inflationary phase. When the comoving horizon is $r_{i}$, the FOPT starts as
$\chi$ tunnels to the true vacuum $\chi_{-}$. In Section II.2, we present one
example of a dynamical mechanism for this transition, achieved by coupling
$\chi$ and the inflaton $\phi$ through the $Z_{2}$-breaking $\mu$ term.

For $\Delta V\ll V_{\rm max}$, the thin-wall approximation is applicable and
the quantum tunneling rate can be calculated in terms of the $O(4)$ symmetric
bounce action Linde (1983); Coleman (1977); Callan and Coleman (1977)

$\displaystyle S_{4}=\frac{27\pi^{2}S_{1}^{4}}{2\Delta
V^{3}}\sim\lambda^{1/2}\left(\frac{m}{\mu}\right)^{3}\gg 1,$ (3)

where $S_{1}\sim m^{3}/\lambda$ is the bounce action along the radial
direction. The quantum phase transition rate per volume is

$\displaystyle\Gamma_{\rm PT}\sim r_{c}^{-4}S_{4}^{2}e^{-S_{4}},$ (4)

where $r_{c}=3S_{1}/\Delta V$ is the size of the critical bubble for the
$O(4)$ symmetric bubble under the thin-wall approximation. Since $\Gamma_{\rm
PT}$ is exponentially sensitive to $S_{4}$, one can easily achieve the
requisite small $\Gamma_{\rm PT}/H_{\rm inf}^{4}$ for a reasonable choice of
parameters. Under these assumptions, and with the additional requirement that
the number of $e$-folds after the tunneling begins is small enough (see
Section III for details), the phase transition does not complete by the end of
inflation. That is, at the end of inflation there are regions where the scalar
field is in the false vacuum, separating bubbles of true vacuum nucleated
during inflation.

After inflation, we assume the Universe instantaneously reheats to a
radiation-dominated state with a temperature $T_{\rm rh}$. We assume
$\Gamma_{\rm PT}$ remains constant after inflation, and that the PT completes
at a temperature $T_{*}$. As the horizon changes after inflation, a constant
$\Gamma_{\rm PT}$ requires that the critical bubble size $r_{c}$ is less than
the physical horizon size at any time. Otherwise, the bounce action will
acquire horizon dependent corrections from the curved spacetime Coleman and De
Luccia (1980). That is, we require

$\displaystyle r_{c}^{-1}\sim\frac{\mu}{\lambda^{1/2}}\gg H_{\rm inf}.$ (5)

This and all other requirements stated in this section can be satisfied by
choosing appropriate values of $m$, $\lambda$, $\mu$, and $H_{\rm inf}$.

###  II.2 Inflaton-Triggered Phase Transition

Our scenario requires a phase transition that is triggered during inflation,
rather than by a coupling to a thermalized Standard Model bath. While more
complex examples might be motivated by other theoretical considerations, a
simple realization of this phenomenon is to couple the $\chi$ field to the
$\phi$ through the small $Z_{2}$-breaking $\mu$ term in Eq. (2):

$\displaystyle
V(\chi,\phi)=-\frac{1}{2}m^{2}\chi^{2}+\frac{\mu(\phi)}{3}\chi^{3}+\frac{\lambda}{4}\chi^{4}.$
(6)

To trigger the FOPT, $\mu(\phi)$ must be some function of $\phi$ that switches
sign as $\phi$ slow-rolls past some critical value of $\phi_{i}$. To simplify
the later discussion, we focus on a specific form of $\mu(\phi)$ that behaves
like a step function when $\phi$ passes $\phi_{i}$. One simple realization of
such a $\mu(\phi)$ function is

$\displaystyle\mu(\phi)=\mu\,\tanh\left(\frac{\phi-\phi_{i}}{\Delta\phi}\right)\,,$
(7)

where the mass scale $\Delta\phi$ is chosen to be much smaller than change of
$\phi$ in a Hubble time. We also require $\mu$ to be small such that the
$\tfrac{1}{3}\mu(\phi)\chi^{3}$ term has a negligible contribution to the
inflaton potential. We note that the isocurvature signals we consider in this
work are insensitive to the specific choice of $\mu(\phi)$, provided that it
behaves like a step function with a quick transition.

For choices of $\mu(\phi)$ that have the necessary temperature-independent
trigger behavior, this model has the novel feature that the phase transition
begins during inflation when $\phi$ rolls past a critical value $\phi_{i}$.
When this happens, $\chi_{+}$ is no longer the global minimum of the potential
and a FOPT can occur as $\chi$ tunnels to the true minimum at $\chi_{-}$ (see
Figure 2). This tunneling occurs stochastically; as bubbles of true vacuum
form at random times and locations, they expand to the horizon size and
inflate along with it (as we will discuss in detail in the next section).
Therefore, the largest bubbles after inflation are those that nucleated
earliest, when $\phi\approx\phi_{i}$ and the comoving horizon size was
$r_{i}$. Therefore, the distribution of bubbles has a characteristic feature
around $r_{i}$. In our simple scenario, the PT rate remains constant and the
phase transition can be fully specified by $\gamma_{\rm PT}\equiv\Gamma_{\rm
PT}/H_{\rm inf}^{4}$ and $r_{i}$.

We emphasize that all the interesting large-scale isocurvature effects
considered in this paper are present as long as the three requirements stated
at the beginning of this section are satisfied. To demonstrate the resulting
constraints from the CMB data, we restrict ourselves to the simple toy model
parameterized by $\gamma_{\rm PT}$ and $r_{i}$.

##  III Stochastic Bubble Distribution from Incomplete Phase Transition during
inflation

In this section, we discuss the dynamics of bubbles nucleated during inflation
and describe their stochastic distribution. From the distribution of bubbles,
we will calculate density perturbations in the $\chi$ field, which will lead
to DR isocurvature with non-Gaussian statistics discussed in Section IV.

For this calculation, we work in the homogeneous inflationary background,
neglecting the effect of metric perturbations on the dynamics of the $\chi$
field during inflation. This approximation is valid if
$\bar{\rho}_{\chi}/(\bar{\rho}_{\phi}+\bar{\rho}_{\chi})\ll 1$ and
$\delta\rho_{\chi}/\bar{\rho}_{\chi}\gg\delta\rho_{\phi}/\bar{\rho}_{\phi}$.
Here $\bar{\rho}$ denotes the averaged energy density of the corresponding
field and $\delta\rho$ its perturbation. These conditions will be satisfied
for all model parameters considered in this work that lead to constraints
driven by DR isocurvature. Under this approximation, each bubble will expand
with spherical symmetry.

The bubble wall dynamics are calculated in detail in Appendix A: for the model
parameters and approximations under which we are working, within a Hubble time
the wall accelerates to a terminal velocity which is close to the speed of
light. Given the short acceleration time, we can take the bubble wall velocity
to be $v_{w}\approx 1$ immediately after nucleation. Under this approximation,
a bubble that nucleates at time $t^{\prime}$ has a comoving radius at time
$t$:

$\begin{split}r(t,t^{\prime})=&\frac{1}{H_{\rm
inf}a(t^{\prime})}\left(1-\frac{a(t^{\prime})}{a(t)}\right)+\frac{r_{c}}{a(t^{\prime})}\\\
\approx&\frac{1}{H_{\rm inf}a(t^{\prime})}.\end{split}$ (8)

Here we drop the second term (as $a(t^{\prime})/a(t)$ decays exponentially
with time) and the third term (as $r_{c}\ll H_{\rm inf}^{-1}$ from Eq. (5)).
After neglecting these terms, the radius only depends on the nucleation time
$t^{\prime}$ and so we will omit the $t$ argument: writing $r(t,t^{\prime})\to
r(t^{\prime})$. This means that the bubble will quickly expand to the size of
the comoving horizon at the time of nucleation. As the comoving horizon
decreases during inflation, the bubble becomes super-horizon soon after
nucleation and is fixed in size during the later evolution of the Universe.

Figure 3: Simulated bubble distribution at different times within a cube with
side-length $5r_{i}$ for $\gamma_{\rm PT}=1\times 10^{-4}$. The simulation
starts at $t=0$ and $r_{i}$ is the hypothetical radius of a bubble nucleated
at $t=0$. The bubble distribution is projected to the $x$-$y$ plane with the
colors representing the $z$ coordinate of the center of each bubble.

As the nucleation of bubbles is a stochastic process, the location and time of
nucleation events are randomly distributed. Analytically, we can calculate the
probability that a point in space remains in the false vacuum at a given time.
Without loss of generality, we set $t=0$ to be the moment when the inflaton
rolls past the trigger value of $\phi_{i}$ (see Section II.2) and bubbles can
start nucleating. For a point $(\mathbf{x},t)$ to be in false vacuum, a bubble
cannot have nucleated at any time $t^{\prime}<t$ in any location within a
sphere of radius $r(t^{\prime})$ centered at $\mathbf{x}$. In terms of the
Hubble-normalized phase transition rate $\gamma_{\rm PT}\equiv\Gamma_{\rm
PT}/H_{\rm inf}^{4}$, the expected number of nucleation events within this
region of space-time is

$\begin{split}I(t)=&\frac{4\pi}{3}\Gamma_{\rm
PT}\int\displaylimits_{0}^{t}dt^{\prime}a(t^{\prime})^{3}r(t^{\prime})^{3}\\\
=&\frac{4\pi}{3}\gamma_{\rm PT}H_{\rm inf}t\,.\end{split}$ (9)

As a result, the probability of false vacuum at $(\mathbf{x},t)$ is

$p_{\rm false}(t)=e^{-I(t)}=e^{-t/\tau_{\rm PT}},$ (10)

where

$\tau_{\rm PT}^{-1}\equiv\frac{4\pi}{3}\gamma_{\rm PT}H_{\rm inf}\,.$ (11)

Outside of the bubbles the $\chi$ field has an energy density of $\Delta V$,
while inside the energy density is zero. Thus, the spatially averaged energy
density of the $\chi$ field at time $t$ is

$\bar{\rho}_{\chi}(t)=\Delta Vp_{\rm false}(t)=\Delta Ve^{-t/\tau_{\rm PT}}.$
(12)

We are interested in the scenario where the phase transition remains
incomplete when inflation ends at time $t_{e}$. This occurs if
$t_{e}/\tau_{\rm PT}<1$. In that case, the $\chi$ field will have a
distribution of mostly non-overlapping bubbles with radii between $r_{e}\equiv
r(t_{e})$ and $r_{i}\equiv r(0)$. In this scenario, we can approximate the
energy density of the $\chi$ field during inflation as

$\rho_{\chi}(\mathbf{x},t)=\Delta
V\left[1-\sum_{I:t_{I}<t}\Theta(r_{I}-|\mathbf{x}-\mathbf{x_{I}}|)\right],$
(13)

where $r_{I}\equiv r(t_{I})$ and $\mathbf{x}_{I}$ are the comoving radius and
center location of bubble $I$, respectively, and the sum runs over all bubbles
nucleated before $t$.

To illustrate the time evolution of the distribution of bubbles, we perform a
simulation of bubble nucleation for an incomplete FOPT during inflation. We
simulate a cubic region of space with comoving volume $(50\,r_{i})^{3}$ and
$\gamma_{\rm PT}=10^{-4}$. We start our simulation at $t=0$, nucleating
bubbles at a rate given by

$\frac{dN}{dt}=\mathcal{V}_{\rm false}(t)a(t)^{3}\Gamma_{\rm PT},$ (14)

where $\mathcal{V}_{\rm false}(t)$ is the volume of space in the false vacuum
at $t$. When a bubble is nucleated, the location of its center is sampled
uniformly from the region of space that is in the false vacuum. We show the
resulting spatial distribution of bubbles at three different times in Figure
3. For illustrative purposes, we show only the $(5r_{i})^{3}$ volume that
contains the first (and largest) bubble nucleated in the original volume.

Next, we calculate the density contrast parameter at the end of inflation,
which is defined as

$\displaystyle\delta_{\chi}(\mathbf{x},t_{e})$ $\displaystyle\equiv$
$\displaystyle\frac{\rho_{\chi}(\mathbf{x},t_{e})-\bar{\rho}_{\chi}(t_{e})}{\bar{\rho}_{\chi}(t_{e})}.$
(15)

In Fourier space this is

$\delta_{\chi}(\mathbf{k},t_{e})=-e^{t_{e}/\tau_{\rm
PT}}\frac{4\pi}{k^{3}}\sum_{I:t_{I}<t_{e}}e^{-i\mathbf{k}\cdot\mathbf{x_{I}}}\mathcal{A}(kr_{I})+{\cal
C}\times\delta^{3}(\mathbf{k}),$ (16)

where $\mathcal{A}(y)\equiv\sin{y}-y\cos{y}$. We ignore the term proportional
to the delta function in Eq. (16) in our analysis, as the $\mathbf{k}=0$
Fourier mode will not contribute to CMB observables with multipoles $\ell>0$.

The two point function of $\delta_{\chi}$ is then given by

$\displaystyle\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle$
$\displaystyle=$ $\displaystyle e^{2t_{e}/\tau_{\rm
PT}}\frac{(4\pi)^{2}}{{k^{\prime}}^{3}k^{3}}$
$\displaystyle\times\sum_{IJ}\langle
e^{-i(\mathbf{k}\cdot\mathbf{x_{I}}+\mathbf{k^{\prime}}\cdot\mathbf{x_{J}})}\mathcal{A}(kr_{I})\mathcal{A}(k^{\prime}r_{J})\rangle.$

The double sum can be split into terms where $I=J$ and those where $I\neq J$.
We will neglect the terms where $I\neq J$, as their sum is suppressed by a
factor $\gamma_{\rm PT}$ relative to the sum over $I=J$ (see Appendix B). We
can replace the ensemble average of the $I=J$ terms with an integral over
bubble coordinates:

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle=e^{2t_{e}/\tau_{PT}}\frac{(4\pi)^{2}}{k^{3}{k^{\prime}}^{3}}N\int
d^{4}x\,p_{1}(x)e^{-i(\mathbf{k}+\mathbf{k^{\prime}})\cdot\bm{x}}\mathcal{A}(kr(t))\mathcal{A}(k^{\prime}r(t))\,.$
(18)

Here, $p_{1}(x)$ is the single-bubble probability distribution function over
nucleation coordinates $x=(\mathbf{x},t)$, and is given by

$p_{1}(x)=\frac{1}{\mathcal{V}}\left(\frac{1}{N}\frac{dN}{dt}\right),$ (19)

where $N$ is the total number of bubbles nucleated in a comoving volume
$\mathcal{V}$ before $t_{e}$ and $dN/dt$ is given by Eq. (14) with

$\mathcal{V}_{\rm false}(t)=\mathcal{V}p_{\rm
false}(t)=\mathcal{V}e^{-t/\tau_{\rm PT}}.$ (20)

Figure 4: Power Spectra of the density contrast of $\chi$ normalized by
$\gamma_{\rm PT}$ for different values of $\gamma_{\rm PT}$ and $H_{\rm
inf}t_{e}=50$. The largest value of $\gamma_{\rm PT}$ is chosen such that the
incomplete PT condition $4\pi\gamma_{\rm PT}H_{\rm inf}t_{e}/3<1$ is barely
satisfied. The dashed line is a broken power law with the same asymptotic
behavior as the power spectrum in the limit of $\gamma_{\rm PT}H_{\rm
inf}t_{e}\ll 1$ (see Eq. (22)) .

After performing the $\mathbf{x}$ integral of Eq. (18), the two point function
can be written in terms of the dimensionless power spectrum $P_{\chi}(k)$ as

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle\equiv
2\pi^{2}\frac{(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k^{\prime}})}{k^{3}}P_{\chi}(k).$
(21)

where

$\displaystyle P_{\chi}(k)$ $\displaystyle=$ $\displaystyle\,8\gamma_{\rm
PT}e^{(8\pi\gamma_{\rm PT}H_{\rm inf}t_{e})/3}(kr_{i})^{-4\pi\gamma_{\rm
PT}/3}$ (22)
$\displaystyle\times\int\displaylimits_{kr_{e}}^{kr_{i}}du\,u^{-4+4\pi\gamma_{\rm
PT}/3}\mathcal{A}(u)^{2}$ $\displaystyle\approx$ $\displaystyle\,\gamma_{\rm
PT}\left\\{\begin{array}[]{cc}\frac{8}{27}(kr_{i})^{3}&\quad kr_{i}\ll 1\\\
4.2&\quad kr_{i}\gg 1\end{array}\right.$ (25)

to leading order of $\gamma_{\rm PT}$ in the limit $\gamma_{\rm PT}H_{\rm
inf}t_{e}\ll 1$. The result for small $kr_{i}$ is calculated analytically by
Taylor-expanding the integrand and the result for large $kr_{i}$ is found by
numerically calculating the limit as $kr_{i}\to\infty$.555We numerically
integrate Eq. (22) for a series of values of $kr_{i}$ until the integral
converges to a constant asymptotic value at 1% precision. As there are no
bubbles larger than $r_{i}$, the $k^{3}$ dependence when $kr_{i}\ll 1$ is a
consequence of the fact that the two-point function of $\delta_{\chi}$ is
constant for modes larger than $r_{i}$. The flat power spectrum when
$kr_{i}\gg 1$ can be understood as follows: when the nucleation rate is low,
the majority of the Universe is in the false vacuum at all times $0<t<t_{e}$,
and thus the rate of bubble nucleation in any physical volume of fixed size is
the same at all times resulting in a scale invariant distribution of bubbles.
In Figure 4, we show numerical results for the power spectrum as a function of
$k$ for various values of $\gamma_{\rm PT}$ along with a broken power law with
the same asymptotic behavior as Eq. (22). The broken power law provides a good
approximation to the power spectrum provided $\gamma_{\rm PT}H_{\rm
inf}t_{e}\ll 1$. Therefore, assuming $H_{\rm inf}t_{e}\sim 50$, the
approximation is valid for $\gamma_{\rm PT}\lesssim 10^{-3}$.

##  IV Imprints of a non-thermal FOPT

After inflation and reheating, we assume the Universe is radiation dominated.
As the Hubble parameter drops as $\propto T^{2}$, $\Gamma_{\rm PT}/H(T)^{4}$
will increase rapidly and reach unity at the nucleation temperature $T_{*}$,
when the phase transition completes. We assume the energy density in the
$\chi$ field will be converted to DR at $T_{*}$ by couplings between DR and
the $\chi$ field. As a result, the DR will inherit the $\chi$ field’s
inhomogeneities, which, on large scales, follow the stochastic distribution of
bubbles nucleated during inflation. This process will generate a relative
entropy perturbation between DR and photons as well as a curvature
perturbation which is distinct from the curvature perturbation produced by
quantum fluctuations of the inflaton. The relative entropy perturbation will
contribute to the DR isocurvature mode, while both curvature perturbations
will contribute to the adiabatic mode. The isocurvature mode also has a non-
negligible 3-point function and therefore has non-Gaussian statistics.

We begin this section by calculating the angular power spectrum of CMB
temperature anisotropies in terms of the adiabatic and isocurvature power
spectra. We then identify the super-horizon initial conditions for the DR
isocurvature and adiabatic modes and calculate the power spectra. We show in
Appendix C that – for the class of models considered in this work – the FOPT-
induced curvature Barir  _et al._ (2023) has a negligible impact on
cosmological observables compared to the DR isocurvature. The DR isocurvature
therefore sets the relevant constraints on these models of FOPTs. Finally, we
calculate the isocurvature bispectrum from the non-Gaussian statistics.

In the synchronous gauge, the perturbed FRW metric is written as

$ds^{2}=a(\tau)^{2}\left[-d\tau^{2}+(\delta_{ij}+h_{ij})dx^{i}dx^{j}\right],$
(26)

where $\tau$ is the conformal time.666While it was natural to use physical
time $t$ during inflation, the initial conditions for CLASS are most
conveniently written in terms of $\tau$. Therefore, all functions that
directly enter the CLASS calculation will be written in terms of $\tau$ rather
than $t$. The metric perturbation $h_{ij}$ can be written in Fourier space as

$h_{ij}(\mathbf{k},\tau)=\left[\hat{k}_{i}\hat{k}_{j}h(\mathbf{k},\tau)+\left(\hat{k}_{i}\hat{k}_{j}-\frac{1}{3}\delta_{ij}\right)6\eta(\mathbf{k},\tau)\right],$
(27)

where $h$ and $\eta$ denote the trace and traceless longitudinal part of
$h_{ij}$ respectively.

In general, we can split a perturbation
$X\in\\{h,\eta,\delta_{a},\theta_{a},\sigma_{a},...\\}$ into a sum of
isocurvature and adiabatic modes. Here $\delta_{a}$, $\theta_{a}$, and
$\sigma_{a}$ correspond to the $\ell=0,1$, and $2$ moments of the phase-space
perturbation of species $a$. In particular,
$\delta\equiv\delta\rho/\bar{\rho}$. For radiation species, moments $\ell>2$
are needed, and it is useful to define $F_{a,\ell}$ as the $\ell^{\rm\,th}$
moment of the phase-space perturbation of relativistic species $a$.

For the FOPT we study in this paper, there are two linearly independent modes:
the adiabatic mode that contains the curvature perturbation, and the
isocurvature mode that contains the relative entropy perturbation. Therefore
each perturbation variable can be written as a sum over these two modes using
the decomposition

$\displaystyle X(\mathbf{k},\tau)=c^{\rm ad}(\mathbf{k})X^{\rm
ad}(k,\tau)+c^{\rm iso}(\mathbf{k})X^{\rm iso}(k,\tau).$ (28)

Here $c^{A}(\mathbf{k})$ is a time-independent coefficient (where $A\in\\{\rm
ad,iso\\}$) which is common for all perturbation variables, and
$X^{A}(k,\tau)$ is perturbation-specific but independent of the direction
$\mathbf{\hat{k}}$. The initial conditions for the modes of $X^{A}(k,\tau)$
are obtained by solving the Boltzmann and Einstein equations in the super-
horizon limit $(k\tau\ll 1)$.

The initial conditions for the adiabatic mode with DR are a straightforward
extension of the adiabatic initial conditions for radiation in the standard
cosmology Ma and Bertschinger (1995): for example, $\delta^{\rm ad}_{\rm
dr}=\delta^{\rm ad}_{\gamma}=\delta^{\rm ad}_{\nu}$. The initial conditions
for the DR isocurvature mode can be derived in analogy to neutrino density
isocurvature Bucher  _et al._ (2000). As in the standard procedure, we choose
a basis such that $\sum_{a}\delta\rho^{\rm iso}_{a}=0$, where the sum runs
over the isocurvature energy density perturbations of all radiation species
($a=\gamma$, $\nu$, ${\rm dr}$). This ensures that the total energy density of
radiation in the isocurvature mode is homogeneous at the initial time $\tau\to
0$ in the radiation-dominated Universe and thus does not source curvature.
Since there is no isocurvature in neutrinos in our models, their density
perturbations are equal to that of photons $(\delta_{\gamma}^{\rm
iso}=\delta_{\nu}^{\rm iso}$). We also choose the normalization such that
$\delta_{\rm dr}^{\rm iso}=1$ when $k\tau\to 0$. These requirements fix the
initial conditions for the DR isocurvature mode:

$\displaystyle\delta_{\rm dr}^{\rm iso}=1~{},~{}\delta_{\gamma}^{\rm
iso}=\delta_{\nu}^{\rm iso}=\frac{-R_{\rm dr}}{1-R_{\rm dr}}~{}~{}~{}(k\tau\to
0),$ (29)

where $R_{\rm dr}\equiv\bar{\rho}_{\rm
dr}/(\bar{\rho}_{\gamma}+\bar{\rho}_{\nu}+\bar{\rho}_{\rm dr})$. The complete
initial conditions for the DR isocurvature mode in the synchronous gauge are
given in Appendix D. Notably, in the limit $R_{\rm dr}\ll 1$, all moments of
the photon and neutrino phase-space perturbations in the isocurvature mode
($F^{\rm iso}_{\gamma,\ell}$ and $F_{\nu,\ell}^{\rm iso}$) are proportional to
$R_{\rm dr}$, which can be seen based on their initial conditions and time
evolution controlled by the coupled Boltzmann hierarchy and Einstein
equations.

The coefficients $c^{A}(\mathbf{k})=\\{c^{\rm ad}(\mathbf{k}),c^{\rm
iso}(\mathbf{k})\\}$ encode the initial random configuration of each mode. For
the FOPT models of interest, the modes are approximately uncorrelated so the
statistics of Gaussian observables are determined by

$\langle
c^{A}(\mathbf{k})c^{A}(\mathbf{k^{\prime}})\rangle=2\pi^{2}\frac{(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k^{\prime}})}{k^{3}}P_{A}(k),$
(30)

where $P_{A}(k)=\\{P_{\rm ad}(k),P_{\rm iso}(k)\\}$ are the power spectra of
the adiabatic and isocurvature modes.

One relevant observable in the CMB is the photon temperature anisotropy
$\Delta\equiv\Delta T/\bar{T}$, which can be written in a Fourier-Legendre
series as

$\Delta(\mathbf{x},\mathbf{\hat{n}},\tau)=\sum_{\ell=0}^{\infty}(-i)^{\ell}(2\ell+1)\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\mathbf{k}\cdot\mathbf{x}}\Delta_{\ell}(\mathbf{k},\tau)P_{\ell}(\mathbf{\hat{k}}\cdot\mathbf{\hat{n}}),$
(31)

where $\mathbf{\hat{n}}$ is the direction of the momentum and $\Delta_{\ell}$
is proportional to the multipole moment of the photon phase-space density
perturbation $\Delta_{\ell}=F_{\gamma,\ell}/4$. Its angular two-point function
is

$\langle\Delta(\mathbf{x},\mathbf{\hat{n}},\tau_{0})\Delta(\mathbf{x},\mathbf{\hat{n}^{\prime}},\tau_{0})\rangle=\frac{1}{4\pi}\sum_{\ell}(2\ell+1)C^{TT}_{\ell}P_{\ell}(\mathbf{\hat{n}}\cdot\mathbf{\hat{n}^{\prime}}).$
(32)

Here $C^{TT}_{\ell}$ is the temperature anisotropy angular power spectrum:

$\begin{split}C^{TT}_{\ell}=4\pi\int d(\ln{k})\left(P_{\rm
ad}(k)|\Delta_{\ell}^{\rm ad}(k,\tau_{0})|^{2}+\right.\\\ \left.P_{\rm
iso}(k)|\Delta_{\ell}^{\rm iso}(k,\tau_{0})|^{2}\right).\end{split}$ (33)

In this expression, $\Delta_{\ell}^{A}(k,\tau)$ is the photon transfer
function, which is obtained by solving the Boltzmann and Einstein equations
for $F_{\gamma,\ell}^{A}(k,\tau)$ with the initial conditions for mode
$A=\\{\rm ad,iso\\}$. The dependence on the particle physics models is
therefore encapsulated in $P_{\rm ad}(k)$ and $P_{\rm iso}(k)$.

We focus on a class of models for which the contribution of the FOPT to
$P_{\rm ad}(k)$ is small compared to the adiabatic modes arising from
inflation (see Appendix C). For this class of models, the dominant effects on
CMB observables is from $P_{\rm iso}(k)$ so we can set $P_{\rm ad}(k)$ to be
the standard $\Lambda$CDM adiabatic power spectrum

$\displaystyle P_{\rm ad}(k)=A_{s}\left(\frac{k}{k_{\rm
pivot}}\right)^{n_{s}-1},$ (34)

where $A_{s}$ and $n_{s}$ are the scalar amplitude and spectral index. $k_{\rm
pivot}$ is the pivot scale, conventionally defined as $k_{\rm pivot}\equiv
0.05\,{\rm Mpc}^{-1}$ Aghanim  _et al._ (2020). The photon transfer function
in the isocurvature mode scales as $\Delta^{\rm iso}_{\ell}=\frac{1}{4}F^{\rm
iso}_{\gamma,\ell}\propto R_{\rm dr}$ for $R_{\rm dr}\ll 1$ (see the
discussion below Eq. (29)). Therefore, the effect of isocurvature on the
angular power spectrum is determined by $R_{\rm dr}^{2}P_{\rm iso}(k)$. This
behavior will be manifest in our results.

###  IV.1 Isocurvature

To obtain $P_{\rm iso}(k)$ in Eq. (33), we consider the gauge-invariant
relative entropy perturbation between the DR and the photons:

$\displaystyle\mathcal{S}_{{\rm{dr}},\gamma}\equiv-3{\cal
H}\left(\frac{\delta\rho_{\rm dr}}{\rho_{\rm
dr}^{\prime}}-\frac{\delta\rho_{\gamma}}{\rho_{\gamma}^{\prime}}\right).$ (35)

Here ′ denotes $\partial/\partial\tau$, and ${\cal H}=aH$ is the conformal
Hubble parameter. With the definition in Eq. (28) and the isocurvature mode
relations in Eq. (29), as $k\tau\to 0$ this reduces to

$\displaystyle\mathcal{S}_{{\rm{dr}},\gamma}=\frac{3}{4}\frac{c^{\rm
iso}(\mathbf{k})}{1-R_{\rm dr}}.$ (36)

Therefore, the isocurvature power spectrum $P_{\rm iso}$ can be directly
related to
$\langle\mathcal{S}_{{\rm{dr}},\gamma}\,\mathcal{S}_{{\rm{dr}},\gamma}\rangle$.

To calculate $\mathcal{S}_{{\rm{dr}},\gamma}$, we first determine the
relations between $\delta_{\rm dr}$ and $\delta_{\chi}$. As the DR we are
interested in was generated by rapid bubble percolation around a temperature
$T_{*}$ (given by $H(T_{*})^{4}=\Gamma_{\rm PT}$), we can assume an
instantaneous transfer of energy from the $\chi$ field into DR when the
temperature of the Universe reaches $T_{*}$. We therefore work in the gauge
where the temperature of the Universe is uniform just before the phase
transition completes. This can be achieved by the following gauge
transformation:

$\delta\tau(\mathbf{x},\tau)=-\frac{\delta\rho_{\rm SM}}{\rho_{\rm
SM}^{\prime}},$ (37)

where SM means a total over particle species in the Standard Model. In this
gauge, the standard model density perturbation (defined as
$\delta\hat{\rho}_{\rm SM}\equiv\delta\rho_{\rm SM}+\rho_{\rm
SM}^{\prime}\,\delta\tau$) vanishes and the density perturbation $\chi$
remains unchanged ($\delta\hat{\rho}_{\chi}=\delta\rho_{\chi}$).

Since the temperature is uniform in the new gauge, the phase transition
completes at the same time $\tau_{*}$ everywhere and the large-scale
perturbations (super-horizon at $T_{*}$) in $\hat{\rho}_{\chi}$ will be
directly converted to those in $\hat{\rho}_{\rm dr}$. The sub-horizon
perturbations of DR around $T_{*}$ are more complicated, and depend on the
dynamics of the FOPT. However, as we are only interested in modes that are
accessible to CMB, we can safely neglect contributions from sub-horizon
dynamics at $T_{*}$, assuming $T_{*}\gg T_{\rm CMB}$. In this paper, we can
therefore approximate

$\displaystyle\hat{\delta}_{\rm
dr}(\mathbf{k},\tau_{*})\approx\delta_{\chi}(\mathbf{k},t_{*})\approx\delta_{\chi}(\mathbf{k},t_{e}),$
(38)

where $t=t_{*}$ when the phase transition completes (at proper time
$\tau_{*}$). The last approximation in Eq. (38) holds because bubbles
nucleated at times $t\in[t_{e},t_{*}]$ have radii that are subhorizon at
$T_{*}$.

The Standard Model density perturbations prior to the end of the FOPT are
generated from inflaton fluctuations and thus are all adiabatic, satisfying

$\frac{\delta\rho_{a}}{\rho_{a}^{\prime}}=\frac{\delta\rho_{\rm SM}}{\rho_{\rm
SM}^{\prime}}$ (39)

for any standard model particle species $a$. Therefore density perturbations
of each individual Standard Model species vanish in the new gauge as well
($\delta\hat{\rho}_{a}=0$) and the relative entropy perturbation is simply

$\displaystyle\mathcal{S}_{{\mathrm{dr}},\gamma}=\frac{3}{4}\hat{\delta}_{\rm
dr}.$ (40)

As $\mathcal{S}_{{\mathrm{dr}},\gamma}$ is gauge invariant, we can directly
relate Eq. (36) and Eq. (40). Together with Eq. (38), we find

$\displaystyle c^{\rm iso}(\mathbf{k})\approx(1-R_{\rm
dr})\delta_{\chi}(\mathbf{k},\tau_{e}).$ (41)

Therefore, $P_{\rm iso}$ is proportional to $P_{\chi}$:

$\displaystyle P_{\rm iso}(k)\approx(1-R_{\rm dr})^{2}P_{\chi}(k).$ (42)

As $P_{\chi}$ exhibits the broken power law behavior in Eq. (22), for the
remainder of our work, we take $P_{\rm iso}$ as

$\displaystyle P_{\rm iso}(k)=f_{\rm
iso}^{2}A_{s}\left\\{\begin{array}[]{cc}(k/k_{i})^{3}&\quad k\leq k_{i}\\\
1&\quad k>k_{i}\end{array}\right.,$ (45)

where $A_{s}$ is the amplitude of the scalar curvature power spectrum and
$f_{\rm iso}$ is a coefficient that sets the amplitude of the isocurvature
power spectrum. In our model, $f_{\rm iso}^{2}A_{s}\approx 4.2\gamma_{\rm
PT}$. We note that we use $k_{i}$ as the scale of the transition of the broken
power law, which is related to $r_{i}^{-1}$ by an ${\cal O}(1)$ factor.

###  IV.2 Non-Gaussianity

The perturbations generated by the FOPT have non-Gaussian statistics that may
be observable in the CMB. The CMB power spectra discussed before are based on
two-point functions, which only contain Gaussian information. To quantify non-
Gaussianity in the CMB, we need to examine higher point functions of the
photon temperature anisotropy $\Delta$. One example is the three-point
function, which comes predominantly from the isocurvature mode for the models
considered, given by

$\begin{split}\langle\Delta(\mathbf{\hat{n}_{1}})\Delta(\mathbf{\hat{n}_{2}})\Delta(\mathbf{\hat{n}_{3}})\rangle=&\left(\frac{20}{3}\right)^{3}\sum_{\ell_{1},\ell_{2},\ell_{3}}(2\ell_{1}+1)(2\ell_{2}+1)(2\ell_{3}+1)(-i)^{\ell_{1}+\ell_{2}+\ell_{3}}\\\
\times&\int\frac{d^{3}k_{1}}{(2\pi)^{3}}\int\frac{d^{3}k_{2}}{(2\pi)^{3}}\int\frac{d^{3}k_{3}}{(2\pi)^{3}}(2\pi)^{3}\delta^{3}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})P_{\ell_{1}}(\mathbf{\hat{k}_{1}}\cdot\mathbf{\hat{n}_{1}})P_{\ell_{2}}(\mathbf{\hat{k}_{2}}\cdot\mathbf{\hat{n}_{2}})P_{\ell_{3}}(\mathbf{\hat{k}_{3}}\cdot\mathbf{\hat{n}_{3}})\\\
\times&B_{\rm iso}(k_{1},k_{2},k_{3})\Delta_{\ell_{1}}^{\rm
iso}(k_{1},\tau_{0})\Delta_{\ell_{2}}^{\rm
iso}(k_{2},\tau_{0})\Delta_{\ell_{3}}^{\rm iso}(k_{3},\tau_{0}).\end{split}$
(46)

We have introduced the isocurvature bispectrum $B_{\rm iso}$, which we define
using the convention Akrami  _et al._ (2020b):

$\begin{split}&\langle c^{\rm iso}(\mathbf{k_{1}})c^{\rm
iso}(\mathbf{k_{2}})c^{\rm iso}(\mathbf{k_{3}})\rangle\\\
&\equiv(2\pi)^{3}\delta^{3}(\mathbf{k_{1}}+\mathbf{k_{2}}+\mathbf{k_{3}})\left(\frac{20}{3}\right)^{3}B_{\rm
iso}(k_{1},k_{2},k_{3}).\end{split}$ (47)

To calculate the bispectrum for our model, we use Eq. (41) to relate $c^{\rm
iso}(\mathbf{k})$ to $\delta_{\chi}(\mathbf{k},t_{e})$ and perform steps
similar to those for the calculation of $P_{\chi}$, leading to

$\begin{split}B_{\rm
iso}(k_{1},k_{2},k_{3})=&-\left(\frac{3\pi}{5}\right)^{3}(1-R_{\rm
dr})^{3}e^{4\pi\gamma_{\rm PT}H_{\rm inf}t_{e}}\frac{\gamma_{\rm
PT}}{k_{1}^{3}k_{2}^{3}}\\\
\times&\int\displaylimits_{k_{3}r_{e}}^{k_{3}r_{i}}duu^{-4}\mathcal{A}(\alpha_{1}u)\mathcal{A}(\alpha_{2}u)\mathcal{A}(u),\end{split}$
(48)

where $\alpha_{1}\equiv k_{1}/k_{3}$, $\alpha_{2}\equiv k_{2}/k_{3}$, and
$\mathcal{A}(u)$ is defined below Eq. (16). In Section VII.3 we will estimate
constraints by comparing this result to the equilateral and local bispectra
which have been studied in the literature Akrami  _et al._ (2020b). To
simplify the comparison, we will write down expressions for $B_{\rm iso}$
evaluated at wavenumbers in the squeezed $k_{1}\ll k_{2}\approx k_{3}$) and
equilateral ($k_{1}\approx k_{2}\approx k_{3}$) configurations. For $R_{\rm
dr}\ll 1$ and $\gamma_{\rm PT}H_{\rm inf}t_{e}\ll 1$ these expressions are

$\begin{split}B_{\rm
iso}(k_{1},k,k)=&-\left(\frac{3\pi}{5}\right)^{3}\frac{\gamma_{\rm
PT}}{k^{6}}\\\
\times&\left\\{\begin{array}[]{cc}\frac{1}{3}\int\displaylimits_{0}^{kr_{i}}duu^{-1}\mathcal{A}(u)^{2}&k_{1}\ll
k\\\
\int\displaylimits_{0}^{kr_{i}}duu^{-4}\mathcal{A}(u)^{3}&k_{1}=k\end{array}\right.\,.\end{split}$
(49)

Restricting to the case where $kr_{i}\gg 1$, the asymptotic behavior of the
integrals in these expressions is

$\begin{split}\frac{1}{3}\int\displaylimits_{0}^{kr_{i}}duu^{-1}\mathcal{A}(u)^{2}\propto(kr_{i})^{2}\\\
\int\displaylimits_{0}^{kr_{i}}duu^{-4}\mathcal{A}(u)^{3}\propto{\rm
constant},\end{split}$ (50)

so

$B_{\rm iso}(k_{1},k,k)\propto\left\\{\begin{array}[]{cr}k^{-4}&\quad k_{1}\ll
k\\\ k^{-6}&\quad k_{1}=k\end{array}\right..$ (51)

##  V Effects on CMB observations

In the previous section, we have calculated the DR isocurvature power spectrum
from non-thermal FOPTs. In this section, we implement the model in CLASS, and
simulate the CMB angular power spectra. As mentioned before, we can
approximate the adiabatic power spectrum as unchanged compared to the standard
$\Lambda$CDM model. The presence of DR (parameterized by $\Delta N_{\rm eff}$)
shifts CMB observables as a result of the change in the background energy
density – such shifts create the standard limit on $\Delta N_{\rm eff}\lesssim
0.3$. In addition, DR isocurvature creates additional shifts in CMB
observables which are the main focus of this work.

As explained in the introduction to Section IV, the effect of DR isocurvature
on the CMB depends on $R_{\rm dr}^{2}P_{\rm iso}(k)$. To simplify comparisons
with limits on DR with adiabatic perturbations, we write $R_{\rm dr}$ in terms
of the equivalent number of effective neutrino species, $\Delta N_{\rm
eff}\equiv 3.044\times(\bar{\rho}_{\rm dr}/\bar{\rho}_{\nu})$. Around photon
decoupling, $R_{\rm dr}\approx 0.13\times\Delta N_{\rm eff}$. The isocurvature
power spectrum (defined in Eq. (45)) depends on two additional parameters:
$f_{\rm iso}$ determines the amplitude, and $k_{i}$ is the wavenumber at which
the broken power law transitions. This transition wavenumber also corresponds
to the comoving scale at the start of the FOPT. Given the form of $P_{\rm
iso}$, the effect of DR isocurvature on the CMB is proportional to $\Delta
N_{\rm eff}^{2}f_{\rm iso}^{2}$ for fixed $k_{i}$. We express $k_{i}$ in terms
of multiples of $k_{\rm pivot}\equiv 0.05\,{\rm Mpc}^{-1}$, the pivot scale of
the adiabatic power spectrum given in Eq. (34).

To show the effects of DR isocurvature on the CMB observables, we simulate the
CMB power spectra with difference choices of $k_{i}$, and compare them to the
case with purely adiabatic initial conditions. In Figure 5, we show the
difference of the temperature angular power spectrum $C_{\ell}^{TT}$ between
$\Lambda$CDM and three models with DR, all with $\Delta N_{\rm eff}=0.1$. In
one case, we consider only the adiabatic effects of DR (setting $f_{\rm
iso}=0$). In the other two we include isocurvature with $k_{i}=k_{\rm pivot}$
and $0.1\,k_{\rm pivot}$ (with $f_{\rm iso}=10$ in both cases). The acoustic
peaks in models with isocurvature shift towards higher $\ell$ as compared to
the adiabatic case with the same $\Delta N_{\rm eff}$. This result is
consistent with the literature (see for example Savelainen  _et al._ (2013)).
In Figure 6, we show the fractional difference of angular power spectra
between the DR isocurvature and adiabatic cases (again holding $\Delta N_{\rm
eff}=0.1$) for the $\rm TT$ (left) and $\rm EE$ (right) modes. For the DR
isocurvature simulations again we fix $f_{\rm iso}=10$ and vary $k_{i}$. The
oscillations seen in both panels are a consequence of the shift in the
acoustic peaks due to isocurvature, as discussed.

Figure 5: The difference of temperature anisotropy angular power spectra
between cases with dark radiation ($\Delta N_{\rm eff}=0.1$) and the
$\Lambda$CDM model. The blue and orange lines denote DR isocurvature with
different $k_{i}$ and fixed $f_{\rm iso}=10$. The green line shows the result
with adiabatic initial conditions. The pivot scale is chosen to be $k_{\rm
pivot}=0.05\,\textrm{Mpc}^{-1}$.

For all DR isocurvature models shown in Figures 5 and 6, differences between
the CMB power spectra disappear at large $\ell$. This occurs because CMB
spherical harmonics with $\ell\gtrsim\ell_{i}\sim k_{i}d_{\rm SLS}$ mostly
receive contributions from wavenumbers $k\gtrsim k_{i}$, where the
isocurvature power spectrum is flat and independent of $k_{i}$. Here $d_{\rm
SLS}\sim 10^{4}\,\rm{Mpc}$ is the comoving distance to the surface of last
scattering.777Thus, for $k_{i}=k_{\rm pivot}$, $\ell_{i}\sim 500$. For
$k_{i}\lesssim 2\times 10^{-3}\,k_{\rm pivot}$, the corresponding transition
spherical harmonic is $\ell_{i}\lesssim 1$ and the entire CMB angular power
spectrum can only probe the plateau of the DR isocurvature power spectrum. For
such $k_{i}$, the effect of DR isocurvature saturates and the CMB becomes
insensitive to further reductions in $k_{i}$.

Figure 6: The fractional difference of CMB power spectra ($C^{TT}_{\ell}$ in
the left panel, $C^{EE}_{\ell}$ in the right panel) between cases with dark
radiation isocurvature and adiabatic initial conditions for different $k_{i}$.
We set $f_{\rm iso}=10$ for DR isocurvature, while $\Delta N_{\rm eff}=0.1$
for both adiabatic and isocurvature cases. The pivot scale is chosen to be
$k_{\rm pivot}=0.05\,\textrm{Mpc}^{-1}$.

##  VI Methodology and Data Sets

To evaluate the cosmological constraints on the FOPT model parameters, we
perform a full likelihood analysis using cosmological datasets. To compute the
observables within these datasets, we employ a modified version of the CLASS
code Ghosh  _et al._ (2022); Blas  _et al._ (2011) along with MontePython
Audren  _et al._ (2013); Brinckmann and Lesgourgues (2019) to perform the MCMC
analysis.

We use the datasets from the following experiments to set limits on the FOPT
model:

  * •

Planck 2018 high-$\ell$ TTTEEE and low-$\ell$ TT, EE anisotropy angular power
spectra and the gravitational lensing data Aghanim  _et al._ (2020).

  * •

The BAO measurements from the Six-degree Field Galaxy Survey (6dFGS) Beutler
_et al._ (2011), the Sloan Digital Sky Survey (SDSS) DR7 MGS Ross  _et al._
(2015), and the LOWZ galaxy samples of BOSS DR12 Alam  _et al._ (2017).

Our baseline cosmology consists of the usual six cosmological parameters of
$\Lambda$CDM. To this, we add the three isocurvature parameters discussed
previously: $\Delta N_{\rm eff}$, $k_{i}$, $f_{\rm iso}$. We assume three
active neutrinos species, one with mass $0.06\,{\rm eV}$ and two massless. As
explained in Section V, the effects of DR isocurvature are sensitive to the
product $\Delta N_{\rm eff}f_{\rm iso}$ while the effects of the background
cosmology are sensitive to $\Delta N_{\rm eff}$. It is also convenient to
express the transition wavenumber $k_{i}$ in terms of the CMB pivot wavenumber
$k_{\rm pivot}$. Therefore we use $\Delta N_{\rm eff}$,
$k_{i}/k_{\mathrm{pivot}}$ and $\Delta N_{\rm eff}f_{\rm iso}$ as the primary
parameters for our MCMC runs with a log-flat prior for each parameter. The
boundaries of these priors are:
$\log_{10}(\Delta{}N_{\mathrm{eff}})\in[-5,-0.3]$,
$\log_{10}(k_{i}/k_{\mathrm{pivot}})\in[-4,3]$ and
$\log_{10}(\Delta{}N_{\mathrm{eff}}f_{\mathrm{iso}})\in[-2.3,2.3]$. The
$\Lambda$CDM parameters have standard wide flat priors. The FOPT isocurvature
power spectrum of our simplified model (calculated in Section IV.1) is a good
approximation provided $\gamma_{\rm PT}\lesssim 10^{-3}$. Therefore, the
particular FOPT model studied in this work maps onto the priors used in the
likelihood analysis for $f_{\rm iso}\lesssim 10^{3}$. We note that our
cosmological model and priors are sufficiently general that other isocurvature
generation mechanisms may lead to cosmologies that map onto subsets of our
parameter space. With this in mind, we show the resulting posteriors for
$f_{\rm iso}$ outside of the strict domain of validity for this class of
FOPTs.

We use the Metropolis-Hastings algorithm with the Gelman-Rubin convergence
criterion $R-1<0.02$ Gelman and Rubin (1992) for the MCMC chains. To analyze
MCMC chains and calculate Bayesian posteriors, we use the GetDist package
Lewis (2019). The best-fit values are obtained using simulated annealing for
the minimization of $\chi^{2}$ values Schöneberg  _et al._ (2022). The
procedure of simulated annealing systematically lowers the temperature during
an MCMC chain run, effectively navigating and optimizing against noisy
likelihood functions with numerous local maxima.

##  VII Results

Having discussed our theoretical and statistical methodology as well as the
datasets that we use in our analysis, we are ready to calculate posteriors for
the parameters of DR isocurvature model. The mean and best-fit values of each
parameter in our fit are shown in Table 1. The posterior distributions ($65\%$
and $95\%$ CL) for these new physics parameters are shown in Figure 7. The
limits on $\Delta N_{\rm eff}f_{\rm iso}$ become approximately independent of
$\Delta N_{\rm eff}$ for $\Delta N_{\rm eff}\ll 0.1$ and are also independent
of $k_{i}$ for $k_{i}\ll k_{\rm pivot}$.

Figure 7:  The two dimensional marginalized posterior distributions for the
FOPT new physics parameters $\Delta N_{\rm eff}$, $\Delta N_{\rm eff}f_{\rm
iso}$ and $k_{i}/k_{\rm pivot}$, when analyzing Planck+BAO (see text for
details). We omit the posteriors of the six $\Lambda$CDM parameters for
clarity.

As seen in Table 1, there is a slight preference (compared to $\Lambda$CDM) in
the data for non-zero isocurvature, with best-fit parameters
$k_{i}/k_{\mathrm{pivot}}=0.36$, $\Delta N_{\rm eff}=7.5\times 10^{-3}$, and
$\Delta N_{\rm eff}f_{\rm iso}=1.58$. We show in Figure 8 the fractional
difference in $C_{\ell}^{TT}$ between this best-fit model and the best-fit
$\Lambda$CDM model compared to the Planck 2018 residuals (also normalized to
the $\Lambda$CDM result). We also show the fractional difference compared to
$\Lambda$CDM of an adiabatic-only model with the same parameters as the DR
isocurvature best-fit other than setting $f_{\rm iso}=0$.

To better understand the features of our posteriors, it is useful to hold one
of the three new physics parameters fixed and consider posteriors over the
other two. We do this in the next two subsections, holding $f_{\rm iso}$ fixed
in Section VII.1 and $k_{i}$ fixed in Section VII.2. Finally, we estimate
limits from non-Gaussianities in Section VII.3, though a dedicated search is
needed to derive a robust bound.

Parameters | FOPT
---|---
$\Omega_{\mathrm{b}}h^{2}$ |  $0.0225$ ($0.0223\pm 0.0002$)
$\Omega_{\mathrm{c}}h^{2}$ |  $0.1199$ ($0.01954\pm 0.0011$)
$100*\theta{}_{s}$ |  $1.0421$ ($1.0420\pm 0.0003$)
$\ln(10^{10}A_{s})$ |  $3.0459$ ($3.0474\pm 0.0146$)
$n_{s}$ |  $0.9678$ ($0.9679\pm 0.0398$)
$\tau_{\mathrm{reio}}$ |  $0.0580$ ($0.0564\pm 0.0073$)
$\log_{10}(\Delta{}N_{\mathrm{eff}})$ |  $-2.1244$ ($<-1.0741$)
$\log_{10}(k_{i}/k_{\mathrm{pivot}})$ |  $-0.4496$ (unconstrained)
$\log_{10}(\Delta{}N_{\mathrm{eff}}f_{\mathrm{iso}})$ |  $0.1974$ ($<1.1453$)
$\Delta\chi^{2}_{\mathrm{tot}}$ | $-0.68$

Table 1: The best-fit parameters of the FOPT isocurvature model resulting from
fits to the Planck+BAO datasets. The mean and $1\sigma$ variation for each
parameter when fit to the data are shown in parenthesis or the 95% CL upper
bound. Figure 8: The change (relative to $\Lambda$CDM) in $C_{\ell}^{TT}$
power spectrum normalized to $C_{\ell,\Lambda{\rm CDM}}^{TT}$ for the FOPT
best-fit model (see Table 1), an adiabatic model and the Planck 2018 data
Aghanim  _et al._ (2020). For the adiabatic case we use the same parameters as
in the FOPT best-fit except we turn off isocurvature ($f_{\rm iso}=0$).

###  VII.1 $\Delta N_{\rm eff}$ constraints with fixed $f_{\rm iso}$

Figure 9: The two dimensional marginalized posterior distributions for the
FOPT new physics parameters $\Delta N_{\rm eff}$ and $k_{i}/k_{\rm pivot}$,
holding $f_{\rm iso}$ fixed.

In this subsection, we fix $f_{\rm iso}$ and perform the scan over $\Delta
N_{\rm eff}$, $k_{i}$ and the other six $\rm\Lambda$CDM parameters. This
corresponds to a scenario where the phase transition rate $\gamma_{\rm PT}\sim
f_{\rm iso}^{2}\,A_{s}$ is held fixed while the horizon size at the beginning
of the transition $r_{i}\sim k_{i}^{-1}$ varies. In Figure 9, we show
posteriors for this scan for three values of $f_{\rm iso}$: 10, 500 and 1000.

For the largest two values of $f_{\rm iso}$ we find that the 2$\sigma$
contours approach the adiabatic constraint $\Delta N_{\rm eff}<0.3$ for large
$k_{i}$. This is expected as the DR isocurvature power spectrum modifies the
CMB power spectrum primarily at $\ell$ above the range currently probed by
Planck when $k_{i}\gtrsim 10^{2}\,k_{\rm pivot}$, making the CMB power
spectrum observationally indistinguishable from the adiabatic-only result. As
$k_{i}$ decreases, the limit on $\Delta N_{\rm eff}$ strengthens rapidly. For
sufficiently small $k_{i}$, CMB data is mainly sensitive to the plateau of the
isocurvature power spectrum (see Eq. (45)) and the constraint approaches an
asymptotic value which depends on $\Delta N_{\rm eff}f_{\rm iso}$.

For $f_{\rm iso}=10$ we find a preference for the range
$k_{i}\sim(0.1-1)\times k_{\rm pivot}$. That is, models with this combination
of $k_{i}$ and $f_{\rm iso}$ fit the data better than models with large
$k_{i}$ (which are themselves indistinguishable from the adiabatic model).

###  VII.2 $\Delta N_{\rm eff}$ constraints with fixed $k_{i}$

Figure 10: The two dimensional marginalized posterior distributions for the
FOPT new physics parameters $\Delta N_{\rm eff}$ and $\Delta N_{\rm eff}f_{\rm
iso}$, holding $k_{i}$ fixed.

Next, we fix $k_{i}$ and calculate the posteriors in $\Delta N_{\rm eff}$ and
$\Delta N_{\rm eff}f_{\rm iso}$. This corresponds to the scenario where the
horizon size at the beginning of the PT is fixed, and the PT rate varies. The
results of our MCMC fit to the data are shown in Figure 10 for three choices
of $(k_{i}/k_{\mathrm{pivot}})$: $10^{3},1$, and $10^{-4}$.

Again, we see that when $k_{i}$ is large ($k_{i}=10^{3}\,k_{\rm pivot}$) the
constraints on $\Delta N_{\rm eff}$ approximate the adiabatic result for all
values of $f_{\rm iso}$ within our priors. For the other two values of
$k_{i}$, isocurvature becomes important and the model can be constrained even
for very small values of $\Delta N_{\rm eff}$. For $\Delta N_{\rm eff}\ll
0.1$, the adiabatic and background effects are negligible and the CMB data
mainly has sensitivity to the combination $\Delta N_{\rm eff}f_{\rm iso}$
which parameterizes the strength of isocurvature effects. In this regime, we
find the $95\%$ CL constraints $\Delta N_{\rm eff}f_{\rm iso}<1.25$ for
$k_{i}=10^{-4}\,k_{\rm pivot}$ and $\Delta N_{\rm eff}f_{\rm iso}<2.09$ for
$k_{i}=k_{\rm pivot}$. We have confirmed through explicit calculation using
linear-flat priors on $\Delta N_{\rm eff}$ that these results are robust
against our choice of priors.

The limits on $\Delta N_{\rm eff}f_{\rm iso}$ are stronger for the smallest
value of $k_{i}$ as the CMB only probes the plateau (the maximum) of the
isocurvature power spectrum. For the larger value, the CMB is sensitive to the
portion of the power spectrum which is proportional to $k^{3}$ and suppressed
relative to the plateau. As $\Delta N_{\rm eff}$ becomes $\mathcal{O}(0.1)$,
the adiabatic effects become important. At this point, the limits deviate from
$\Delta N_{\rm eff}f_{\rm iso}=\mathrm{const.}$ and we (approximately) recover
the adiabatic result: $\Delta N_{\rm eff}<0.3$.

The constraint on $\Delta N_{\rm eff}f_{\rm iso}$ for $\Delta N_{\rm eff}\ll
0.1$ can be rewritten as a constraint on $\Delta N_{\rm eff}$ in terms of
other quantities directly related to FOPTs. As shown below Eq. (45),
$A_{s}f_{\rm iso}^{2}\approx 4.2\gamma_{\rm PT}$ and $\gamma_{\rm
PT}=\Gamma_{\rm PT}/H_{\rm inf}^{4}$. Using the definitions of nucleation
temperature $T_{*}$ ($\Gamma_{\rm PT}=H^{4}(T_{*})$) and reheating temperature
$T_{\rm rh}$ ($H_{\rm inf}=H(T_{\rm rh})$) for instantaneous reheating, we can
rewrite $\gamma_{\rm PT}=(T_{*}/T_{\rm rh})^{8}$ in the radiation-dominated
Universe. Therefore, the constraint $\Delta N_{\rm eff}f_{\rm
iso}<\beta(k_{i})$ is equivalent to

$\displaystyle\Delta N_{\rm eff}<2.8\times 10^{-5}\left(\frac{T_{\rm
rh}}{T_{*}}\right)^{4}\left(\frac{\beta(k_{i})}{1.25}\right).$ (52)

for $T_{*}/T_{\rm rh}\gtrsim 0.1$. $\beta(k_{i})$ is mildly sensitive to
$k_{i}$ for $k_{i}\lesssim k_{\rm pivot}$ as demonstrated by our limits on
$\Delta N_{\rm eff}f_{\rm iso}$ for small $\Delta N_{\rm eff}$ discussed
earlier this subsection.

For the two smallest values of $k_{i}$ shown in Figure 10, there is a
preference for models with isocurvature ($\Delta N_{\rm eff}f_{\rm iso}\sim
0.1-1$) compared to the adiabatic-only model which is well approximated by the
smallest values of $\Delta N_{\rm eff}f_{\rm iso}$ allowed by our prior. This
effect is more prominent for $k_{i}=k_{\rm pivot}$ which is close to the best-
fit value of $k_{i}$.

###  VII.3 $\Delta N_{\rm eff}$ constraints from non-Gaussianity

In this section, we will estimate constraints on this class of FOPTs from non-
Gaussianity in the CMB. Currently, there is no dedicated search for non-
Gaussianity with DR isocurvature produced in a FOPT. The most relevant study
is on non-Gaussianity with neutrino density isocurvature Akrami  _et al._
(2020b).

Based on Eq. (46), the effect of non-Gaussianity in our study will be
comparable in size to that of neutrino density isocurvature when

$\begin{split}&B_{\rm iso}(k_{1},k_{2},k_{3})\Delta_{\ell_{1}}^{\rm
iso}(k_{1},\tau_{0})\Delta_{\ell_{2}}^{\rm
iso}(k_{2},\tau_{0})\Delta_{\ell_{3}}^{\rm iso}(k_{3},\tau_{0})\sim\\\
&B_{\nu,\rm iso}(k_{1},k_{2},k_{3})\Delta_{\ell_{1}}^{\nu,\rm
iso}(k_{1},\tau_{0})\Delta_{\ell_{2}}^{\nu,\rm
iso}(k_{2},\tau_{0})\Delta_{\ell_{3}}^{\nu,\rm
iso}(k_{3},\tau_{0}).\end{split}$ (53)

where $B_{\nu,\rm iso}$ and $\Delta_{\ell}^{\nu,\rm iso}$ are the bispectrum
and transfer function for the neutrino density isocurvature model. This
transfer function can be simply related to the transfer function for DR
isocurvature by accounting for the different energy fraction of DR and
neutrinos. As discussed in Section IV, the transfer function for DR
isocurvature is proportional to $R_{\rm dr}$. Therefore, the transfer function
for neutrino isocurvature can be well-approximated by

$\Delta^{\nu,\rm iso}_{\ell}(k,\tau)\sim\frac{R_{\nu}}{R_{\rm dr}}\Delta^{\rm
iso}_{\ell}(k,\tau),$ (54)

where
$R_{i}\equiv\bar{\rho}_{i}/(\bar{\rho}_{\gamma}+\bar{\rho}_{\nu}+\bar{\rho}_{\rm
dr})$. Subsituting this result into Eq. (53), we find that the effect of non-
Gaussianity will be similar in the two models when

$B_{\rm iso}(k_{1},k_{2},k_{3})\sim\left(\frac{R_{\nu}}{R_{\rm
dr}}\right)^{3}B_{\nu,\rm iso}(k_{1},k_{2},k_{3}).$ (55)

The available constraints on neutrino isocurvature non-Gaussianity Akrami  _et
al._ (2020b) assume a local bispectrum which peaks in the squeeze
configuration:

$B_{\rm local}(k_{1},k,k)\propto\frac{f_{\rm NL}^{\rm
local}}{k_{1}^{4-n_{s}}k^{4-n_{s}}}\quad k_{1}\ll k$ (56)

where $f_{\rm NL}^{\rm local}$ is a dimensionless parameter that quantifies
the amount of non-Gaussianity present. The FOPT bispectrum, on the other hand,
is independent of $k_{1}$ in the squeeze configuration (as shown in Eq. (51))
so the bound on $f^{\rm local}_{\rm NL}$ does not apply to our case.

The bispectrum from the FOPT is similar to the equilateral template for which
the dominant effects come from the equilateral configuration Akrami  _et al._
(2020b)

$B_{\rm equil}(k,k,k)=6\left(\frac{2\pi^{2}}{25}\right)^{2}A_{s}^{2}f_{\rm
NL}^{\rm equil}\frac{1}{k^{6}}\left(\frac{k_{\rm
pivot}}{k}\right)^{2(1-n_{s})}.$ (57)

However, measurements of $f_{\rm NL}^{\rm equil}$ for isocurvature non-
Gaussianity do not exist in the literature. Based on measurements of
equilateral and local non-Gaussianity for curvature Akrami  _et al._ (2020b),
we assume that a search for equilateral non-Gaussianity for isocurvature would
find a value of $f_{\rm NL}^{\rm equil}$ within an order of magnitude of
$f_{\rm NL}^{\rm local}$. Under this assumption and given the 2$\sigma$ upper
bound on $|f_{\rm NL}^{\rm local}|$ for neutrino density isocurvature Akrami
_et al._ (2020b), we can estimate an approximate bound $|f_{\rm NL}^{\rm
equil}|\lesssim 10^{3}$ for the same cosmological model. This result can be
translated into a bound on FOPT non-Gaussianity by evaluating Eq. (55) in the
equilateral configuration ($k_{1}\approx k_{2}\approx k_{3}$) near the pivot
scale $k_{\rm pivot}$ and taking $B_{\rm\nu,iso}=B_{\rm equil}$. After
numerical integration of the second line of Eq. (49), the FOPT bispectrum can
be written as

$B_{\rm iso}(k,k,k)=-0.97A_{s}f_{\rm iso}^{2}\frac{1}{k^{6}}.$ (58)

Combining Eqs. (55), (57) and (58) leads to the approximate constraint

$\Delta N_{\rm eff}^{3}f_{\rm iso}^{2}\lesssim 2\times 10^{-4}.$ (59)

In Section VII.2 we found $\Delta N_{\rm eff}f_{\rm iso}\lesssim 1$ in the
regime $k_{i}\lesssim k_{\rm pivot}$ and $\Delta N_{\rm eff}\ll 0.3$.
Therefore the dominant constraint on FOPT may come from non-Gaussianity for
$\Delta N_{\rm eff}\gtrsim 2\times 10^{-4}$, though we again emphasize that a
dedicated study is required to set a robust bound.

##  VIII Conclusions

We have showed that for a broad class of non-thermal FOPTs, the scalar field
undergoes a phase transition that nucleates bubbles during inflation. If the
PT rate is small during inflation, the bubbles do not collide and their size
is determined by the horizon size at the time of nucleation. The FOPT will
eventually complete after reheating, at which point we assume most of the
remaining energy density in the scalar field converts into dark radiation. The
DR inherits the inhomogeneity of the scalar field, leading to a DR
isocurvature mode.

If the PT starts sufficiently early during inflation, isocurvature in the DR
can create measurable perturbations in the CMB. In this work, we performed a
detailed calculation of the DR isocurvature power spectrum and implemented it
in CLASS code. Using Planck 2018 and BAO data, we calculated the constraint on
the energy density of DR (in terms of the effective relativistic degrees of
freedom $\Delta N_{\rm eff}$). We demonstrated that the constraint on $\Delta
N_{\rm eff}$ can be much stronger than the limit derived assuming adiabatic
initial conditions. In particular, when the temperature $T_{*}$ at which the
phase transition completes, is close to the reheating temperature $T_{\rm
rh}$, the isocurvature perturbations set a limit of $\Delta N_{\rm
eff}\lesssim 10^{-5}(T_{*}/T_{\rm rh})^{-4}$. This limit weakens as $T_{*}$
decreases and approaches the adiabatic constraint of $\Delta N_{\rm eff}<0.3$
when $T_{*}$ becomes much less than $T_{\rm rh}$.

Since the distribution of the perturbation caused by bubbles is intrinsically
non-Gaussian, it could leave sizable non-Gaussianity signals in the CMB. We
have also estimated the strength of the non-Gaussianity constraint on FOPTs
from the CMB by comparing to neutrino density isocurvature. However, we showed
that the bispectrum calculated for this class of FOPTs is different than those
studied previously. Therefore, a dedicated search is needed to derive a robust
bound on FOPTs from non-Gaussianity in the CMB.

Moreover, our study can also put a lower bound on the scale of inflation if we
can determine $\Delta N_{\rm eff}$ and $T_{*}$. These quantities can be
derived if we observe the direct gravitational wave spectrum from the FOPT. We
will show in a follow-up paper that the recent nano-Hz gravitational wave data
could already set a lower limit on the scale of the reheating temperature. In
the future, it will be interesting to explore the correlation of gravitational
wave signals and cosmological data to further probe FOPTs.

##  IX Acknowledgements

We thank Subhajit Ghosh, Soubhik Kumar and Chen Sun for useful discussions.
This work was supported by DOE grant DOE-SC0010008. Key parts of this work
were motivated by discussions at the Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-2210452.

##  Appendix A Bubble Wall Dynamics

As the $\chi$ field tunnels to the true vacuum in a region of spacetime, a
wall forms as $\chi$ smoothly varies from $\chi_{-}$ inside the region to
$\chi_{+}$ outside. We follow the approach of Ref. Mégevand and Membiela
(2023) to calculate the field profile and dynamics inside the wall, working in
the thin wall approximation which assumes that the wall thickness is smaller
than any other length scale in the problem.

For this calculation, it is convenient to work in Gaussian-normal coordinates
adapted to a hypersurface $\Sigma$ along the wall. On the hypersurface, the
field takes on a constant value
$\chi(x^{\mu})=\chi_{\Sigma}\in(\chi_{-},\chi_{+})$. Points on $\Sigma$ are
parameterized as $x^{\mu}=X^{\mu}(\xi^{a})$ for three coordinates
$\xi^{a}=(\xi^{0},\xi^{1},\xi^{2})$. The fourth coordinate $n$ measures the
proper distance along a geodesic that originates from $X^{\mu}(\xi^{a})$ with
a tangent vector that is orthogonal to $\Sigma$. In this coordinate system, an
arbitrary point near $\Sigma$ has Gaussian-normal coordinates which can be
related to the original coordinates by

$x^{\mu}=X^{\mu}(\xi^{a})+nN^{\mu}(\xi^{a})+\mathcal{O}(n^{2})$ (60)

where $N^{\mu}(\xi^{a})$ is a unit normal vector satisfying

$\begin{split}N_{\mu}N^{\mu}=&1\\\ N_{\mu}\partial_{a}X^{\mu}=&0.\\\
\end{split}$ (61)

On the hypersurface $\Sigma$, $n=0$.

The scalar field $\chi$ obeys the standard equation of motion:

$g^{\mu\nu}\left[\partial_{\mu}\partial_{\nu}\chi-\Gamma^{\rho}_{\mu\nu}\partial_{\rho}\chi\right]=V^{\prime}(\chi).$
(62)

In the thin-wall approximation, $\partial_{n}\chi\gg\partial_{a}\chi$,
allowing us to drop derivatives other than those with respect to $n$. Thus, in
Gaussian-normal coordinates, the equations of motion are approximately888This
equation is exact at $n=0$ as, in the Gaussian normal coordinates,
$\left.\partial_{a}\chi\right|_{n=0}=0$.

$\partial^{2}_{n}\chi-K(n)\partial_{n}\chi=V^{\prime}(\chi(n)),$ (63)

where $K(n)$ is the mean extrinsic curvature of the hypersurface at $n$. In
the thin wall approximation, $K(n)$ does not vary considerably over the width
of the wall, so we can replace $K(n)\to K(0)\equiv K$. Multiplying Eq. (63) by
$\partial_{n}\chi$ and integrating, we obtain

$\frac{1}{2}(\partial_{n}\chi)^{2}-K\int\displaylimits_{-\infty}^{n}dn(\partial_{n}\chi)^{2}=V(\chi(n)),$
(64)

where $V(\chi)$ is the potential defined in Eq. (2) with an appropriate
constant added such that $V$ is positive definite. Here we have defined $n$
such that $\chi(n)\to\chi_{-}$ as $n\to-\infty$. Taking $n\to\infty$ in Eq.
(64) gives

$-K\sigma=\Delta V$ (65)

where $\sigma$ is the surface tension, defined as

$\sigma\equiv\int dn(\partial_{n}\chi)^{2}.$ (66)

This is equivalent to the radial bounce action $S_{1}$, introduced in Section
II. Since $K$ can be written in terms of the shape and dynamics of the wall,
Eq. (65) can be used as a dynamic equation for the wall after we calculate
$\sigma$.

To determine the surface tension $\sigma$, we solve Eq. (64). To do so, it is
useful to consider the order of magnitude of the terms in the equation for
values of $n$ inside the wall:

$\begin{split}K\int\displaylimits_{-\infty}^{n}dn(\partial_{n}\chi)^{2}\sim&\mathcal{O}\left(\Delta
V\right)\\\ V(\chi)\sim&\mathcal{O}\left(V_{\rm max}\right),\end{split}$ (67)

where $V_{\rm max}$ is the maximum potential value between the two minima. The
thin-wall approximation requires that $\Delta V\ll V_{\rm max}$. Therefore, an
approximate solution for the profile of the bubble wall can be found by
setting $K=0$ in Eq. (64). However, from Eq. (65) this solution is only
consistent if we set $\Delta V=0$ as well. We find a self-consistent
approximate solution by splitting the potential into two parts:

$\begin{split}V=&V_{0}+V_{1}\\\
V_{0}\equiv&\frac{\lambda}{4}\left(\chi^{2}-\frac{m^{2}}{\lambda}\right)^{2}\\\
V_{1}\equiv&\frac{\mu}{3}\chi^{3}+\rm const.\end{split}$ (68)

where $V_{0}$ is the $Z_{2}$ symmetric part of the potential and the constant
on the last line makes $V_{1}(\chi_{-})\approx 0$.

Therefore, an approximate solution for the bubble wall profile can be found by
solving

$\frac{1}{2}(\partial_{n}\chi)^{2}\approx V_{0}(\chi(n)).$ (69)

subject to boundary conditions $\chi(\pm\infty)=\pm m/\sqrt{\lambda}$. The
solution is

$\chi=\frac{m}{\sqrt{\lambda}}\tanh\left(\frac{n}{\ell}\right),$ (70)

where $\ell\equiv\sqrt{2}/m$. This allows us to calculate $\sigma$ explicitly:

$\sigma=\frac{2\sqrt{2}}{3}\frac{m^{3}}{\lambda}\,.$ (71)

Having calculated the surface tension, we next consider the mean extrinsic
curvature of the surface $\Sigma$:

$K=-(g^{\mu\nu}\nabla_{\mu}\tilde{N}_{\nu})|_{n=0}.$ (72)

Here $\tilde{N}_{\nu}(n,\xi^{a})$ is an arbitrary extension of $N_{\nu}$ to a
unit vector field that agrees with $N_{\nu}$ on $\Sigma$. In our case,
$\Sigma$ is the surface of a spherical bubble nucleated at some time $t_{I}$
with a radius

$r(t,t_{I})=\int\displaylimits_{t_{I}}^{t}\frac{dt^{\prime}v_{w}(t^{\prime})}{a(t^{\prime})}+\frac{r_{c}}{a(t_{I})}$
(73)

where the critical radius $r_{c}$ is the initial physical radius of the bubble
as mentioned in Section II. Points $x^{\mu}=(t,r,\theta,\phi)$ on $\Sigma$
satisfy

$0=F(x^{\mu})\equiv r-r(t,t_{I})$ (74)

The unit normal vector $N_{\mu}$ is then given by

$N_{\mu}=\frac{\partial_{\mu}F}{\sqrt{|\partial_{\mu}F\partial^{\mu}F|}}=\gamma_{w}(-v_{w},a(t)\mathbf{\hat{r}}).$
(75)

and

$n=(x^{\mu}-X^{\mu})N_{\mu}=\frac{F(x^{\mu})}{\sqrt{|\partial_{\mu}F\partial^{\mu}F|}}=a(t)\gamma_{w}(r-r(t,t_{I})).$
(76)

Using Eqs. (72) and (75), we calculate the mean curvature of $\Sigma$. Since
$N_{\mu}$ does not depend explicitly on $r$ but $n$ depends on $t$ and $r$, we
can extend the domain of $N_{\mu}$ away from $n=0$ and take
$\tilde{N}_{\mu}\equiv N_{\mu}$. Then the mean curvature can be written as

$K=-\partial_{t}(\gamma_{w}v_{w})-3H\gamma_{w}v_{w}-\frac{2\gamma_{w}}{a(t)r(t,t_{I})}.$
(77)

Substituting this result into Eq. (65) leads to a dynamic equation for the
expansion of the bubble wall

$\partial_{t}(\gamma_{w}v_{w})+3H\gamma_{w}v_{w}+\frac{2\gamma_{w}}{a(t)r(t,t_{I})}=\frac{\Delta
V}{\sigma}.$ (78)

When a bubble first nucleates its physical size is the critical radius
$r_{c}\equiv 3S_{1}/\Delta V=3\sigma/\Delta V$ and its wall is at rest Coleman
(1977); Linde (1983). Putting $a(t_{I})r(t_{I},t_{I})=r_{c}$ and $v_{w}=0$ in
Eq. (78), we can see that, the bubble wall will accelerate and expand until it
reaches the terminal velocity given by

$\left.(\gamma_{w}v_{w})\right|_{\infty}=\frac{\Delta V}{3\sigma H_{\rm
inf}}=\frac{1}{r_{c}H_{\rm inf}}\gg 1,$ (79)

where we have used the condition in Eq. (5). This implies the terminal wall
velocity is $v_{w}\approx 1$, with $\gamma_{w}\approx 1/(r_{c}H_{\rm inf})$,
after a time $t-t_{I}\sim H_{\rm inf}^{-1}$.

##  Appendix B Two-Bubble Terms

Here we quantify the errors from our calculation of the power spectrum of
$\delta_{\rm\chi}$ in Section III. The power spectrum is determined by the
two-point function which can be written as the double sum over bubbles given
in Eq. (III). The double sum can be split up as

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle=\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle^{(1)}+\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle^{(2)}$
(80)

where the first term on the right is the contribution from terms where $I=J$
and the second term is the contribution from terms where $I\neq J$. In Section
III, we calculated the power spectrum from the $I=J$ terms only and here we
will approximate the correction from including the $I\neq J$ terms.

The correction to the 2-point function from the $I\neq J$ terms can be written
as

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle^{(2)}=e^{2t_{e}/\tau_{\rm
PT}}\frac{(4\pi)^{2}}{k^{3}{k^{\prime}}^{3}}(N^{2}-N)\int
d^{4}xd^{4}x^{\prime}p_{2}(x,x^{\prime})e^{-i(\mathbf{k}\cdot\mathbf{x}+\mathbf{k^{\prime}}\cdot\mathbf{x^{\prime}})}\mathcal{A}(kr(t))\mathcal{A}(k^{\prime}r(t^{\prime}))$
(81)

where $p_{2}(x,x^{\prime})$ is the joint probability density function for
bubbles nucleated at space-time coordinates $x$ and $x^{\prime}$. For
$t<t^{\prime}$, we will write the joint probability density function as

$p_{2}(x,x^{\prime})=p(x|x^{\prime})p_{1}(x^{\prime})$ (82)

where

$p_{1}(x^{\prime})=\frac{1}{N}\Gamma a(t^{\prime})^{3}p_{\rm
false}(t^{\prime})$ (83)

and the conditional probability can be written as

$p(x|x^{\prime})=\frac{1}{N}\Gamma a(t)^{3}\Theta(\Delta x-r(t))p_{\rm
false}(t|t^{\prime},\Delta x)\quad\quad t<t^{\prime}.$ (84)

Here, $\Delta x\equiv|\mathbf{x}-\mathbf{x^{\prime}}|$ and

$\begin{split}p_{\rm false}(t|t^{\prime},\Delta x)=&e^{-J(t,\Delta
x)}\quad\quad t<t^{\prime}\\\ J(t,\Delta
x)=&\int\displaylimits_{0}^{t}dt^{\prime\prime}\Gamma
a(t^{\prime\prime})^{3}\left[\mathcal{V}_{\mathcal{H}}(t^{\prime\prime})-\mathcal{V}_{\mathcal{O}}(t^{\prime\prime},\Delta
x)\right]\end{split}$ (85)

where $\mathcal{V_{H}}(t^{\prime\prime})$ is the comoving Hubble volume at
$t^{\prime\prime}$ and $\mathcal{V_{O}}(t^{\prime\prime},\Delta x)$ is the
overlap volume between two spheres of comoving radius $r(t^{\prime\prime})$ at
$\mathbf{x}$ and $\mathbf{x^{\prime}}$. $J(t,\Delta x)\leq t/\tau_{\rm PT}$
where equality is acheived if $\mathcal{V_{O}}(t^{\prime\prime},\Delta x)=0$
for $t^{\prime\prime}\in[0,t]$. The $\Theta$ function in Eq. (84) ensures that
the probability goes to zero for a bubble nucleated at coordinates
$\mathbf{x},t$ that encloses the point $\mathbf{x^{\prime}}$. For
$t>t^{\prime}$ we use $p_{2}(x,x^{\prime})=p(x^{\prime}|x)p_{1}(x)$ rather
than Eq. (82) and exchange $x\leftrightarrow x^{\prime}$ in Eq. (84).

Since $p_{2}(x,x^{\prime})$ only depends on $\mathbf{x}$ and
$\mathbf{x^{\prime}}$ through $\Delta x$, we change integration variables to

$\begin{split}\mathbf{X}\equiv&\frac{\mathbf{x}+\mathbf{x^{\prime}}}{2}\\\
\mathbf{\Delta x}\equiv&\mathbf{x}-\mathbf{x^{\prime}}\end{split}$ (86)

and write
$\mathbf{k}\cdot\mathbf{x}+\mathbf{k^{\prime}}\cdot\mathbf{x^{\prime}}=(\mathbf{k}+\mathbf{k^{\prime}})\cdot\mathbf{X}+(\mathbf{k}-\mathbf{k^{\prime}})\cdot\mathbf{\Delta
x}/2$. Integration over $\mathbf{X}$ is then trivial. If $t_{e}/\tau_{\rm
PT}\ll 1$, we can safely set $\exp{\left[-J(t,\Delta x)\right]}\to 1$ and
$\exp{(-t/\tau_{\rm PT})}\to 1$, leading to the expression

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle^{(2)}=2\frac{(4\pi)^{2}}{k^{6}}(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k^{\prime}})\Gamma^{2}\int\displaylimits_{0}^{t_{e}}dt\int\displaylimits_{t}^{t_{e}}dt^{\prime}a(t)^{3}a(t^{\prime})^{3}\mathcal{A}(kr(t))\mathcal{A}(kr(t^{\prime}))\int
d^{3}(\Delta x)e^{-i\mathbf{k}\cdot\mathbf{\Delta x}}\Theta(\Delta x-r(t))$
(87)

The $\mathbf{\Delta x}$ integral leads to a term proportional to
$\delta^{3}(\mathbf{k})$ and a finite term. We ignore the term proportional to
the delta function and change to dimensionless integration variables $u\equiv
kr(t)$ and $u^{\prime}\equiv kr(t^{\prime})$ to obtain

$\langle\delta_{\chi}(\mathbf{k})\delta_{\chi}(\mathbf{k^{\prime}})\rangle^{(2)}=-\frac{128\pi^{3}}{k^{3}}(2\pi)^{3}\delta^{3}(\mathbf{k}+\mathbf{k^{\prime}})\gamma_{\rm
PT}^{2}\int\displaylimits_{kr_{e}}^{kr_{i}}duu^{-4}\mathcal{A}(u)^{2}\int\displaylimits_{kr_{e}}^{u}du^{\prime}{u^{\prime}}^{-4}\mathcal{A}(u^{\prime})+(\propto\delta^{3}(\mathbf{k}))$
(88)

This contribution to the two point function is proportional to $\gamma_{\rm
PT}^{2}$ so it is suppressed by an extra factor of $\gamma_{\rm PT}$ compared
to the result found in Section III. We checked numerically that this term can
be neglected for all values of $\gamma_{\rm PT}$ considered in this work.

##  Appendix C Curvature Produced by FOPT

In the main text we calculated the contribution of the FOPT to the
isocurvature power spectrum. In this section we will show that the FOPT can
also source curvature perturbations during inflation due to non-adiabatic
pressure in bubble walls. We will start by calculating the components of the
stress-energy tensor for the $\chi$ field and extract the non-adiabatic
pressure. We will then use our expression for the non-adiabatic pressure to
calculate the evolution of the comoving curvature perturbation during
inflation.

###  C.1 Stress-Energy Tensor

The stress-energy tensor for the $\chi$ field is given by the standard
expression for a scalar field

${T^{\mu}}_{\nu}=g^{\mu\rho}\partial_{\rho}\chi\partial_{\nu}\chi-\delta^{\mu}_{\nu}\left(\frac{1}{2}g^{\rho\sigma}\partial_{\rho}\chi\partial_{\sigma}\chi+V(\chi)\right)$
(89)

As discussed in Appendix A, when a bubble is nucleated, $\chi$ only depends on
the coordinate $n$ orthogonal to the wall (with unit vector $N_{\nu}$). As a
result, the stress-energy tensor for a single bubble can be written as

${T^{\mu}}_{\nu}=\left[N^{\mu}N_{\nu}-\frac{1}{2}\delta^{\mu}_{\nu}\right](\partial_{n}\chi)^{2}-\delta^{\mu}_{\nu}V(\chi).$
(90)

The definition of $n$ for a bubble nucleated at the origin at time $t_{I}$ is
given in Eq. (76). We split the potential into its two components
$V=V_{0}+V_{1}$ and substitute Eq. (69) for $V_{0}$. Then the stress energy
tensor can be written as

${T^{\mu}}_{\nu}=\left[N^{\mu}N_{\nu}-\delta^{\mu}_{\nu}\right](\partial_{n}\chi)^{2}-\delta^{\mu}_{\nu}V_{1}(\chi)$
(91)

The term proportional to $(\partial_{n}\chi)^{2}$ vanishes everywhere except
on the wall. The term proportional to $V_{1}$ is negligible inside the wall
compared to the first term but switches from $0$ for $n\lesssim 0$ to $+\Delta
V$ for $n\gtrsim 0$. This motivates us to define “wall” and “bulk”
contributions to the stress energy tensor

$\begin{split}{{(T_{\rm
wall})}^{\mu}}_{\nu}=&\left[N^{\mu}N_{\nu}-\delta^{\mu}_{\nu}\right](\partial_{n}\chi)^{2}\\\
{{(T_{\rm bulk})}^{\mu}}_{\nu}=&V_{1}(\chi(n))\end{split}$ (92)

where the full tensor is the sum of the two:

${T^{\mu}}_{\nu}={{(T_{\rm wall})}^{\mu}}_{\nu}+{{(T_{\rm
bulk})}^{\mu}}_{\nu}$ (93)

In this work, the most important properties of the stress tensor are the
energy density $\rho$ and the pressure $p$. Using Eq. (75), we can write the
“wall” and “bulk” components of the single bubble energy density as

$\begin{split}\rho_{\rm wall}=&-{{(T_{\rm
wall})}^{0}}_{0}=\gamma_{w}^{2}(\partial_{n}\chi)^{2}\\\ \rho_{\rm
bulk}=&-{{(T_{\rm bulk})}^{0}}_{0}=V_{1}(\chi(n))\end{split}$ (94)

and of the pressure as

$\begin{split}p_{\rm wall}=&\frac{1}{3}{{(T_{\rm
wall})}^{i}}_{i}=\gamma_{w}^{2}(v_{w}^{2}-2/3)(\partial_{n}\chi)^{2}=(v_{w}^{2}-2/3)\rho_{\rm
wall}\\\ p_{\rm bulk}=&\frac{1}{3}{{(T_{\rm
bulk})}^{i}}_{i}=-V_{1}(\chi(n))=-\rho_{\rm bulk}\end{split}$ (95)

####  C.1.1 Wall Energy Density

We next calculate $\bar{\rho}_{\rm wall}$ and $\delta\rho_{\rm
wall}(\mathbf{k},t)$ for the whole Universe at time $t$. To do so, we will sum
the first line of Eq. (94) over bubbles nucleated before $t$

$\rho_{\rm
wall}(\mathbf{x},t)=\gamma_{w}^{2}\sum_{I=1}^{N}\Theta(t-t_{I})\left.(\partial_{n}\chi)^{2}\right|_{n=n_{I}}$
(96)

where $n_{I}\equiv a(t)\gamma_{w}(|\mathbf{x}-\mathbf{x_{I}}|-r(t,t_{I}))$ is
the normal coordinate for bubble $I$, $\mathbf{x_{I}}$ is the center position
of the bubble, and for $\gamma_{w}$ we take the terminal value for all bubbles
given in Eq. (79). Using the solution for the profile of the bubble wall in
Eq. (70), the wall energy density can be written as

$\rho_{\rm wall}=\frac{\gamma_{w}\Delta
V}{4\ell}\sum_{I}\Theta(t-t_{I})\operatorname{sech}^{4}{(n_{I}/\ell)}$ (97)

The spatial average of the wall energy is

$\begin{split}\bar{\rho}_{\rm
wall}=&\frac{1}{\mathcal{V}}\int\displaylimits_{\mathcal{V}}d^{3}x\rho_{g}\\\
=&\frac{4\pi}{3}\frac{\Delta
V}{\mathcal{V}a(t)}\sum_{I=1}^{N}\Theta(t-t_{I})r(t,t_{I})^{2}.\end{split}$
(98)

We substitute Eq. (8) for $r(t,t_{I})$ and use the probability density in Eq.
(19) to convert the sum over bubbles into an integral over time

$\bar{\rho}_{\rm wall}=\frac{4\pi}{3}\gamma_{\rm PT}\Delta Ve^{-t/\tau_{PT}}.$
(99)

The perturbation $\rho_{\rm wall}-\bar{\rho}_{\rm wall}$ of the wall’s energy
density in Fourier space is

$\delta\rho_{\rm wall}=\int d^{3}x\rho_{\rm
wall}(\mathbf{k})e^{-i\mathbf{k}\cdot\mathbf{x}}+{\rm
const.}\times\delta^{3}(\mathbf{k}),$ (100)

where the second term is only non-zero for $\mathbf{k}=0$. This term can be
dropped, as only the non-zero wavenumbers enter our analysis. Assuming that
$r_{I}\gg\ell/(\gamma_{w}a)$ for all relevant bubbles, and that we are
interested in modes such that $k\ell/(\gamma_{w}a)\ll 1$, the integral over
$\mathbf{x}$ results in the closed form solution

$\delta\rho_{\rm wall}=\frac{4\pi}{3}\frac{\Delta
V}{Ha(t)v_{w}k^{2}}\sum_{I}\Theta(t-t_{I})\mathcal{B}(kr(t_{I}))e^{-i\mathbf{k}\cdot\mathbf{x_{I}}},$
(101)

where

$\mathcal{B}(y)\equiv y\sin{y}.$ (102)

####  C.1.2 Bulk Energy Density

The bulk energy density is given by the point-wise function of $\chi$ in the
second line of Eq. (94). In a region of space where bubble $I$ is the only
bubble, the behavior of $\chi$ is given by Eq. (70) with $n\to n_{I}$. Since
we only care about the bulk energy density on scales much larger than the
thickness of the wall, we can approximate

$\begin{split}\chi(n_{I})\approx&\frac{2m}{\sqrt{\lambda}}\left[\Theta(n_{I})-\frac{1}{2}\right]\\\
=&\frac{2m}{\sqrt{\lambda}}\left[\Theta(|\mathbf{x}-\mathbf{x_{I}}|-r(t,t_{I}))-\frac{1}{2}\right],\end{split}$
(103)

leading to

$\rho_{\mathrm{bulk},I}=V_{1}(\chi(n_{I}))\approx\Delta
V\left[1-\Theta(r(t,t_{I})-|\mathbf{x}-\mathbf{x_{I}}|)\right].$ (104)

where $\Delta V\equiv 2\mu m^{3}/(3\lambda^{3/2})$. The fraction of space in
true vacuum by the end of inflation is $\sim\gamma_{\rm PT}H_{\rm inf}t_{e}\ll
1$. As a result, most bubbles do not overlap and the bulk energy density for
the whole Universe can be approximated as

$\rho_{\rm bulk}=\Delta
V\left[1-\sum_{I}\Theta(t-t_{I})\Theta(r(t,t_{I})-|\mathbf{x}-\mathbf{x_{I}}|)\right].$
(105)

Spatially averaging this bulk density results in

$\bar{\rho}_{\rm bulk}=\Delta Vp_{\rm false}(t)=\Delta Ve^{-t/\tau_{PT}},$
(106)

which dominates over the wall contribution, Eq. (106). The density
perturbation for the bulk contribution is

$\begin{split}\delta\rho_{\rm bulk}(\mathbf{k},t)=&-\Delta
V\sum_{I}\Theta(t-t_{I})\\\ \times&\int
d^{3}xe^{-i\mathbf{k}\cdot\mathbf{x}}\Theta(r(t,t_{I})-|\mathbf{x}-\mathbf{x_{I}}|)\\\
+&{\rm const.}\times\delta^{3}(\mathbf{k}).\end{split}$ (107)

Again dropping the second term and performing the integral over $\mathbf{x}$

$\delta\rho_{\rm bulk}(\mathbf{k},t)=\frac{4\pi\Delta
V}{k^{3}}\sum_{I}\Theta(t-t_{I})\mathcal{A}(kr(t,t_{I}))e^{-i\mathbf{k}\cdot\mathbf{x_{I}}},$
(108)

where we have defined

$\mathcal{A}(y)\equiv y\cos{y}-\sin{y}.$ (109)

###  C.2 Curvature Power Spectrum

In the standard inflationary model, the inflaton field produces subhorizon
quantum fluctuations in the comoving curvature perturbation (denoted as
$\mathcal{R}$) with a power spectrum given by Malik and Wands (2009)

$P_{\mathcal{R_{\phi}}}(k)\equiv
P_{\mathcal{R}}(k,t_{k})=\left.\frac{H^{4}}{4\pi^{2}\dot{\bar{\phi}}^{2}}\right|_{t_{k}},$
(110)

where $t_{k}$ is the time that wave-number $k$ exits the horizon
($a(t_{k})H_{\rm inf}\sim k$). In the class of models of FOPT considered in
this paper, bubbles from the FOPT quickly become larger than the horizon. As a
result, their effect on the subhorizon dynamics is mostly from modifying the
background energy density $\bar{\rho}_{\chi}$, which can be neglected since we
work in the region
$\bar{\rho}_{\chi}\ll(\bar{\rho}_{\chi}+\bar{\rho}_{\phi})$. Thus, the
subhorizon dynamics are almost unaltered from standard inflation.

Outside of the horizon, the curvature perturbation evolves according to Gordon
_et al._ (2000); Malik and Wands (2005)

$\dot{\mathcal{R}}=\frac{H}{\bar{\rho}+\bar{p}}\delta p_{\rm
nad}+\mathcal{O}\left(\frac{k^{2}}{a^{2}H^{2}}\right),$ (111)

where $\delta p_{\rm nad}$ is the non-adiabatic pressure perturbation

$\delta p_{\rm nad}\equiv\delta
p-\frac{\dot{\bar{p}}}{\dot{\bar{\rho}}}\delta\rho.$ (112)

The non-adiabatic pressure of the inflaton is suppressed outside the horizon
Wands (2008), while for the $\chi$ field it is sourced by the energy of the
bubble wall. Setting $v_{w}=1$ and using results in Appendix C.1, the
superhorizon evolution of $\mathcal{R}$ is

$\dot{\mathcal{R}}\approx\frac{4}{3}\frac{H}{\dot{\bar{\phi}}^{2}}\delta\rho_{\rm
wall}~{}~{}~{}~{}(t>t_{k}),$ (113)

where $\delta\rho_{\rm wall}$ is defined in Eq. (100). The curvature
$\mathcal{R}$ at the end of inflation ($t=t_{e}$) is then

$\begin{split}&\mathcal{R}(\mathbf{k},t_{e})=\mathcal{R}(\mathbf{k},t_{k})+\frac{16\pi}{9}\frac{\Delta
V}{\dot{\bar{\phi}}^{2}}\frac{1}{k^{2}}\\\
\times&\sum_{I}kr(t_{I})\sin{(kr(t_{I}))}e^{-i\mathbf{k}\cdot\mathbf{x_{I}}}\int_{t_{k}}^{t_{e}}dt\frac{1}{a(t)}\Theta(t-t_{I}),\end{split}$
(114)

where $I$ runs over all nucleated bubbles.

To calculate the power spectrum of $\mathcal{R}$ at the end of inflation, we
follow the same procedure as in Eq. (III): we ignore cross-terms in the double
sum, write the ensemble average as an integral, integrate over spatial
coordinates, and change the $t$ integration variable to $u\equiv kr(t)$. The
curvature power spectrum is then given as:

$\begin{split}&P_{\cal R}(k,t_{e})-P_{\cal
R}(k,t_{k})=\frac{128}{81}\left(\frac{\Delta
V}{\dot{\bar{\phi}}^{2}}\right)^{2}\gamma_{\rm PT}\\\
\times&\left[\frac{1}{2}\left(kr(t_{k})-\frac{1}{2}\sin{(2kr(t_{k}))}\right)+\int\displaylimits_{kr(t_{k})}^{kr_{i}}du\frac{\sin^{2}{u}}{u^{2}}\right].\\\
\quad&\end{split}$ (115)

The first term in the square brackets comes from bubbles nucleated after
$t_{k}$ while the second term comes from bubbles nucleated before $t_{k}$. If
FOPT starts after $t_{k}$ ($t_{i}>t_{k}$), the second term in the brackets of
Eq. (115) is absent since there are no bubbles nucleated before $t_{k}$. The
power spectrum in this case is simply given by the first term with $t_{k}\to
t_{i}$. Since $r(t)$ is the radius of the horizon at $t$, then the definition
of $t_{k}$ implies $r(t_{k})\sim 1/k$ so the first term in the brackets is
independent of $k$.

Given our definition of $P_{\rm ad}(k)$ in Eq. (30), we find that $P_{\rm
ad}(k)=P_{\cal R}(k,t_{e})$. There is also a contribution to the curvature
power spectrum from bubble percolation when the phase transition completes
around $T_{*}$ Freese and Winkler (2023); Liu  _et al._ (2023); Elor  _et al._
(2023). As mentioned in section IV.1, we can neglect those contributions when
studying CMB observables since we consider scenarios where $T_{*}\gg T_{\rm
CMB}$. Bringing all of this together, we can obtain $P_{\rm ad}(k)$:

$P_{\rm ad}(k)\approx P_{\cal R_{\phi}}(k)+\frac{128}{81}\left(\frac{\Delta
V}{\dot{\bar{\phi}}^{2}}\right)^{2}\gamma_{\rm
PT}\begin{dcases}\frac{1}{3}(kr_{i})^{3}&kr_{i}\ll 1\\\ 0.95&kr_{i}\gg
1\end{dcases}.$ (116)

Here the inflaton contribution $P_{\cal{R}_{\phi}}(k)$ is usually
parameterized as $P_{\cal R_{\phi}}(k)=A_{s}(k/k_{\rm pivot})^{n_{s}-1}$,
where $A_{s}$ and $n_{s}$ are the amplitude and spectral index respectively.
The contribution from FOPT has a distinct $k$ dependence: for $kr_{i}\ll 1$ it
goes as $(kr_{i})^{3}$ and for $kr_{i}\to\infty$ it approaches a constant.

Since we focus on isocurvature effects of the FOPT in this work, we choose the
model parameters such that the contribution of FOPT to $P_{\rm ad}(k)$ (second
term of Eq. (116)) is negligible. This translates to the requirement that:

$\left(\frac{\Delta V}{\dot{\bar{\phi}}^{2}}\right)^{2}\gamma_{\rm PT}\ll
A_{s}.$ (117)

Moreover, we also require the effect of FOPT on CMB angular power spectrum
(e.g., $C_{\ell}^{TT}$) is dominantly from isocurvature. Given the scaling of
the isocurvature effects at the end of Section IV.1, this condition is given
as:

$\left(\frac{\Delta V}{\dot{\bar{\phi}}^{2}}\right)^{2}\ll R_{\rm dr}^{2}.$
(118)

##  Appendix D Initial Conditions for the Dark Radiation isocurvature mode

In this section, we will present the full set of initial conditions for
perturbations in the DR isocurvature mode in the synchronous gauge and in
terms of the conformal time $\tau$ (see also Ghosh  _et al._ (2022)). We
assume DR is free-streaming, analogous to Standard Model neutrinos. We keep
terms up to ${\cal O}((k\tau)^{2})$. For terms that are zero up to this order,
we retain the leading non-vanishing term. The initial conditions are

$\displaystyle\eta^{\rm iso}=\frac{-R_{\rm dr}+R_{\rm dr}^{2}+R_{\rm
dr}R_{\nu}}{6(1-R_{\rm dr})(15+4R_{\rm dr}+4R_{\nu})}(k\tau)^{2}$
$\displaystyle h^{\rm iso}=\frac{R_{\rm dr}}{40(1-R_{\rm
dr})}\omega_{b}k^{2}\tau^{3}$ $\displaystyle\delta_{\gamma}^{\rm
iso}=\delta_{\nu}^{\rm iso}=\frac{-R_{\rm dr}}{1-R_{\rm
dr}}\left(1-\frac{1}{6}(k\tau)^{2}\right)$ $\displaystyle\delta_{b}^{\rm
iso}=\frac{R_{\rm dr}}{8(1-R_{\rm dr})}(k\tau)^{2}$
$\displaystyle\delta_{c}^{\rm iso}=\frac{-R_{\rm dr}}{80(1-R_{\rm
dr})}\omega_{b}k^{2}\tau^{3}$ (119) $\displaystyle\delta_{\rm dr}^{\rm
iso}=1-\frac{1}{6}(k\tau)^{2}$ $\displaystyle\theta_{\gamma}^{\rm
iso}=\theta_{b}^{\rm iso}=\frac{-R_{\rm dr}}{4(1-R_{\rm dr})}k^{2}\tau$
$\displaystyle\theta_{\rm dr}^{\rm iso}=\frac{1}{4}k^{2}\tau$
$\displaystyle\sigma_{\nu}^{\rm iso}=\frac{-19R_{\rm dr}}{30(1-R_{\rm
dr})(15+4R_{\rm dr}+4R_{\nu})}k^{2}\tau^{2}$ $\displaystyle\sigma_{\rm
dr}^{\rm iso}=\frac{15-15R_{\rm dr}+4R_{\nu}}{30(1-R_{\rm dr})(15+4R_{\rm
dr}+4R_{\nu})}k^{2}\tau^{2},$

where $\omega_{b}\equiv\sqrt{8\pi G/3}\,a(\tau_{\rm
ini})\bar{\rho}_{b}(\tau_{\rm ini})/\sqrt{\bar{\rho}_{r}(\tau_{\rm ini})}$ and
$\bar{\rho}_{b}$,$\bar{\rho}_{r}$ is the background baryon and radiation
density respectively. At the time initial conditions are set,
$\omega_{b}\tau\ll 1$ can be treated as an expansion parameter.

## References

  * Kuzmin  _et al._ (1985) V. Kuzmin, V. Rubakov, and M. Shaposhnikov, Physics Letters B 155, 36 (1985).
  * Shaposhnikov (1986) M. E. Shaposhnikov, JETP Lett. 44, 465 (1986).
  * Shaposhnikov (1987) M. E. Shaposhnikov, Nucl. Phys. B 287, 757 (1987).
  * Cohen  _et al._ (1993) A. G. Cohen, D. B. Kaplan, and A. E. Nelson, Ann. Rev. Nucl. Part. Sci. 43, 27 (1993), arXiv:hep-ph/9302210 .
  * Trodden (1999) M. Trodden, Rev. Mod. Phys. 71, 1463 (1999), arXiv:hep-ph/9803479 .
  * Riotto and Trodden (1999) A. Riotto and M. Trodden, Ann. Rev. Nucl. Part. Sci. 49, 35 (1999), arXiv:hep-ph/9901362 .
  * Dine and Kusenko (2003) M. Dine and A. Kusenko, Rev. Mod. Phys. 76, 1 (2003), arXiv:hep-ph/0303065 .
  * Morrissey and Ramsey-Musolf (2012) D. E. Morrissey and M. J. Ramsey-Musolf, New J. Phys. 14, 125003 (2012), arXiv:1206.2942 [hep-ph] .
  * Creminelli  _et al._ (2002) P. Creminelli, A. Nicolis, and R. Rattazzi, JHEP 03, 051 (2002), arXiv:hep-th/0107141 .
  * Randall and Servant (2007) L. Randall and G. Servant, JHEP 05, 054 (2007), arXiv:hep-ph/0607158 .
  * Nardini  _et al._ (2007) G. Nardini, M. Quiros, and A. Wulzer, JHEP 09, 077 (2007), arXiv:0706.3388 [hep-ph] .
  * Konstandin  _et al._ (2010) T. Konstandin, G. Nardini, and M. Quiros, Phys. Rev. D 82, 083513 (2010), arXiv:1007.1468 [hep-ph] .
  * Konstandin and Servant (2011) T. Konstandin and G. Servant, JCAP 12, 009 (2011), arXiv:1104.4791 [hep-ph] .
  * Bunk  _et al._ (2018) D. Bunk, J. Hubisz, and B. Jain, Eur. Phys. J. C 78, 78 (2018), arXiv:1705.00001 [hep-ph] .
  * Baratella  _et al._ (2019) P. Baratella, A. Pomarol, and F. Rompineve, JHEP 03, 100 (2019), arXiv:1812.06996 [hep-ph] .
  * Bruggisser  _et al._ (2018) S. Bruggisser, B. Von Harling, O. Matsedonskyi, and G. Servant, JHEP 12, 099 (2018), arXiv:1804.07314 [hep-ph] .
  * Megías  _et al._ (2018) E. Megías, G. Nardini, and M. Quirós, JHEP 09, 095 (2018), arXiv:1806.04877 [hep-ph] .
  * Agashe  _et al._ (2020) K. Agashe, P. Du, M. Ekhterachian, S. Kumar, and R. Sundrum, JHEP 05, 086 (2020), arXiv:1910.06238 [hep-ph] .
  * Fujikura  _et al._ (2020) K. Fujikura, Y. Nakai, and M. Yamada, JHEP 02, 111 (2020), arXiv:1910.07546 [hep-ph] .
  * Agashe  _et al._ (2021) K. Agashe, P. Du, M. Ekhterachian, S. Kumar, and R. Sundrum, JHEP 02, 051 (2021), arXiv:2010.04083 [hep-th] .
  * Ares  _et al._ (2020) F. R. Ares, M. Hindmarsh, C. Hoyos, and N. Jokela, JHEP 21, 100 (2020), arXiv:2011.12878 [hep-th] .
  * Agrawal and Nee (2021) P. Agrawal and M. Nee, JHEP 10, 105 (2021), arXiv:2103.05646 [hep-ph] .
  * Levi  _et al._ (2023) N. Levi, T. Opferkuch, and D. Redigolo, JHEP 02, 125 (2023), arXiv:2212.08085 [hep-ph] .
  * Csáki  _et al._ (2023) C. Csáki, M. Geller, Z. Heller-Algazi, and A. Ismail, JHEP 06, 202 (2023), arXiv:2301.10247 [hep-ph] .
  * Eröncel  _et al._ (2023) C. Eröncel, J. Hubisz, S. J. Lee, G. Rigo, and B. Sambasivam, (2023), arXiv:2305.03773 [hep-ph] .
  * Mishra and Randall (2023) R. K. Mishra and L. Randall, (2023), arXiv:2309.10090 [hep-ph] .
  * Mishra and Randall (2024) R. K. Mishra and L. Randall, (2024), arXiv:2401.09633 [hep-ph] .
  * Bai and Schwaller (2014) Y. Bai and P. Schwaller, Phys. Rev. D 89, 063522 (2014), arXiv:1306.4676 [hep-ph] .
  * Kribs and Neil (2016) G. D. Kribs and E. T. Neil, Int. J. Mod. Phys. A 31, 1643004 (2016), arXiv:1604.04627 [hep-ph] .
  * Cline (2022) J. M. Cline, SciPost Phys. Lect. Notes 52, 1 (2022), arXiv:2108.10314 [hep-ph] .
  * Shelton and Zurek (2010) J. Shelton and K. M. Zurek, Phys. Rev. D 82, 123512 (2010), arXiv:1008.1997 [hep-ph] .
  * Dutta and Kumar (2011) B. Dutta and J. Kumar, Physics Letters B 699, 364 (2011).
  * Holthausen  _et al._ (2013) M. Holthausen, J. Kubo, K. S. Lim, and M. Lindner, JHEP 12, 076 (2013), arXiv:1310.4423 [hep-ph] .
  * Schwaller (2015) P. Schwaller, Phys. Rev. Lett. 115, 181101 (2015), arXiv:1504.07263 [hep-ph] .
  * Hall  _et al._ (2023) E. Hall, T. Konstandin, R. McGehee, and H. Murayama, Phys. Rev. D 107, 055011 (2023), arXiv:1911.12342 [hep-ph] .
  * Bai and Korwar (2022) Y. Bai and M. Korwar, Phys. Rev. D 105, 095015 (2022), arXiv:2109.14765 [hep-ph] .
  * Bottaro  _et al._ (2021) S. Bottaro, M. Costa, and O. Popov, JHEP 11, 055 (2021), arXiv:2104.14244 [hep-ph] .
  * Pasechnik  _et al._ (2023) R. Pasechnik, M. Reichert, F. Sannino, and Z.-W. Wang, (2023), arXiv:2309.16755 [hep-ph] .
  * Karwal and Kamionkowski (2016) T. Karwal and M. Kamionkowski, Phys. Rev. D 94, 103523 (2016), arXiv:1608.01309 [astro-ph.CO] .
  * Mörtsell and Dhawan (2018) E. Mörtsell and S. Dhawan, JCAP 09, 025 (2018), arXiv:1801.07260 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Niedermann and Sloth (2021) F. Niedermann and M. S. Sloth, Phys. Rev. D 103, L041303 (2021), arXiv:1910.10739 [astro-ph.CO] .
  * Niedermann and Sloth (2020) F. Niedermann and M. S. Sloth, Phys. Rev. D 102, 063527 (2020), arXiv:2006.06686 [astro-ph.CO] .
  * Niedermann and Sloth (2022) F. Niedermann and M. S. Sloth, Phys. Rev. D 105, 063509 (2022), arXiv:2112.00770 [hep-ph] .
  * Kosowsky  _et al._ (1992a) A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. Lett. 69, 2026 (1992a).
  * Kosowsky  _et al._ (1992b) A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. D 45, 4514 (1992b).
  * Kosowsky and Turner (1993) A. Kosowsky and M. S. Turner, Phys. Rev. D 47, 4372 (1993), arXiv:astro-ph/9211004 .
  * Kamionkowski  _et al._ (1994) M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D 49, 2837 (1994), arXiv:astro-ph/9310044 .
  * Caprini  _et al._ (2016) C. Caprini _et al._ , JCAP 04, 001 (2016), arXiv:1512.06239 [astro-ph.CO] .
  * Caprini  _et al._ (2020) C. Caprini _et al._ , JCAP 03, 024 (2020), arXiv:1910.13125 [astro-ph.CO] .
  * Caldwell  _et al._ (2022) R. Caldwell _et al._ , Gen. Rel. Grav. 54, 156 (2022), arXiv:2203.07972 [gr-qc] .
  * Agazie  _et al._ (2023) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Afzal  _et al._ (2023) A. Afzal _et al._ (NANOGrav), Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023) J. Antoniadis _et al._ (EPTA), (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Reardon  _et al._ (2023) D. J. Reardon _et al._ , Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Xu  _et al._ (2023) H. Xu _et al._ , Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Auclair  _et al._ (2023) P. Auclair _et al._ (LISA Cosmology Working Group), Living Rev. Rel. 26, 5 (2023), arXiv:2204.05434 [astro-ph.CO] .
  * Kawamura  _et al._ (2006) S. Kawamura _et al._ , Class. Quant. Grav. 23, S125 (2006).
  * Kawamura  _et al._ (2021) S. Kawamura _et al._ , PTEP 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Harry  _et al._ (2006) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner, and E. S. Phinney, Class. Quant. Grav. 23, 4887 (2006), [Erratum: Class.Quant.Grav. 23, 7361 (2006)].
  * Froustey  _et al._ (2020) J. Froustey, C. Pitrou, and M. C. Volpe, JCAP 12, 015 (2020), 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, JCAP 04, 073 (2021), arXiv:2012.02726 [hep-ph] .
  * Akita and Yamaguchi (2020) K. Akita and M. Yamaguchi, JCAP 08, 012 (2020), arXiv:2005.07047 [hep-ph] .
  * Bucher  _et al._ (2000) M. Bucher, K. Moodley, and N. Turok, Phys. Rev. D 62, 083508 (2000), arXiv:astro-ph/9904231 .
  * Wands  _et al._ (2000) D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, Phys. Rev. D 62, 043527 (2000), arXiv:astro-ph/0003278 .
  * Gordon  _et al._ (2000) C. Gordon, D. Wands, B. A. Bassett, and R. Maartens, Phys. Rev. D 63, 023506 (2000), arXiv:astro-ph/0009131 .
  * Lyth  _et al._ (2003) D. H. Lyth, C. Ungarelli, and D. Wands, Phys. Rev. D 67, 023503 (2003), arXiv:astro-ph/0208055 .
  * Malik and Wands (2005) K. A. Malik and D. Wands, JCAP 02, 007 (2005), arXiv:astro-ph/0411703 .
  * Wands (2008) D. Wands, Lect. Notes Phys. 738, 275 (2008), arXiv:astro-ph/0702187 .
  * Freese and Winkler (2023) K. Freese and M. W. Winkler, Phys. Rev. D 107, 083522 (2023), arXiv:2302.11579 [astro-ph.CO] .
  * Elor  _et al._ (2023) G. Elor, R. Jinno, S. Kumar, R. McGehee, and Y. Tsai, (2023), arXiv:2311.16222 [hep-ph] .
  * Akrami  _et al._ (2020a) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020a), arXiv:1807.06211 [astro-ph.CO] .
  * Ghosh  _et al._ (2022) S. Ghosh, S. Kumar, and Y. Tsai, JCAP 05, 014 (2022), arXiv:2107.09076 [astro-ph.CO] .
  * Adshead  _et al._ (2020) P. Adshead, G. Holder, and P. Ralegankar, JCAP 09, 016 (2020), arXiv:2006.01165 [astro-ph.CO] .
  * Guth and Weinberg (1983) A. H. Guth and E. J. Weinberg, Nuclear Physics B 212, 321 (1983).
  * Liddle and Wands (1991) A. R. Liddle and D. Wands, Monthly Notices of the Royal Astronomical Society 253, 637 (1991).
  * Turner  _et al._ (1992) M. S. Turner, E. J. Weinberg, and L. M. Widrow, Phys. Rev. D 46, 2384 (1992).
  * Copeland  _et al._ (1994) E. J. Copeland, A. R. Liddle, D. H. Lyth, E. D. Stewart, and D. Wands, Phys. Rev. D 49, 6410 (1994), arXiv:astro-ph/9401011 .
  * Baccigalupi and Perrotta (2000) C. Baccigalupi and F. Perrotta, Mon. Not. Roy. Astron. Soc. 314, 1 (2000), arXiv:astro-ph/9911530 .
  * Barir  _et al._ (2023) J. Barir, M. Geller, C. Sun, and T. Volansky, Phys. Rev. D 108, 115016 (2023), arXiv:2203.00693 [hep-ph] .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 07, 034 (2011), arXiv:1104.2933 [astro-ph.CO] .
  * Lesgourgues (2011) J. Lesgourgues, (2011), arXiv:1104.2932 [astro-ph.IM] .
  * Lesgourgues and Tram (2011) J. Lesgourgues and T. Tram, JCAP 09, 032 (2011), arXiv:1104.2935 [astro-ph.CO] .
  * Akrami  _et al._ (2020b) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A9 (2020b), arXiv:1905.05697 [astro-ph.CO] .
  * Montandon  _et al._ (2021) T. Montandon, G. Patanchon, and B. van Tent, JCAP 01, 004 (2021), arXiv:2007.05457 [astro-ph.CO] .
  * Linde (1983) A. D. Linde, Nucl. Phys. B216, 421 (1983), [Erratum: Nucl. Phys.B223,544(1983)].
  * Coleman (1977) S. R. Coleman, Phys. Rev. D15, 2929 (1977), [Erratum: Phys. Rev.D16,1248(1977)].
  * Callan and Coleman (1977) C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D 16, 1762 (1977).
  * Coleman and De Luccia (1980) S. Coleman and F. De Luccia, Phys. Rev. D 21, 3305 (1980).
  * Ma and Bertschinger (1995) C.-P. Ma and E. Bertschinger, Astrophys. J. 455, 7 (1995), arXiv:astro-ph/9506072 .
  * Savelainen  _et al._ (2013) M. Savelainen, J. Valiviita, P. Walia, S. Rusak, and H. Kurki-Suonio, Phys. Rev. D 88, 063010 (2013), arXiv:1307.4398 [astro-ph.CO] .
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 02, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Brinckmann and Lesgourgues (2019) T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 [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, 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, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statist. Sci. 7, 457 (1992).
  * Lewis (2019) A. Lewis, (2019), arXiv:1910.13970 [astro-ph.IM] .
  * 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, Phys. Rept. 984, 1 (2022), arXiv:2107.10291 [astro-ph.CO] .
  * Mégevand and Membiela (2023) A. Mégevand and F. A. Membiela, JCAP 06, 007 (2023), arXiv:2302.13349 [gr-qc] .
  * Malik and Wands (2009) K. A. Malik and D. Wands, Phys. Rept. 475, 1 (2009), arXiv:0809.4944 [astro-ph] .
  * Liu  _et al._ (2023) J. Liu, L. Bian, R.-G. Cai, Z.-K. Guo, and S.-J. Wang, Phys. Rev. Lett. 130, 051001 (2023), arXiv:2208.14086 [astro-ph.CO] .