# Finite Bubble Statistics Constrain Late Cosmological Phase Transitions

Gilly Elor  Weinberg Institute, Department of Physics, University of Texas at
Austin, Austin, TX 78712, USA  Ryusuke Jinno  Research Center for the Early
Universe (RESCEU), Graduate School of Science, The University of Tokyo, Tokyo
113-0033, Japan  Soubhik Kumar  Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA  Berkeley
Center for Theoretical Physics, Department of Physics, University of
California, Berkeley, CA 94720, USA  Theoretical Physics Group, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USA  Robert McGehee  William
I. Fine Theoretical Physics Institute, School of Physics and Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  Leinweber Center for
Theoretical Physics, Department of Physics,
University of Michigan, Ann Arbor, MI 48109, USA  Yuhsin Tsai  Department of
Physics, University of Notre Dame, IN 46556, USA

###### Abstract

We consider first order cosmological phase transitions (PT) happening at late
times, below Standard Model (SM) temperatures $T_{\rm PT}\lesssim$ GeV. The
inherently stochastic nature of bubble nucleation and the finite number of
bubbles associated with a late-time PT lead to superhorizon fluctuations in
the PT completion time. We compute how such fluctuations eventually source
curvature fluctuations with universal properties, independent of the
microphysics of the PT dynamics. Using Cosmic Microwave Background (CMB) and
Large Scale Structure (LSS) measurements, we constrain the energy released in
a dark-sector PT. For 0.1 eV $\lesssim T_{\rm PT}\lesssim$ keV this constraint
is stronger than both the current bound from additional neutrino species
$\Delta N_{\rm eff}$, and in some cases, even CMB-S4 projections. Future
measurements of CMB spectral distortions and pulsar timing arrays will also
provide competitive sensitivity for keV $\lesssim T_{\rm PT}\lesssim$ GeV.

††preprint: FTPI-MINN-23-20††preprint: UTWI-39-2023

PTs have been studied extensively for decades in models of baryogenesis Cohen
_et al._ (1990); Fromme  _et al._ (2006); Hall  _et al._ (2020); Elor  _et
al._ (2022), (asymmetric) dark matter Cohen  _et al._ (2008); Zurek (2014);
Baker and Kopp (2017); Hall  _et al._ (2023); Asadi  _et al._ (2021a, b); Hall
_et al._ (2022); Elor  _et al._ (2023); Asadi  _et al._ (2022), extended Higgs
sectors Profumo  _et al._ (2007); Chacko  _et al._ (2006); Schwaller (2015);
Ivanov (2017), and spontaneously broken conformal symmetry Creminelli  _et
al._ (2002); Randall and Servant (2007); Nardini  _et al._ (2007); Konstandin
_et al._ (2010); Konstandin and Servant (2011); Baratella  _et al._ (2019);
Agashe  _et al._ (2020, 2021); Ares  _et al._ (2020); Levi  _et al._ (2023);
Mishra and Randall (2023), among others. PTs may also generate gravitational
waves (GW) Kosowsky  _et al._ (1992a, b); Kosowsky and Turner (1993);
Kamionkowski  _et al._ (1994); Caprini  _et al._ (2016, 2020); Caldwell  _et
al._ (2022); Auclair  _et al._ (2023) that can be observed in the near future.

_Late-time_ PTs have garnered attention due to their possible connections to
puzzling observations from PTAs Agazie  _et al._ (2023); Antoniadis  _et al._
(2023a) and the ‘$H_{0}$ tension’: a discrepancy between the direct
measurements of the Hubble constant $H_{0}$ Riess  _et al._ (2022) and its
value inferred from the CMB Aghanim  _et al._ (2020) (for a review see
Schöneberg  _et al._ (2022)). These PTs occur below SM temperatures of a GeV
(redshift $z\sim 10^{13}$) and before matter-radiation equality ($z\approx
3400$). For instance, a PT around $z\approx 10^{4}$ is motivated by the
proposed New Early Dark Energy (NEDE) solution to the $H_{0}$ problem
Wetterich (2004); Doran and Robbers (2006); Poulin  _et al._ (2019);
Niedermann and Sloth (2021, 2020). Since the Hubble tension favors such
“early-time” solutions Kamionkowski and Riess (2022), other ideas also use
such PTs Niedermann and Sloth (2022); Freese and Winkler (2021). A PT at
$z\sim 10^{10}$ has also been proposed as a source of the observed stochastic
GW background measured by PTAs Agazie  _et al._ (2023); Afzal  _et al._
(2023); Antoniadis  _et al._ (2023b); Reardon  _et al._ (2023); Xu  _et al._
(2023). Even later PTs may ameliorate the cosmological constant problem Bloch
_et al._ (2020).

Due to constraints from big bang nucleosynthesis (BBN) and the CMB, late-time
PTs that occur entirely in a dark sector with no significant reheating to SM
particles are favored Bai and Korwar (2022). Thus, we focus on PTs that only
release GWs and other forms of dark radiation (DR). The gravitational
backreaction on the SM sector is the only way to identify and constrain such
dark-sector PTs. For larger couplings between the dark sector and SM,
constraints stronger than our model-independent ones may exist. A well-known
constraint on post-BBN dark PTs is the bound on the number of additional
neutrinos, $\Delta N_{\rm eff}<0.29$ at $95\%$ CL Aghanim  _et al._ (2020);
Cielo  _et al._ (2023), derived from Baryon Acoustic Oscillations (BAO) and
CMB measurements. This places an upper bound on the fraction of DR energy
density compared to the total radiation energy density $f_{\rm
DR}\equiv\rho_{\rm DR}/\rho_{\text{tot}}\lesssim 0.04$.

A PT proceeds via nucleation of bubbles of true vacuum inside the metastable
phase. To estimate the typical number of bubbles, consider a comoving volume
corresponding to multipole $\ell\sim 10^{3}$. If $T_{\text{PT}}\sim\text{
TeV}$, there are an enormous number of bubbles inside that comoving volume:
$N_{b}\sim\left[a({\rm TeV})H({\rm
TeV})/(10^{3}a_{0}H_{0})\right]^{3}\sim\left[\tau_{0}/(10^{3}\tau({\rm
TeV}))\right]^{3}\sim 10^{34}$. Here $\tau({\rm TeV})$ and $\tau_{0}$ are the
conformal times at $T={\rm TeV}$ and today, respectively, with the
corresponding scale factors denoted by $a({\rm TeV})$ and $a_{0}$. However, if
$T_{\text{PT}}\sim\text{ keV}$, for example, the number of bubbles is much
less, $N_{b}\sim 10^{6}$.

In this Letter, we demonstrate that the finite number of bubbles involved in a
PT gives rise to super-horizon perturbations in PT completion time $\propto
1/\sqrt{N_{b}}$ which eventually contribute to curvature perturbations. We
present the first calculation of these perturbations using gauge-invariant
observables. The curvature power spectrum on large length scales follows a
universal power-law scaling in the PT parameters, independent of its
microscopic details.111While previous studies have addressed the effect of a
PT on curvature perturbations at a parametric level Freese and Winkler (2023,
2022); Liu  _et al._ (2023), including the _big bubble problem_ in
inflationary models Copeland  _et al._ (1994), they have not precisely
determined the magnitude of the power spectrum. Even for a dark PT with no
direct coupling to the SM, CMB and LSS measurements constrain the resulting
curvature perturbations. This in turn constrains $f_{\rm DR}$ even more than
the current and projected $\Delta N_{\rm eff}$ limits when
$T_{\text{PT}}\lesssim 1\leavevmode\nobreak\ \text{keV}$. Additionally, the
power-law scale dependence of this new contribution to curvature perturbations
can create distinct signatures in the CMB and matter power spectrum, enabling
us to identify the origin of the perturbation as a late-time PT.

#### Super-horizon fluctuations in percolation time from number of bubbles.

—A PT proceeds through bubble nucleation and expansion; a point in
space transitions to the true vacuum when a bubble engulfs it. However, this
process is inherently stochastic, and therefore, not all points in space
transition to the true vacuum at the same time. To quantify this, the bubble
nucleation rate per unit time and volume is Hogan (1983); Enqvist  _et al._
(1992); Hindmarsh  _et al._ (2015)

$\begin{split}\Gamma=\Gamma_{0}e^{-S(t)}\approx\Gamma_{0}e^{-S(t_{f})}e^{\beta(t-t_{f})},\end{split}$
(1)

with $S$ the bounce action for nucleating a critical bubble and
$\beta\approx-{\rm d}S(t_{f})/{\rm d}t$. Here $t_{f}$ is a reference time to
measure the progression of the PT. If the nucleation rate increases with time,
the fraction of the Universe in the metastable phase can be determined as
Enqvist  _et al._ (1992),

$\begin{split}h(t)=\exp\left(-e^{\beta(t-t_{f})}\right).\end{split}$ (2)

This implies that for $\beta\gg H_{\rm PT}$, the Hubble scale at $t_{f}$, the
PT completes very rapidly around the time $t_{f}$, within a small fraction of
Hubble time. However, as mentioned above, the PT does not complete exactly at
the same time everywhere. Given the factor of $\beta$ in Eqs. (1) and (2), we
expect the variance of the PT completion time to be inversely proportional to
$\beta$, i.e., a slower PT with a small value of $\beta$ will exhibit a larger
variance and vice versa.

To characterize the variation in the PT completion time, we denote the time at
which a point $\vec{x}$ transitions to true vacuum by $t_{c}(\vec{x})$, and
compute the two-point function $H_{\rm PT}^{2}\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$. $\delta
t_{c}(\vec{x})=t_{c}(\vec{x})-\bar{t}_{c}$ with $\bar{t}_{c}$ the average time
of conversion (see the Supplementary Material for a detailed definition). We
include factors of $H_{\rm PT}$ so that $H_{\rm PT}\delta t_{c}$ is
dimensionless. Practically, it is easier to write $H_{\rm PT}^{2}\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle=(H_{\rm
PT}/\beta)^{2}\times\beta^{2}\langle\delta t_{c}(\vec{x})\delta
t_{c}(\vec{y})\rangle$ and calculate the last factor, thus separating the
effect of cosmic expansion from the PT dynamics. It is more convenient to work
with the dimensionless Fourier transformed two-point function denoted by:

$\displaystyle\mathcal{P}_{\delta
t}(k)=\frac{k^{3}}{2\pi^{2}}\left(\frac{H_{\rm PT}}{\beta}\right)^{2}\int{\rm
d}^{3}r\leavevmode\nobreak\ e^{i\vec{k}\cdot\vec{r}}\beta^{2}\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle,$ (3)

where $\vec{r}=\vec{x}-\vec{y}$. Physically, this characterizes the
correlation of $t_{c}$ between any two points separated by a distance $\sim
1/k$. $\mathcal{P}_{\delta t}(k)$ has two distinct contributions: single-
bubble and double-bubble, corresponding to $t_{c}(\vec{x})$ and
$t_{c}(\vec{y})$ being set by the same bubble or two bubbles (see the
Supplementary Material).

$\mathcal{P}_{\delta t}(k)$ changes qualitatively for modes smaller or larger
than the typical bubble size $v_{w}/\beta$, where $v_{w}$ is the bubble wall
velocity. For modes smaller than or comparable to the bubble size, we
analytically compute $\mathcal{P}_{\delta t}(k)$ in the Supplementary Material
assuming constant wall velocity. However, due to intricate fluid dynamics and
magnetohydrodynamics effects, a translation between $\mathcal{P}_{\delta
t}(k)$ and sourced curvature perturbations is involved and model-dependent.
Therefore, the dependence of the curvature power spectrum on $k$ for
$k_{p}\equiv k/a_{\rm PT}\gtrsim\beta/v_{w}$ is less universal and varies as
the properties of the PT change. (Here $k_{p}$ is a physical wavenumber and
$a_{\rm PT}$ is the scale factor at $t_{f}$.) On the other hand, scales
$k_{p}\ll\beta/v_{w}$ have many bubbles contributing to the correlation
function within a spatial volume of linear size $1/k_{p}$. Thus,
$\mathcal{P}_{\delta t}$ is more universal and less sensitive to the details
of the PT thanks to the central limit theorem.

We can understand the behavior of $\mathcal{P}_{\delta t}$ for
$k_{p}\ll\beta/v_{w}$ as follows. The average separation between bubbles is
Enqvist  _et al._ (1992); Hindmarsh  _et al._ (2015),
$d_{b}\approx(8\pi)^{1/3}v_{w}/\beta$. In a given volume $V$, there are $N\sim
V/d_{b}^{3}$ independent regions where bubble nucleation can take place in an
uncorrelated fashion. As a result, the standard deviation in PT completion
time, when averaged over this entire volume, scales as $1/\sqrt{N}$. Thus,
given $N\sim 1/(k_{p}d_{b})^{3}$ for a scale $k_{p}$, we expect ${\cal
P}_{\delta t}\propto(k_{p}d_{b})^{3}$. Also, the combination $\beta\times t$
is what appears in the nucleation rates in Eqs. (1) and (2), so ${\cal
P}_{\delta t}\propto 1/\beta^{2}$. Thus, for $k_{p}\ll\beta/v_{w}$ and $H_{\rm
PT}\ll\beta$, the dimensionless power spectrum scales as

$\begin{split}{\cal P}_{\delta t}&\sim(8\pi)\left(\frac{H_{\rm
PT}}{\beta}\right)^{2}\left(\frac{v_{w}k_{p}}{\beta}\right)^{3}\\\
&=c(8\pi)v_{w}^{3}(k\tau_{\rm PT})^{3}\left(\frac{H_{\rm
PT}}{\beta}\right)^{5}.\end{split}$ (4)

We have noted that $\tau_{\rm PT}=1/(a_{\rm PT}H_{\rm PT})$ for PTs during
radiation domination and introduced a constant prefactor $c$. From the
qualitative arguments above we expect $c\sim 1$; a detailed calculation in the
Supplementary Material gives $c\approx 2.8$. Another way to understand the
$k^{3}$ scaling of ${\cal P}_{\delta t}$ for small $k$ is through Eq. (3). We
show in the Supplementary Material that the correlation $\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle\sim e^{-\beta ra_{\rm PT}/2}$ for
$\beta ra_{\rm PT}\gg 1$. As a result for $k_{p}\ll\beta/v_{w}$, the
exponential phase in Eq. (3) does not contribute, and the $k$-dependence of
${\cal P}_{\delta t}$ comes solely from the $k^{3}$ prefactor.

The result for ${\cal P}_{\delta t}$ is in Fig. 1, where $\xi\equiv
k_{p}d_{b}=(8\pi)^{1/3}v_{w}(k\tau_{\rm PT})(H_{\rm PT}/\beta)$. For $\xi\ll
1$, ${\cal P}_{\delta t}\sim\xi^{3}(H_{\rm PT}/\beta)^{2}$, as expected from
Eq. (4). However, close to $\xi\sim 1$, we see a deviation from that scaling.
Below, we only use the result for ${\cal P}_{\delta t}$ for $\xi\leq 1$ since
the regime $\xi\gg 1$ is sensitive to turbulence and magnetohydrodynamics
effects. However, these sub-horizon inhomogeneities also give rise to density
perturbations. In dark sectors in which sound waves dominantly source the GWs,
the resulting constraints may even be stronger Ramberg  _et al._ (2023).

We have described how the PT completion time fluctuates on superhorizon and
‘superbubble’ scales. While such fluctuations are of ‘isocurvature’ type
initially (since they do not induce a change in energy density), eventually
they source curvature perturbations.

Figure 1: Dimensionless power spectrum of phase transition time fluctuation,
re-scaled by $(\beta/H_{\rm PT})^{2}$ and plotted against comoving wavenumber
ratio $\xi$ representing the perturbation mode relative to typical bubble
separation. The PT spectrum (red) derived in the Supplementary Material is
independent of $(\beta,\tau_{\rm PT})$, unlike adiabatic perturbations
(purple, orange).

#### Curvature perturbations from fluctuations in percolation time.

—Consider two acausal patches, $A$ and $B$, where the PT takes place.
Since our analysis relies on $\xi\sim 0.1-1$ and $\beta/H_{\rm PT}\lesssim
10^{3}$, we are in a regime where ${\cal P}_{\delta t}\gg A_{s}$, where
$A_{s}=2.1\times 10^{-9}$ Akrami  _et al._ (2020) is the magnitude of the
inflationary scalar power spectrum (we express ${\cal P}_{\delta t}$ in terms
of gauge-invariant observables below). Thus, we can ignore effects due to
$A_{s}$ and assume the PT takes place in a Universe that is a priori
homogeneous in different patches. We will see how ${\cal P}_{\delta t}$ leads
to inhomogeneities in the dark sector and how they then feed back into the SM
sector, weighted by factors of $f_{\rm DR}$.

Take the two patches $A$ and $B$ to each have size $v_{w}/\beta$ and equal
energy density. $\rho_{F}$ is their (equal) false vacuum energy density and
$t_{c}^{A}$ and $t_{c}^{B}$ their respective PT completion times.
$t_{c}^{A}\neq t_{c}^{B}$ in general and we define the difference
$t_{c}^{B}-t_{c}^{A}\equiv\delta t_{c}\ll t_{c}^{A,B}\sim 1/H_{\rm PT}$, with
$\delta t_{c}>0$. When $A$ and $B$ undergo the PT, $\rho_{F}$ is converted
into DR with an energy density $\rho_{\rm DR}$. Right at $t_{c}^{A}$, the
energy densities in $A$ and $B$ are identical and the curvature perturbation
is still zero. However, there is a non-zero isocurvature perturbation in DR at
this time. This subsequently induces curvature perturbations as time evolves
since DR and vacuum energy redshift differently. In other words, the equation
of state of the Universe is not barotropic, i.e., the total pressure is not a
definite function of the total energy, $p\neq p(\rho)$. As a result, the
curvature perturbation is not constant (see e.g. Garcia-Bellido and Wands
(1996); Wands  _et al._ (2000)) and evolves with time after the PT occurs.

We assume that (i) the PT takes place when the dark-sector energy density is
dominated by the false vacuum and (ii) $\rho_{F}$ is entirely converted to
$\rho_{\rm DR}$ after the PT. We can then write the DR energy density in the
two patches at a later time $t_{\rm fin}$ as

$\begin{split}\rho_{\rm DR}^{A,B}(t_{\rm
fin})=\rho_{F}\left(\frac{t_{c}^{A,B}}{t_{\rm fin}}\right)^{2}\,.\end{split}$
(5)

This shows that the energy densities of DR in the two patches are different
and a nonzero DR density perturbation has been sourced by the DR isocurvature
perturbation.222The different values of $\rho_{\rm DR}$ in $A$ and $B$ changes
Hubble in the two patches, altering the energy-density redshift, but this
correction is $\mathcal{O}(\delta\rho_{\rm DR}/\rho_{\rm SM})$ and negligible
for our leading-order analysis. We can compute this density perturbation using
Eq. (5), $\delta\rho_{\rm DR}/\rho_{\rm DR}=2\delta t_{c}/t_{c}$. Since we are
working to leading order in perturbations, $\rho_{\text{DR}}\equiv\rho_{\rm
DR}^{A,B}(t_{\rm fin})$ and likewise for $t_{c}$ in the denominators of the
previous expression.

To compute the associated curvature perturbation, we can use the spatially
flat gauge (for a review, see Malik and Wands (2009)), which amounts to
comparing the energy densities in patches $A$ and $B$ at a common time $t_{\rm
fin}$ when the scale factors are identical. Then the curvature perturbation
(on uniform-density hypersurfaces) is

$\begin{split}\zeta=-\frac{H_{\rm PT}\delta\rho_{\rm DR}}{(\dot{\rho}_{\rm
DR}+\dot{\rho}_{\rm SM})}&=\frac{1}{4}f_{\rm DR}\frac{\delta\rho_{\rm
DR}}{\rho_{\rm DR}}=\frac{1}{2}f_{\rm DR}\frac{\delta
t_{c}}{t_{c}}\,.\end{split}$ (6)

So far, we have been assuming the dark-sector PT takes place when the dark
sector is vacuum-energy dominated, i.e., the latent heat released in the PT is
much larger than the temperature of the dark radiation bath. If not, we can
include the standard parameter $\alpha_{\rm PT}=\rho_{F}/\rho_{\rm
DR*}\lesssim 1$ where $\rho_{\rm DR*}$ is the energy density in the dark
radiation bath at the percolation time. Including this factor, we arrive at
the curvature power spectrum,

$\mathcal{P}_{\zeta}(k)=\frac{1}{4}\left(\frac{\alpha_{\rm PT}}{1+\alpha_{\rm
PT}}\right)^{2}f_{\rm DR}^{2}\mathcal{P}_{\delta t}(k)+\mathcal{P}_{\rm
ad}(k)\,.$ (7)

In the last term, we have reintroduced the uncorrelated adiabatic perturbation
$\mathcal{P}_{\rm ad}$. We take the pivot scale $k_{*}=0.05$ Mpc-1,
$A_{s}=2.1\times 10^{-9}$, and tilt $n_{s}=0.966$ Aghanim  _et al._ (2020)
when calculating $\mathcal{P}_{\rm ad}$. The constraints on
$\mathcal{P}_{\zeta}(k)$ can then be used to bound $f_{\rm DR}$ for different
dark PT parameters. For an alternate derivation that relates ${\cal P}_{\delta
t}$ to ${\cal P}_{\zeta}$ without relying on the separate Universe approach
followed here, together with a derivation using the $\delta N$-formalism
Langlois  _et al._ (2008), see the Supplementary Material.

Since we ignored the presence of inflationary, adiabatic perturbations while
analyzing $\mathcal{P}_{\delta t}$, Eq. (7) is valid only for
$\mathcal{P}_{\delta t}\gg{\cal P}_{\rm ad}$. In practice, given the current
precision $\Delta{\cal P}_{\zeta}/{\cal P}_{\zeta}\sim 5\%$ on CMB scales, the
above restriction puts an upper bound $f_{\rm DR}\lesssim 0.4$, above which
the effects of $A_{s}$ would be relevant for determining ${\cal P}_{\delta
t}$. Once the PT ends, all the dominant energy densities are in radiation, and
superhorizon $\zeta$-modes remain constant until they enter the horizon. We
note that PT also generates DR isocurvature, with a size roughly given by
${\cal P}_{\delta t}$, implying that isocurvature vanishes in the limit of
${\cal P}_{\delta t}\rightarrow 0$. The Planck constraint on DR isocurvature
Ghosh  _et al._ (2022) is similar to the constraints on curvature
perturbation. Therefore, we will not consider the effect of isocurvature
perturbations separately, but rather study their effects via the constraints
on ${\cal P}_{\zeta}$.

Figure 2: $2\sigma$ exclusion bounds on DR energy density fraction from
current observations, as derived in this work. We show bounds using the CMB
(Planck 2018 Akrami  _et al._ (2020)) and Lyman-$\alpha$ Bird  _et al._ (2011)
(pumpkin orange), as well as the FIRAS constraint Fixsen and Mather (2002)
(candy corn yellow). The light grey region represents the existing bound
$\Delta N_{\rm eff}\geq 0.29$ Aghanim  _et al._ (2020), and the dotted grey
lines, the projected bounds from the Simons Observatory (SO) $\Delta N_{\rm
eff}\geq 0.1$ Ade  _et al._ (2019) and CMB-S4 $\Delta N_{\rm eff}\geq 0.03$
Abazajian  _et al._ (2016). Future projections from SuperPIXIE Chluba  _et
al._ (2012a) (red), assuming sensitivity of
$\Delta\rho_{\gamma}/\rho_{\gamma}\sim 10^{-8}$, and PTA Lee  _et al._ (2021)
(maroon) are also depicted. We display the NEDE model’s preferred region
Niedermann and Sloth (2020) (darker grey) and the PT generating the potential
stochastic GW background Afzal  _et al._ (2023). To assess existing NEDE model
bounds, focus on $\beta/H_{\rm PT}$ values within the indicated bounds and
disregard the grey region representing the $\Delta N_{\rm eff}$ bound.

#### Time evolution of curvature perturbations.

—Perturbation modes with $k\tau_{\rm PT}\ll 1$ are outside the horizon
when the PT takes place and we can characterize their subsequent cosmological
evolution by just specifying ${\cal P}_{\zeta}(k)$. However, modes with
$\xi\approx 1$ correspond to $k\tau_{\rm PT}\approx(8\pi)^{-1/3}(\beta/H_{\rm
PT})v_{w}^{-1}\gtrsim 1$ for $\beta/H_{\rm PT}\gtrsim 10$, implying such modes
are already inside the horizon when the PT takes place. To derive constraints
based on ${\cal P}_{\zeta}(k)$ for a sub-horizon $k$-mode, we need to take
into account that there is no sub-horizon evolution for a time
$\Delta\tau\sim\tau_{\rm PT}-k^{-1}$, between mode reentry and the PT.

CMB temperature perturbations undergo diffusion damping while inside the
horizon. A delay in subhorizon evolution by an amount $\Delta\tau$ implies PT-
induced perturbations undergo less damping for a given $k$ compared to
$\Lambda$CDM expectations. Starting with the same value of ${\cal
P}_{\zeta}(k)$, the CMB anisotropies are larger in the PT scenario compared to
$\Lambda$CDM. In this work, we take a conservative approach by not including
this enhancement and leave a more precise computation for future work.

Perturbations in dark matter experience logarithmic growth in the radiation-
dominated era upon horizon reentry. In $\Lambda$CDM cosmology, the power in a
$k$-mode at the time of matter-radiation equality ($\tau_{\rm eq}\approx 110$
Mpc) is $\Delta_{\rm DM}(1/k,\tau_{\rm eq})^{2}{\cal P}_{\zeta}(k)$, where
$\Delta_{\rm DM}$ denotes the matter transfer function: $\Delta_{\rm
DM}(\tau_{i},\tau)\approx 6.4\ln(0.44\tau/\tau_{i})$ Hu and Sugiyama (1996);
Dodelson (2003). For the PT scenario, the analogous expression is $\Delta_{\rm
DM}(\tau_{\rm PT},\tau_{\rm eq})^{2}{\cal P}_{\zeta}(k)$. We use a rescaled
and weaker constraint ${\cal P}_{\zeta}(k)\times(\Delta_{\rm DM}(\tau_{\rm
PT},\tau_{\rm eq})/\Delta_{\rm DM}(1/k,\tau_{\rm eq}))^{2}$ for
$k>\tau^{-1}_{\rm PT}$ to take into account dark matter clustering bounds such
as from Lyman-$\alpha$ and future PTA constraint on DM clustering.

#### Cosmological Constraints.

—In Fig. 2, we present $2\sigma$ exclusion bounds on $f_{\text{DR}}$
using current constraints on $\mathcal{P}_{\zeta}$ and projected future
sensitivities. We translate the comoving time $\tau_{\rm PT}$ in Eq. (4) into
the SM temperature and redshift at the time of the PT. CMB Akrami  _et al._
(2020) and Lyman-$\alpha$ Bird  _et al._ (2011) measurements set upper bounds
on $\mathcal{P}_{\zeta}$ for $k$-modes up to $k\lesssim 3$ Mpc-1. Our analysis
excludes the pumpkin orange regions for various $\beta/H_{\rm PT}$ since in
those regions, the PT contribution to ${\cal P}_{\zeta}$ is too large. Other
constraints from ultracompact minihalos impacting PTAs may be relevant for
$T_{\text{PT}}\gtrsim 1\text{ MeV}$, but have unknown uncertainties related to
the time of DM collapse Clark  _et al._ (2016a, b).

The $T_{\rm PT}$ dependence of the constraints in Fig. 2 can be understood as
follows. During a radiation-dominated epoch, $k_{\rm peak}\propto T_{\rm PT}$,
where $k_{\rm peak}$ is the comoving wavenumber of the peak in ${\cal
P}_{\delta t}(k)$. The constraint on $f_{\rm DR}$ for a given $T_{\rm PT}$
then depends on whether $k_{\rm peak}$ or the IR tail of ${\cal P}_{\delta
t}(k)$ lies within the range probed by a given observable. Suppose that for a
range of $T_{\rm PT}$, the corresponding range of $k_{\rm peak}$ is directly
constrained by an observable. Then if the constraint on ${\cal P}_{\zeta}(k)$
over that range of $k_{\rm peak}$ is flat, the associated constraint on
$f_{\rm DR}$ is also flat with respect to $T_{\rm PT}$, resulting in the
plateaus in Fig. 2. This is because ${\cal P}_{\delta t}(k_{\rm peak})$ does
not change as $T_{\rm PT}$ varies. This is what happens for the CMB bound for
$T_{\rm PT}\lesssim 1$ eV for $\beta/H_{\rm PT}=10$. For larger $T_{\rm PT}$,
$k_{\rm peak}$ lies outside the region probed by ${\cal P}_{\zeta}(k)$
constraints; constraints are only sensitive to the tail of the ${\cal
P}_{\delta t}$ distribution and $\propto f_{\rm DR}^{2}k^{3}$ (from (4) and
(7)). In those regions, as $T_{\rm PT}$ is increased, the bound on $f_{\rm
DR}$ goes as $1/T_{\rm PT}^{3/2}$ (since $T_{\rm PT}\propto k$). A similar
transition from a plateau behavior is also seen at $\sim$100 eV for
Ly-$\alpha$ and at $\sim$$100$ MeV for PTA.

Notably, for $\beta/H_{\rm PT}\lesssim 400$, the bounds we derive from DR
inhomogeneities are stronger than current $\Delta N_{\rm eff}$ constraints
that track the homogeneous abundance of DR. For PTs that occur before BBN,
there is a stricter bound of $\Delta N_{\rm eff}\lesssim 0.23$ when applying
more observational constraints; for those PTs after, the constraint is
slightly weaker at $\Delta N_{\rm eff}\lesssim 0.31$ Yeh  _et al._ (2022).
Since our new CMB$+$Ly-$\alpha$ bounds are in the latter range, and we want to
show analogous $\Delta N_{\rm eff}$ constraints for CMB-S4, we plot Fig. 2
using the well-known $\Delta N_{\rm eff}<0.29$ Aghanim  _et al._ (2020). For
$\beta/H_{\rm PT}\sim 10$ and $T_{\rm PT}\lesssim$ keV, our analysis using
CMB+Ly-$\alpha$ constrains such PTs as much as or better than the future
Simons Observatory (SO) Ade  _et al._ (2019) and CMB-S4 projections on $\Delta
N_{\rm eff}$ Abazajian  _et al._ (2016).

The NEDE model in Niedermann and Sloth (2021, 2020) favors $\alpha_{\rm
PT}f_{\rm DR}$ in the upper (lower) dark gray region ($\pm 1\sigma$) from the
Planck18+BAO+LSS fit with (without) SH0ES data. While large values of
$\beta/H_{\rm PT}$ generally require extra model-building, the model in
Niedermann and Sloth (2020) assumes $\beta/H_{\rm PT}\gtrsim 100$ and permits
$\beta/H_{\rm PT}$ as large as $\sim 10^{3}$ by including a field to trigger
the PT. Still, our $\mathcal{P}_{\zeta}$ bound effectively disfavors the
preferred NEDE region in $\left(T_{\text{PT}},f_{\text{DR}}\right)$ for all
values of $\beta/H_{\rm PT}\lesssim$ 320 (230) with (without) SH0ES data.

For a large $\mathcal{P}_{\zeta}(k)$ with $k\lesssim 5400$ Mpc-1 the PT can
impact the dissipation of acoustic modes in photon-baryon perturbations,
altering the photon’s blackbody spectrum and inducing $\mu$\- and
$y$-distortions Chluba  _et al._ (2012b); Hooper  _et al._ (2023):

$\displaystyle X\\!\\!\simeq$ $\displaystyle\\!\\!A\displaystyle{\int_{k_{\rm
min}}^{\infty}}\frac{{\rm
d}k}{k}\mathcal{P}_{\zeta}(k)\left[Be^{-\frac{k}{5400/{\rm
Mpc}}}-Ce^{-(\frac{k}{31.6/{\rm Mpc}})^{2}}\right]$ (8)

where $k_{\rm min}=1$ Mpc-1, $(A,B,C)_{X}=(2.2,1,1)_{\mu}$ and
$(0.4,0,-1)_{y}$. Comparing this to the FIRAS bound of $|\mu|<9.0\times
10^{-5}$ and $|y|<1.5\times 10^{-5}$ Mather  _et al._ (1994); Fixsen  _et al._
(1996), we derive the exclusion bound labelled as ‘FIRAS’. When lowering
$T_{\rm PT}$, the $y$-distortion bound takes over the $\mu$-bound around
$T_{\rm PT}=10^{3}$ $(10^{2})$ eV for the $\beta/H_{\rm PT}=10$ $(100)$ case.
In contrast to Ref. Liu  _et al._ (2023), our findings indicate that the FIRAS
constraint is less stringent than the $\Delta N_{\rm eff}$ constraint, even
for PT with small $\beta=10H_{\rm PT}$.333Ref. Liu  _et al._ (2023) assumes
the spatial energy density perturbation directly equals $\alpha_{\rm PT}f_{\rm
DR}(\beta/H_{\rm PT})^{-5/2}$, omitting the numerical prefactor derived in the
Supplementary Material. Their procedure in deriving $\mathcal{P}_{\zeta}(k)$
is also sensitive to the choice of a window function, while ours is not.
Instead, our analysis derives $\mathcal{P}_{\zeta}(k)$ by starting with the
primordial fluctuations and following the energy-momentum conservation
equation for linear perturbations Malik and Wands (2009). This could explain
why our bounds differ.

Current $\mathcal{P}_{\zeta}$ measurements are less sensitive to PTs than the
$\Delta N_{\rm eff}$ constraint for $T_{\rm PT}\gtrsim\text{ keV}$, but
several proposed searches can constrain $\mathcal{P}_{\zeta}$ more powerfully
and constrain weaker dark PTs. Super-PIXIE aims to measure the CMB with a
sensitivity of $\Delta\rho_{\gamma}/\rho_{\gamma}\sim 10^{-8}$ Chluba  _et
al._ (2019) and the associated constraint is shown in red. PTAs can also probe
$\mathcal{P}_{\zeta}$ by observing the phase shift in periodic pulsar signals
mainly caused by the Doppler effect induced by an enhanced dark matter
structure that accelerates Earth or a pulsar. The PTA sensitivity curves
(maroon) use $\mathcal{P}_{\zeta}$ sensitivity derived in Lee  _et al._ (2021)
that assumes 20 years of observations of 200 pulsars. This future sensitivity
to PTs with $T_{\text{PT}}\gtrsim\text{ MeV}$ may test exotic DM models which
rely on them Elor  _et al._ (2023). Also shown is the $2\sigma$-preferred PT
region for the GW background hinted at by NANOGrav Agazie  _et al._ (2023);
Afzal  _et al._ (2023) (darker grey, see, e.g., Franciolini  _et al._ (2023)
for alternative GW spectrum assumptions). At face value, this region conflicts
with the $\Delta N_{\rm eff}$ constraint, but this prominent GW signal could
largely originate from supermassive black hole mergers. With enhanced PTA
measurements, we might still detect the PT signal within a comparable $T_{\rm
PT}$ range. Then PTA measurements of $\mathcal{P}_{\zeta}$ could complement
the GW detection.

#### Discussion.

—We have demonstrated that finite bubble statistics can lead to
superhorizon fluctuations in the PT completion time, regardless of the PT
details. These fluctuations source curvature perturbations that affect the
CMB, LSS, and other observables. Utilizing these, we find our constraints are
in tension with some of the best fit regions of the NEDE models proposed to
ameliorate the Hubble tension. At superhorizon scales, the (dimensionless)
power spectrum of these fluctuations has a $k^{3}$-model-independent scaling
since it is just determined by Poisson statistics. This contribution makes the
total curvature perturbation scale non-invariant. Thus, the associated CMB
phenomenology shares some similarities with the scale non-invariant effects
due to ‘primordial features’ Chluba  _et al._ (2015) produced during inflation
and models with ‘blue-tilted’ curvature perturbation Kasuya and Kawasaki
(2009); Kawasaki  _et al._ (2013); Chung and Tadepalli (2022); Ebadi  _et al._
(2023).

In our analysis, we have kept the $\Lambda$CDM parameters fixed. However,
given the model-independent shape, one can do a joint analysis where both
dark-sector and $\Lambda$CDM parameters are varied. We have also not
considered constraints from modes that are smaller than typical bubbles as
those are more model dependent. However, in the context of specific models one
can obtain stronger constraints from such modes. We leave these for future
work.

_Acknowledgments._ We thank Matthew Buckley, Zackaria Chacko, Peizhi Du, Anson
Hook, Gordan Krnjaic, Toby Opferkuch, Davide Racco, Albert Stebbins, Gustavo
Marques-Tavares, and Neal Weiner for helpful comments on the draft and
discussions. The research of GE is supported by the National Science
Foundation (NSF) Grant Number PHY-2210562, by a grant from University of Texas
at Austin, and by a grant from the Simons Foundation. The work of RJ is
supported by JSPS KAKENHI Grant Number 23K19048. SK is supported partially by
the NSF grant PHY-2210498 and the Simons Foundation. RM was supported in part
by DOE grant DE-SC0007859. YT is supported by the NSF grant PHY-2112540. RJ,
SK, RM, and YT thank the Mainz Institute for Theoretical Physics (MITP) of the
Cluster of Excellence PRISMA+ (project ID 39083149) for their hospitality
while a portion of this work was completed. GE, SK, and YT thank Aspen Center
for Physics (supported by NSF grant PHY-2210452) for their hospitality while
this work was in progress.

## References

  * Cohen  _et al._ (1990) A. G. Cohen, D. B. Kaplan, and A. E. Nelson, Phys. Lett. B 245, 561 (1990).
  * Fromme  _et al._ (2006) L. Fromme, S. J. Huber, and M. Seniuch, JHEP 11, 038 (2006), arXiv:hep-ph/0605242 .
  * Hall  _et al._ (2020) E. Hall, T. Konstandin, R. McGehee, H. Murayama, and G. Servant, JHEP 04, 042 (2020), arXiv:1910.08068 [hep-ph] .
  * Elor  _et al._ (2022) G. Elor _et al._ , in  _Snowmass 2021_ (2022) arXiv:2203.05010 [hep-ph] .
  * Cohen  _et al._ (2008) T. Cohen, D. E. Morrissey, and A. Pierce, Phys. Rev. D 78, 111701 (2008), arXiv:0808.3994 [hep-ph] .
  * Zurek (2014) K. M. Zurek, Phys. Rept. 537, 91 (2014), arXiv:1308.0338 [hep-ph] .
  * Baker and Kopp (2017) M. J. Baker and J. Kopp, Phys. Rev. Lett. 119, 061801 (2017), arXiv:1608.07578 [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] .
  * Asadi  _et al._ (2021a) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, Phys. Rev. Lett. 127, 211101 (2021a), arXiv:2103.09822 [hep-ph] .
  * Asadi  _et al._ (2021b) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, Phys. Rev. D 104, 095013 (2021b), arXiv:2103.09827 [hep-ph] .
  * Hall  _et al._ (2022) E. Hall, R. McGehee, H. Murayama, and B. Suter, Phys. Rev. D 106, 075008 (2022), arXiv:2107.03398 [hep-ph] .
  * Elor  _et al._ (2023) G. Elor, R. McGehee, and A. Pierce, Phys. Rev. Lett. 130, 031803 (2023), arXiv:2112.03920 [hep-ph] .
  * Asadi  _et al._ (2022) P. Asadi _et al._ , (2022), arXiv:2203.06680 [hep-ph] .
  * Profumo  _et al._ (2007) S. Profumo, M. J. Ramsey-Musolf, and G. Shaughnessy, JHEP 08, 010 (2007), arXiv:0705.2425 [hep-ph] .
  * Chacko  _et al._ (2006) Z. Chacko, H.-S. Goh, and R. Harnik, Phys. Rev. Lett. 96, 231802 (2006), arXiv:hep-ph/0506256 .
  * Schwaller (2015) P. Schwaller, Phys. Rev. Lett. 115, 181101 (2015), arXiv:1504.07263 [hep-ph] .
  * Ivanov (2017) I. P. Ivanov, Prog. Part. Nucl. Phys. 95, 160 (2017), arXiv:1702.03776 [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] .
  * Baratella  _et al._ (2019) P. Baratella, A. Pomarol, and F. Rompineve, JHEP 03, 100 (2019), arXiv:1812.06996 [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] .
  * 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] .
  * Levi  _et al._ (2023) N. Levi, T. Opferkuch, and D. Redigolo, JHEP 02, 125 (2023), arXiv:2212.08085 [hep-ph] .
  * Mishra and Randall (2023) R. K. Mishra and L. Randall, (2023), arXiv:2309.10090 [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] .
  * Auclair  _et al._ (2023) P. Auclair _et al._ (LISA Cosmology Working Group), Living Rev. Rel. 26, 5 (2023), arXiv:2204.05434 [astro-ph.CO] .
  * Agazie  _et al._ (2023) G. Agazie _et al._ (NANOGrav), Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023a) J. Antoniadis _et al._ , (2023a), arXiv:2306.16214 [astro-ph.HE] .
  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [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] .
  * 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] .
  * Wetterich (2004) C. Wetterich, Phys. Lett. B 594, 17 (2004), arXiv:astro-ph/0403289 .
  * Doran and Robbers (2006) M. Doran and G. Robbers, JCAP 06, 026 (2006), arXiv:astro-ph/0601544 .
  * 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] .
  * 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] .
  * Kamionkowski and Riess (2022) M. Kamionkowski and A. G. Riess, (2022), arXiv:2211.04492 [astro-ph.CO] .
  * Niedermann and Sloth (2022) F. Niedermann and M. S. Sloth, Phys. Rev. D 105, 063509 (2022), arXiv:2112.00770 [hep-ph] .
  * Freese and Winkler (2021) K. Freese and M. W. Winkler, Phys. Rev. D 104, 083533 (2021), arXiv:2102.13655 [astro-ph.CO] .
  * Afzal  _et al._ (2023) A. Afzal _et al._ (NANOGrav), Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023b) J. Antoniadis _et al._ (EPTA), (2023b), 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] .
  * Bloch  _et al._ (2020) I. M. Bloch, C. Csáki, M. Geller, and T. Volansky, JHEP 12, 191 (2020), arXiv:1912.08840 [hep-ph] .
  * Bai and Korwar (2022) Y. Bai and M. Korwar, Phys. Rev. D 105, 095015 (2022), arXiv:2109.14765 [hep-ph] .
  * Cielo  _et al._ (2023) M. Cielo, M. Escudero, G. Mangano, and O. Pisanti, (2023), arXiv:2306.05460 [hep-ph] .
  * Freese and Winkler (2023) K. Freese and M. W. Winkler, Phys. Rev. D 107, 083522 (2023), arXiv:2302.11579 [astro-ph.CO] .
  * Freese and Winkler (2022) K. Freese and M. W. Winkler, Phys. Rev. D 106, 103523 (2022), arXiv:2208.03330 [astro-ph.CO] .
  * 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] .
  * 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 .
  * Hogan (1983) C. J. Hogan, Physics Letters B 133, 172 (1983).
  * Enqvist  _et al._ (1992) K. Enqvist, J. Ignatius, K. Kajantie, and K. Rummukainen, Phys. Rev. D 45, 3415 (1992).
  * Hindmarsh  _et al._ (2015) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 92, 123009 (2015), arXiv:1504.03291 [astro-ph.CO] .
  * Ramberg  _et al._ (2023) N. Ramberg, W. Ratzinger, and P. Schwaller, JCAP 02, 039 (2023), arXiv:2209.14313 [hep-ph] .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Garcia-Bellido and Wands (1996) J. Garcia-Bellido and D. Wands, Phys. Rev. D 53, 5437 (1996), arXiv:astro-ph/9511029 .
  * 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 .
  * Malik and Wands (2009) K. A. Malik and D. Wands, Phys. Rept. 475, 1 (2009), arXiv:0809.4944 [astro-ph] .
  * Langlois  _et al._ (2008) D. Langlois, F. Vernizzi, and D. Wands, JCAP 12, 004 (2008), arXiv:0809.4646 [astro-ph] .
  * Ghosh  _et al._ (2022) S. Ghosh, S. Kumar, and Y. Tsai, JCAP 05, 014 (2022), arXiv:2107.09076 [astro-ph.CO] .
  * Bird  _et al._ (2011) S. Bird, H. V. Peiris, M. Viel, and L. Verde, Mon. Not. Roy. Astron. Soc. 413, 1717 (2011), arXiv:1010.1519 [astro-ph.CO] .
  * Fixsen and Mather (2002) D. J. Fixsen and J. C. Mather, Astrophys. J. 581, 817 (2002).
  * Ade  _et al._ (2019) P. Ade _et al._ (Simons Observatory), JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Chluba  _et al._ (2012a) J. Chluba, R. Khatri, and R. A. Sunyaev, Mon. Not. Roy. Astron. Soc. 425, 1129 (2012a), arXiv:1202.0057 [astro-ph.CO] .
  * Lee  _et al._ (2021) V. S. H. Lee, A. Mitridate, T. Trickle, and K. M. Zurek, JHEP 06, 028 (2021), arXiv:2012.09857 [astro-ph.CO] .
  * Hu and Sugiyama (1996) W. Hu and N. Sugiyama, Astrophys. J. 471, 542 (1996), arXiv:astro-ph/9510117 .
  * Dodelson (2003) S. Dodelson,  _Modern Cosmology_ (Academic Press, Elsevier Science, 2003).
  * Clark  _et al._ (2016a) H. A. Clark, G. F. Lewis, and P. Scott, Mon. Not. Roy. Astron. Soc. 456, 1394 (2016a), [Erratum: Mon.Not.Roy.Astron.Soc. 464, 2468 (2017)], arXiv:1509.02938 [astro-ph.CO] .
  * Clark  _et al._ (2016b) H. A. Clark, G. F. Lewis, and P. Scott, Mon. Not. Roy. Astron. Soc. 456, 1402 (2016b), [Erratum: Mon.Not.Roy.Astron.Soc. 464, 955–956 (2017)], arXiv:1509.02941 [astro-ph.CO] .
  * Yeh  _et al._ (2022) T.-H. Yeh, J. Shelton, K. A. Olive, and B. D. Fields, JCAP 10, 046 (2022), arXiv:2207.13133 [astro-ph.CO] .
  * Chluba  _et al._ (2012b) J. Chluba, A. L. Erickcek, and I. Ben-Dayan, Astrophys. J. 758, 76 (2012b), arXiv:1203.2681 [astro-ph.CO] .
  * Hooper  _et al._ (2023) D. Hooper, A. Ireland, G. Krnjaic, and A. Stebbins, (2023), arXiv:2308.00756 [astro-ph.CO] .
  * Mather  _et al._ (1994) J. C. Mather, E. S. Cheng, D. A. Cottingham, J. Eplee, R. E., D. J. Fixsen, T. Hewagama, R. B. Isaacman, K. A. Jensen, S. S. Meyer, P. D. Noerdlinger, S. M. Read, L. P. Rosen, R. A. Shafer, E. L. Wright, C. L. Bennett, N. W. Boggess, M. G. Hauser, T. Kelsall, J. Moseley, S. H., R. F. Silverberg, G. F. Smoot, R. Weiss, and D. T. Wilkinson, Astrophys. J. 420, 439 (1994).
  * Fixsen  _et al._ (1996) D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright, Astrophys. J. 473, 576 (1996), arXiv:astro-ph/9605054 .
  * Chluba  _et al._ (2019) J. Chluba _et al._ , Bull. Am. Astron. Soc. 51, 184 (2019), arXiv:1903.04218 [astro-ph.CO] .
  * Franciolini  _et al._ (2023) G. Franciolini, D. Racco, and F. Rompineve, (2023), arXiv:2306.17136 [astro-ph.CO] .
  * Chluba  _et al._ (2015) J. Chluba, J. Hamann, and S. P. Patil, Int. J. Mod. Phys. D 24, 1530023 (2015), arXiv:1505.01834 [astro-ph.CO] .
  * Kasuya and Kawasaki (2009) S. Kasuya and M. Kawasaki, Phys. Rev. D 80, 023516 (2009), arXiv:0904.3800 [astro-ph.CO] .
  * Kawasaki  _et al._ (2013) M. Kawasaki, N. Kitajima, and T. T. Yanagida, Phys. Rev. D 87, 063519 (2013), arXiv:1207.2550 [hep-ph] .
  * Chung and Tadepalli (2022) D. J. H. Chung and S. C. Tadepalli, Phys. Rev. D 105, 123511 (2022), arXiv:2110.02272 [astro-ph.CO] .
  * Ebadi  _et al._ (2023) R. Ebadi, S. Kumar, A. McCune, H. Tai, and L.-T. Wang, (2023), arXiv:2307.01248 [astro-ph.CO] .
  * Jinno and Takimoto (2019) R. Jinno and M. Takimoto, JCAP 01, 060 (2019), arXiv:1707.03111 [hep-ph] .
  * Jinno and Takimoto (2017) R. Jinno and M. Takimoto, Phys. Rev. D 95, 024009 (2017), arXiv:1605.01403 [astro-ph.CO] .

Supplementary Material for Finite Bubble Statistics Constrain Late
Cosmological Phase Transitions
Gilly Elor, Ryusuke Jinno, Soubhik Kumar, Robert McGehee, and Yuhsin Tsai

##  I Power spectrum of the transition time

In this Supplementary Material we derive analytic expressions for the power
spectrum of the transition time $t_{c}(\vec{x})$. The calculation is similar
to Refs. Jinno and Takimoto (2019, 2017). We subtract the average transition
time $\bar{t}_{c}$ and define $\delta t_{c}(\vec{x})\equiv
t_{c}(\vec{x})-\bar{t}_{c}$, and calculate the power spectrum of the
dimensionless quantity $\beta\delta t_{c}(\vec{x})$. It is handy to introduce
the dimensionless power spectrum as usual

$\displaystyle{\cal P}_{\delta t}$ $\displaystyle\equiv{\cal P}_{\beta\delta
t_{c}}(k)=\frac{k^{3}}{2\pi^{2}}P_{\beta\delta t_{c}}(k),$ (S1) $\displaystyle
P_{\beta\delta t_{c}}(k)$ $\displaystyle=\int d^{3}r\leavevmode\nobreak\
e^{i\vec{k}\cdot\vec{r}}\beta^{2}\langle\delta t_{c}(\vec{x})\delta
t_{c}(\vec{y})\rangle=\int_{0}^{\infty}4\pi r^{2}dr\leavevmode\nobreak\
\frac{\sin(kr)}{kr}\beta^{2}\langle\delta t_{c}\delta t_{c}\rangle(r).$ (S2)

Though $\langle\delta t_{c}\delta t_{c}\rangle(r)$ depends only on
$r\equiv|\vec{x}-\vec{y}|$, we sometimes use the notation $\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$.

Since we consider the two-point function of the transition time, we have two
distinct contributions in $\langle\delta t_{c}(\vec{x})\delta
t_{c}(\vec{y})\rangle$. They are illustrated in Fig. S1. We fix the evaluation
points $\vec{x}$ and $\vec{y}$, and consider different realizations of bubble
nucleation history. Since the transition time for each spatial point is
determined by the bubble whose wall passes that point for the first time, the
possibilities are that the transition times for the two points are determined
by one single bubble, or that they are determined by two different bubbles.
These two cases correspond to the left and right panels of Fig. S1.

Let us consider the probability for $t_{c}(\vec{x})$ and $t_{c}(\vec{y})$ to
be in the infinitesimal intervals $t_{x}<t_{c}(\vec{x})<t_{x}+dt_{x}$ and
$t_{y}<t_{c}(\vec{y})<t_{y}+dt_{y}$, respectively. For this to happen, note
that $\vec{x}$ and $\vec{y}$ must be in the false vacuum at $t_{x}$ and
$t_{y}$. Thus we need (1) and “(2-s) or (2-d)” in the following:

  * (1)

$\vec{x}$ and $\vec{y}$ remain in the false vacuum at $t_{x}$ and $t_{y}$,
respectively.

  * (2)

    * (2-s)

One bubble nucleates at the intersection of the two past cones (red diamond in
the right panel of Fig. S2).

    * (2-d)

Two bubbles nucleate on the surface of the past cones, one being in the
$\vec{x}$ side and the other in the $\vec{y}$ side (blue bands in the right
panel of Fig. S2).

Hereafter we assume luminal walls $v_{w}=1$ for simplicity, and thus the cones
are light-cones. Once we obtain the probability for these to happen, the
correlator $\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$ reduces
to $\sum_{t_{x}}\sum_{t_{y}}$ (probability) $\times$ (value of $\delta
t_{x}\delta t_{y}$). Since the two contributions are distinct, we may
decompose the correlator into two terms

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$
$\displaystyle=\langle\delta t_{c}(\vec{x})\delta
t_{c}(\vec{y})\rangle^{(s)}+\langle\delta t_{c}(\vec{x})\delta
t_{c}(\vec{y})\rangle^{(d)},$ (S3)

with each being called the single- and double-bubble contribution. Note that
the condition (1) is common to (2-s) and (2-d). In the following we use the
nucleation rate

$\displaystyle\Gamma(t_{n})$ $\displaystyle=\Gamma_{*}e^{\beta(t_{n}-t_{*})},$
(S4)

with $\beta$ being constant and $t_{*}$ being the typical transition time.
Hereafter we adopt $\beta=1$ unit. We sometimes use a shorthand notation
$x\equiv(t_{x},\vec{x})$ and $y\equiv(t_{y},\vec{y})$.

Now let us move on to the detailed calculation. We first calculate the average
transition time $\bar{t}_{c}$. For a spatial point $\vec{x}$ to experience the
transition between $[t_{x},t_{x},dt_{x}]$, it must be in the false vacuum at
$t_{x}$. Such a probability, which we call “survival probability”, is given by

$\displaystyle P_{\rm surv}(x)$ $\displaystyle=e^{-\int_{z\in({\rm
inside\leavevmode\nobreak\ the\leavevmode\nobreak\ past\leavevmode\nobreak\
cone\leavevmode\nobreak\ of\leavevmode\nobreak\ }x)}d^{4}z\leavevmode\nobreak\
\Gamma(t_{z})}$
$\displaystyle=\exp\left[-\int_{-\infty}^{t_{x}}dt_{n}\leavevmode\nobreak\
4\pi(t_{x}-t_{n})^{2}\Gamma_{*}e^{-(t_{n}-t_{*})}\right]$
$\displaystyle=\exp\left[-8\pi\Gamma_{*}e^{t_{x}-t_{*}}\right].$ (S5)

Hence, the average transition time is obtained as

$\displaystyle\bar{t}_{c}$
$\displaystyle=\int_{-\infty}^{\infty}dt_{x}\leavevmode\nobreak\
t_{x}\int_{-\infty}^{t_{x}}dt_{n}\leavevmode\nobreak\
4\pi(t_{x}-t_{n})^{2}\Gamma(t_{n})P_{\rm surv}(x)$
$\displaystyle=t_{*}-\ln(8\pi\Gamma_{*})-\gamma,$ (S6)

with $\gamma$ being the Euler-Mascheroni constant.

We next calculate the probability for (1) to happen. Following a similar
procedure as above, we get

$\displaystyle P_{\rm surv}(x,y)$
$\displaystyle=e^{-I(x,y)}=e^{-\int_{z\in({\rm inside\leavevmode\nobreak\
the\leavevmode\nobreak\ past\leavevmode\nobreak\ cones\leavevmode\nobreak\
of\leavevmode\nobreak\ }\vec{x}{\rm\leavevmode\nobreak\
and\leavevmode\nobreak\ }\vec{y})}d^{4}z\leavevmode\nobreak\
\Gamma(t_{z})}=e^{-\Gamma_{*}{\cal I}(x,y)e^{t_{\langle x,y\rangle}-t_{*}}},$
(S7) $\displaystyle{\cal I}(x,y)$
$\displaystyle=8\pi\left[e^{t_{x,y}/2}+e^{-t_{x,y}/2}+\frac{t_{x,y}^{2}-(r^{2}+4r)}{4r}e^{-r/2}\right].$
(S8)

Here we defined $t_{\langle x,y\rangle}\equiv(t_{x}+t_{y})/2$, $t_{x,y}\equiv
t_{x}-t_{y}$, $t_{\langle x,y\rangle,n}\equiv t_{\langle x,y\rangle}-t_{n}$,
$t_{\langle x,y\rangle,xn}\equiv t_{\langle x,y\rangle}-t_{xn}$, and
$t_{\langle x,y\rangle,yn}\equiv t_{\langle x,y\rangle}-t_{yn}$. In the
following we calculate the probability for “(1) and (2-s)” and that for “(1)
and (2-d)”, and estimate the correlator $\beta^{2}\langle\delta
t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$. The final expression is the sum
of Eqs. (S11) and (S16).

Figure S1: Two-dimensional illustration of the single-bubble (left) and
double-bubble (right) contributions. The black dots are the evaluation points
$\vec{x}$ and $\vec{y}$, while the stars denote the nucleation points of the
bubbles. The color corresponds to the transition time for each spatial point.
For the left case, the transition time at both points is determined by one
single bubble nucleating around the center, while for the right case it is
determined by two different bubbles nucleating at different locations.

Figure S2: Past-cone geometry in 3D spacetime. Embedded are 2D spatial slices
of past cones. Left: Past cone of $\vec{x}$. This geometry is used to
calculate the average transition time $\bar{t}_{c}$. The survival probability
$P(x)$ is defined as the probability for no bubble to nucleate inside the
cone. A bubble nucleating on the surface of the cone passes $\vec{x}$ at the
evaluation time $t_{x}$. Right: Past cones of $\vec{x}$ and $\vec{y}$. This
geometry is used to calculate the two-point correlation function
$\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$. For the single-
bubble case one single bubble nucleates at the intersection of the cone, while
for the double-bubble two different bubbles nucleate on the surface of the
each past cone.

Single-bubble contribution. The single-bubble contribution corresponds to the
cases in which one single bubble nucleates at the red diamond in the right
panel of Fig. S2. Noting that the red diamond forms a circle along the $\phi$
direction in three-dimensional space, we obtain

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle^{(s)}$
$\displaystyle=\int_{|t_{x}-t_{y}|<r}dt_{x}dt_{y}\leavevmode\nobreak\
(t_{x}-\bar{t}_{c})(t_{y}-\bar{t}_{c})\int_{-\infty}^{(t_{x}+t_{y}-r)/2}dt_{n}\int_{0}^{2\pi}d\phi\leavevmode\nobreak\
\frac{(t_{x}-t_{n})(t_{y}-t_{n})}{r}P_{\rm surv}(x,y)\Gamma(t_{n})$
$\displaystyle=\int_{-\infty}^{\infty}dt_{\langle
x,y\rangle}\int_{-r}^{r}dt_{x,y}\int_{r/2}^{\infty}dt_{\langle
x,y\rangle,n}\leavevmode\nobreak\ \frac{2\pi}{r}e^{-\Gamma_{*}{\cal
I}(x,y)e^{t_{\langle x,y\rangle}-t_{*}}}\Gamma_{*}e^{t_{\langle
x,y\rangle}-t_{*}}e^{-t_{\langle x,y\rangle,n}}$
$\displaystyle\qquad\times\left(t_{\langle
x,y\rangle}+\frac{t_{x,y}}{2}-\bar{t}_{c}\right)\left(t_{\langle
x,y\rangle}-\frac{t_{x,y}}{2}-\bar{t}_{c}\right)\left(t_{\langle
x,y\rangle,n}+\frac{t_{x,y}}{2}\right)\left(t_{\langle
x,y\rangle,n}-\frac{t_{x,y}}{2}\right).$ (S9)

We first integrate out $t_{\langle x,y\rangle}$ and get

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle^{(s)}$
$\displaystyle=\int_{-r}^{r}dt_{x,y}\int_{r/2}^{\infty}dt_{\langle
x,y\rangle,n}\leavevmode\nobreak\ \frac{2\pi e^{-t_{\langle
x,y\rangle,n}}}{r\,{\cal I}(x,y)}\left(t_{\langle
x,y\rangle,n}^{2}-\frac{t_{x,y}^{2}}{4}\right)\left[\left(\ln\left(\frac{{\cal
I}(x,y)}{8\pi}\right)\right)^{2}-\frac{t_{x,y}^{2}}{4}+\frac{\pi^{2}}{6}\right],$
(S10)

and then integrate out $t_{\langle x,y\rangle,n}$ and get

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle^{(s)}$
$\displaystyle=\int_{-r}^{r}dt_{x,y}\leavevmode\nobreak\ \frac{2\pi
e^{-r/2}}{r\,{\cal
I}(x,y)}\left(\frac{r^{2}}{4}+r+2-\frac{t_{x,y}^{2}}{4}\right)\left[\left(\ln\left(\frac{{\cal
I}(x,y)}{8\pi}\right)\right)^{2}-\frac{t_{x,y}^{2}}{4}+\frac{\pi^{2}}{6}\right].$
(S11)

Double-bubble contribution. The single-bubble contribution corresponds to the
cases in which two different bubbles nucleate on the surface of the past cones
in Fig. S2.

$\displaystyle\langle\Delta t_{c}(\vec{x})\Delta t_{c}(\vec{y})\rangle^{(d)}$
$\displaystyle=\int_{|t_{x}-t_{y}|<r}dt_{x}dt_{y}\leavevmode\nobreak\
(t_{x}-\bar{t}_{c})(t_{y}-\bar{t}_{c})\int_{-\infty}^{t_{x}}dt_{xn}\int_{-\infty}^{t_{y}}dt_{yn}\int
dc_{x}\int
dc_{y}\int_{0}^{2\pi}d\phi_{x}\int_{0}^{2\pi}d\phi_{y}\leavevmode\nobreak\ $
$\displaystyle\qquad(t_{x}-t_{xn})^{2}(t_{y}-t_{yn})^{2}P_{\rm
surv}(x,y)\Gamma(t_{xn})\Gamma(t_{yn})$
$\displaystyle=\int_{-\infty}^{\infty}dt_{\langle
x,y\rangle}\int_{-r}^{r}dt_{x,y}\int_{0}^{\infty}dt_{x,xn}\int_{0}^{\infty}dt_{y,yn}\leavevmode\nobreak\
16\pi^{2}e^{-\Gamma_{*}{\cal I}(x,y)e^{t_{\langle
x,y\rangle}-t_{*}}}\Gamma_{*}^{2}e^{2t_{\langle
x,y\rangle}-2t_{*}}e^{-(t_{x,xn}+t_{y,yn})}$ $\displaystyle\qquad\times
t_{x,xn}^{2}t_{y,yn}^{2}\left(t_{\langle
x,y\rangle}+\frac{t_{x,y}}{2}-\bar{t}_{c}\right)\left(t_{\langle
x,y\rangle}-\frac{t_{x,y}}{2}-\bar{t}_{c}\right)f_{x}f_{y}.$ (S12)

Note that the integration ranges for $c_{x}$ and $c_{y}$ can be nontrivial
since the blue bands in the right panel of Fig. S2 do not necessarily form
complete circles (i.e., complete spheres in three dimensions). The fractions
$f_{x}$ and $f_{y}$ in the last line take account of this

$\displaystyle f_{x}$
$\displaystyle=1\quad\left(t_{x,xn}<\frac{t_{x,y}+r}{2}\right),\quad\frac{(t_{x,xn}+r)^{2}-(t_{x,xn}-t_{x,y})^{2}}{4rt_{x,xn}}\quad\left(t_{x,xn}>\frac{t_{x,y}+r}{2}\right),$
(S13) $\displaystyle f_{y}$
$\displaystyle=1\quad\left(t_{y,yn}<\frac{-t_{x,y}+r}{2}\right),\quad\frac{(t_{y,yn}+r)^{2}-(t_{y,yn}+t_{x,y})^{2}}{4rt_{y,yn}}\quad\left(t_{y,yn}>\frac{-t_{x,y}+r}{2}\right).$
(S14)

We first integrate out $t_{\langle x,y\rangle}$ and get

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle^{(d)}$
$\displaystyle=\int_{-r}^{r}dt_{x,y}\int_{0}^{\infty}dt_{x,xn}\int_{0}^{\infty}dt_{y,yn}\leavevmode\nobreak\
\frac{16\pi^{2}e^{-(t_{x,xn}+t_{y,yn})}}{{\cal I}^{2}(x,y)}$
$\displaystyle\qquad\times
t_{x,xn}^{2}t_{y,yn}^{2}\left[\left(\ln\left(\frac{{\cal
I}(x,y)}{8\pi}\right)-1\right)^{2}-\frac{t_{x,y}^{2}}{4}+\frac{\pi^{2}}{6}-1\right]f(x)f(y),$
(S15)

and then integrate out $t_{x,xn}$ and $t_{y,yn}$ and get

$\displaystyle\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle^{(d)}$
$\displaystyle=\int_{-r}^{r}dt_{x,y}\leavevmode\nobreak\
\frac{16\pi^{2}}{{\cal I}^{2}(x,y)}$
$\displaystyle\qquad\times\left[4-\frac{e^{-t_{x,y}/2-r/2}}{2r}(r+t_{x,y}+4)(r-t_{x,y})-\frac{e^{t_{x,y}/2-r/2}}{2r}(r-t_{x,y}+4)(r+t_{x,y})\right.$
$\displaystyle\qquad\left.+\frac{e^{-r}}{16r^{2}}((r+4)^{2}-t_{x,y}^{2})(r^{2}-t_{x,y}^{2})\right]\left[\left(\ln\left(\frac{{\cal
I}(x,y)}{8\pi}\right)-1\right)^{2}-\frac{t_{x,y}^{2}}{4}+\frac{\pi^{2}}{6}-1\right].$
(S16)

Final expressions. The final expression for the correlator
$\beta^{2}\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$ is Eqs.
(S11) and (S16) substituted into Eq. (S3). These are one-dimensional integrals
and easy to evaluate. Fig. S3 shows the two contributions and their sum. As
expected, both decay exponentially at large $r$. Now it is straightforward to
evaluate $P_{\beta\delta t_{c}}$. Fig. S4 shows $P_{\beta\delta t_{c}}$
calculated from Eq. (S2). For small or large $k$, one may approximate the
spectrum with ${\cal P}_{\beta\delta t_{c}}\simeq 70(k/\beta)^{3}$ or ${\cal
P}_{\beta\delta t_{c}}\simeq 0.7(k/\beta)^{-3}$, respectively.

Figure S3: Single-bubble (blue) and double-bubble (red) contributions to the
correlator $\beta^{2}\langle\delta t_{c}(\vec{x})\delta t_{c}(\vec{y})\rangle$
and their sum (yellow). Figure S4: Power spectrum ${\cal P}_{\delta
t}\equiv{\cal P}_{\beta\delta t_{c}}$. The dashed lines are ${\cal
P}_{\beta\delta t_{c}}=100(k/\beta)^{3}$ and ${\cal P}_{\beta\delta
t_{c}}=(k/\beta)^{-3}$ for comparison.

##  II Alternative Derivation of Curvature Perturbation

Here we provide a derivation of (7) using the superhorizon evolution equation
for the curvature perturbation $\zeta$. On large scales, where spatial
gradients can be ignored, $\zeta$ follows the equation Garcia-Bellido and
Wands (1996); Wands  _et al._ (2000)

$\begin{split}\dot{\zeta}=-\frac{H}{\rho+p}\delta p_{\rm nad},\end{split}$
(S17)

where

$\displaystyle\delta p_{\rm nad}=\dot{p}\left(\frac{\delta
p}{\dot{p}}-\frac{\delta\rho}{\dot{\rho}}\right),$ (S18)

and $H$ the Hubble rate, $\rho$ the total energy density and $p$ the total
pressure. Now consider a dark sector PT that takes place at time $t_{c}$ that
instantaneously converts the false vacuum energy into DR. This assumption of
instantaneousness is justified since percolation takes place within a time
window much narrower than a Hubble time, for $\beta\gg H_{\rm PT}$, as we
consider.

We can then write the total energy density as,

$\displaystyle\rho=\rho_{F}\Theta(t_{c}-t)+\rho_{F}(a_{c}/a)^{4}\Theta(t-t_{c})+\rho_{\rm
SM},$ (S19)

and

$\displaystyle
p=-\rho_{F}\Theta(t_{c}-t)+(\rho_{F}/3)(a_{c}/a)^{4}\Theta(t-t_{c})+\rho_{\rm
SM}/3.$ (S20)

Here $\Theta(x)$ is the Heaviside theta function and $a_{c}$ is the scale
factor at $t_{c}$. From these we can compute the time derivatives:

$\displaystyle\dot{\rho}$
$\displaystyle=-\rho_{F}\delta(t_{c}-t)+(-4H)\rho_{F}(a_{c}/a)^{4}\Theta(t-t_{c})+\rho_{F}(a_{c}/a)^{4}\delta(t-t_{c})+\dot{\rho}_{\rm
SM},$ (S21)
$\displaystyle=(-4H)\rho_{F}(a_{c}/a)^{4}\Theta(t-t_{c})+\dot{\rho}_{\rm SM},$
(S22) $\displaystyle\equiv B+\dot{\rho}_{\rm SM},$ (S23)

and

$\displaystyle\dot{p}$
$\displaystyle=\rho_{F}\delta(t_{c}-t)+(-4H)(\rho_{F}/3)(a_{c}/a)^{4}\Theta(t-t_{c})+(\rho_{F}/3)(a_{c}/a)^{4}\delta(t-t_{c})+\dot{\rho}_{\rm
SM}/3,$ (S24)
$\displaystyle=(4/3)\rho_{F}\delta(t_{c}-t)+(-4H)(\rho_{F}/3)(a_{c}/a)^{4}\Theta(t-t_{c})+\dot{\rho}_{\rm
SM}/3,$ (S25) $\displaystyle=(4/3)\rho_{F}\delta(t_{c}-t)+B/3+\dot{\rho}_{\rm
SM}/3.$ (S26)

We note $B\ll\dot{\rho}_{\rm SM}$ since $f_{\rm DR}\ll 1$ for our analysis.
Using this we can approximate,

$\displaystyle\frac{\dot{p}}{\dot{\rho}}=\frac{(4/3)\rho_{F}\delta(t_{c}-t)+B/3+\dot{\rho}_{\rm
SM}/3}{B+\dot{\rho}_{\rm
SM}}\approx\frac{1}{3}+\frac{4\rho_{F}\delta(t-t_{c})}{3\dot{\rho}_{\rm SM}}.$
(S27)

So far we have computed the time derivatives of the homogeneous quantities
$\rho$ and $p$ with time. To evaluate $\delta p_{\rm nad}$, we also compute
the variations $\delta\rho$ and $\delta p$, which are variations in the total
energy density and pressure, respectively, as we compare different acausal
regions. As explained in the main text, we neglect inflationary fluctuations
and therefore $\delta\rho_{\rm SM}=0$. However, since the PT does not complete
everywhere at the same time, both $\delta t_{c}\neq 0$ and $\delta a_{c}\neq
0$. It is these variations that eventually source a non-zero curvature
perturbation $\zeta$. The variation in the total energy density is given by,

$\displaystyle\delta\rho$ $\displaystyle=\rho_{F}\delta(t_{c}-t)\delta
t_{c}+(4\delta
a_{c}/a_{c})\rho_{F}(a_{c}/a)^{4}\Theta(t-t_{c})+\rho_{F}(a_{c}/a)^{4}\delta(t-t_{c})(-\delta
t_{c})$ (S28) $\displaystyle=(4\delta
a_{c}/a_{c})\rho_{F}(a_{c}/a)^{4}\Theta(t-t_{c}),$ (S29)

and in the total pressure by,

$\displaystyle\delta p$ $\displaystyle=-\rho_{F}\delta(t_{c}-t)\delta
t_{c}+4(\delta
a_{c}/a)(\rho_{F}/3)(a_{c}/a)^{4}\Theta(t-t_{c})+(\rho_{F}/3)(a_{c}/a)^{4}\delta(t-t_{c})(-\delta
t_{c}),$ (S30) $\displaystyle=-(4/3)\rho_{F}\delta(t_{c}-t)\delta
t_{c}+4(\delta a_{c}/a)(\rho_{F}/3)(a_{c}/a)^{4}\Theta(t-t_{c}).$ (S31)

We now evaluate, to leading order in $f_{\rm DR}$,

$\displaystyle\delta
p-\frac{\dot{p}}{\dot{\rho}}\delta\rho\approx-\frac{4}{3}\rho_{F}\delta(t_{c}-t)\delta
t_{c}.$ (S32)

Hence we can finally evaluate the change in $\zeta$ using (S17),

$\displaystyle\dot{\zeta}$ $\displaystyle\approx-\frac{3H}{4\rho_{\rm
SM}}\left(-\frac{4}{3}\rho_{F}\delta(t_{c}-t)\delta t_{c}\right)$ (S33)
$\displaystyle\approx H\frac{\rho_{F}}{\rho_{\rm SM}}\delta(t_{c}-t)\delta
t_{c}$ (S34) $\displaystyle\approx\frac{\rho_{F}}{\rho_{\rm
SM}}\delta(t_{c}-t)\frac{\delta t_{c}}{2t_{c}}.$ (S35)

Since the curvature perturbation before the PT is zero (since we assume
inflationary fluctuations to be much smaller than the fluctuations induced by
the dark PT), we can integrate the above to arrive at

$\displaystyle\zeta=\frac{\rho_{F}}{\rho_{\rm SM}}\frac{\delta
t_{c}}{2t_{c}}=f_{\rm DR}\frac{\delta t_{c}}{2t_{c}}\leavevmode\nobreak\
\text{(following the PT)}.$ (S36)

This matches with the expression in (6), obtained using the separate Universe
approach.

For subhorizon modes, additional gradient terms would appear in Eq. (S17)
Malik and Wands (2009). However, those terms are homogeneous in $\zeta$ and
the Newtonian potential, and therefore do not source $\zeta$ via $t_{c}$
fluctuations. Hence, we do not consider the gradient terms when evaluating how
$\delta t_{c}$ sources $\zeta$. Once $\zeta$ is generated, however, we
estimate the subhorizon evolution as mentioned in the main text.

###  II.1 Derivation using the $\delta N$ formalism

Figure S5: An illustration of the energy density evolution of dark PT.
$\rho_{r,{\rm SM}}$ is the energy density of SM particles dominated by
radiation. False vacuum energy of the two patches $A$ and $B$ reheats into DR
at $t_{c}^{A,B}$. The total (DR and SM) energy density of the two patches are
equal at $t_{f}^{A,B}$.

We can also derive the form of $\zeta$ by comparing the energy densities in
the two patches $A$ and $B$, as discussed in the main text, and using the
$\delta N$ formalism. The evolution of the energy densities is shown
schematically in Fig. S5. The points $\alpha$ and $\beta$ denote the time of
PT for patch $A$ and $B$, respectively. At the point $\delta$, patch $A$ has a
smaller energy density compared to patch $B$, as DR redshifts from point
$\alpha$ to $\delta$. This energy difference can be thought to arise from a
hypothetical curvaton field. In particular, we can translate the energy
difference between $\beta$ and $\delta$, to an energy density difference at
$t_{c}^{A}$, between $\alpha$ and $\gamma$. Extrapolating to an even earlier
time, we can imagine this difference to arise from primordial fluctuations of
a curvaton field $\chi$, on a spatially flat hypersurface. Given Eq. (5), we
find $\delta\rho_{\chi}=(2\delta t_{c}/t_{c})\bar{\rho}_{\chi}$, where
$\bar{\rho}_{\chi}$ is the homogeneous energy density in $\chi$. The advantage
of this approach is that we can now directly utilize the standard expressions
derived in the curvaton scenario, see, e.g., Langlois  _et al._ (2008).

We denote the energy densities of the SM radiation and $\chi$ by $\rho_{r,{\rm
SM}}$ and $\rho_{\chi}$, respetively, and write them in terms of curvature
perturbation on uniform density hypersurfaces for each species Langlois  _et
al._ (2008):

$\rho_{r,{\rm SM}}=\bar{\rho}_{r,{\rm SM}}e^{4(\zeta_{\rm inf}-\delta
N)}\,,\quad\rho_{\chi}=\bar{\rho}_{\chi}e^{3(\zeta_{\chi}-\delta N)}\,.$ (S37)

Here $\zeta_{\rm inf}$ is the curvature perturbation of the SM radiation
(sourced by inflaton), and $\zeta_{\chi}$ is the curvature perturbation of
$\chi$ right before its decay (sourced by PT time fluctuation). On the uniform
density hypersurface, the total radiation energy density, $\rho_{r,{\rm
SM}}+\rho_{\rm DR}\equiv\rho_{r}=\bar{\rho}_{r}\exp[4(\zeta_{r}-\delta N)]$,
right after the curvaton decay equals $\bar{\rho}_{r}$, giving
$\zeta_{r}=\delta N$. Given that $\rho_{\rm DR}=\rho_{\chi}$ right after the
decay, we have

$f_{\rm DR}e^{3(\zeta_{\chi}-\zeta_{r})}+(1-f_{\rm DR})e^{4(\zeta_{\rm
inf}-\zeta_{r})}=1\,,\quad f_{\rm DR}=\frac{\bar{\rho}_{\rm
DR}}{\bar{\rho}_{r,{\rm SM}}+\bar{\rho}_{\rm DR}}\,.$ (S38)

The leading order expansion of the above equation in $\zeta_{\chi}$ and
$\zeta_{\rm inf}$ gives

$\zeta_{r}=\frac{3f_{\rm DR}}{4-f_{\rm DR}}\zeta_{\chi}+\frac{4(1-f_{\rm
DR})}{4-f_{\rm DR}}\zeta_{\rm inf}\,.$ (S39)

If $f_{\rm DR}\to 0$ ($f_{\rm DR}=1$), $\zeta_{r}=\zeta_{\rm inf}$
($\zeta_{r}=\zeta_{\chi}$) as expected. In the post-inflation era, where the
curvaton energy density is very subdominant, the uniform energy density
hypersurface is characterized by $\delta N=\zeta_{\rm inf}$, and the local
density of $\chi$ is

$\bar{\rho}_{\chi}e^{3(\zeta_{\chi}-\delta
N)}=\bar{\rho}_{\chi}e^{3(\zeta_{\chi}-\zeta_{\rm
inf})}=\bar{\rho}_{\chi}+\delta\rho_{\chi}\,.$ (S40)

The leading order expansion gives

$3\zeta_{\chi}=\frac{\delta\rho_{\chi}}{\bar{\rho}_{\chi}}+3\zeta_{\rm
inf}=2\frac{\delta t_{c}}{\bar{t}_{c}}+3\zeta_{\rm inf}\,,$ (S41)

and the total curvature perturbation

$\zeta_{r}=\frac{2f_{\rm DR}}{4-f_{\rm DR}}\frac{\delta
t_{c}}{\bar{t}_{c}}+\zeta_{\rm inf}\,.$ (S42)

Since the PT perturbation is uncorrelated to the inflaton perturbation, we
have for $f_{\rm DR}\ll 1$,

$\mathcal{P}_{\zeta_{r}}\approx\frac{1}{4}f_{\rm DR}^{2}\mathcal{P}_{\delta
t}+\mathcal{P}_{\rm ad}\,,$ (S43)

where we have denoted the power spectrum of $\zeta_{\rm inf}$ by
$\mathcal{P}_{\rm ad}$. This reproduces Eq. (7) for $\alpha_{\rm PT}\gtrsim
1$.


# Short proofs on $k$-extendible graphs

Shenwei Huang College of Computer Science, Nankai University, Tianjin 300350,
China. Partially supported by the National Natural Science Foundation of China
(11801284) and Natural Science Foundation of Tianjin (20JCYBJC01190). Yongtang
Shi Center for Combinatorics and LPMC, Nankai University, Tianjin 300071,
China. Partially supported by the National Natural Science Foundation of China
(No. 11922112) and Natural Science Foundation of Tianjin (Nos. 20JCZDJC00840,
20JCJQJC00090).

(October 7, 2021)

###### Abstract

In this note, we give short inductive proofs of two known results on
$k$-extendible graphs based on a property proved in [Qinglin Yu, A note on
$n$-extendable graphs. Journal of Graph Theory, 16:349-353, 1992].

##  1 Introduction

A graph $G$ is $k$-extendible if it satisfies the following conditions:

  * •

$|G|\geq 2k+2$;

  * •

$G$ is connected;

  * •

$G$ has a perfect matching;

  * •

for every matching $M_{k}$ of $G$ of size $k$, there is a perfect matching of
$G$ containing $M_{k}$.

The notion of $k$-extendible graphs was first defined and studied by Plummer
[5]. In particular, 2-extendible bipartite graphs play an important role in
the study of Pólya’s permanent problem [6] whose solution was obtained by
Robertsen, Seymour and Tomas [7] and independently by McCuaig [4]. We refer to
the monograph of Lovász and Plummer [3] for a detailed account of 1-extendible
graphs.

Our Contribution. Based on a property of $k$-extendible graphs proved by Yu
[8], we give short inductive proofs of two known results on $k$-extendible
graphs. Our proofs are much simpler than the existing proofs due to the fact
that the property allows us to apply the inductive hypothesis on subgraphs of
the given $k$-extendible graph rather than on the given graph itself.

##  2 Preliminaries

Let $G=(V,E)$ be a graph. For $S\subseteq V$ and a subgraph $H$ of $G$, the
neighborhood of $S$ in $H$, denoted by $N_{H}(S)$, is the set of vertices in
$H$ that are adjacent to some vertex in $S$. The size of $N_{H}(S)$ is denoted
by $d_{H}(S)$. If $S=\\{v\\}$, we simply write $N_{H}(v)$ and $d_{H}(v)$
instead of $N_{H}(\\{v\\})$ and $d_{H}(\\{v\\})$, respectively. For a subset
$S\subseteq V$, we denote by $G[S]$ the subgraph induced by $S$. Given a
matching $M$ of $G$, we denote by $V(M)$ the set of vertices that are
endvertices of edges in $M$. The minimum degree and the matching number of $G$
are denoted by $\delta(G)$ and $\alpha^{\prime}(G)$, respectively. For other
standard terminology we refer to [1].

We start with two simple propositions of $k$-extendible graphs that were
obtained by Plummer [5]. Since the proofs are short, we include them here for
the sake of completeness.

######  Proposition 1 ([5]).

Every $k$-extendible graph is $(k-1)$-extendible.

###### Proof.

Let $G$ be a $k$-extendible graph. By contradiction, let
$M=\\{a_{1}b_{1},\ldots,a_{k-1}b_{k-1}\\}$ be a matching of size $k-1$ that is
not contained in a perfect matching of $G$. Since $G$ is $k$-extendible, $M$
is a maximal matching of $G$. This implies that

$S=V(G)\setminus\\{a_{1},b_{1},\ldots,a_{k-1},b_{k-1}\\}$

is independent. Since $|G|\geq 2k+2$, it follows that $|S|\geq 4$. Since $M$
is not a maximum matching of $G$, it follows from Berge’s Theorem (see [1])
that there exists an $M$-augmenting path, that is, an $u$-$v$ path $P$ such
that $u,v\notin V(M)$ and edges of $P$ are alternating between $E(G)\setminus
M$ and $M$, starting with an edge not in $M$. Then $M^{\prime}=(E(P)\setminus
M)\cup(M\setminus E(P))$ is a matching of size $k$ with
$V(M^{\prime})=V(M)\cup\\{u,v\\}$. Since $G-V(M^{\prime})=S\setminus\\{u,v\\}$
is an independent set of size at least 2, $M^{\prime}$ cannot be extended to a
perfect matching of $G$. This is a contradiction. ∎

######  Proposition 2 ([5]).

Every 1-extendible graph is 2-connected.

###### Proof.

Suppose by contradiction that $G$ is a 1-extendible graph but not 2-connected.
Then there exists a cut vertex $v$ such that $G-v$ has components
$C_{1},\ldots,C_{t}$ for some $t\geq 2$. Since $G$ is connected, $v$ has a
neighbor $u_{i}\in C_{i}$ for each $1\leq i\leq t$. Since $G$ is 1-extendible,
there is a perfect matching containing $vu_{1}$ and this implies that
$|C_{1}|$ is odd. On the other hand, there is a perfect matching containing
$vu_{2}$ and this implies that $|C_{1}|$ is even. This is a contradiction. ∎

The following property of $k$-extendible graphs was proved by Yu [8] whose
proof used Theorem 1 below. Here we give a new proof that avoids the use of
Theorem 1.

######  Proposition 3 ([8]).

Let $G$ be a $k$-extendible graph with $k\geq 2$. Then for every edge $uv\in
E(G)$, $G-\\{u,v\\}$ is $(k-1)$-extendible.

###### Proof.

Let $e=uv\in E(G)$ and $G^{\prime}=G-\\{u,v\\}$. By Propostion 1, $G$ is
1-extendible and so 2-connected by Propostion 2. Since $G$ is $k$-extendible,
every matching of size $k-1$ of $G^{\prime}$ can be extended to a perfect
matching of $G^{\prime}$. So it remains to show that $G^{\prime}$ is
connected. Suppose by contradiction that $G^{\prime}$ has components
$C_{1},C_{2},\ldots,C_{t}$ for $t\geq 2$. Since $G$ is 2-connected, each of
$u$ and $v$ has a neighbor in each component $C_{i}$. Let $s\in C_{1}$ be a
neighbor of $u$ and $t\in C_{2}$ be a neighbor of $v$. Since $G$ is
$k$-extendible with $k\geq 2$, there is a perfect matching of $G$ containing
$\\{us,vt\\}$ by Propostion 1. This implies that $|C_{1}|$ is odd. Then there
is no perfect matching of $G$ containing $uv$. This contradicts that $G$ is
1-extendible. Therefore, $G-\\{u,v\\}$ is $(k-1)$-extendible. ∎

##  3 New Proofs

In this section, we present our new proofs of two known results on
$k$-extendible graphs. The first result was proved by Plummer [5] on the
connectivity of $k$-extendible graphs. The overall strategy of Plummer [5] was
to apply the inductive hypothesis on the input graph (due to Propostion 1) and
then used a variation of Menger’s Theorem. Our proof below, on the other hand,
is simpler due to the fact that we were able to apply the inductive hypothesis
on subgraphs of the input graph due to Propostion 3.

######  Theorem 1 ([5]).

Every $k$-extendible graph is $(k+1)$-connected.

###### Our Proof.

Let $G$ be a $k$-extendible graph. We prove by induction on $k$. The base case
is Propostion 2. Now suppose that $k\geq 2$ and the statement is true for
$(k-1)$-extendible graphs. By Propostion 3, $G-\\{u,v\\}$ is
$(k-1)$-extendible for every edge $uv\in E(G)$ and so is $k$-connected by the
inductive hypothesis. By Propostion 1 and Propostion 2, it follows that
$\delta(G)\geq 2$. For any vertex $v\in V(G)$, let $u$ be a neighbor of $v$.
Since $d(u)\geq 2$, $u$ has a neighbor $w$ other than $v$. Since
$H=G-\\{u,w\\}$ is $k$-connected, $d_{H}(v)\geq k$ and thus $d_{G}(v)\geq
k+1$. This shows that $\delta(G)\geq k+1$.

Now let $S\subseteq V(G)$ be an arbitrary set with $|S|=k$. Let $s\in S$ and
$t$ be a neighbor of $s$. We show that $G-S$ is connected.

Case 1. $t\in S$. Then

$G-S=(G-\\{s,t\\})-(S\setminus\\{s,t\\})$

is connected, since $G-\\{s,t\\}$ is $k$-connected.

Case 2. $t\notin S$. Let $G^{\prime}=G-(S\cup\\{t\\})$. Note that

$G^{\prime}=(G-\\{s,t\\})-(S\setminus\\{s\\}).$

Since $G-\\{s,t\\}$ is $k$-connected, $G^{\prime}$ is connected. Since
$\delta(G)\geq k+1$, $t$ has a neighbor in $G^{\prime}$. Therefore,

$G-S=G[V(G^{\prime})\cup\\{t\\}]$

is connected. ∎

The second result is on $k$-extendible bipartite graphs. The celebrated Hall’s
Theorem gives a necessary and sufficient condition for a balanced bipartite
graph to have a perfect matching. It turns out that $k$-extendible bipartite
graphs have a similar characterization.

######  Theorem 2 ([2]).

Let $G=(X,Y)$ be a connected bipartite graph with a perfect matching and
$|G|\geq 2k+2$. Then $G$ is $k$-extendible if and only if $|N(A)|\geq|A|+k$
for every subset $A\subseteq X$ with $1\leq|A|\leq|X|-k$.

Theorem 2 was first stated and proved by Brualdi and Perfect [2] in the
language of matrices (Theorem 2.1 in [2]). Here we give two graph-theoretical
proofs. The first one relies on Propostion 3 while the second one is based on
the König-Ore Formula.

######  Theorem 3 (The König-Ore Formula).

Let $G=(X,Y)$ be a bipartite graph. Then

$\alpha^{\prime}(G)=|X|-\max_{S\subseteq X}(|S|-|N(S)|).$

###### Our First Proof of Theorem 2.

We first prove the sufficiency. Take a matching

$M=\\{x_{1}y_{1},\ldots,x_{k}y_{k}\\}$

of size $k$. Let $X^{\prime}=X\setminus\\{x_{1},\ldots,x_{k}\\}$ and
$Y^{\prime}=Y\setminus\\{y_{1},\ldots,y_{k}\\}$. Denote by $H$ the subgraph of
$G$ induced by $X^{\prime}\cup Y^{\prime}$. Note that every nonempty subset
$A$ of $X^{\prime}$ has $1\leq|A|\leq|X|-k$. It follows from the assumption
that $|N_{G}(A)|\geq|A|+k$. This implies that $|N_{H}(A)|\geq|A|$. By the
König-Ore Formula, $H$ has a perfect matching $M_{H}$. It follows that
$M_{H}\cup M$ is a perfect matching of $G$ containing $M$. This shows that $G$
is $k$-extendible.

We now prove the necessity by induction on $k$.

Base Case: $k=1$. By contradiction, let $A$ be a subset of $X$ with
$1\leq|A|\leq|X|-1$ such that $|N(A)|<|A|+1$. Since $G$ has a perfect
matching, $|N(A)|\geq|A|$. It then follows that $|N(A)|=|A|$. Since $G$ is
connected, there is an edge $e=xy$ between $N(A)$ and $X\setminus A$. So there
is no perfect matching of $G$ containing $e$, simply because there are not
enough vertices in $G-\\{x,y\\}$ to match vertices in $A$.

Inductive Step: We assume that $k\geq 2$ and the statement is true for $k-1$.
let $A$ be an arbitrary subset of $X$ with $1\leq|A|\leq|X|-k$. If $N(A)=Y$,
then

$|N(A)|=|Y|=|X|\geq|A|+k.$

So we may assume that $N(A)\neq Y$. Since $G$ is connected, there is an edge
$xy\in E(G)$ such that $y\in N(A)$ and $x\in X\setminus A$. Let
$G^{\prime}=G-\\{x,y\\}=(X\setminus\\{x\\},Y\setminus\\{y\\})$. By Propostion
3, $G^{\prime}$ is $(k-1)$-extendible. On the other hand, $A\subseteq
X\setminus\\{x\\}$ has

$1\leq|A|\leq|X|-k=(|X|-1)-(k-1).$

By the inductive hypothesis, $|N_{G^{\prime}}(A)|\geq|A|+(k-1)$. Since
$N_{G}(A)=N_{G^{\prime}}(A)\cup\\{y\\}$, it follows that
$|N_{G}(A)|\geq|A|+k$. ∎

###### Our Second Proof of Theorem 2.

The difference lies in the inductive step of the necessity. We assume that
$k\geq 2$ and the statement is true for $k-1$. By contradiction, let $A$ be a
subset of $X$ with $1\leq|A|\leq|X|-k$ such that $|N(A)|<|A|+k$. By Propostion
1, $G$ is $(k-1)$-extendible. By the inductive hypothesis,
$|N(A)|\geq|A|+(k-1)$. It follows that

$|N(A)|=|A|+(k-1).$ (1)

Let $X^{\prime}=X\setminus A$ and write $B=N(A)$. Note that $|B|\geq k$ for
otherwise Equation 1 would be contradicted. Denote by $H$ the subgraph induced
by $B\cup X^{\prime}$. If $H$ has a matching of size $k$, then it cannot be
extended to a perfect matching of $G$ (because there are not enough vertices
to match vertices in $A$). So the matching number of $H$ is at most $k-1$. By
the König-Ore Formula,

$\alpha^{\prime}(H)=|B|-\max_{S\subseteq B}(|S|-|N_{H}(S)|)\leq k-1.$

So there exists a subset $S\subseteq B$ such that
$|S|-|N_{H}(S)|\geq|B|-(k-1)$. Since $|B|\geq k$, $|S|\geq 1$. Moreover,
$|N_{H}(S)|\leq|S|+(k-1)-|B|$. Therefore,

$|N_{G}(S)|\leq|A|+|N_{H}(S)|\leq|A|+|S|+(k-1)-|B|=|S|,$

where the last equality follows from Equation 1. Since $|A|\leq|X|-k$, it
follows that $|B|=|N(A)|=|A|+(k-1)\leq|X|-1$. Hence, $S\subseteq B$ violates
the condition for $G$ to be 1-extendible. ∎

##  4 Concluding Remarks

The fact that our proof of Propostion 3 does not use Theorem 1 makes our new
proof of Theorem 1 self-contained. To the best of our knowledge, our first
proof of Theorem 2 is new and self-contained. The second proof, in essence, is
the graph counterpart of the proof given in [2] stated in matrix language.
That proof used the Frobenius-König Theorem which is the matrix counterpart of
the the König-Ore Formula. However, we feel that it may be convenient for
graph theorists to have a graph-theoretical proof. So we include our second
proof here as well.

## References

  * [1] J. A. Bondy and U. S. R. Murty.  Graph Theory.  Springer, 2008.
  * [2] R. A. Brualdi and H. Perfect.  Extension of partial diagonals of matrices I.  Monatshefte für Mathematik, 75:385–379, 1971.
  * [3] L. Lovász and M. D. Plummer.  Matching Theory.  North–Holland, 1986.
  * [4] W. McCuaig.  Pólya’s permanent problem.  The Electronic Journal of Combinatorics, 11:R79, 2004.
  * [5] M. D. Plummer.  On $n$-extendible graphs.  Discrete Mathematics, 31:201–210, 1980.
  * [6] G. Pólya.  Aufgabe 424.  Arch. Math. Phys., 20:271, 1913.
  * [7] N. Robertsen, P. Seymour, and R. Thomas.  Permanents, Pfaffian orientations, and even directed circuits.  Annals of Mathematics, 150:929–975, 1999.
  * [8] Q. Yu.  A note on $n$-extendable graphs.  Journal of Graph Theory, 16:349–353, 1992.