aainstitutetext: Institute of Theoretical Physics, Faculty of Physics,
University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Polandbbinstitutetext:
CERN, Theoretical Physics Department, Geneva 23 CH-1211,
Switzerlandccinstitutetext: KTH Royal Institute of Technology, Department of
Physics, SE-10691 Stockholm, Swedenddinstitutetext: The Oskar Klein Centre for
Cosmoparticle Physics, AlbaNova University Centre, SE-10691 Stockholm, Sweden.

# Phases of Pseudo-Nambu-Goldstone Bosons

Fotis Koutroulis b Matthew McCullough a,c,d Marco Merchand a Stefan Pokorski a
and Kazuki Sakurai<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the vacuum dynamics of pseudo-Nambu-Goldstone bosons (pNGBs) for
$\text{SO}(N+1)\rightarrow\text{SO}(N)$ spontaneous and explicit symmetry
breaking. We determine the magnitude of explicit symmetry breaking consistent
with an EFT description of the effective potential at zero and finite
temperatures. We expose and clarify novel additional vacuum transitions that
can arise for generic pNGBs below the initial scale of
$\text{SO}(N+1)\rightarrow\text{SO}(N)$ spontaneous symmetry breaking, which
may have phenomenological relevance. In this respect, two phenomenological
scenarios are analyzed: thermal and supercooled dark sector pNGBs. In the
thermal scenario the vacuum transition is first-order but very weak. For a
supercooled dark sector we find that, depending on the sign of the explicit
symmetry breaking, one can have a symmetry-restoring vacuum transition
$\text{SO}(N-1)\rightarrow\text{SO}(N)$ which can be strongly first-order,
with a detectable stochastic gravitational wave background signal.

††preprint: CERN-TH-2023-172

##  1 Introduction

PNGBs Nambu:1960tm ; Goldstone:1961eq  arise in nature, as phonons, magnons,
pions and in a broad range of theoretical scenarios. It is no surprise that
they are abundant. It is a theorem that whenever a continuous global symmetry
is spontaneously broken that NGBs will arise Goldstone:1961eq . Furthermore,
it is widely believed that there can be no exact continuous global symmetries
in nature (more precisely, in gravitational theories Hawking:1975vcx ;
Zeldovich:1976vq ; Zeldovich:1977be ; Banks:1988yz ; Banks:2010zn ), in which
case any NGB will, in reality, be a pNGB. Thus, while the effective field
theory (EFT) description of the low-energy behaviour of exact NGBs is an
interesting object for theoretical study, it is likely that in nature the
physics below the scale of spontaneous symmetry breaking is dominated by the
scalar potential generated for pNGBs, since it contains the most relevant
operators.

Since the structure of the pNGB potential determines the vacuum dynamics it is
well-motivated to map the connections between explicit symmetry breaking
sources in a UV theory and the vacuum structure and dynamics in the IR, since
this aspect is physically relevant for pNGBs that are realised in nature. Once
this map is firmly established one can then determine and/or classify the
plausible phases of pNGB vacua and their dynamics.

Ref. Durieux:2021riy  established the first part of this programme for an
$\text{SO}(N+1)\to\text{SO}(N)$ spontaneous and explicit symmetry breaking
pattern. The fundamental building blocks of explicit symmetry breaking were
found to be the irrep spurions of $\text{SO}(N+1)$ which preserve an
$\text{SO}(N)$ subgroup. Each such spurion gives rise, in the IR, to a unique
Gegenbauer scalar potential which is an eigenfunction of the Laplacian on the
$N$-sphere. Any general pNGB potential for $\text{SO}(N+1)\to\text{SO}(N)$ can
thus be decomposed as a sum of Gegenbauer polynomials. Note that this is
strongly analogous to the solution of the Hydrogen wavefunction in quantum
mechanics. The angular momentum $|j,0\rangle$ eigenstates correspond to a non-
zero expectation value for the spin-$j$ irrep of $\text{SO}(3)$ which gives
rise to the $j^{\text{th}}$ Legendre polynomial, which is simply an
$\text{SO}(3)\to\text{SO}(2)$ Gegenbauer polynomial. Any wavefunction which is
a superposition of angular momentum eigenstates may be written as a sum of
Legendre polynomials. Thus what we are familiar with for angular momentum in
Hydrogen maps to the pNGBs of $\text{SO}(N+1)\to\text{SO}(N)$ breaking, where
the spatial rotation global symmetry becomes an internal global symmetry.

With this organisation of pNGB potentials complete the next logical step,
which is to understand the vacuum dynamics, is the focus of this work.
Throughout we are concerned with the same $\text{SO}(N+1)\to\text{SO}(N)$
spontaneous and explicit symmetry breaking pattern. We focus for the most
part, as a benchmark, on a single Gegenbauer pNGB potential, in the
understanding that the lessons learned will map, in a straightforward way,
into a sum of Gegenbauer potentials for any form of pNGB potential.

We begin by ascertaining the conditions under which the EFT description of the
potential is valid, both at zero and finite temperature (specifically in the
region of an interesting vacuum transition). This effectively places a
quantitative constraint on the magnitude of the explicit symmetry breaking
tolerable. Violation of this constraint implies a potential for which one does
not have a controlled series expansion in the explicit symmetry breaking,
whether at tree-level or at higher loop orders.

Subject to this constraint we then explore the vacuum dynamics for pNGBs,
which we find to be rich and varied. It should be noted that throughout there
is explicit $\text{SO}(N+1)\to\text{SO}(N)$ breaking thus, in terms of exact
global symmetries, there is no formal phase transition, since only
$\text{SO}(N)$ is an exact symmetry of the Lagrangian. However, since this
explicit symmetry breaking is small, one does have a sense in which the
fields, which play the role of order parameters, undergo vacuum transitions.

In this work we find that below the scale of spontaneous
$\text{SO}(N+1)\to\text{SO}(N)$ breaking, which is driven by the development
of a non-zero value for the $\text{SO}(N+1)$ radial mode, there are
generically additional pNGB vacuum transitions. There is an additional
critical temperature at which the pNGBs themselves develop a vacuum
expectation value, triggering a further stage of spontaneous
$\text{SO}(N)\to\text{SO}(N-1)$ breaking. This breaking is due to the explicit
symmetry breaking, but the change in order parameter is independent of the
magnitude of the explicit symmetry breaking. The reverse can also occur, with
a pattern of $\text{SO}(N+1)\to\text{SO}(N-1)$ breaking followed by a further
stage of $\text{SO}(N-1)\to\text{SO}(N)$ symmetry restoration at lower
temperatures.

It follows to determine the nature of these pNGB vacuum transitions. There are
two classes to consider, namely thermal and supercooled. In the thermal case
we find that the transition is generically weakly first-order. On the other
hand, when the pNGB sector is supercooled we find that the vacuum transition,
leading to symmetry restoration, can be strong enough to generate detectable
GW signatures. We finish with conclusions and future speculations.

##  2 pNGB Potential Regime of Validity

We consider an EFT containing the pNGBs $\psi$ arising from the spontaneous
breaking of an approximate global symmetry at the scale $f$. We define the
action at zero temperature as

$\mathcal{L}=\tfrac{1}{2}g_{ij}(\psi)\partial_{\mu}\psi^{i}\partial^{\mu}\psi^{j}+\mathcal{O}(\partial^{4})-\varepsilon
V_{\varepsilon}(\psi)-\varepsilon^{2}V_{\varepsilon^{2}}(\psi)-\mathcal{O}(\varepsilon^{3})+...+\mathcal{L}^{\text{CT}}~{}~{},$
(1)

where we have Taylor expanded in derivatives and in $\varepsilon$, which is,
by assumption for pNGBs, a small parameter associated with a source of
explicit symmetry breaking. $\mathcal{L}^{\text{CT}}$ represents the
counterterms required for renormalisation.

Before commencing with any concrete calculations some considerations are in
order concerning the validity of this EFT. To be effective, it must be valid
for some range of energies and field scales. For the former, scattering
amplitudes involving derivatives will scale as $(p^{2}/M^{2})^{j}$, where $j$
is some integer and $M$ is the cutoff energy of the EFT, often associated with
the mass of the radial mode of spontaneous symmetry breaking or some other UV
scale such as the mass scale of intermediate vector resonances. In any case,
the EFT description breaks down, by assumption, whenever $|p^{2}|\sim M^{2}$.

Equally important is the parameter $\varepsilon$. In order to be considered
pNGBs there must be some range of field values over which there is some
sensible notion of perturbative calculability within the EFT and of a scale
separation with the UV. For pNGBs the field range is periodic in the
spontaneous symmetry-breaking scale $\sim 2\pi f$. Due to this periodicity we
will require that the EFT description is valid and affords a degree of
perturbative calculability over all pNGB field values.

To determine the potential limits on the magnitude of $\varepsilon$ it is
helpful to consider the case of pions. Were the quark masses to be comparable
to the QCD scale, or the QED gauge coupling to be $e\sim 4\pi$ in the vicinity
of the QCD scale, there would be no sense in which one would have had light
pions at all, as they would naturally have mass at the QCD scale. Following
this, it is tempting to diagnose EFT validity using the pNGB masses. However,
mass-scale separation alone seems insufficient. For instance, in a scenario
with two large sources of explicit symmetry breaking
$\varepsilon_{1},\varepsilon_{2}\sim 1$ one could in principle fine-tune their
independent contributions to a pNGB potential to give a small mass-squared in
the global vacuum, generating a scale separation $m^{2}_{\psi}\ll M^{2}$.
However, one would have no control over perturbative corrections to the form
of the pNGB potential, either at tree-level at the matching scale or in the IR
at higher loops, due to the underlying magnitude of explicit symmetry
breaking. We must therefore be more pragmatic in determining the requirement
on $\varepsilon$ for the EFT description to be valid. The condition cannot
simply be that $m^{2}_{\psi}\ll M^{2}$, which is seemingly necessary but not
sufficient. Therefore we opt for the imprecise, but practical, condition that
the pNGB potential at $\mathcal{O}(\varepsilon)$ must be a good approximation
to the full potential with all quantum corrections included. In other words,
while $\mathcal{O}(\varepsilon^{2})$ and higher terms will exist, they must
not qualitatively alter the form of the pNGB potential.

The one-loop Coleman-Weinberg potential provides a useful diagnostic in this
respect. For pNGBs this is given by Binetruy:1984yx ; Alonso:2015fsp ;
Durieux:2021riy

$V^{\rm
CW}=\frac{1}{2}\mathrm{Tr}\int\frac{d^{4}p}{(2\pi)^{4}}\log\left[p^{2}+\varepsilon
g^{-1}\left(\frac{\delta^{2}V_{\varepsilon}}{\delta\psi^{2}}-\frac{\delta
V_{\varepsilon}}{\delta\psi}\,\Gamma\right)\right]~{}~{},$ (2)

where $\Gamma$ are the Christoffel symbols. The field-dependent curvature (or
mass-squared) entering this expression is

$\mathcal{M}_{\varepsilon}^{2}(\psi)=\varepsilon
g^{-1}\left(\frac{\delta^{2}V_{\varepsilon}}{\delta\psi^{2}}-\frac{\delta
V_{\varepsilon}}{\delta\psi}\,\Gamma\right)~{}~{},$ (3)

whose trace is simply the Laplace-Beltrami operator acting on the space
spanned by the pNGBs. Notably, this depends on the geometry of the manifold on
which the pNGBs live. In all of our applications we will be interested in the
scenarios in which the spontaneous symmetry breaking pattern is

$\frac{\text{SO}(N+1)}{\text{SO}(N)}\cong\mathcal{S}^{N}~{}~{},$ (4)

which we recall consists of the set of points a fixed distance from the origin
in $\mathcal{R}^{N+1}$. For the sake of illustration, we focus on scenarios in
which the explicit symmetry breaking follows the same pattern, preserving the
$\text{SO}(N)$ subgroup. As a result, we may parameterise the $N$ Goldstone
bosons on this manifold through the unit vector living in $\mathcal{R}^{N+1}$
as

$\bm{\phi}=f\sin\frac{\Pi}{f}\begin{pmatrix}\text{n}_{1}\\\ \text{n}_{2}\\\
\vdots\\\ \text{n}_{N}\\\ \cot\frac{\Pi}{f}\end{pmatrix}\;,$ (5)

where $\mathbf{n}\cdot\mathbf{n}=1$. Thus, in this picture, $\Pi/f$
essentially corresponds to the angle between the Goldstone boson direction and
a given arbitrarily chosen axes in $\mathcal{R}^{N+1}$.

In these coordinates we have that the relevant mass-squared matrix is

$\mathcal{M}_{\varepsilon}^{2}(\Pi)=\varepsilon\begin{pmatrix}\frac{\cot(\frac{\Pi}{f})}{f}V^{\prime}_{\varepsilon}\mathbb{1}_{N-1}&\mathbb{0}\\\
\mathbb{0}&V_{\varepsilon}^{\prime\prime}\end{pmatrix}\;.$ (6)

where

$V^{\prime}_{\varepsilon}\equiv\frac{\partial
V_{\varepsilon}}{\partial\Pi}\,\,\,~{}~{}~{}~{}{\rm
and}\,\,\,~{}~{}~{}~{}V_{\varepsilon}^{\prime\prime}\equiv\frac{\partial^{2}V_{\varepsilon}}{\partial\Pi^{2}}~{}~{}.$
(7)

Thus, considering the traces of products of this matrix which will arise in
perturbative calculations, it suffices to consider the Laplace-Beltrami
operator

$\displaystyle\Delta_{\mathcal{S}^{N}}V_{\varepsilon}=V_{\varepsilon}^{\prime\prime}+(N-1)\cot\frac{\Pi}{f}\frac{V^{\prime}_{\varepsilon}}{f}~{}~{}.$
(8)

As a result, truncating the momentum integral at the UV-cutoff, the zero-
temperature effective potential at one-loop is

$\displaystyle V$ $\displaystyle=$ $\displaystyle
V^{(0)}+V^{\text{CW}}+V^{\text{CT}}$ $\displaystyle=$
$\displaystyle\varepsilon\left[V_{\varepsilon}+\frac{M^{2}}{32\pi^{2}}\Delta_{\mathcal{S}^{N}}V_{\varepsilon}+V^{\text{CT}}_{\varepsilon}\right]+$
$\displaystyle\varepsilon^{2}\bigg{[}V_{\varepsilon^{2}}+\frac{1}{64\pi^{2}}\bigg{\\{}\left(V_{\varepsilon}^{\prime\prime}\right)^{2}\bigg{(}\log\left(\frac{\varepsilon}{M^{2}}V_{\varepsilon}^{\prime\prime}\right)-\frac{1}{2}\bigg{)}$
$\displaystyle+(N-1)\left(\frac{\cot\frac{\Pi}{f}}{f}V^{\prime}_{\varepsilon}\right)^{2}\bigg{(}\log\left(\frac{\varepsilon}{M^{2}}\frac{\cot\frac{\Pi}{f}}{f}V^{\prime}_{\varepsilon}\right)-\frac{1}{2}\bigg{)}\bigg{\\}}$
$\displaystyle+\frac{M^{2}}{32\pi^{2}}\Delta_{\mathcal{S}^{N}}V_{\varepsilon^{2}}+V^{\text{CT}}_{\varepsilon^{2}}\bigg{]}+\mathcal{O}(\varepsilon^{3})+...$

Here the terms denoted $V^{\text{CT}}$ represent the counterterms required to
renormalise the pNGB potential and $V^{(0)}$ is the tree-level scalar
potential. Thus we see that if $\Delta_{\mathcal{S}^{N}}V_{\varepsilon}$ has a
very different functional form to $V_{\varepsilon}$, the counterterm potential
cannot be similar in form to $V_{\varepsilon}$, implying some level of fine-
tuning between UV/threshold corrections, which must exist, and the bare
potential in order to realise the form of $V_{\varepsilon}$. If, however, they
are of a similar functional form then the $\mathcal{O}(\varepsilon)$
corrections will not destabilise the pNGB potential at that order. We will
return to this possibility in due course.

More immediately relevant is that the $\mathcal{O}(\varepsilon^{2})$ effective
potential corrections may significantly modify the qualitative nature of the
potential. This would signify the breakdown of the effective description of
the pNGB potential. Thus we will only work with EFTs for the pNGBs in which
$\varepsilon$ is sufficiently small that the physics of the zero-temperature
potential is well described at leading order in $\varepsilon$, hence

$V\approx\varepsilon\left(V_{\varepsilon}+\frac{M^{2}}{32\pi^{2}}\Delta_{\mathcal{S}^{N}}V_{\varepsilon}+V^{\text{CT}}_{\varepsilon}\right)~{}~{},$
(10)

is a reasonable approximation to the pNGB potential at zero temperature. This
can only be diagnosed on a case-by-case basis, and so we leave further
discussion of this aspect until a specific model has been chosen.

Now moving to finite temperature and following by analogy with the Coleman-
Weinberg potential, under the same set of assumptions, the full finite-
temperature potential at one-loop is, to a leading approximation,

$V(T)=V^{(0)}+V^{\text{CW}}+V^{\text{CT}}+V^{\text{T}}~{}~{},$ (11)

where Dolan:1973qd

$\displaystyle V^{\text{T}}$ $\displaystyle=$
$\displaystyle\frac{T^{4}}{2\pi^{2}}\mathrm{Tr}J_{B}\left(\frac{\mathcal{M}^{2}(\Pi)}{T^{2}}\right)~{}~{},$
(12) $\displaystyle=$
$\displaystyle\frac{T^{4}}{2\pi^{2}}\left(J_{B}\left(\frac{\varepsilon
V_{\varepsilon}^{\prime\prime}}{T^{2}}\right)+(N-1)J_{B}\left(\frac{\varepsilon\cot\left(\frac{\Pi}{f}\right)V_{\varepsilon}^{\prime}}{fT^{2}}\right)\right),$
(13)

and the function $J_{B}$ is

$J_{B}(x)=\int^{\infty}_{0}dyy^{2}\log\left(1-\exp^{-\sqrt{y^{2}+x}}\right)~{}~{}.$
(14)

Since we now have a new energy scale in the theory, $T$, we ought to
reconsider the conditions under which one has an appropriate description of
the physics. For $T\to 0$ we have that $V^{\text{T}}\to 0$, as expected, thus
at very low temperatures we may simply use the zero-temperature effective
potential already described.

At high temperatures we may also perform an expansion, in which case

$V^{\text{T}}\approx-N\frac{\pi^{2}}{90}T^{4}+\varepsilon\frac{T^{2}}{24}\Delta_{\mathcal{S}^{N}}V_{\varepsilon}-\frac{T}{12\pi}\left(\varepsilon
V_{\varepsilon}^{\prime\prime}\right)^{3/2}-(N-1)\frac{T}{12\pi}\left(\varepsilon\cot\frac{\Pi}{f}\frac{V^{\prime}_{\varepsilon}}{f}\right)^{3/2}+...~{}~{}.$
(15)

The validity of this expansion rests on two separate aspects. The first is
that the high-temperature expansion should be convergent, hence when the
system lies at high enough temperatures we require that the physics is, to a
good approximation, described by the second term alone, with the third
remaining a subleading correction. The second aspect concerns the non-
analyticity of the $J_{B}$ function, and hence of the third term of eq. (15).
This non-analyticity generates imaginary terms in the effective potential in
regions where $\partial^{2}V^{(0)}(\Pi)/\partial\Pi^{2}<0$. Since the
effective potential is, by definition, a real scalar quantity this signals a
breakdown in the effective description of the physics.

Without committing to a specific model in which one can calculate the
magnitude of the various terms this is as far as we may proceed, thus we now
commit to a specific class of scenarios.

##  3 Gegenbauer Goldstones

Experience with many physical systems, including electrostatics and
thermodynamics, suggests that when one encounters the Laplacian the natural
functions to work with are the eigenfunctions, satisfying an equation of the
form $\Delta_{\mathcal{S}^{N}}V_{\varepsilon}(\Pi)\propto
V_{\varepsilon}(\Pi)$. This is an eigenfunction problem and the solutions
which are analytic in $\Pi$ are the well-known Gegenbauer polynomials
Durieux:2021riy

$\Delta_{\mathcal{S}^{N}}G^{\frac{N-1}{2}}_{n}(\cos\Pi/f)=-\frac{n(n+N-1)}{f^{2}}G^{\frac{N-1}{2}}_{n}(\cos\Pi/f)~{}~{},$
(16)

where the eigenvalues and eigenfunctions are characterised by the two
integers, $N\geq 1$ and $n\geq 0$. In the application to the pNGB potential,
these integers are related to the explicit symmetry breaking pattern
$\text{SO}(N+1)\to\text{SO}(N)$ realised by a symmetry-breaking spurion in the
$n$-index symmetric irrep of $\text{SO}(N+1)$ Durieux:2021riy .

Motivated by this we will thus consider a zero-temperature pNGB potential of
the form

$\displaystyle V(\Pi,0)$ $\displaystyle\approx$
$\displaystyle\varepsilon_{n}V_{\varepsilon_{n}}+\mathcal{O}(\varepsilon^{2})$
(17) $\displaystyle\approx$
$\displaystyle\varepsilon_{n}f^{2}M^{2}G^{\frac{N-1}{2}}_{n}(\cos\Pi/f)+\mathcal{O}(\varepsilon^{2})+...~{}~{}.$

where note that from now on $\varepsilon$ would carry the subscript $n$ to
distinguish the above choice from the general pNGB case of eq. (10). No
summation over the index $n$ is implied. The typical shape of the Gegenbauer
potential at zero temperature ($T=0$) is shown in the left
($\varepsilon_{n}<0$) and right ($\varepsilon_{n}>0$) panels of Fig. 1. Note
that for positive $\varepsilon_{n}$ the global minimum is at a scale
$\langle\Pi\rangle\sim 5.1f/n$ Durieux:2021riy , whereas for negative
$\varepsilon_{n}$ the global minimum is at the origin. Importantly, this
potential is radiatively stable, since at leading order in this spurion only
this term can arise irrespective of the UV physics.

Figure 1: A cartoon picture showing the functional form of the Gegenbauer
thermal effective potential given by eq. (22), for the temperature asymptotics
$T=0$ and $T\gg 0$, when the symmetry-breaking parameter, $\varepsilon_{n}$,
is either positive or negative. The high-temperature limit terminates below
the radial mode mass $M$, otherwise the original, approximate, symmetry is
restored and the effective description of the model in terms of pNGBs is lost.
Cooling down will lead to $SO(N)$ symmetry restoration or breaking depending
on the sign of $\varepsilon_{n}$.

Since any general potential may be constructed from a linear sum of Gegenbauer
polynomials the lessons learnt from studying the single polynomial case will,
in generic cases, extend to more general pNGB potentials that can arise for
the $\text{SO}(N+1)\to\text{SO}(N)$ case.

###  3.1 pNGB Potentials at Zero Temperature

With this model we may now return to our general requirement of eq. (10). We
consider the zero-temperature potential at one-loop

$\displaystyle V^{(1)}(\Pi,0)$ $\displaystyle\approx$
$\displaystyle\varepsilon_{n}\left[\left(1-\frac{n(n+(N-1))M^{2}}{32\pi^{2}f^{2}}\right)V_{\varepsilon_{n}}+V_{\varepsilon_{n}}^{\rm
CT}\right]$
$\displaystyle+\varepsilon_{n}^{2}\left[V_{\varepsilon_{n}^{2}}+V_{\varepsilon_{n}^{2}}^{\rm
CT}-\frac{1}{128\pi^{2}}\left(\Bigl{(}V_{\varepsilon_{n}}^{\prime\prime}(\Pi)\Bigr{)}^{2}+(N-1)\frac{\cot^{2}\frac{\Pi}{f}}{f^{2}}\Bigl{(}V_{\varepsilon_{n}}^{\prime}(\Pi)\Bigr{)}^{2}\right)+...\right]~{}~{},$

where the ellipses denote the logarithmic terms. We see that at
$\mathcal{O}(\varepsilon_{n})$ the quadratic divergence may be absorbed into a
counterterm of the same functional form as the initial potential, reflecting
the radiative stability of this potential. However, we also see that,
regardless of the form of the potential at $\mathcal{O}(\varepsilon_{n}^{2})$,
there are calculable terms proportional to $\varepsilon_{n}^{2}$. In order for
the EFT to be valid it is necessary that these terms are subdominant to the
leading one.

Since it is the point at which the second derivative of the potential is
maximal in magnitude, to establish the maximal permitted value of
$\varepsilon_{n}$ we now focus our discussion around the origin of field
space. The Gegenbauer potential and its derivatives scale there as

$\displaystyle V_{\varepsilon_{n}}(0)$ $\displaystyle=$ $\displaystyle
f^{2}M^{2}\frac{(n+N-2)!}{n!(N-2)!}~{},$ $\displaystyle
V_{\varepsilon_{n}}^{\prime\prime}(\Pi)\Big{|}_{\Pi=0}$ $\displaystyle=$
$\displaystyle\cot\frac{\Pi}{f}\,\frac{V_{\varepsilon_{n}}^{\prime}(\Pi)}{f}\Big{|}_{\Pi=0}=-M^{2}(N-1)\frac{(n+N-1)!}{(n-1)!N!}~{}~{}.$
(19)

Thus we find the condition

$\displaystyle\frac{\varepsilon_{n}^{2}}{128\pi^{2}}N\Bigl{(}V_{\varepsilon_{n}}^{\prime\prime}(0)\Bigr{)}^{2}$
$\displaystyle\ll$ $\displaystyle\varepsilon_{n}V_{\varepsilon_{n}}(0)~{},$
(20)

which, under eq. (3.1), is reduced to

$|\varepsilon_{n}|\ll
128\pi^{2}\,\frac{f^{2}}{M^{2}}\frac{N!(n-1)!}{(n+N-1)!n(N-1)(n+N-1)}\equiv\varepsilon^{0}_{n,\rm
max}~{},$ (21)

as a necessary condition for the EFT expansion to be valid at zero
temperature, hence the upper-script $0$ in $\varepsilon_{n,\rm max}$ refers to
the zero-temperature case.

###  3.2 pNGB Potentials at Finite Temperature

After renormalization, for this class of potentials the high (enough)
temperature form is approximately

$V(\Pi,T)\approx\varepsilon_{n}f^{2}M^{2}\left(1-\frac{n(n+N-1)}{24}\frac{T^{2}}{f^{2}}\right)G^{\frac{N-1}{2}}_{n}(\cos\Pi/f)+\mathcal{O}(\varepsilon^{2})+...~{}~{}.$
(22)

Thus, for temperatures satisfying

$T^{2}\gtrsim T_{F}^{2}=\frac{24}{n(n+N-1)}f^{2}~{}~{},$ (23)

where we refer to $T_{F}$ as the “Flipping Temperature”, the overall sign of
the scalar potential has changed, indicating a transition in the position of
the global minimum relative to the zero-temperature potential, see Fig. 1. The
functional form of the scalar potential remains unchanged up to the overall
factor. We must, however, determine whether we may trust the EFT expansion at
this temperature by checking the magnitude of the next term in the finite-
temperature expansion.

We proceed as for the zero-temperature case, but now using the thermal
potential in eq. (15). The effective potential becomes

$\displaystyle V(\Pi,T)$ $\displaystyle\approx$
$\displaystyle-N\frac{\pi^{2}T^{4}}{90}+\varepsilon_{n}\Bigl{[}1-\frac{n(n+N-1)T^{2}}{24f^{2}}\Bigr{]}V_{\varepsilon_{n}}(\Pi)$
(24) $\displaystyle-$
$\displaystyle\frac{T\,(\varepsilon_{n}V^{\prime\prime}_{\varepsilon_{n}}(\Pi))^{\frac{3}{2}}}{12\pi}-(N-1)\frac{T\,(\cot{\frac{\Pi}{f}}\,\varepsilon_{n}V^{\prime}_{\varepsilon_{n}}(\Pi))^{\frac{3}{2}}}{12\pi
f^{3/2}}+{\cal O}({\varepsilon_{n}}^{2})~{}~{}.$

Focusing around the origin of the field space and noting that the second
derivative of the Gegenbauer polynomial is negative there, the relevant
constraint reads

$\Big{|}\frac{T^{2}}{24}\varepsilon_{n}\Delta_{\mathcal{S}^{N}}V_{\varepsilon_{n}}(0)\Big{|}\gg\Big{|}N\frac{T}{12\pi}\left(\varepsilon_{n}V_{\varepsilon_{n}}^{\prime\prime}(0)\right)^{3/2}\Big{|}~{}.$
(25)

This is a necessary condition for the validity of the EFT expansion at a given
temperature. For $T\approx T_{F}$ we get

$|\varepsilon_{n}|\ll
6\pi^{2}\frac{f^{2}}{M^{2}}\frac{N!(n-1)!}{(n+N-1)!n(N-1)(n+N-1)}\equiv\varepsilon^{T_{F}}_{n,\rm
max}~{}~{}.$ (26)

This is a stronger bound than at zero temperature, since

$\varepsilon^{T_{F}}_{n,\rm max}=\frac{3}{64}\varepsilon^{0}_{n,\rm
max}~{}~{}.$ (27)

The condition eq. (25) is necessary for validity at any temperature but not
sufficient. A stronger bound is obtained for $T=T_{\rm Crit}$, the ‘Critical
Temperature’, at which the vacuum transition is initiated. In general $T_{\rm
Crit}>T_{F}$, with the former defined as the temperature where the potential
energy of the two relevant phases becomes degenerate (or the two phases have
equal free energy density)

$V(0,T_{\rm Crit})=V(\braket{\Pi},T_{\rm Crit})~{},$ (28)

where $\braket{\Pi}$ is the pNGB value at the degenerate vacuum. From Fig. 1
note that no matter which cooling-down picture we consider, the potential
admits one global minimum around the field-space origin justifying our choice
of $V(0,T_{\rm Crit})$ as the free energy of one of the degenerate phases.

Using the effective potential of eq. (24), assuming for now
$\varepsilon_{n}>0$, the above equality gives

$T_{\rm Crit}^{2}+[B_{\varepsilon}T^{2}_{F}]\frac{T_{\rm
Crit}}{f}-T_{F}^{2}=0~{}~{},$ (29)

with the solution

$T_{\rm
Crit}=\frac{1}{2}\Biggl{[}-B_{\varepsilon}+\sqrt{\frac{4f^{2}}{T_{F}^{2}}+B^{2}_{\varepsilon}}\Biggr{]}\frac{T^{2}_{F}}{f}~{}~{}.$
(30)

$B_{\varepsilon}$ is a dimensionless parameter defined as

$\displaystyle B_{\varepsilon}$ $\displaystyle=$ $\displaystyle\frac{f\Delta
V_{\varepsilon,3/2}}{12\pi\Delta
V_{\varepsilon}}\approx\frac{f}{T_{F}}\Bigg{\\{}\frac{T_{F}N\big{(}\varepsilon_{n}V_{\varepsilon_{n}}^{\prime\prime}(0)\big{)}^{\frac{3}{2}}}{12\pi\,\varepsilon_{n}V_{\varepsilon_{n}}(0)}\Bigg{\\}}-\frac{f\big{(}\varepsilon_{n}V_{\varepsilon_{n}}^{\prime\prime}(\braket{\Pi})\big{)}^{\frac{3}{2}}}{12\pi\varepsilon_{n}V_{\varepsilon_{n}}(0)}$
(31)

where we have defined

$\Delta
V_{\varepsilon,3/2}=N\big{(}\varepsilon_{n}V_{\varepsilon_{n}}^{\prime\prime}(0)\big{)}^{\frac{3}{2}}-\big{(}\varepsilon_{n}V_{\varepsilon_{n}}^{\prime\prime}(\braket{\Pi})\big{)}^{\frac{3}{2}}~{}~{},$
(32)

and

$\Delta
V_{\varepsilon}=\varepsilon_{n}V_{\varepsilon_{n}}(0)\Bigg{(}1-\frac{V_{\varepsilon_{n}}(\braket{\Pi})}{V_{\varepsilon_{n}}(0)}\Bigg{)}\approx\varepsilon_{n}V_{\varepsilon_{n}}(0)>0~{}~{}.$
(33)

The notion of $T_{\rm Crit}$ and the validity of the EFT breaks down if
$B_{\varepsilon}$ has large imaginary part. Note that the term included in
$\\{\cdots\\}$ above, which is purely imaginary, has been used in eq. (25) to
derive the $\varepsilon$ bound of eq. (26). However, that bound is not
sufficient to make the left hand side of eq. (28) (and as a consequence
$B_{\varepsilon}$) to a good approximation real. It is found that only for an
$|\varepsilon_{n}|$ which is at least $\mathcal{O}(10^{-2})$ smaller than
$\varepsilon^{T_{F}}_{n,\rm max}$ the $\frac{f}{T_{F}}\\{\cdot\cdot\cdot\\}$
term can safely be neglected from $B_{\varepsilon}$ and the latter then
becomes

$B_{\varepsilon}\approx-\frac{f\big{(}\varepsilon_{n}V_{\varepsilon}^{\prime\prime}(\braket{\Pi})\big{)}^{\frac{3}{2}}}{12\pi\varepsilon_{n}V_{\varepsilon_{n}}(0)}<0~{}~{}~{}{\rm
and}~{}~{}~{}|B_{\varepsilon}|\ll 1~{}~{},$ (34)

and is real so we can safely evaluate the critical temperature. This stronger
bound is used in this paper as the sufficient condition for the validity of
the EFT in the whole relevant range of temperatures. Under that condition we
obtain that $T_{\rm Crit}\gtrsim T_{F}$ within a few percent. The two
temperatures are sometimes identified in our qualitative discussion but kept
distinct in the numerical calculations.

To summarise, we see that for this class of pNGB potentials there are
hierarchies of vacuum transitions. Starting from zero temperature as the
temperature is raised there will be a vacuum transition in the vicinity of the
flipping temperature. Depending on the sign of the spurion this will be from
zero pNGB vev to a non-vanishing one, with $\langle\Pi\rangle\propto f/n$, or
vice-versa. The nature of this transition is not yet clear from this analysis,
yet its existence is clear. Going to even higher temperatures, above the mass
scale of the radial mode in the UV completion the standard symmetry-restoring
transition occurs. These scenarios are illustrated in Fig. 2.

Figure 2:  Schematic phase diagram for radiatively and thermally stable pNGB
potentials, for $\varepsilon_{n}>0$ (left) and $\varepsilon_{n}<0$ (right).
Throughout there is explicit breaking $\text{SO}(N+1)\to\text{SO}(N)$. At high
temperatures, above the mass of the radial mode, an approximate
$\text{SO}(N+1)$ is restored. For $\varepsilon_{n}>0$ at lower temperatures,
$\text{SO}(N+1)$ is spontaneously broken and at some lower temperature the
exact $\text{SO}(N)$ is also spontaneously broken. Whereas for
$\varepsilon_{n}<0$ at lower temperatures, $\text{SO}(N+1)$ is spontaneously
broken to $\text{SO}(N-1)$ and at some lower temperature the exact
$\text{SO}(N)$ is restored.

It is surprising and rather non-trivial that for a single spontaneous symmetry
breaking scenario, with a single explicit symmetry-breaking spurion in a
symmetric irrep one has a hierarchy of vacuum transitions at hierarchical
scales. It remains to determine the nature of this new vacuum transition.

##  4 Cosmological Gegenbauer Phases

Having outlined the general phase structure of pNGB potentials it remains to
determine any potential observable consequences of the additional pNGB vacuum
transitions. We consider a dark sector (DS) containing pNGBs with two initial
conditions after the end of inflation; thermal and supercooled, however in
both cases colder than the visible sector. Given the natural origins and
ubiquity of light pNGBs in quantum field theories, and given the clear
evidence for the existence of dark matter, a DS scenario is well motivated and
plausible. In both cases we also investigate potential stochastic GW
Background signatures arising from the vacuum transitions.

###  4.1 Hot Dark Sector

We assume that the early universe dynamics is governed by the inflaton which,
at the end of inflation, starts to oscillate about the minimum of its
potential thus, due to its coupling to the Standard Model fields, the universe
enters the reheating period. At the same time we consider a DS of pNGBs which
is completely decoupled from (or may have an extremely small coupling to) the
SM, such that it will not thermalize with the SM fields. The DS temperature,
$T_{h}$, could be above or below the visible one, $T_{v}$, depending on how
strongly each sector couples to the inflaton. The ratio of temperatures after
reheating, $\xi_{\rm DS}=T_{h}/T_{v}$, is heavily constrained by Big Bang
Nucleosynthesis (BBN) and Cosmic Microwave Background (CMB) measurements
Planck:2018vyg ; Yeh:2022heq .

As noted, we assume $\xi_{\rm DS}<1$. This type of scenario has been
investigated in Breitbach:2018ddu ; Fairbairn . The case of $\xi_{\rm DS}>1$
is more delicate since it requires an out-of-equilibrium mechanism to inject
entropy back into the SM before BBN, see e.g. Ertas:2021xeh . For model-
independent studies regarding the constraints on DS vacuum transition
parameters see also PhysRevD.105.095015 ; Bringmann:2023opz .111Here we will
not deal with the case where $\xi_{\rm DS}=1$, which could happen either by
thermalization of the DS with the SM thermal bath or due to specific initial
conditions where the inflaton couples democratically to both sectors. We
escape the former by assuming the DS has a negligible interaction or never
comes into contact with SM and the latter by considering a different evolution
of the two sectors during reheating.

A general investigation of the nature of the transition is challenging and
essentially beyond the reach of standard computations. However, subject to the
requirement of small enough $\varepsilon_{n}$, discussed in the previous
section, we may have some control in the vicinity of the flipping temperature.

To proceed let us recall that the scalar potential in the DS is a Gegenbauer
polynomial. The vacuum structure of such a potential is non-trivial given that
different local minima coexist for a wide range of temperatures (see Fig. 1).
Analysing its thermal history in the following, a vacuum transition is
expected to occur. Particularly, for $\varepsilon_{n}>0$, $\Pi$ obtains a non-
zero vacuum expectation value and spontaneously breaks the SO(N) symmetry.

Before getting into a description of the phase transition details let us
present an analytic estimate for the transition strength $\alpha$, assuming it
takes place around $T\approx T_{F}$. To quantify $\alpha$ we use the latent
heat released normalized to the radiation energy density, which can be written
as

$\alpha(T)\equiv\frac{1}{\rho_{R}}\left(\Delta
V(\Pi,T)-\frac{T}{4}\Delta\frac{\partial V(\Pi,T)}{\partial T}\right)~{},$
(35)

where the difference between the false and true vacuum is taken. The energy
density is

$\displaystyle\rho_{R}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}\,g_{\Pi}^{*}T_{h}^{4}}{30}+\frac{\pi^{2}\,g_{\rm
SM}^{*}(T_{v})T_{v}^{4}}{30}$ (36) $\displaystyle=$
$\displaystyle\frac{\pi^{2}\,T_{h}^{4}}{30}\Big{(}N+\frac{g^{*}_{\rm
SM}(T_{v})}{\xi_{DS}^{4}}\Big{)}$

Since we consider a phase transition within the DS, the Hubble rate and the
other relevant parameters are functions of $T_{h}$. We keep $T_{v}$ as a fixed
initial parameter and the number of degrees of freedom in the DS corresponds
to the number of pNGBs, i.e., $g_{\Pi}^{*}=N$. We evaluate the radiation
degrees of freedom of the SM, $g_{\rm SM}^{*}$, from tabulated data in
Saikawa:2018rcs  and we keep them constant for temperatures in the vicinity of
the phase transition.

Making use of the high-temperature expansion we have that the potential energy
difference between false and true vacua is

$\displaystyle\Delta V(\Pi,T)$ $\displaystyle\approx$ $\displaystyle
V(0,T)-V(\braket{\Pi},T)=\left[1-\frac{T^{2}}{T_{F}^{2}}\right]\Delta
V_{\varepsilon}-\frac{T\,\Delta V_{\varepsilon,3/2}}{12\pi}~{},$ (37)

while the partial derivative with respect to temperature becomes

$\displaystyle\frac{T}{4}\Delta\frac{\partial V(\Pi,T)}{\partial T}$
$\displaystyle=$ $\displaystyle\frac{T}{4}\left[\frac{\partial
V(\Pi,T)}{\partial T}\Big{|}_{\Pi=0}-\frac{\partial V(\Pi,T)}{\partial
T}\Big{|}_{\Pi=\braket{\Pi}}\right]$ (38) $\displaystyle\approx$
$\displaystyle-\frac{1}{2}\frac{T^{2}}{T_{F}^{2}}\Delta
V_{\varepsilon}-\frac{1}{4}\frac{T\,\Delta V_{\varepsilon,3/2}}{12\pi}~{},$

thus

$\displaystyle\alpha(T)$ $\displaystyle\approx$ $\displaystyle\frac{\Delta
V_{\varepsilon}}{\rho_{R}}\Biggl{(}\left[1-\frac{1}{2}\frac{T^{2}}{T_{F}^{2}}\right]-\frac{3}{4}\frac{T\,\Delta
V_{\varepsilon,3/2}}{12\pi\Delta V_{\varepsilon}}\Biggr{)}~{}.$ (39)

Focusing around $T_{F}$, which is used as a proxy for the nucleation
temperature $T_{n}$ since we have verified they are very close numerically,
the second term in the above equation reduces to the $\\{\cdot\cdot\cdot\\}$
term of eq. (31) which, as follows from the discussion above eq. (34), has to
be very small for the validity of the EFT. Thus, the transition strength
becomes

$\alpha(T_{F})\approx\frac{\Delta
V_{\varepsilon}}{\rho_{R}}\left[1-\frac{1}{2}\frac{T_{F}^{2}}{T_{F}^{2}}\right]=\frac{\Delta
V_{\varepsilon}}{2\rho_{R}}~{}~{}.$ (40)

By setting $\varepsilon_{n}=10^{-2}\varepsilon^{T_{F}}_{n,\rm max}$ we obtain

$\alpha(T_{F})\lesssim\frac{0.002}{\left(1+\frac{g_{\rm
SM}^{*}(T_{v})}{\xi_{\rm DS}^{4}N}\right)}~{}~{}.$ (41)

The phase transition is weak because of the strong upper bound on
$\varepsilon_{n}$, which also controls the magnitude of the explicit breaking
of the original symmetry. This value of $\alpha$ corresponds to, at most, a
very weakly first-order transition and suppressed gravitatonal wave spectrum.

Since the phase transition occurs at finite temperature under the presence of
a non-negligible thermal plasma formed out of a system of pNGBs, the expanding
bubble walls transmit a substantial energy density and pressure to the
surrounding plasma. Hence, the dominant source of GW production is the motion
of the plasma itself, expressed in the form of sound waves. As described in
greater detail in an app. (A), for the GW spectrum, under the assumption of
small $\alpha$, the peak of the spectrum is Hindmarsh:2013xza ;
Hindmarsh:2015qta ; Hindmarsh:2017gnf ; Hindmarsh:2020hop

$\Omega_{\rm sw}(\text{Peak})h^{2}\approx 4\times
10^{-7}\,\left(R_{*}H_{*}\right)^{2}\left(\kappa_{\rm
sw}\,\alpha\right)^{\frac{3}{2}}\,,$ (42)

where $\kappa_{\rm sw}$ encodes kinetic energy normalized to vacuum energy. We
evaluate the efficiency factor $\kappa_{\rm sw}$ using the numerical fits of
Espinosa:2010hh . $R_{*}$ is the average bubble size at collision. As
described in app. (A), we find that numerically, at the time of the
transition, one has $R_{*}H_{*}\sim 10^{-6}$. Hence we expect at most to have
a spectral peak of magnitude

$\Omega_{\rm sw}(\text{Peak})h^{2}\lesssim 4\times 10^{-23}\,,$ (43)

well below the expected reach of future gravitational wave detectors.

###  4.2 Supercooled Dark Sector

Let us now explore the extreme possibility that our pNGB DS is supercoooled,
parameterized as $\xi_{\rm DS}\approx 0$. This may occur if, for instance, the
DS is very weakly coupled to the inflaton. We also discuss the role of
$\varepsilon_{n}$’s sign. In the previous section we have assumed that
$\varepsilon_{n}>0$. However, in principle, $\varepsilon_{n}$ can be either
positive or negative and, as we explain in the following, the choice of sign
impacts the cosmology of the DS.

It is possible that the expansion rate of the universe is initially much
faster than the bubble nucleation rate in a supercooled DS.222Since the DS is
almost decoupled from the SM it will evolve independently, so we consider that
the visible sector is “frozen” to a given temperature $T_{v}$. As a
consequence the DS can enter a period of supercooling, remaining in a local
minimum until quantum tunneling towards another local or a global minimum
takes place.

For $\varepsilon_{n}\gtrsim 0$ the vacuum dynamics of a supercooled DS is
governed by the zero-temperature potential of eq. (17).

Figure 3: Successive tunneling towards the true vacuum for the benchmark
scenario $n=15,N=4$. The colouring shows that we move from a higher
$\braket{\Pi}$ ( purple dot) down to smaller values until the DS reaches the
deepest minimum ( red dot).

Such a scenario has interesting phenomenology as the associated potential
possesses various local minima and as a consequence the supercooled DS could
in principle exhibit successive vacuum transitions, depicted on Fig. 3, via
tunneling. For an indicative example we consider the case when the DS is
initially in the minimum depicted by the  purple dot in Fig. 3 with associated
vev $\braket{\Pi_{\rm purple}}$. We calculate the probability of tunneling
towards its nearest neighbor  blue dot with associated vev $\braket{\Pi_{\rm
blue}}$. For this transition it is clear that the barrier between the vacua is
large compared to the energy difference between them, therefore the thin wall
approximation Coleman4  is a well motivated analytic approach. According to
this approximation and following Fairbairn , the probability of nucleating a
critical bubble via quantum tunneling is

$\Gamma_{4}=A_{4}~{}e^{-S_{4}}\equiv\frac{1}{R_{0}^{4}}\left(\frac{S_{4}}{2\pi}\right)^{2}e^{-S_{4}}$
(44)

where $S_{4}$ is the $O(4)$-symmetric bounce solution and $R_{0}$ is the size
of the nucleating bubble.

Moreover, following the cosine-like approximation to the Gegenbauer potential
provided in Eq. (2.12) of Durieux:2021riy , and employing the triangle
approximation to the cosine potential, for which an anlytic expression was
derived in Duncan:1992ai , in the thin wall approximation the bounce action
$S_{4}$ scales as

$\displaystyle S_{4}\approx\frac{32\pi^{2}}{3}\frac{(\Delta
V_{\text{Max}}(\Pi))^{2}(\Delta\Pi)^{4}}{(\Delta V(\Pi))^{3}}~{}~{},$ (45)

where $\Delta\Pi$ is the leading order change in vev between vacua, $\Delta
V(\Pi)$ is the change in vacuum energy between the two vacua and $\Delta
V_{\text{Max}}(\Pi)$ is the change in vacuum energy between the vacuum and the
top of the barrier between them.

The resulting expression for the bounce, in the large $n$ limit, is

$S_{4}\sim\frac{2^{3-n-N}n^{2}\pi^{5}\Gamma(n+N)}{3(N-1)^{4}\Gamma\left(\frac{n+1}{2}\right)\Gamma\left(\frac{N}{2}\right)\Gamma\left(\frac{n+N-1}{2}\right)}\times\frac{\varepsilon^{0}_{n,\rm
max}}{\varepsilon_{n}}~{}~{},$ (46)

which ultimately scales proportional to $n!/((n/2)!)^{2}$, quickly becoming
very large for large $n$. We also have that

$R_{0}^{4}\approx\frac{S_{4}}{\pi^{2}\Delta V(\Pi)}$ (47)

so substituting the above relations back to eq. (44) it becomes clear that for
$\varepsilon_{n,\rm max}^{0}/\varepsilon_{n}$ satisfying the criteria for a
controlled EFT expansion the exponential becomes extremely small. The
condition for a successful completion of the vacuum transition is

$\Gamma_{4}\gtrsim H^{4}~{}~{},$ (48)

which is difficult to fulfil. In conclusion, if the DS is for some reason
localized at the  purple dot then it will face an extremely slow decay rate,
compared with the expansion of the universe, such that it will never
completely tunnel to the  blue dot in a time scale which is relevant, leading
to an eternally-inflating DS.

Figure 4: The bounce solution $S_{4}$ evaluated numerically as a function of
$\varepsilon$ for the  green dot $\to$ red dot transition as they are
represented in Fig. 3.

Naturally one is led to consider the other tunneling possibilities. Naïvely
for transitions closer to the true global minimum one, such as the  green to
red dot vacuum transition (see Fig. 3), one does not expect a dramatic change
since the difference in vacuum energy and the height of the barrier grow in a
correlated manner, however eq. (45) suggests that the change in vacuum energy
may ultimately dominate such that faster tunnelling may be possible. In such
transitions the energy difference is comparable to the barrier height, hence
the thin wall approximation cannot be trusted and a numerical analysis of the
bounce action is required. To this end we rely again on a modified version of
CosmoTransitions Wainwright:2011kj  code. The numerical analysis of the bounce
solution as a function of $\varepsilon_{n}$, for the benchmark scenario
studied here, is shown in Fig. 4, demonstrating that only a case of a large
$\varepsilon_{n}$, well above the upper value for an effective description of
the pNGB potential, admits values of $S_{4}$ which could allow the vacuum
transition to complete.

To conclude, we find that a supercooled vacuum transition in a DS with a
single Gegenbauer potential and $\varepsilon_{n}>0$, is highly unlikely to
successfully complete unless $\varepsilon_{n}$ violates the EFT bound, in
which case calculability is called into question.

### PT from a flipped potential

Now consider the case with $\varepsilon_{n}<0$, as displayed in Fig. 5. We
focus on the transition from the second minimum to the origin. Notice that
this process corresponds to a symmetry-restoring phase transition since the
pNGB order parameter $\Pi$ has a zero vev in the true vacuum. This transition
is outside the validity of the thin-wall approximation thus we compute the
constant decay rate, eq. (44), numerically. To estimate the bubble radius at
nucleation, $R_{0}$, we use the value at which the field profile function is
halfway between the two minima.

Figure 5:  Inverted tree-level Gegenbauer potential. With the transition from
the green dot to the red dot considered.

The Hubble rate is written as

$H^{2}\equiv\frac{\pi^{2}g^{*}_{SM}(T_{v})T_{v}^{4}}{90M_{\text{Pl}}^{2}}+\frac{\Delta
V(\Pi,0)}{3M_{\rm Pl}^{2}}~{},$ (49)

where the first term comes from the standard radiation degrees of freedom. The
second term above is the vacuum contribution and we have assumed that the DS
temperature remains negligibly small. For simplicity, we fix the value of
$\varepsilon_{n}=10^{-2}\varepsilon^{0}_{n,\rm max}$ and the resonance mass
scale to $M=4\pi f$. Thus only $N$, $n$ and the symmetry breaking scale $f$
are free parameters.

The tunnelling rate $\Gamma_{4}$ is independent of the visible sector
temperature and instead all the temperature dependence is encoded in eq. (49).
We also find that the polynomial order $n$ has a negligible impact on the
decay rate. Once one fixes $N$, $n$ and $f$, one has that
$\Gamma_{4}/H^{4}\propto 1/T_{v}^{8}$ for large $T_{v}$. As the temperature
drops the vacuum contribution starts dominating the Hubble rate and
$\Gamma_{4}/H^{4}\approx\text{const}$. This behavior is displayed in Fig. 6
for $N=10$ and $n=20$ and several values of symmetry breaking scale $f$. One
can observe from this figure that the nucleation temperature is directly
proportional to the compositeness scale $f$, as expected on dimensional
grounds.

Figure 6:  Ratio of nucleation rate to Hubble volume as a function of visible
sector temperature for different values of the compositeness scale. The
horizontal line marks the nucleation condition while the vertical lines help
visualize the intersection point. At high temperatures
$\Gamma_{4}/H^{4}\propto 1/T_{v}^{8}$ while as the temperature drops the
vacuum contribution begins dominating the Hubble rate and
$\Gamma_{4}/H^{4}\approx\text{const}$.

Notice that if a transition is too slow to occur at $T_{v}=0$ then it cannot
start for any $T_{v}$. In addition, since the potential is effectively
temperature-independent, the strength parameter of the phase transition is
approximately

$\alpha(T_{v})\approx\frac{\Delta V(\Pi,0)}{\rho_{R}}~{}.$ (50)

In Fig. 7 we show this transition strength (colorbar) alongside the behavior
of the nucleation temperature as a function of symmetry breaking scale for two
benchmark values of $N$. The number of pNGBs, $N$, significantly impacts the
possible range of nucleation temperature due to the fact that, in our chosen
parametrization, $N$ affects the barrier height and thus, through the bounce
action, impacts the tunneling rate exponentially. The lines terminate at the
symmetry breaking scale $f$ for which the nucleation rate matches the minimum
value $\Gamma_{4}\approx H^{4}$, as can be inferred from Fig. 6. Close to this
point, the nucleation condition becomes numerically ambiguous. For smaller
values of $f$ the lines are truncated at values with extremely weak vacuum
transitions. It can be observed that the strongest phase transitions are
associated with the largest possible symmetry breaking scale and can attain
values $\alpha\approx\mathcal{O}(1)$.

Figure 7:  Nucleation temperature as a function of symmetry breaking scale $f$
with the colorbar displaying the strength parameter $\alpha$ at the
percolation temperature.

For very strong phase transitions the latent heat released accelerates the
wall to relativistic velocities and the effects of the thermal plasma are
suppressed. Thus the DS plasma of pNGBs exerts negligible friction on the wall
and one has $v_{w}\approx 1$. In this case the GW signal is sourced by the
collision of the walls and not by the sound waves, thus the treatment differs
from sect. 4.1. To estimate the time scale of the transition we consider the
bubble number density, which for a constant decay rate reads Cutting:2018tjt
333In this expression, the gamma function $\Gamma(x)$ should not be confused
with the decay rate $\Gamma_{4}$.

$\frac{1}{R_{*}^{3}}=\frac{1}{4}\left(\frac{\Gamma_{4}}{v_{w}}\right)^{3/4}\Gamma\left(\frac{1}{4}\right)\left(\frac{3}{\pi}\right)^{1/4}=\frac{1}{8\pi}\frac{\beta^{3}}{v_{w}^{3}}.$
(51)

The GW spectrum from bubble collisions is estimated as Bringmann:2023opz

$\Omega_{\rm GW}(f)h^{2}=\tilde{\Omega}\times
S\left(\frac{f_{g}}{f_{\text{col}}}\right),$ (52)

Figure 8:  GW spectrum from bubble collisions for the strongest signals found.
The red contours are the violin curves for the NANOGrav $15$ yr data obtained
from NANOGrav:2023hvm  using the public tool Rohatgi2022 . The integrated
sensitivity curves for LISA and BBO were obtained using Thrane:2013oya .

where we write the amplitude in terms of mean bubble separation as

$\tilde{\Omega}\approx 1.7\times 10^{-5}\ \tilde{\Omega}_{\text{bw}}(H_{\rm
min}R_{*})^{2}(8\pi)^{-2/3}\left(\frac{\kappa_{\phi}\alpha(T_{p})}{1+\alpha(T_{p})}\right)^{2}\left(\frac{g_{*}(T_{p})}{100}\right)^{-1/3},$
(53)

where $H_{\rm min}^{2}=\Delta V/3M_{\rm Pl}^{2}$, the coefficient
$\kappa_{\phi}$ is obtained from the detonation approximation from
Espinosa:2010hh  and the spectral function is given by

$S(x)=\frac{19x^{14/5}}{5+14x^{19/5}}~{}.$ (54)

After red-shifting the peak amplitude we have that

$f_{\text{col}}=1.7\times
10^{-5}(R_{*}H_{\text{min}})^{-1}(8\pi)^{1/3}\left(\frac{T_{p}}{100\
\text{GeV}}\right)\left(\frac{g_{*}(T_{p})}{100}\right)^{1/6}\left(\frac{f_{\text{peak}}}{\beta}\right)~{}~{}\text{Hz},$
(55)

with $f_{\text{peak}}/\beta\approx 0.2$, $\tilde{\Omega}_{\rm bw}\approx
0.08$. In the expressions above we have used a slightly more precise
percolation temperature, at which the probability to find a region of space-
time still in the false vacuum has decreased to about $P(T_{p})\sim e^{-1}$.

We show, in Fig. 8, the predicted GW spectrum from bubble collisions for three
benchmark values of $N$ where in each case we select the value of $f$ which
maximises the strength of the phase transition. We display the sensitivities
of the future detectors LISA LISACosmologyWorkingGroup:2022jok ;
Caprini:2019pxz  and BBO Corbin:2005ny . As we can observe from this figure,
the case $N=9$ could potentially explain the recently observed common-red
spectrum from the NANOGrav $15$ yr data NANOGrav:2023hvm  which is shown as
the gray curves.

##  5 Summary and conclusions

The vacuum structure and dynamics of theories possessing pNGB fields in the IR
is of theoretical interest and physical importance. Indeed, the vacuum
structure of QCD itself is a rich subject rendered tractable by studying the
vacuum structure of the pNGBs Witten:1980sp ; Creutz:1995wf ; Smilga:1998dh ;
March-Russell:2020lkq . In this work we have explored a complementary facet of
pNGB vacua which arises if explicit symmetry breaking occurs due to a spurion
in a non-minimal representation. Here, again, there are metastable vacua,
however they exist for different field values, as described in Durieux:2021riy
. In this work, we have focused on the same $\text{SO}(N+1)\to\text{SO}(N)$
symmetry breaking pattern and investigated the resulting vacuum dynamics,
which are found to be much richer than one might naïvely expect.

Our main result is that the ‘primary’ phase transition associated with
spontaneous $\text{SO}(N+1)\to\text{SO}(N)$ breaking when the radial mode
obtains a vacuum expectation value is not the end of the story. Below this
scale the pNGBs will typically undergo additional vacuum transitions unless
the sources of explicit symmetry breaking take the most minimal form.

These vacuum transitions may occur in two ways. Thermally, there is a second
critical temperature scale, the ‘Flipping Temperature’, which scales
proportional to $T_{F}\propto f/n$ and can thus naturally be well below the
spontaneous symmetry breaking scale $f$. Crucially, at this temperature the
functional form of the pNGB potential remains the same, to leading order in
the spurion. However, the overall sign flips, such that the higher temperature
minimum becomes the lower temperature maximum, and vice-versa for the higher
temperature maximum. As a result, in the vicinity of the flipping temperature
an additional vacuum transition occurs. We find this is likely weakly first-
order, at least for parameters consistent with a controlled EFT.

The second possibility arises non-thermally, if the pNGB sector becomes
supercooled in a metastable state, which is not implausible given the
existence of $\sim n$ different metastable vacua. In this case multiple vacuum
transitions can occur, with the most likely being to a nearest neighbour. As
the field approaches the global minimum the final vacuum transition can be
strong enough to generate observable GWs.

The vacuum structure of our universe is of prime importance and interest in
physics. It determines the ultimate fate of the observable universe and may
carry lessons about the deep UV and quantum gravity itself Obied:2018sgi .
Spontaneous symmetry breaking is ubiquitous in nature, for which Nambu-
Goldstone bosons are the physical manifestation of the vacuum structure.
Similarly, pNGBs manifest, through their vacuum structure, patterns of
explicit symmetry breaking. As a result, physically relevant lessons
concerning the vacuum structure and cosmological dynamics of nature may be
learned by studying pNGBs, perhaps even the case in which the Higgs boson is a
pNGB; a case we leave to further study.

## Acknowledgements

The authors would like to thank Marek Lewicki and Andreas Mantziris for useful
discussions. The research of F.K., M.M., S.P. and K.S. leading to these
results has received funding from the Norwegian Financial Mechanism for years
2014-2021, grant nr DEC-2019/34/H/ST2/00707. M.M. also acknowledges support
from the Polish National Science Center grant 2018/31/D/ST2/02048. K.S. is
partially supported by the National Science Centre, Poland, under research
grant 2017/26/E/ST2/00135.

##  Appendix A Hot Sector Calculations

We now detail a numerical investigation of the phase transition for a hot DS.
The theory of the vacuum decay from a local false minima to the true global
minima at zero and finite temperature has been studied extensively
Coleman:1977py ; Callan:1977pt ; Coleman:1980aw ; Linde:1980tt ; Linde:1981zj
. When the temperature is non-negligible the transition proceeds through
thermal fluctuations by the nucleation of true vacuum bubbles within the space
filled with false vacuum energy. The probability of decay per unit time and
volume is

$\Gamma_{3}(T)=\left(\frac{S_{3}(T)}{2\pi
T}\right)^{3/2}T^{4}e^{-S_{3}(T)/T}~{},$ (56)

where $S_{3}(T)/T$ is the finite temperature Euclidean action of our pNGB
model and is less than the zero-temperature one, $S_{4}$, around $T_{F}$.

The true vacuum bubble nucleates when the decay rate becomes comparable to the
expansion rate of the universe. Namely, we define the bubble nucleation
temperature by

$\Gamma_{3}\approx H^{4}\big{|}_{T\equiv T_{n}}~{},$ (57)

where the Hubble rate is given by

$H^{2}=\frac{\rho_{R}}{3M_{\text{Pl}}^{2}}+\frac{\Delta
V(\Pi,T)}{3M_{\text{Pl}}^{2}}~{},$ (58)

which includes the contribution from the potential energy difference between
false and true minima and $M_{\text{Pl}}=2.4\times 10^{18}$ GeV is the reduced
Planck mass.

As mentioned earlier, the hidden and visible sectors have independent
temperatures and cool at different rates. From eq. (36) above we can read off
the total effective number of degrees of freedom as

$g_{*}=\left(N+\frac{g^{*}_{\rm SM}(T_{v})}{\xi_{DS}^{4}}\right)~{}.$ (59)

The time scale of the transition is given by

$\frac{\beta}{H}\equiv
T\frac{d}{dT}\left(\frac{S_{3}(T)}{T}\right)\bigg{|}_{T\to T_{n}}\,.$ (60)

To compute the action we solve the equation of motion for the system, also
known as the bounce solution. This can be considerably simplified by
considering the parametrization of eq. (5) and allowing for a vev only in the
$\Pi$ direction such that

$\Box\Pi-\frac{\partial V(\Pi,T)}{\partial\Pi}=0~{}~{}.$ (61)

We use a modified version of the publicly available code CosmoTransitions
Wainwright:2011kj  to compute the Euclidean action.

Finally, it is necessary to have an estimate for the bubble wall velocity.
This requires an out-of-equilibrium computation of the deviation from
equilibrium of all the particle distribution functions. While this is still a
very active area of research Moore:2000wx ; Dorsch:2021nje ;
Friedlander:2020tnq ; Cline:2020jre ; Cline:2021iff ; Laurent:2020gpg ;
Megevand:2013hwa ; Dorsch:2018pat ; Kozaczuk:2015owa ; Konstandin:2014zta ;
Moore:1995si ; Ai:2021kak ; Lewicki:2021pgr ; Ellis:2022lft ;
BarrosoMancha:2020fay ; Liu:1992tn ; Laurent:2022jrs , here we will adopt the
analytic estimate of Lewicki:2021pgr ; Ellis:2022lft

$v_{w}=\begin{cases}\sqrt{\frac{\Delta V}{\alpha\rho_{R}}}\quad\quad{\rm
for}\quad\sqrt{\frac{\Delta V}{\alpha\rho_{R}}}<v_{J}(\alpha)~{},\\\
1\quad\quad\quad\quad{\rm for}\quad\sqrt{\frac{\Delta V}{\alpha\rho_{R}}}\geq
v_{J}(\alpha)\,,\end{cases}$ (62)

where $\alpha$ is the transition strength given in eq. (35) and
$v_{J}=\frac{1}{\sqrt{3}}\frac{1+\sqrt{3\alpha^{2}+2\alpha}}{1+\alpha}$ the
Chapman-Jouguet velocity which defines the upper limit for which hydrodynamic
solutions can be found. Although this result is valid for simple extensions of
the SM, in our case, we expect it to give us a realistic estimate. The reason
is that we expect the friction force on the bubble wall to become significant
due to the mass of the pNGBs at the metastable vacuum.

The sound wave source template reads444We notice that there are several
templates for the GW which derive from fits to different numerical
simulations. In particular the template we use do not match those of, e.g.
Bringmann:2023opz  but we nevertheless expect that our conclusion remain
qualitatively the same regardless of which template is used.
Hindmarsh:2013xza ; Hindmarsh:2015qta ; Hindmarsh:2017gnf ; Hindmarsh:2020hop
as a function of the frequency, $f_{g}$,

$\Omega_{\rm sw}(f_{g})h^{2}=4.13\times
10^{-7}\,\left(R_{*}H_{*}\right)\left(1-\frac{1}{\sqrt{1+2\tau_{\rm
sw}H_{*}}}\right)\left(\frac{\kappa_{\rm
sw}\,\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}S_{\rm
sw}(f_{g})\,,$ (63)

where $R_{*}$ is the average bubble size at collision and the spectral
function is

$S_{\rm sw}(f_{g})=\left(\frac{f_{g}}{f_{\rm
sw}}\right)^{3}\left[\frac{4}{7}+\frac{3}{7}\left(\frac{f_{g}}{f_{\rm
sw}}\right)^{2}\right]^{-\frac{7}{2}},$ (64)

$g_{*}$ is given in eq. (59) and all the quantities of the GW spectrum are
evaluated at the nucleation temperature $T_{*}=T_{n}\approx T_{F}$. The
frequency at the peak of the spectrum is given by

$f_{\rm sw}\,=2.6\times 10^{-5}\,{\rm
Hz}\left(R_{*}H_{*}\right)^{-1}\left(\frac{T_{*}}{100{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,,$ (65)

while the duration of the sound wave source reads Ellis:2018mja ;
Ellis:2019oqb ; Hindmarsh:2017gnf ; Ellis:2020awk

$\tau_{\rm sw}H_{*}=\frac{H_{*}R_{*}}{U_{f}}\,,\quad
U_{f}\approx\sqrt{\frac{3}{4}\frac{\alpha}{1+\alpha}\kappa_{\rm sw}}\,.$ (66)

For the mean bubble separation we use

$H_{*}R_{*}\approx(8\pi)^{\frac{1}{3}}\left(\frac{\beta}{H}\right)^{-1}.$ (67)

For all the computations that follow in this subsection we have fixed the
explicit symmetry breaking parameter to
$\varepsilon_{n}=10^{-2}\times\varepsilon_{n,\rm max}^{T_{F}}$. For the UV
scale at which we expect new resonances to appear we have fixed $M=4\pi f$.
Lastly, since the details of the dynamics of the finite temperature phase
transition are to a good approximation controlled by the DS flipping
temperature $T_{F}$, our only free parameters for this analysis are $n$, $N$,
$f$, $\xi_{\rm DS}$ and $T_{v}$ where the visible sector temperature is fixed
above $T_{F}$.

Figure 9: The GW spectrum for sound waves with $n=10-58$ (orange curves). The
explicit symmetry breaking parameter has been set to
$\varepsilon=10^{-2}\times\varepsilon_{n,max}^{T_{F}}$. The symmetry breaking
scale was fixed to $f=1$ TeV and the number of pNGBs to $N=4$. The temperature
of the visible sector was fixed to $T_{v}=2T_{F}$ for each benchmark. Figure
10: The GW spectrum from sound waves for $n=20$. The explicit symmetry
breaking parameter has been set to
$\varepsilon_{n}=10^{-2}\times\varepsilon_{n,max}^{T_{F}}$. The symmetry
breaking scale was fixed to $f=1$ TeV and the number of pNGBs to $N=4$. The
temperature of the visible sector was fixed as specified on the plot legends.

We present the predictions of the GW spectrum in Fig. 9 and Fig. 10 below. In
Fig. 9 we display the variation of the signal as a function of the Gegenabuer
polynomial order, $n$, while fixing $N=4$, $f=1$ TeV and $T_{v}=2T_{F}$. We
notice that the amplitude of the signal is very small compared with the
expected experimental sensitivities, in particular we find $\alpha\approx
0.002$ (in agreement with our analytic prediction for the transition strength
given in eq. (41)), $\beta/H\approx 10^{6}$ and $v_{w}\approx 0.06$. We do not
observe strong dependence on the polynomial order $n$. Recall that
$T_{n}\approx T_{F}\sim f/n$, hence the flipping temperature is numerically
very close to the critical and the nucleation temperature.

In Fig. 10, we instead vary the ratio of hidden to visible temperatures by
choosing different values of $T_{v}/T_{F}$ while setting $N=4$, $n=20$ and
$f=1$ TeV. In this case we notice a substantial reduction in the amplitude as
we increase the temperature hierarchy. This is expected as the amplitude
formula eq. (63) is inversely proportional to the total number of degrees of
freedom, in agreement with the results of Breitbach:2018ddu . Furthermore we
have verified numerically that varying other parameters of the potential do
not substantially change the amplitude of the signal and, irrespectively of
the adopted benchmark, we obtain a strength parameter of about $\alpha\approx
0.002$ while for the inverse timescale $\beta/H\approx 10^{6}$ and
$v_{w}\approx 0.06$. These numerical values are indicative of a very weak and
quick transition, if not a crossover, motivating our initial choice of using
$T_{n}$ in the GW template formula rather than the percolation temperature.

## References

  * (1) Y. Nambu, Quasiparticles and Gauge Invariance in the Theory of Superconductivity, Phys. Rev. 117 (1960) 648–663.
  * (2) J. Goldstone, Field Theories with Superconductor Solutions, Nuovo Cim. 19 (1961) 154–164.
  * (3) S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43 (1975) 199–220. [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * (4) Y. B. Zeldovich, A New Type of Radioactive Decay: Gravitational Annihilation of Baryons, Phys. Lett. A 59 (1976) 254.
  * (5) Y. B. Zeldovich, A Novel Type of Radioactive Decay: Gravitational Baryon Annihilation, Zh. Eksp. Teor. Fiz. 72 (1977) 18–21.
  * (6) T. Banks and L. J. Dixon, Constraints on String Vacua with Space-Time Supersymmetry, Nucl. Phys. B 307 (1988) 93–108.
  * (7) T. Banks and N. Seiberg, Symmetries and Strings in Field Theory and Gravity, Phys. Rev. D 83 (2011) 084019, [arXiv:1011.5120].
  * (8) G. Durieux, M. McCullough, and E. Salvioni, Gegenbauer Goldstones, JHEP 01 (2022) 076, [arXiv:2110.06941].
  * (9) P. Binetruy and M. K. Gaillard, Temperature Corrections in the Case of Derivative Interactions, Phys. Rev. D 32 (1985) 931–937.
  * (10) R. Alonso, E. E. Jenkins, and A. V. Manohar, A Geometric Formulation of Higgs Effective Field Theory: Measuring the Curvature of Scalar Field Space, Phys. Lett. B 754 (2016) 335–342, [arXiv:1511.00724].
  * (11) L. Dolan and R. Jackiw, Symmetry Behavior at Finite Temperature, Phys. Rev. D 9 (1974) 3320–3341.
  * (12) Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * (13) T.-H. Yeh, J. Shelton, K. A. Olive, and B. D. Fields, Probing physics beyond the standard model: limits from BBN and the CMB independently and combined, JCAP 10 (2022) 046, [arXiv:2207.13133].
  * (14) M. Breitbach, J. Kopp, E. Madge, T. Opferkuch, and P. Schwaller, Dark, Cold, and Noisy: Constraining Secluded Hidden Sectors with Gravitational Waves, JCAP 07 (2019) 007, [arXiv:1811.11175].
  * (15) M. Fairbairn, E. Hardy, and A. Wickens, Hearing without seeing: gravitational waves from hot and cold hidden sectors, JHEP 07 (2019) 044, [arXiv:1901.11038].
  * (16) F. Ertas, F. Kahlhoefer, and C. Tasillo, Turn up the volume: listening to phase transitions in hot dark sectors, JCAP 02 (2022), no. 02 014, [arXiv:2109.06208].
  * (17) Y. Bai and M. Korwar, Cosmological constraints on first-order phase transitions, Phys. Rev. D 105 (May, 2022) 095015.
  * (18) T. Bringmann, P. F. Depta, T. Konstandin, K. Schmidt-Hoberg, and C. Tasillo, Does NANOGrav observe a dark sector phase transition?, arXiv:2306.09411.
  * (19) K. Saikawa and S. Shirai, Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model, JCAP 05 (2018) 035, [arXiv:1803.01038].
  * (20) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Gravitational waves from the sound of a first order phase transition, Phys. Rev. Lett. 112 (2014) 041301, [arXiv:1304.2433].
  * (21) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Numerical simulations of acoustically generated gravitational waves at a first order phase transition, Phys. Rev. D 92 (2015), no. 12 123009, [arXiv:1504.03291].
  * (22) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Shape of the acoustic gravitational wave power spectrum from a first order phase transition, Phys. Rev. D 96 (2017), no. 10 103520, [arXiv:1704.05871]. [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * (23) M. B. Hindmarsh, M. Lüben, J. Lumma, and M. Pauly, Phase transitions in the early universe, SciPost Phys. Lect. Notes 24 (2021) 1, [arXiv:2008.09136].
  * (24) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, Energy Budget of Cosmological First-order Phase Transitions, JCAP 06 (2010) 028, [arXiv:1004.4187].
  * (25) S. Coleman, Fate of the false vacuum: Semiclassical theory, Phys. Rev. D 15 (May, 1977) 2929–2936.
  * (26) M. J. Duncan and L. G. Jensen, Exact tunneling solutions in scalar field theory, Phys. Lett. B 291 (1992) 109–114.
  * (27) C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields, Comput. Phys. Commun. 183 (2012) 2006–2013, [arXiv:1109.4189].
  * (28) D. Cutting, M. Hindmarsh, and D. J. Weir, Gravitational waves from vacuum first-order phase transitions: from the envelope to the lattice, Phys. Rev. D 97 (2018), no. 12 123513, [arXiv:1802.05712].
  * (29) NANOGrav Collaboration, A. Afzal et al., The NANOGrav 15 yr Data Set: Search for Signals from New Physics, Astrophys. J. Lett. 951 (2023), no. 1 L11, [arXiv:2306.16219].
  * (30) A. Rohatgi, Webplotdigitizer: Version 4.6, 2022.
  * (31) E. Thrane and J. D. Romano, Sensitivity curves for searches for gravitational-wave backgrounds, Phys. Rev. D 88 (2013), no. 12 124032, [arXiv:1310.5300].
  * (32) LISA Cosmology Working Group Collaboration, P. Auclair et al., Cosmology with the Laser Interferometer Space Antenna, arXiv:2204.05434.
  * (33) C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, Reconstructing the spectral shape of a stochastic gravitational wave background with LISA, JCAP 11 (2019) 017, [arXiv:1906.09244].
  * (34) V. Corbin and N. J. Cornish, Detecting the cosmic gravitational wave background with the big bang observer, Class. Quant. Grav. 23 (2006) 2435–2446, [gr-qc/0512039].
  * (35) E. Witten, Large N Chiral Dynamics, Annals Phys. 128 (1980) 363.
  * (36) M. Creutz, Quark masses and chiral symmetry, Phys. Rev. D 52 (1995) 2951–2959, [hep-th/9505112].
  * (37) A. V. Smilga, QCD at theta similar to pi, Phys. Rev. D 59 (1999) 114021, [hep-ph/9805214].
  * (38) J. March-Russell and R. Petrossian-Byrne, QCD, Flavor, and the de Sitter Swampland, arXiv:2006.01144.
  * (39) G. Obied, H. Ooguri, L. Spodyneiko, and C. Vafa, De Sitter Space and the Swampland, arXiv:1806.08362.
  * (40) S. R. Coleman, The Fate of the False Vacuum. 1. Semiclassical Theory, Phys. Rev. D 15 (1977) 2929–2936. [Erratum: Phys.Rev.D 16, 1248 (1977)].
  * (41) C. G. Callan, Jr. and S. R. Coleman, The Fate of the False Vacuum. 2. First Quantum Corrections, Phys. Rev. D 16 (1977) 1762–1768.
  * (42) S. R. Coleman and F. De Luccia, Gravitational Effects on and of Vacuum Decay, Phys. Rev. D 21 (1980) 3305.
  * (43) A. D. Linde, Fate of the False Vacuum at Finite Temperature: Theory and Applications, Phys. Lett. B 100 (1981) 37–40.
  * (44) A. D. Linde, Decay of the False Vacuum at Finite Temperature, Nucl. Phys. B 216 (1983) 421. [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * (45) G. D. Moore, Electroweak bubble wall friction: Analytic results, JHEP 03 (2000) 006, [hep-ph/0001274].
  * (46) G. C. Dorsch, S. J. Huber, and T. Konstandin, A sonic boom in bubble wall friction, JCAP 04 (2022), no. 04 010, [arXiv:2112.12548].
  * (47) A. Friedlander, I. Banta, J. M. Cline, and D. Tucker-Smith, Wall speed and shape in singlet-assisted strong electroweak phase transitions, Phys. Rev. D 103 (2021), no. 5 055020, [arXiv:2009.14295].
  * (48) J. M. Cline and K. Kainulainen, Electroweak baryogenesis at high bubble wall velocities, Phys. Rev. D 101 (2020), no. 6 063525, [arXiv:2001.00568].
  * (49) J. M. Cline, A. Friedlander, D.-M. He, K. Kainulainen, B. Laurent, and D. Tucker-Smith, Baryogenesis and gravity waves from a UV-completed electroweak phase transition, Phys. Rev. D 103 (2021), no. 12 123529, [arXiv:2102.12490].
  * (50) B. Laurent and J. M. Cline, Fluid equations for fast-moving electroweak bubble walls, Phys. Rev. D 102 (2020), no. 6 063516, [arXiv:2007.10935].
  * (51) A. Mégevand, Friction forces on phase transition fronts, JCAP 07 (2013) 045, [arXiv:1303.4233].
  * (52) G. C. Dorsch, S. J. Huber, and T. Konstandin, Bubble wall velocities in the Standard Model and beyond, JCAP 12 (2018) 034, [arXiv:1809.04907].
  * (53) J. Kozaczuk, Bubble Expansion and the Viability of Singlet-Driven Electroweak Baryogenesis, JHEP 10 (2015) 135, [arXiv:1506.04741].
  * (54) T. Konstandin, G. Nardini, and I. Rues, From Boltzmann equations to steady wall velocities, JCAP 09 (2014) 028, [arXiv:1407.3132].
  * (55) G. D. Moore and T. Prokopec, How fast can the wall move? A Study of the electroweak phase transition dynamics, Phys. Rev. D 52 (1995) 7182–7204, [hep-ph/9506475].
  * (56) W.-Y. Ai, B. Garbrecht, and C. Tamarit, Bubble wall velocities in local equilibrium, JCAP 03 (2022), no. 03 015, [arXiv:2109.13710].
  * (57) M. Lewicki, M. Merchand, and M. Zych, Electroweak bubble wall expansion: gravitational waves and baryogenesis in Standard Model-like thermal plasma, JHEP 02 (2022) 017, [arXiv:2111.02393].
  * (58) J. Ellis, M. Lewicki, M. Merchand, J. M. No, and M. Zych, The scalar singlet extension of the Standard Model: gravitational waves versus baryogenesis, JHEP 01 (2023) 093, [arXiv:2210.16305].
  * (59) M. Barroso Mancha, T. Prokopec, and B. Swiezewska, Field-theoretic derivation of bubble-wall force, JHEP 01 (2021) 070, [arXiv:2005.10875].
  * (60) B.-H. Liu, L. D. McLerran, and N. Turok, Bubble nucleation and growth at a baryon number producing electroweak phase transition, Phys. Rev. D 46 (1992) 2668–2688.
  * (61) B. Laurent and J. M. Cline, First principles determination of bubble wall velocity, Phys. Rev. D 106 (2022), no. 2 023501, [arXiv:2204.13120].
  * (62) J. Ellis, M. Lewicki, and J. M. No, On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal, JCAP 04 (2019) 003, [arXiv:1809.08242].
  * (63) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, Gravitational wave energy budget in strongly supercooled phase transitions, JCAP 06 (2019) 024, [arXiv:1903.09642].
  * (64) J. Ellis, M. Lewicki, and J. M. No, Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source, JCAP 07 (2020) 050, [arXiv:2003.07360].