# How robust are gravitational wave predictions from cosmological phase
transitions?

Peter Athron<EMAIL_ADDRESS>Department of Physics and Institute of
Theoretical Physics, Nanjing Normal University, Nanjing, 210023, China
Lachlan Morris<EMAIL_ADDRESS>School of Physics and Astronomy,
Monash University, Melbourne, Victoria 3800, Australia  Zhongxiu Xu
<EMAIL_ADDRESS>Department of Physics and Institute of Theoretical
Physics, Nanjing Normal University, Nanjing, 210023, China

###### Abstract

Gravitational wave (GW) predictions of cosmological phase transitions are
almost invariably evaluated at either the nucleation or percolation
temperature. We investigate the effect of the transition temperature choice on
GW predictions, for phase transitions with weak, intermediate and strong
supercooling. We find that the peak amplitude of the GW signal varies by a
factor of a few for weakly supercooled phase transitions, and by an order of
magnitude for strongly supercooled phase transitions. The variation in
amplitude for even weakly supercooled phase transitions can be several orders
of magnitude if one uses the mean bubble separation, while the variation is
milder if one uses the mean bubble radius instead. We also investigate the
impact of various approximations used in GW predictions. Many of these
approximations introduce at least a 10% error in the GW signal, with others
introducing an error of over an order of magnitude.

##  I Introduction

We are now in an era where existing gravitational wave (GW) data can have an
impact on our understanding of physics beyond the Standard Model (BSM) of
particle physics. Very recently pulsar timing array experiments have detected
a stochastic GW background (SGWB) [1, 2, 3, 4] and find that new physics
explanations have a slight preference over less exotic sources [5]. Existing
data on GWs from the LIGO/VIRGO network [6] is also constraining well-
motivated Pati-Salam models that can lead to gauge coupling unification [7] as
well as models of the dark sector [8].

However, with this exciting progress also comes significant challenges. It is
now essential that we have reliable calculations of the GW spectra for BSM
models where we understand the uncertainties involved and the effects of
various approximations and assumptions that are commonly used. There are many
challenging calculations involved in going from a particular BSM scenario to a
predicted GW spectrum; see Ref. [9] for a review. Quantities derived from the
effective potential can strongly depend on the method used [10] and
uncertainties in the GW spectra from effective potential computations have
been investigated in Ref. [11]. Here we show that even if the effective
potential calculation was under full control, there are many other challenges
for reliable predictions of GW spectra.

Since the first direct detection of GWs [12] in 2015, there has been
substantial progress in understanding how to characterise phase transitions
and extract GW predictions. Here we mention a few important points. Sound
waves are expected to be the largest source of GWs following Ref. [13] which
showed that sound waves source last long after the bubbles have merged.
However, more recently it has been shown that in many cases the lifetime is
nonetheless significantly shorter than the Hubble time [14, 15] and
suppression factors were introduced [16, 17] to account for the finite
lifetime of the source. These suppression factors were subsequently refined to
address issues stemming from the derivation of the Hubble time as the maximum
lifetime of the source [18]. Furthermore, the modelling of GWs from sound
waves has improved considerably from simulations [19, 20] and the construction
of the sound shell model [21] and its further development [22, 23, 24].
Significant improvements have also been made in determining the kinetic energy
fraction that is available to source GWs. New parameterisations have been
developed that go beyond simplified models such as the bag model, first for
the case where bubbles expand as supersonic detonations [25] and later
generalised to cover subsonic deflagrations and hybrids [26]. These advances
have both improved predictions and raised questions about our previous and
current understanding of how sensitive GW experiments can be to first-order
phase transitions.

In particular, strongly supercooled phase transitions present significant
challenges for calculations and may lead to erroneous explanations of GW
signals [27]. We therefore treat the extent of supercooling as an important
parameter when considering the uncertainties and compare scenarios with weak,
intermediate, and strong supercooling. Previously, we have shown that in the
presence of supercooling various possible choices of transition temperature
decouple [28] and it has been argued that the percolation temperature should
be used [17, 29, 30, 28]. Here we show explicitly that the peak amplitude and
frequency of the GW spectrum — and thus the resulting signal-to-noise ratio
(SNR) at a detector — are sensitive to the choice of transition temperature.
This is especially true for strongly supercooled phase transitions as one
might expect, but is also true for weakly supercooled phase transitions. We
show that if one chooses the nucleation temperature as the transition
temperature (as is very common practice), then the peak amplitude, peak
frequency, and SNR can change by orders of magnitude compared to when using
the percolation temperature. This has a huge impact on the the prospects for
detection. However, such a drastic change only arises when using the mean
bubble separation as the characteristic length scale. If one is more careful
about the choice of length scale, the discrepancy can potentially be reduced
to a factor of a few.

Additionally, we investigate how the predictions can be affected by different
estimates of the thermal parameters which determine the GW spectrum. We
compare various parameterisations of the kinetic energy fraction, which
determines the energy available for sourcing GWs. Another important factor
that determines the peak GW amplitude and frequency is the timescale during
which the source is active, which is usually replaced by a characteristic
length scale. The mean bubble separation is used as this length scale in
lattice simulations. We compare the impact different estimates of this have on
GW signals, and we qualitatively explore the consequences of using the mean
bubble radius instead. Finally, because the turbulence contribution to the
overall GW signal is not well modelled, but could be significant, we also
compare several different choices for the energy available for sourcing GWs
from turbulence and show the impact that this can have on the SNR.

In section II we describe first-order phase transitions and supercooling in
more detail, and we define important milestone temperatures. In section III we
describe how properties of the phase transition and the thermal parameters are
computed in particle physics models. We also discuss various estimates for
these thermal parameters that are made in the literature. We briefly describe
how we use these thermal parameters to predict GW spectra in section IV. We
then introduce the model we use to obtain a first-order phase transition in
section V. Finally, we present our results in section VI and provide
concluding remarks in section VII.

##  II First-order phase transitions and supercooling

As the Universe cools down the shape of the effective potential changes such
that minima (or phases) can appear and disappear and cosmological phase
transitions take place. These cosmological phase transitions play an important
role in particle physics, such as breaking the electroweak symmetry and
thereby generating masses for the fundamental particles via the Higgs
mechanism. Further, if a phase transition is of first order (i.e. there is a
potential barrier separating the phases), GWs are produced in the process.

A potential barrier between the phases prevents an instantaneous transition
from the local metastable minimum to the deeper minimum on the other side of
the barrier. Instead, the phase transition must proceed via either tunnelling
through the barrier or fluctuating over it. This first becomes possible when
the Universe cools below the critical temperature, $T_{c}$, where the free
energy densities of the two minima are degenerate. Below $T_{c}$ the
transition begins through a stochastic process where the tunnelling or
fluctuations occur at localised points in space-time, and when this happens
bubbles of the new phase can form and grow in a process known as bubble
nucleation. The phase transition completes if the bubbles of the new phase
fill the whole universe. More precisely, because it is a stochastic process we
define the the completion temperature, $T_{f}$, to be the temperature when the
fraction of the universe left in the false vacuum (i.e. the old phase) is less
then $1\%$, $P_{f}(T_{f})<0.01$.

When this process takes a long time to complete $T_{f}$ may be much smaller
than the critical temperature $T_{c}$ at which the new minimum first becomes
energetically favoured. This is known as supercooling in analogy with the
phenomenon where liquids are supercooled well below their freezing point. All
first-order cosmological phase transitions exhibit some degree of supercooling
because they do not happen instantly. However, the temperature change can vary
from $T_{f}$ being within 1% of $T_{c}$ to being orders of magnitude smaller.

The degree of supercooling can have a significant impact on a phase transition
and is an important characteristic when comparing phase transitions.
Increasing supercooling may boost the energy released in the phase transition
and the amplitude of resultant GWs, but too much supercooling can lead to the
transition failing to complete.

Strongly supercooled phase transitions admit qualitatively different behaviour
compared to weakly supercooled phase transitions. Because the nucleation rate
is lower, the smaller number of bubbles that are nucleated grow to much larger
sizes. This means that the number of bubbles per Hubble volume, $N$, can be
less than one during the period where most of the bubbles are colliding or
even by the time the phase transition has completed [28]. This can be
expressed more precisely as follows. The nucleation temperature $T_{n}$ is
defined by the condition $N(T_{n})=1$. Usually $T_{n}$ is higher than the
percolation temperature $T_{p}$, defined by the moment when the false vacuum
fraction, $P_{f}$, is roughly 71%: $P_{f}(T_{p})=0.71$. Roughly speaking,
$T_{p}$ is where the bubbles should be in contact with each other (see section
4.7.2 of Ref. [9] for more details). In strongly supercooled scenarios the
nucleation temperature can be reached some time after most of the the bubble
collisions have taken place. In more extreme cases the phase transition may
complete, reaching $P_{f}(T_{f})<0.01$, before $N(T)=1$. In such cases there
is no nucleation temperature. However, strongly supercooled scenarios can also
have enough bubble nucleation such that $N(T)=1$ is reached relatively early
in the phase transition but the transition is still slow leading to a
substantial gap between $T_{n}$ and $T_{p}$ or $T_{f}$. Thus, the nucleation
temperature is not coupled with the actual progress of the phase transition
and the production of GWs.

##  III Determining properties of the phase transition

###  III.1 Temperatures and length scales

The rate of a phase transition depends strongly on the size and persistence of
the potential barrier. In fast transitions the barrier disappears fairly
quickly. The nucleation rate is initially zero at $T_{c}$ and then increases
rapidly as the barrier dissolves, giving an exponential nucleation rate of the
form

$\displaystyle\Gamma(t)=\Gamma(t_{*})\exp(\beta(t-t_{*})),$ (1)

where $t_{*}$ is some relevant time in the transition (often taken to
correspond to $T_{n}$). In contrast, if the barrier persists at low
temperatures or even at $T=0$, the nucleation rate can instead reach a maximum
at some temperature $T_{\Gamma}$ because lower temperature reduces the
likelihood of thermal fluctuations over the barrier.

The nucleation rate is given by [31]

$\displaystyle\Gamma(T)=T^{4}\left(\frac{S(T)}{2\pi}\right)\exp(-S(T)),$ (2)

where $S(T)$ is the bounce action which we obtain from a modified version of
CosmoTransitions [32].111See appendix F of Ref. [28] for details of the
modifications. If one expresses $S$ as a function of time and Taylor expands
about $t_{*}$,

$\displaystyle S(t)\approx S(t_{*})$
$\displaystyle+\left.\derivative{S}{t}\right|_{t=t_{*}}\\!\\!\\!\\!(t-t_{*})$
(3)
$\displaystyle{}+{}\frac{1}{2}\\!\left.\derivative[2]{S}{t}\right|_{t=t_{*}}\\!\\!\\!\\!(t-t_{*})^{2}+\cdots,$
(4)

then truncating at first order gives the exponential nucleation rate given in
eq. 1, and we can identify

$\displaystyle\beta=-\\!\left.\derivative{S}{t}\right|_{t=t_{*}}.$ (5)

This can be useful because $\beta$ is related to the mean separation of
bubbles, $R_{\text{sep}}$, through [33]

$R_{\text{sep}}=(8\pi)^{\frac{1}{3}}\frac{v_{w}}{\beta}.$ (6)

The mean bubble separation is an important quantity for GW predictions.
Equation 6 should hold when evaluated at the temperature where $P_{f}$ has
decreased to $1/e$, denoted by $T_{e}$. Computing $\beta$ directly from the
bounce action and using eq. 6 to estimate $R_{\text{sep}}$ can simply
calculations significantly.

However, while an exponential nucleation rate is a common assumption and eq. 6
is widely used, these approximations can be problematic in strongly
supercooled scenarios. We will demonstrate the potential consequences of this
in section VI. Note that if the transition temperature $T_{*}$ used to
evaluate $\beta$ is close to the temperature where nucleation rate is
maximised, $T_{\Gamma}$, then $\beta\approx 0$. Further, $\beta$ is negative
when $T_{*}<T_{\Gamma}$. Therefore, the use of $\beta$ entirely breaks down in
these cases. However, because $\beta$ vanishes one can truncate eq. 4 at
second order and obtain a Gaussian nucleation rate,

$\displaystyle\Gamma(t)=\Gamma(t_{*})\exp(\frac{\beta_{V}^{2}}{2}(t-t_{*})^{2}),$
(7)

where

$\displaystyle\beta_{V}=\sqrt{\left.\derivative[2]{S}{t}\right|_{t=t_{\Gamma}}}.$
(8)

We can relate $\beta_{V}$ to $R_{\text{sep}}$ through [14]

$R_{\text{sep}}=\left(\sqrt{2\pi}\frac{\Gamma(T_{\Gamma})}{\beta_{V}}\right)^{\\!\\!-\frac{1}{3}}.$
(9)

It is unclear how well the approximations eq. 6 and eq. 9 perform, so we
include this investigation in our study. We note that we use temperature
rather than time in our analysis, so we employ the usual time-temperature
relation [28]

$\derivative{t}{T}=\frac{-1}{TH(T)}.$ (10)

Thus, $\beta$ and $\beta_{V}$ are in fact calculated from
$\text{d}S/\text{d}T$. The Hubble rate is given by

$\displaystyle H(T)=\sqrt{\frac{8\pi G}{3}\rho_{\text{tot}}(T)},$ (11)

where $\rho_{\text{tot}}$ is the total energy density. We use energy
conservation such that $\rho_{\text{tot}}=\rho_{f}-\rho_{\text{gs}}$, where
$\rho_{f}$ is the false vacuum energy density and $\rho_{\text{gs}}$ is the
ground state energy density. We renormalise the free energy density such that
$\rho_{\text{gs}}=0$, leaving $\rho_{\text{tot}}=\rho_{f}$.

Returning to the full treatment, the nucleation rate in eq. 2 can be used
directly to compute the false vacuum fraction $P_{f}$ as a function of
temperature, given by

$P_{f}(T)=\exp\\!\left[-\frac{4\pi}{3}\\!\int_{T}^{T_{c}}\\!\frac{dT^{\prime}}{T^{\prime
4}}\frac{\Gamma(T^{\prime})}{H(T^{\prime})}\\!\left(\int_{T}^{T^{\prime}}\\!\\!\\!dT^{\prime\prime}\frac{v_{w}(T^{\prime\prime})}{H(T^{\prime\prime})}\right)^{\\!\\!3}\right].$
(12)

Here we have assumed that the Universe is expanding adiabatically and we
neglect the initial radius of the bubble at formation. See Ref. [9] for more
details on the derivations and assumptions. The last undetermined quantity in
eq. 12 is the bubble wall velocity, $v_{w}$. We discuss our treatment of
$v_{w}$ in section III.2.

The number of bubbles nucleated at any given temperature can also be computed
from eq. 2. In the literature it is standard to calculate the nucleation
temperature from an approximation for the number of bubbles per Hubble volume,

$\displaystyle N(T)$
$\displaystyle=\int_{T}^{T_{c}}\\!\\!dT^{\prime}\,\frac{\Gamma(T^{\prime})}{T^{\prime}H^{4}(T^{\prime})}.$
(13)

This implicitly assumes a fast transition so that one can assume $P_{f}=1$
before $T_{n}$, and thus omit $P_{f}$ from the integrand [28].222A factor of
$4\pi/3$ from the spherical Hubble volume is also neglected in this treatment.
In this study we only use $T_{n}$ to show the impact of approximations made in
the literature, so we use the expression in eq. 13 to calculate $T_{n}$ for
consistency.

In contrast, to compute the mean bubble separation we determine the bubble
number density with $P_{f}(T)$ included to account for the fact that true
vacuum bubbles can only nucleate in regions that are still in the false
vacuum. The mean bubble separation is given by

$R_{\text{sep}}(T)=(n_{B}(T))^{-\frac{1}{3}},$ (14)

where

$n_{B}(T)=T^{3}\\!\\!\int_{T}^{T_{c}}\\!dT^{\prime}\frac{\Gamma(T^{\prime})P_{f}(T^{\prime})}{T^{\prime
4}H(T^{\prime})}$ (15)

is the bubble number density. Finally, there are possibly other choices for
the characteristic length scale in GW predictions [34, 35, 14, 29, 9].
However, fits for GW predictions are determined in terms of $R_{\text{sep}}$,
and one cannot directly replace $R_{\text{sep}}$ with alternative length
scales in those fits. Still, we seek to investigate (among other things) the
impact of the choice of $T_{*}$ on the GW predictions (see section VI), so it
is important to understand the impact of $T_{*}$ on various length scales.
Thus, we also consider the mean bubble radius,

$\bar{R}(T)=\frac{T^{2}}{n_{B}(T)}\int_{T}^{T_{c}}\\!\\!dT^{\prime}\frac{\Gamma(T^{\prime})P_{f}(T^{\prime})}{T^{\prime
4}H(T^{\prime})}\\!\int_{T}^{T^{\prime}}\\!\\!dT^{\prime\prime}\frac{v_{w}(T^{\prime\prime})}{H(T^{\prime\prime})}.$
(16)

For more details see section 5.5 of Ref. [9] and references therein.

We can compute the milestone temperatures $T_{n}$, $T_{p}$, $T_{e}$ and
$T_{f}$ using eqs. 13 and 12, and we can similarly use eqs. 14 and 16 to
compute $R_{\text{sep}}$ and $\bar{R}$ at these milestone temperatures or at
arbitrary temperatures. We use PhaseTracer [36] to map the phase structure of
the potential and TransitionSolver [37] to analyse the phase history,
including all relevant phase transitions,333There is another high-temperature
phase transition with $T_{c}\sim 180$ GeV in the intermediate and strong
supercooling scenarios considered in section V. The phase transition is very
fast and is not relevant to our analysis. as well as determine the milestone
temperatures and relevant parameters for GW predictions. The GW fits are
parameterised in terms of thermal parameters, which — in addition to the
transition temperature and the characteristic length scale — also include
hydrodynamic parameters such as the kinetic energy fraction and the bubble
wall velocity.

###  III.2 Hydrodynamic parameters

Here we discuss the hydrodynamic parameters used in GW fits. First we discuss
our best treatment for these parameters, then we introduce several common
variations to this treatment used in the literature. We will investigate the
impact of these variations on the GW signature in section VI.2. All of these
parameters — and all of the quantities that they depend on — should be
evaluated at the transition temperature, $T_{*}$.

In our best treatment, we take $T_{*}=T_{p}$, and we determine the kinetic
energy fraction using the pseudotrace difference between the phases,
corresponding to M2 in Ref. [25]:

$K=\frac{\bar{\theta}_{f}(T_{*})-\bar{\theta}_{t}(T_{*})}{\rho_{\text{tot}}(T_{*})}\kappa_{\bar{\theta}}(\alpha_{\bar{\theta}}(T_{*}),c_{s,f}(T_{*})).$
(17)

Here, $c_{s,f}$ is the speed of sound in the false vacuum and
$\alpha_{\bar{\theta}}$ is the transition strength parameter. The speed of
sound in phase $i$ is given by

$c_{s,i}^{2}(T)=\left.\frac{\partial_{T}V}{T\partial_{T}^{2}V}\right|_{\boldsymbol{\phi}_{i}(T)},$
(18)

where $V$ is the effective potential, or free energy density, and
$\boldsymbol{\phi}_{i}$ is the field configuration for phase $i$. The
transition strength parameter is defined as

$\alpha_{x}(T)=\frac{4\\!\left(x_{f}(T)-x_{t}(T)\right)}{3w_{f}(T)},$ (19)

where $x$ is a hydrodynamic quantity for which various choices exist in the
literature, and $w_{f}$ is the enthalpy density in the false vacuum. We use
the pseudotrace for $x$ in our best treatment, given by [25]

$\bar{\theta}_{i}(T)=\frac{1}{4}\\!\left(\rho_{i}(T)-\frac{p_{i}(T)}{c_{s,t}^{2}(T)}\right)$
(20)

in phase $i$, where $\rho$ and $p$ are the energy density and pressure,
respectively. The pseudotrace generalises the trace anomaly to models where
the speed of sound deviates from $1/\sqrt{3}$. We use the code snippet
provided in the appendix of Ref. [25] to determine the efficiency coefficient
$\kappa_{\bar{\theta}}$.

Turbulence from cosmological phase transitions is not well understood because
current hydrodynamic simulations cannot probe the turbulent regime. Hence, it
is difficult to estimate the efficiency coefficient for turbulence,
$\kappa_{\text{turb}}$, which is needed for turbulence contributions to the
production of GWs. However, it is expected that stronger phase transitions
(with larger $\alpha$) could result in more turbulence developing sooner and
could reduce the lifetime of the sound wave source. Lacking sufficient
modelling of the turbulence source, we consider the efficiency coefficient as
a fraction of $\kappa_{\bar{\theta}}$,

$\kappa_{\text{turb}}=\epsilon\kappa_{\bar{\theta}},$ (21)

and we take $\epsilon=0.05$ as our default treatment.

Finally, for our treatment of the bubble wall velocity, we assume bubbles grow
as supersonic detonations regardless of the extent of supercooling for
simplicity. General friction estimates are beyond the scope of this study, and
neither the ultra-relativistic or non-relativistic limits of friction are
applicable for all benchmark points in our study. We assume the bubbles expand
at the Chapman-Jouguet velocity,

$v_{w}=v_{\text{CJ}}=\frac{1+\sqrt{3\alpha_{\bar{\theta}}(1+c_{s,f}^{2}(3\alpha_{\bar{\theta}}-1))}}{c_{s,f}^{-1}+3\alpha_{\bar{\theta}}c_{s,f}},$
(22)

where temperature dependence has been suppressed. The Chapman-Jouguet velocity
is by no means the most likely supersonic detonation solution, however it does
capture dependence on the transition temperature and ensures a supersonic
detonation regardless of the extent of supercooling. The same cannot be said
for any fixed choice of $v_{w}$.

We now turn to the variations on our best treatment. First, we consider the
effect of setting $T_{*}$ to the other milestone temperatures: $T_{n}$,
$T_{e}$ and $T_{f}$. This involves using our best treatment (e.g. calculating
$K$ using eq. 17) but evaluating all quantities at, for example, $T_{n}$
instead of $T_{p}$. As a reminder, $T_{n}$ can be obtained by the condition
$N(T_{n})=1$ (see eq. 13), while $T_{p}$, $T_{e}$ and $T_{f}$ all come from
conditions on the false vacuum fraction (see eq. 12); specifically,
$P_{f}(T_{p})=0.71$, $P_{f}(T_{e})=1/e$ and $P_{f}(T_{f})=0.01$.

The approach we use for estimating $K$ was developed only recently in Refs.
[25, 26], so it is not yet widely adopted. More approximate treatments are
widespread, which we enumerate here. It is very common to determine $K$
through

$K_{x}=\frac{\kappa_{x}\alpha_{x}}{1+\alpha_{x}},$ (23)

with various choice of $x$ often being made. This parameterisation alone
introduces error in the determination of $K$, regardless of the choice of $x$
(see appendix A for details). The trace anomaly,

$\theta(T)=\frac{1}{4}(\rho(T)-3p(T)),$ (24)

is the closest alternative to $\bar{\theta}$, in fact exactly matching
$\bar{\theta}$ when $c_{s,t}=1/\sqrt{3}$ like in the bag model. The other
common choices for $x$ are the pressure $p$ and the energy density $\rho$. The
efficiency coefficient used for these choices of $x$ was derived in the bag
model, and is given by [38]

$\kappa=\frac{\sqrt{\alpha_{x}}}{0.135+\sqrt{0.98+\alpha_{x}}}$ (25)

for $v_{w}=v_{\text{CJ}}$, which is implicitly dependent on temperature.

In these more approximate treatments of $K$, the enthalpy density in the
denominator of eq. 19 is usually replaced with $w_{f}=\frac{4}{3}\rho_{R}$,
where $\rho_{R}=\frac{\pi^{2}}{30}g_{\text{eff}}T^{4}$ is the radiation energy
density and $g_{\text{eff}}$ is the effective number of relativistic degrees
of freedom. We find the replacement of the enthalpy density in this way (which
comes from the bag model) to be a very good approximation. This replacement
leads to less than 1% error in the GW predictions. Therefore our
$\alpha_{\rho}$ effectively corresponds to the latent heat definition
frequently found in the literature, see eq. 5.35 of Ref. [9]. Similarly
$\alpha_{\theta}$ also effectively corresponds to eq. 5.36 of Ref. [9], which
also frequently appears in the literature, though here one also needs to
substitute $\theta=\frac{1}{4}(\rho-3p)$. One could also replace
$\bar{\theta}$ with $\theta$ in eq. 17 and use eq. 25 for $\kappa$,
corresponding to M3 in Ref. [25]. However, we find this introduces at most 1%
difference in the GW predictions compared to using eq. 23 with $x=\theta$, so
we do not consider this variation in our results.

As described in section III.1, one can approximate the mean bubble separation
$R_{\text{sep}}$ through the often-used thermal parameter $\beta$, or through
$\beta_{V}$. We investigate the error in these approximations for
$R_{\text{sep}}$ and the corresponding effect on GW predictions. We also
demonstrate the impact of using $\bar{R}$ instead of $R_{\text{sep}}$, but we
do not treat this as a variation of the treatment because mapping $\bar{R}$
onto existing GW fits is currently ambiguous.

We also consider alternative treatments of the turbulence efficiency
coefficient. The most obvious variation is to simply choose another arbitrary,
fixed value. We choose $\epsilon_{2}=0.1$, where the subscript ‘2’ denotes the
index of this variation for $\epsilon$. We also consider
$\epsilon_{3}=\left(1-\min(H(T_{*})\tau_{\text{sw}},1\right))^{2/3}$, which
comes from assuming that a reduction in the lifetime of the sound wave source
$\tau_{\text{sw}}$ could boost the turbulence contribution to GW production
[39, 16]. However, the effective lifetime of the sound wave source is more
accurately suppressed by the factor
$\Upsilon=1-1/\sqrt{1+2H(T_{*})\tau_{\text{sw}}}$ derived in Ref. [18]. This
motivates a slightly modified choice: $\epsilon_{4}=(1-\Upsilon)^{2/3}$.

There are of course many other variations to the treatment that could be
considered, but we restrict our study to the variations mentioned thus far.
Changes to the bubble wall velocity could significantly impact the GW
predictions and even the phase transition properties, particularly if the
expansion mode of the bubbles changes from a supersonic detonation.
TransitionSolver currently does not use a full hydrodynamic treatment of
bubble profiles and therefore only provides accurate predictions for
supersonic detonations.444Reheating in the false vacuum for other bubble
expansion modes affects both bubble nucleation and growth [40, 41, 42]. Thus,
we currently cannot explore effect of $v_{w}$ on GW predictions. We explored
the impact of approximations made for the reheating temperature and GW
redshifting factors in Ref. [27], and found that their effects were small. We
do not reconsider these approximations here due to their accuracy. Also, we
explored accuracy of various approximations for $T_{n}$ as a function of
supercooling in Ref. [28]. Here we only calculate $T_{n}$ using eq. 13, but we
note that rougher approximations for $T_{n}$ are unreliable in strongly
supercooled scenarios, and would thus lead to significant errors in GW
predictions.

##  IV Gravitational waves

We consider only the sound wave and turbulence sources of GWs in this study.
The collision source is expected to contribute negligibly due to friction with
the plasma. Even though some of the benchmark points listed in section V admit
strong supercooling, the bubbles nucleate at temperatures where the plasma
still imposes significant friction on the expanding bubble walls. Thus, we do
not expect runaway bubble walls and consequently neglect the collision source
altogether.

The general scaling of the GW equations is predominantly governed by two key
parameters: the kinetic energy fraction $K$ and the characteristic length
scale $L_{*}$. We set $L_{*}=R_{\text{sep}}(T_{p})$ in our best treatment. The
scaling of the peak amplitude $\Omega_{\text{peak}}$ and the peak frequency
$f_{\text{peak}}$ is roughly

$\displaystyle\Omega_{\text{peak}}$ $\displaystyle\propto K^{n}L_{*},$ (26)
$\displaystyle f_{\text{peak}}$ $\displaystyle\propto L_{*}^{-1},$ (27)

where $n=2$ for sound waves and $n=3/2$ for turbulence.

The details of the GW equations we use can be found in appendix A.5 of Ref.
[27]. In addition to the turbulence fit [43] and the sound shell model [21,
22] used for the sound wave source, we also consider another fit for the sound
wave source provided in Ref. [19]. We will refer to this fit as the ‘lattice
fit’ for the sound wave source, for lack of a better name. In this fit, the
redshifted peak amplitude is

$h^{2}\Omega_{\text{sw}}^{\text{lat}}(f)=5.3\\!\times\\!10^{-2}\,\mathcal{R}_{\Omega}K^{2}\\!\left(\\!\frac{HL_{*}}{c_{s,f}}\\!\right)\\!\Upsilon(\tau_{\text{sw}})S_{\text{sw}}(f),$
(28)

the redshifted peak frequency is

$f_{\text{sw}}^{\text{lat}}=1.58\,\mathcal{R}_{f}\\!\left(\frac{1}{L_{*}}\right)\\!\left(\frac{z_{p}}{10}\right),$
(29)

matching one of the key frequencies in the sound shell model, and the spectral
shape is

$S_{\text{sw}}(f)=\left(\frac{f}{f_{\text{sw}}^{\text{lat}}}\right)^{\\!\\!3}\\!\left(\frac{7}{4+3(f/f_{\text{sw}}^{\text{lat}})^{2}}\right)^{\\!\\!\frac{7}{2}}.$
(30)

See Ref. [9] and the appendices of Ref. [27] for details of the redshifting
factors $\mathcal{R}_{f}$ and $\mathcal{R}_{\Omega}$, the lifetime suppression
factor $\Upsilon$, and the simulation-derived factor $z_{p}$ (which is taken
to be $z_{p}=10$). All quantities in the fit are evaluated at $T_{*}$, except
for the redshifting factors. These are instead evaluated at the reheating
temperature, which itself depends on $T_{*}$. Just as in Ref. [27], we do not
include a suppression factor coming from bubbles not reaching their asymptotic
hydrodynamic profiles in the simulations from which the GW fits are obtained.
This suppression factor would likely depend on $T_{*}$ and the extent of
supercooling, however further modelling is required.

We also compute the SNR for the planned space-based GW detector LISA [44].
LISA has a peak sensitivity at the frequency scale $f_{\text{LISA}}\sim
10^{-3}$ Hz, which is the expected scale of GW signals from a first-order
electroweak phase transition [45]. We use the projected sensitivity curve
$\Omega_{\text{LISA}}$ from Refs. [46, 47], plotted in fig. 5. We calculate
the SNR as [47]

$\text{SNR}=\sqrt{\mathcal{T}\\!\int_{0}^{\infty}\\!\\!\\!df\\!\left(\frac{\Omega_{\text{GW}}(f)}{\Omega_{\text{LISA}}(f)}\right)^{\\!\\!2}},$
(31)

where $\Omega_{\text{GW}}$ is the total GW signal from the sound wave and
turbulence sources, and assume an effective observation time $\mathcal{T}$ of
three years, coming from a mission duration of four years and 75% data taking
uptime.

##  V Model

We use the real scalar singlet model — which is a simple yet realistic
extension to the Standard Model — to realise a first-order electroweak phase
transition. Details of this model and our treatment of one-loop corrections
are available in section 4.2 of Ref. [28]. We improve the treatment by adding
extra fermions (including all quarks and the muon and tau), and adding
Boltzmann suppression factors to the Debye corrections. We also appropriately
adjust the radiation degrees of freedom to $g_{*}^{\prime}=22.25$. A similar
treatment in a simpler single-field model was used in Ref. [27].

We consider four benchmark points (BPs) in this study, each with a different
extent of supercooling. All BPs come from a narrow slice of the total
parameter space of the model. We start with M2-BP2 of Ref. [28] and vary only
the mixing angle $\theta_{m}$ to vary the extent of supercooling. The other
input parameters are fixed as $\kappa_{hhs}=-1259.83$ GeV,
$\kappa_{sss}=-272.907$ GeV, $v_{s}=663.745$ GeV and $m_{s}=351.183$ GeV. The
mixing angles and the milestone temperatures for the BPs are listed in table
1. The supercooling increases with the BP index. BP1 represents a typical
weakly supercooled phase transition with only 1 GeV difference between the
onset of bubble nucleation and percolation, and $\alpha_{\bar{\theta}}\approx
0.01$. BP2 represents a moderately supercooled phase transition with
$\alpha_{\bar{\theta}}\approx 0.05$. Both of these BPs have an exponential
nucleation rate, thus we do not calculate $T_{\Gamma}$ for them. BP3
represents a very strongly supercooled phase transition, where the physical
volume of the false vacuum only begins to decrease just below $T_{p}$. While
BP3 has a completion temperature, percolation is questionable [48, 14, 28].
The transition strength parameter is $\alpha_{\bar{\theta}}\approx 1.7$,
beyond the reach of current hydrodynamic simulations of GWs [20]. Thus, one
must be cautious when interpreting GW predictions from BP3, and indeed BP4.
BP4 has even stronger supercooling, so much so that the phase transition does
not complete. The transition strength parameter in BP4 is
$\alpha_{\bar{\theta}}\approx 177$.

| $\theta_{m}$ | $T_{c}$ | $T_{n}$ | $T_{p}$ | $T_{e}$ | $T_{f}$ | $T_{\Gamma}$ | $\log_{10}(\alpha_{\bar{\theta}})$
---|---|---|---|---|---|---|---|---
BP1 | 0.24 | 117.0 | 106.0 | 104.8 | 104.7 | 104.6 | N/A | $-1.938$
BP2 | 0.258 | 108.3 | 78.10 | 74.17 | 73.80 | 73.24 | N/A | $-1.264$
BP3 | 0.262 | 106.2 | N/A | 32.46 | 25.65 | 12.69 | 59.47 | $\;\;\,0.2178$
BP4 | 0.2623 | 106.1 | N/A | 10.09 | N/A | N/A | 59.57 | $\;\;\,2.248$
Table 1: Benchmark points and their corresponding milestone temperatures. The
mixing angle is expressed in radians, and the temperatures have units of GeV.
The transition strength parameter $\alpha_{\bar{\theta}}$ is evaluated at
$T_{p}$.

##  VI Results

###  VI.1 Dependence on the transition temperature

In this section we discuss the impact on GW predictions when varying the
transition temperature, $T_{*}$. The SNR as a function of $T_{*}$ is shown in
fig. 1 for each BP. The SNR varies by orders of magnitude over the duration of
the phase transition. However, GWs are not produced until the phase transition
is well underway, so we restrict our attention to the temperature range
$T\in[T_{f},\max(T_{n},T_{\Gamma})]$.

There are two sets of curves — solid and dashed — which have starkly different
forms in the temperature domain. The solid curves use $L_{*}=R_{\text{sep}}$
while the dashed curves use $L_{*}=\bar{R}$, with everything else in the
treatment being the same between the two sets of curves. The most immediate
difference between the two sets is that the SNR increases with $T_{*}$ when
using $R_{\text{sep}}$, and decreases with $T_{*}$ when using $\bar{R}$. In
fig. 2(a,b) we see that the peak amplitude of GWs follows a similar behaviour:
the amplitude increases (decreases) with $T_{*}$ when using $R_{\text{sep}}$
($\bar{R}$). Inversely, in fig. 3(a,b) we see that the peak frequency of GWs
decreases with $T_{*}$ when using $R_{\text{sep}}$, and increases considerably
slower with $T_{*}$ when using $\bar{R}$.

These observations can be easily explained by investigating the behaviour of
$R_{\text{sep}}$ and $\bar{R}$ as function of $T_{*}$ (see fig. 4). In fact,
we find that the dominant thermal parameter when varying $T_{*}$ is $L_{*}$,
not $K$. In fig. 4(a) we plot choices of the length scale as a function of
$T_{*}$ for BP2 (intermediate supercooling). The mean bubble separation is
large near the start of the phase transition (at higher $T_{*}$) because there
are few bubbles so their separation is large. The separation decreases over
the course of the phase transition (with decreasing $T_{*}$) because new
bubbles are nucleated. The mean bubble radius, on the other hand, begins very
small because the first bubbles to nucleate have not had time to grow
significantly. As the phase transition continues, pre-existing bubbles grow,
but more small bubbles are nucleated, suppressing an increase in the mean
radius. Thus, the mean bubble radius increases over time (i.e. with decreasing
$T_{*}$) but varies less than the mean bubble separation. We also see that the
mean bubble separation estimated using $\beta$ actually emulates the mean
bubble radius. This is unsurprising, because $R_{\text{sep}}$ is supposedly
inversely proportional to $\beta$, and $\beta$ is much higher at the start of
a phase transition with the bounce action diverging at $T_{c}$. Thus,
$R_{\text{sep}}$ estimated using $\beta$ is small at high $T_{*}$, in line
with $\bar{R}$, whereas the true $R_{\text{sep}}$ is large at high $T_{*}$.

The behaviour of $R_{\text{sep}}$ in BP3 (see fig. 4(b)) is more complicated
due to strong supercooling. The expansion of space dilutes the bubble number
density and increases the separation between bubbles. Additionally, bubble
nucleation is negligible well below $T_{\Gamma}$ so new bubbles are not
nucleated to reduce the mean separation. With even stronger supercooling in
BP4 (not shown), $R_{\text{sep}}$ begins to increase rapidly as $T_{*}$ drops
below $T_{p}$. We also see that $\beta$ cannot be used to estimate
$R_{\text{sep}}$ in BP3 (at least below $T_{\Gamma}$). However, one can
instead use $\beta_{V}$ under the Gaussian nucleation rate approximation,
which is seen to reproduce both $R_{\text{sep}}$ and $\bar{R}$ quite well at
$T_{p}$ in this example.

Now that the temperature scaling of the length scales is clear, we can return
to effects on the GW signal. First, the peak frequency for all sources is
inversely proportional to $L_{*}$ and is largely unaffected by any other
thermal parameters. Only the frequency corresponding to the sound shell
thickness scale (in the sound shell model) is directly affected by the
hydrodynamic parameters $K$, $v_{w}$ and $c_{s}$. The two key frequencies in
the sound shell model are less separated with increased supercooling due to
thickening of the sound shell. Otherwise, the behaviour of the peak
frequencies in fig. 3 can be explained purely by the behaviour of the length
scales in fig. 4. If one uses $\bar{R}$, the change in frequency with $T_{*}$
is milder than when using $R_{\text{sep}}$. In general, stronger supercooling
lowers the peak frequency at $T_{p}$.

Next, the peak amplitude for all sources is proportional to $L_{*}$. However,
the amplitude also depends on $K$ and $c_{s}$, as well as $v_{w}$ indirectly
through $K$. Nevertheless, $L_{*}$ typically has a dominant effect on the
amplitude. In the absence of strong supercooling, $R_{\text{sep}}$ changes
considerably with $T_{*}$ while $\bar{R}$ does not. Yet, $K$ and the other
hydrodynamic parameters change very little, so $L_{*}$ still has a dominant
effect even when using $\bar{R}$. With strong supercooling, $K$ and the other
hydrodynamic parameters can vary considerably between $T_{\Gamma}$ and
$T_{p}$. So too can $\bar{R}$, while $R_{\text{sep}}$ remains approximately
constant, and is in fact minimised near $T_{\Gamma}$. The peak amplitude
increases rapidly below $T_{p}$ in BP4 not because of $K$ (which is roughly
unity even at $T_{p}$), but because of the expansion of space causing a rapid
increase in $L_{*}$.555This also causes a rapid decrease in the peak frequency
at low temperature, consistent with the findings in Ref. [27]. These are all
generic features of strongly supercooled phase transitions so the results and
analysis presented here should apply to other BPs and other models.

Combining the peak amplitudes and frequencies of the GW sources, one can then
compare the GW signal to the sensitivity of a GW detector to obtain the SNR.
We consider LISA in this study, but in principle the SNR at any detector could
be calculated. Although we now have a clear picture of the behaviour of the
peak amplitude and frequency, the behaviour of the SNR is complicated by the
sensitivity window of LISA. The SNR is enhanced when the peak frequency
matches the frequency range where LISA is most sensitive; that is, near
$f_{\text{LISA}}\sim 10^{-3}$ Hz. If by varying $T_{*}$ one would obtain a
higher peak amplitude but shift the peak frequency further from LISA’s optimal
frequency range, the SNR could decrease. Thus, investigating the peak
amplitude or peak frequency in isolation will not give a clear indication of
detectability.

In fig. 5 we plot the peak of the GW signal in the amplitude-frequency plane
as a function of $T_{*}$ for BP3 to provide further insight into these
competing effects. When using $R_{\text{sep}}$ for a strongly supercooled
phase transition, the peak frequency (amplitude) increases (decreases) with
decreasing $T_{*}$, until a reversal at $T_{\Gamma}$. However, between
$T_{\Gamma}$ and $T_{p}$ the amplitude increases faster than the frequency
decreases, increasing the SNR at LISA. Meanwhile, if one uses $\bar{R}$ for a
strongly supercooled phase transition, the peak frequency (amplitude)
decreases (increases) with decreasing $T_{*}$. In the example of BP3, the peak
of the GW signal slides across the boundary of LISA’s sensitivity curve,
leading to an almost constant SNR between $T_{\Gamma}$ and $T_{f}$. One can
imagine that a slightly different BP could alter the GW peak scaling slightly,
leading to a substantially different scaling of SNR with $T_{*}$. Naturally,
the curves for $R_{\text{sep}}$ and $\bar{R}$ meet near $T_{f}$ because the
two length scales are very similar near the end of the phase transition (as
was also demonstrated in Ref. [34]).

The GW signal is formed from the sound wave and turbulence contributions,
noting again that we have neglected the collision contribution. We consider
one GW fit for the turbulence throughout, but we present results for two GW
fits for sound waves: the sound shell model and lattice fits. First we compare
the two fits for the sound wave source. Based on the SNR alone (see fig. 1) we
find a significant discrepancy between the two fits at $T_{p}$ in BP1 and BP2.
The fits agree quite well for BP3 and BP4 when using $R_{\text{sep}}$ but this
is a coincidence due to LISA’s sensitivity window. Looking at the peak
amplitudes and frequencies separately for BP3 and BP4 (see fig. 2(c,d) and
fig. 3(c,d)), we see that the predicted GW signals are still different. When
using $\bar{R}$ instead, the SNR of the sound shell model is consistently
smaller in BP1 and BP2 for all $T_{*}$ because the peak frequency is always
above LISA’s optimal frequency, $f_{\text{LISA}}$. The situation is more
complicated in BP3 and BP4 because the peak frequency crosses
$f_{\text{LISA}}$ as $T_{*}$ is varied.

The ratio of peak amplitudes in the two sound wave fits is
$\Omega_{\text{sw}}^{\text{ss}}/\Omega_{\text{sw}}^{\text{lat}}\approx 0.20$
for $v_{w}\sim 1$ and $c_{s,f}\sim 1/\sqrt{3}$, where the superscripts ‘ss’
and ‘lat’ denote the sound shell and lattice fits, respectively. This ratio is
approximately independent of $T_{*}$ and is similar for all BPs. The ratio of
peak frequencies is
$f_{\text{sw}}^{\text{ss}}/f_{\text{sw}}^{\text{lat}}\approx 2.4$ for
$v_{w}\sim 1$ and $c_{s,f}\sim 1/\sqrt{3}$ as in BP3, but increases to roughly
$8.1$ in BP1 where $v_{\text{CJ}}\approx 0.65$. The ratio of peak frequencies
has a slight dependence on $T_{*}$ due to our choice $v_{w}=v_{\text{CJ}}$,
with $v_{\text{CJ}}$ implicitly depending on $T_{*}$ through $\alpha$. The
large frequency ratio in BP1 and BP2 leads to a large difference in the SNR at
LISA between the two sound wave fits. The choice $v_{w}=v_{\text{CJ}}$ results
in a large separation in length scales — $L_{*}$ and $L_{*}\Delta_{w}$ — when
$v_{\text{CJ}}\sim c_{s,f}$, which occurs when $\alpha\ll 1$. Here,
$\Delta_{w}=(v_{w}-c_{s,f})/v_{w}$ is a multiplier for the sound shell
thickness, and can be applied to either $R_{\text{sep}}$ or $\bar{R}$.

Next we compare the sound wave source to the turbulence source. In general,
$\Omega_{\text{turb}}$ decreases faster than $\Omega_{\text{sw}}$ with
decreasing $T_{*}$ when using $R_{\text{sep}}$. This is because both
amplitudes are proportional to the decreasing $L_{*}$, but
$\Omega_{\text{sw}}$ is proportional to the increasing $K^{2}$ while
$\Omega_{\text{turb}}$ is proportional to $K^{3/2}$. Thus, the fractional
contribution of turbulence to the total GW signal decreases with $T_{*}$.
However, when $K\sim 1$, as in BP4 below $T_{p}$, the scaling with $K$ is
equivalent between the two GW sources. The comparison of the two sources does
not change when instead using $\bar{R}$, although the amplitudes now
monotonically increase with decreasing $T_{*}$. There is little insight to
gain when comparing the peak frequencies of the GW sources because they
largely differ by a constant factor. The peak frequency for the turbulence
contribution is between the peak frequencies of the two sound wave fits; it is
larger than that of the lattice fit and smaller than that of the sound shell
model fit. However, because the sound shell thickens with supercooling (at
least when choosing $v_{w}=v_{\text{CJ}}$), we find that the peak frequency of
turbulence closely matches the peak frequency in the sound shell model in
strongly supercooled scenarios. Though, the GW fits were obtained in weak and
intermediate supercooling scenarios, so their use in scenarios with strong
supercooling requires extrapolation and should be interpreted with care.

Finally, one can compare the contribution to the SNR from the sound wave and
turbulence sources. This information cannot be inferred from the results shown
in fig. 1. Instead, we will discuss the turbulence contribution — and the
impact on the SNR when increasing it — in the next section, where we consider
variations of our best treatment.

Figure 1: SNR at LISA as a function of $T_{*}$. From top to bottom: (a) BP1,
(b) BP2, (c) BP3, (d) BP4. The vertical dashed lines correspond to key
temperatures: $T_{\Gamma}$ (magenta), $T_{n}$ (red), $T_{p}$ (green), $T_{e}$
(blue) and $T_{f}$ (black). Completion occurs at the left border of each plot,
except for BP4 where there is no completion. The solid curves correspond to
$L_{*}=R_{\text{sep}}$ and the dashed curves correspond to $L_{*}=\bar{R}$.

Figure 2: Peak GW amplitude as a function of transition temperature. See the
caption of fig. 1 for further details.

Figure 3: Peak GW frequency as a function of transition temperature. See the
caption of fig. 1 for further details.

Figure 4: Characteristic length scale as a function of transition temperature.
From top to bottom: (a) BP2, (b) BP3. The qualitative features of BP1 and BP4
are respectively very similar to those of BP2 and BP3, although
$R_{\text{sep}}$ and $\bar{R}$ increase rapidly near $T_{f}$ in BP4. The
vertical dashed lines correspond to key temperatures: $T_{\Gamma}$ (magenta),
$T_{n}$ (red), $T_{p}$ (green), $T_{e}$ (blue) and $T_{f}$ (black). Completion
occurs at the left border of each plot Figure 5: The peak amplitude and
frequency of the GW signal for BP3 as a function of temperature. Here we show
only the sound shell model for the sound wave source. The noise curve for LISA
is shown in blue.

###  VI.2 Variations of the treatment

We now discuss the impact of individual variations to our best treatment for
GW prediction. These variations involve estimating $R_{\text{sep}}$ using
$\beta$ and $\beta_{V}$, estimating $K$ using other hydrodynamic quantities,
and changing the efficiency coefficient for turbulence, as discussed in
section III.2. The numerical results are stated in tables 2, 3 and 4 for
BP1-3. We do not consider BP4 here because the phase transition does not
complete; besides the results should qualitatively match those of BP3. Note
that studies typically do not vary from our best treatment by one small
change. Usually many approximations are made for all thermal parameters used
in GW predictions. Our investigation here does not encompass such treatments;
instead we point the reader to Ref. [30] where they compare low and high
diligence treatments. However, one cannot easily determine from their results
the effects of individual variations to indicate whether an approximation is
appropriate.

First, we briefly discuss the impact of the varying the transition
temperature, which is otherwise treated in more detail in section VI.1. The
two main parameters affecting the GW predictions are $K$ and $L_{*}$. We see
that $K$ changes by at most a factor of a few between $T_{n}$ and $T_{f}$ even
in the strongly supercooled scenario, BP3.666Evaluating the GW signal at
$T_{f}$ (defined by $P_{f}(T_{f})=0.01$) is not a standard treatment. We show
this variation to demonstrate the limiting behaviour of quantities near the
end of the phase transition. Yet the peak amplitudes and frequencies change by
several orders of magnitude. This is because $R_{\text{sep}}$ changes by
several orders of magnitude between $T_{n}$ and $T_{f}$. Whether the SNR is
higher or lower for some choice of $T_{*}$ depends on where the peak frequency
lies with respect to LISA’s peak sensitivity, $f_{\text{LISA}}$. Because of
this, there is no consistent trend in the effect of $T_{*}$ on the SNR across
the BPs, even though there is a consistent trend in the peak amplitudes and
frequencies.

Next, we find that using $\beta(T_{p})$ to estimate $R_{\text{sep}}(T_{p})$
results in roughly a factor of two error in peak amplitudes and frequencies in
BP1 and BP2. A similar error is present when using $\beta_{V}$ to estimate
$R_{\text{sep}}(T_{p})$ in BP3. However, it is common practice to evaluate
$\beta$ at $T_{n}$ rather than at $T_{p}$, which introduces a larger error as
seen in fig. 4(a). Yet using $\beta(T_{n})$ is more appropriate than using
$R_{\text{sep}}(T_{n})$ simply because the bubble number density changes
faster than $\beta$ between $T_{n}$ and $T_{p}$. We do not consider the
variation $L_{*}=\bar{R}$ here because GW fits are derived in terms of
$R_{\text{sep}}$ rather than $\bar{R}$. An appropriate mapping would need to
be applied to use $\bar{R}$ in the fits, such as multiplying $L_{*}$ by an
unknown constant factor in the fits.

Varying the hydrodynamic quantity $x$ in eq. 23 has a significant impact on
the prediction of $K$ in BP1 and BP2. The effect is considerably smaller in
BP3. This can be understood as follows. The pressure difference $\Delta p$ and
energy density difference $\Delta\rho$ are starkly different at high
temperature, with $\Delta p=0$ and $\Delta\rho\neq 0$ at $T_{c}$. We always
have $\alpha_{p}<\alpha_{\theta}<\alpha_{\rho}$ [25]. Using the pressure
difference underestimates $K$, while using the energy density difference
overestimates $K$. Our results match the findings of Refs. [25, 26]. With
increased supercooling (i.e. at lower temperature), the energy density
approaches the pressure such that $\alpha_{p}\approx\alpha_{\rho}$, and
$c_{s,t}^{2}\approx 1/3$ such that $\bar{\theta}\approx\theta$. Thus, for
strong supercooling we find that all methods to estimate $K$ lead to similar
results, while significant discrepancies arise for weak and intermediate
supercooling.

Lastly, we consider the impact of varying the turbulence efficiency
coefficient, $\kappa_{\text{turb}}$, through variation of $\epsilon$ (see eq.
21). Increasing $\kappa_{\text{turb}}$ can have a large impact on the SNR,
particularly if the peak frequency of turbulence better matches the detector’s
sensitivity window than the peak frequency of sound waves does. The variations
$\epsilon_{3}$ and $\epsilon_{4}$ increase the amplitude of the turbulence
source by two orders of magnitude because $\epsilon$ approaches unity, and
$(1/0.05)^{3/2}\approx 90$. However, $\epsilon_{3}$ predicts zero turbulence
in BP3 because $H(T_{*})\tau_{\text{sw}}>1$. Increasing the turbulence
contribution increases the SNR significantly in BP1 when using the sound shell
model but has little effect when using the lattice fit for sound waves. The
effect is small in BP2 with up to a 50% increase in SNR when using the sound
shell model. The effect is significant in BP3 when using either sound wave
fit.

Variation |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-17})$ |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-18})$ |  $f_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-5})$ |  $f_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-4})$ |  $h^{2}\Omega_{\mathrm{turb}}$ $(\\!\times\\!10^{-20})$ |  $f_{\mathrm{turb}}$ $(\\!\times\\!10^{-5})$ |  $\mathrm{SNR_{lat}}$ $(\\!\times\\!10^{-5})$ |  $\mathrm{SNR_{ss}}$ $(\\!\times\\!10^{-7})$ |  $\alpha$ $(\\!\times\\!10^{-3})$ |  $\kappa$ $(\\!\times\\!10^{-2})$ |  $K$ $(\\!\times\\!10^{-4})$
---|---|---|---|---|---|---|---|---|---|---|---
None | 22.57 | 31.49 | 1422 | 1157 | 21.28 | 3150 | 156.2 | 39.60 | 11.52 | 9.900 | 11.20
$T_{*}=T_{e}$ | 13.97 | 19.50 | 1833 | 1490 | 12.90 | 4061 | 56.44 | 11.24 | 11.57 | 9.921 | 11.27
$T_{*}=T_{f}$ | 11.10 | 15.50 | 2080 | 1685 | 10.16 | 4607 | 33.82 | 6.105 | 11.66 | 9.955 | 11.39
$T_{*}=T_{n}$ | 147000 | 204300 | 2.611 | 2.187 | 5448000 | 5.785 | 10230 | 5026000 | 10.74 | 9.565 | 10.09
$R_{\text{sep}}(\beta)$ | 11.04 | 15.41 | 2062 | 1678 | 10.12 | 4567 | 34.32 | 6.216 |  |  |
$K(\alpha(\theta))$ | 21.09 | 29.44 |  |  | 19.92 |  | 146.0 | 37.03 | 11.46 | 9.466 | 10.72
$K(\alpha(p))$ | 1.403 | 1.957 |  |  | 1.489 |  | 9.711 | 2.509 | 3.590 | 5.317 | 1.902
$K(\alpha(\rho))$ | 261.9 | 365.5 |  |  | 234.7 |  | 1813 | 456.2 | 35.05 | 16.39 | 55.50
$\epsilon_{2}$ |  |  |  |  | 60.18 |  | 156.4 | 54.06 |  |  |
$\epsilon_{3}$ |  |  |  |  | 1776 |  | 166.0 | 1035 |  |  |
$\epsilon_{4}$ |  |  |  |  | 1787 |  | 166.0 | 1041 |  |  |

Table 2: GW predictions and hydrodynamic parameters for BP1. Each row
corresponds to a different variation of our best treatment. Blank cells match
the result of our best treatment (i.e. the top row). Frequencies are stated in
units of GeV, with all other quantities being dimensionless. The scripts ‘ss’
and ‘lat’ respectively denote the sound shell model fit and the lattice fit
for the sound wave source of GWs.

Variation |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-13})$ |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-14})$ |  $f_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-5})$ |  $f_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-4})$ |  $h^{2}\Omega_{\mathrm{turb}}$ $(\\!\times\\!10^{-16})$ |  $f_{\mathrm{turb}}$ $(\\!\times\\!10^{-5})$ | $\mathrm{SNR_{lat}}$ | $\mathrm{SNR_{ss}}$ |  $\alpha$ $(\\!\times\\!10^{-2})$ | $\kappa$ |  $K$ $(\\!\times\\!10^{-3})$
---|---|---|---|---|---|---|---|---|---|---|---
None | 3.590 | 5.673 | 129.6 | 60.20 | 3.898 | 287.0 | 10.08 | 2.031 | 5.450 | 0.2074 | 10.64
$T_{*}=T_{e}$ | 2.552 | 4.042 | 159.9 | 73.75 | 2.662 | 354.2 | 8.763 | 1.204 | 5.575 | 0.2096 | 10.99
$T_{*}=T_{f}$ | 2.146 | 3.410 | 181.7 | 82.91 | 2.187 | 402.5 | 8.110 | 0.8892 | 5.771 | 0.2129 | 11.54
$T_{*}=T_{n}$ | 676.5 | 1046 | 2.189 | 1.098 | 8968 | 4.849 | 1.310 | 5.142 | 4.297 | 0.1857 | 7.597
$R_{\text{sep}}(\beta)$ | 2.019 | 3.191 | 177.5 | 82.45 | 2.078 | 393.1 | 7.510 | 0.8449 |  |  |
$K(\alpha(\theta))$ | 3.372 | 5.329 |  |  | 3.676 |  | 9.469 | 1.908 | 5.362 | 0.2011 | 10.23
$K(\alpha(p))$ | 1.428 | 2.256 |  |  | 1.649 |  | 4.010 | 0.8081 | 3.698 | 0.1682 | 5.997
$K(\alpha(\rho))$ | 14.45 | 22.84 |  |  | 14.61 |  | 40.59 | 8.172 | 10.35 | 0.2736 | 25.68
$\epsilon_{2}$ |  |  |  |  | 11.03 |  | 10.11 | 2.064 |  |  |
$\epsilon_{3}$ |  |  |  |  | 290.2 |  | 11.21 | 3.406 |  |  |
$\epsilon_{4}$ |  |  |  |  | 301.7 |  | 11.26 | 3.462 |  |  |

Table 3: The same as table 2 but for BP2.

Variation |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-7})$ |  $h^{2}\Omega_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-8})$ |  $f_{\mathrm{sw}}^{\mathrm{lat}}$ $(\\!\times\\!10^{-6})$ |  $f_{\mathrm{sw}}^{\mathrm{ss}}$ $(\\!\times\\!10^{-6})$ |  $h^{2}\Omega_{\mathrm{turb}}$ $(\\!\times\\!10^{-10})$ |  $f_{\mathrm{turb}}$ $(\\!\times\\!10^{-6})$ | $\mathrm{SNR_{lat}}$ | $\mathrm{SNR_{ss}}$ | $\alpha$ | $\kappa$ | $K$
---|---|---|---|---|---|---|---|---|---|---|---
None | 1.861 | 3.748 | 9.345 | 23.48 | 6.348 | 20.70 | 249.6 | 307.7 | 1.651 | 0.7175 | 0.4536
$T_{*}=T_{e}$ | 4.318 | 8.872 | 7.908 | 19.12 | 14.74 | 17.52 | 443.7 | 498.2 | 4.257 | 0.8422 | 0.6950
$T_{*}=T_{f}$ | 17.04 | 35.42 | 4.111 | 9.722 | 81.84 | 9.106 | 864.5 | 876.4 | 71.06 | 0.9831 | 0.9803
$R_{\text{sep}}(\beta_{V})$ | 1.193 | 2.402 | 12.80 | 32.17 | 3.394 | 28.36 | 222.6 | 356.9 |  |  |
$K(\alpha(\theta))$ | 1.819 | 3.663 |  |  | 6.227 |  | 244.9 | 301.5 | 1.605 | 0.7269 | 0.4478
$K(\alpha(p))$ | 1.768 | 3.560 |  |  | 6.083 |  | 239.2 | 294.2 | 1.564 | 0.7269 | 0.4409
$K(\alpha(\rho))$ | 1.967 | 3.962 |  |  | 6.646 |  | 261.4 | 323.0 | 1.728 | 0.7383 | 0.4677
$\epsilon_{2}$ |  |  |  |  | 17.95 |  | 700.0 | 742.2 |  |  |
$\epsilon_{3}$ |  |  |  |  | 0 |  | 18.36 | 130.9 |  |  |
$\epsilon_{4}$ |  |  |  |  | 288.4 |  | 11210 | 11230 |  |  |

Table 4: The same as table 2 but for BP3. There is no row for $T_{*}=T_{n}$
because there is no nucleation temperature for BP3. This time there is a row
for $R_{\text{sep}}(\beta_{V})$ instead of $R_{\text{sep}}(\beta)$ because the
bubble nucleation rate is Gaussian rather than exponential. In fact, $\beta$
is negative and leads to invalid predictions.

##  VII Discussion

In this study we have investigated several ambiguities and approximations made
in predictions of GWs from cosmological phase transitions. We considered each
approximation in isolation to provide a clear indication of their individual
effects on the GW signal. We recommend our results be used in conjunction with
the results of Ref. [30] to determine whether a particular set of
approximations can lead to reliable GW predictions. Alternatively, one could
use our best treatment described in section III.2 if feasible, and even
improve on it with a proper treatment of hydrodynamic profile around bubble
walls and a method for estimating friction on the bubble wall.

To our knowledge, our investigation is the first to explicitly determine the
effect of varying the transition temperature, $T_{*}$. We note that our
investigation is fundamentally different from studies that vary thermal
parameters (including $T_{*}$) separately, treating them as independent
quantities. We account for the implicit interdependence of all thermal
parameters.

The correct choice of the transition temperature is still unknown because the
hydrodynamic simulations from which GW fits are obtained hold the temperature
fixed. In fact, evaluating GW predictions at a single temperature may fall out
of favour once modelling of GW production is improved further. We have
demonstrated that using the current set of thermal parameters (in particular
$R_{\text{sep}}$), the GW signal can change by several orders of magnitude
between commonly chosen transition temperatures: $T_{n}$ and $T_{p}$. If a
more appropriate choice of transition temperature lies somewhere between
$T_{n}$ and $T_{p}$, then new GW predictions would significantly differ from
those obtained using the current best treatments which use $T_{*}=T_{p}$.

We argued in section VI.1 that evaluating the GW signal at temperatures above
$T_{n}$ is not meaningful because bubble collisions would not have occurred to
source GWs at that stage in the phase transition. This same reasoning can also
be used to discard $T_{n}$ as a reasonable transition temperature. The only
case where the nucleation temperature reflects a time when collisions are
occurring is in some strongly supercooled phase transitions — where in extreme
cases $T_{n}\sim T_{p}$, counter-intuitively [28]. However, using $T_{n}$ in
strongly supercooled phase transitions is not recommended. For one, it
decouples from the progress of the phase transition, so it does not represent
a consistent stage in the phase transition. Further, the existence of a
nucleation temperature does not indicate whether a phase transition occurs or
completes, as discussed in Ref. [28]. Thus, one must be careful when using
$T_{n}$, and ensure that the phase transition is in fact weakly supercooled.

It is commonly assumed that the GW signal should be similar at $T_{n}$ and
$T_{p}$ for weakly supercooled phase transitions. This is not consistent with
our findings. Calculating the mean bubble separation properly (from the bubble
number density) would suggest orders of magnitude difference in the GW signal
between $T_{n}$ and $T_{p}$. Using the mean bubble radius or $\beta$ instead
still suggests a factor of a few difference in the GW signal between $T_{n}$
and $T_{p}$. The hydrodynamic parameters like the kinetic energy fraction,
however, are similar at the two temperatures.

The mean bubble radius varies much slower with temperature than the mean
bubble separation. Thus, studies evaluating GWs at $T_{n}$ should use the mean
bubble radius or $\beta$ instead of calculating the mean bubble separation
directly from the bubble number density. However, we note that if one could
calculate the bubble number density, then one could calculate $T_{p}$ and use
the recommended treatment outlined in section III.2.

In general, we find that variations of the treatment of GW predictions can
lead to sizeable deviations in the SNR and peak amplitudes and frequencies;
potentially deviations of many orders of magnitude. In the context of GW
predictions from cosmological phase transitions, mild deviation is of the
order of 10%, suggesting that constraints on particle physics models from GW
observations will be hard to apply reliably at this stage. Nevertheless, the
recent emergence of successful GW astronomy offers hope for constraining
particle physics models at scales beyond the reach of particle physics
experiments.

##  VIII Acknowledgements

LM thanks Thomas Konstandin for assistance with numerical accuracy in
calculating $\kappa_{\bar{\theta}}$. LM was supported by an Australian
Government Research Training Program (RTP) Scholarship and a Monash Graduate
Excellence Scholarship (MGES). The work of PA is supported by the National
Natural Science Foundation of China (NNSFC) under grant No. 12150610460 and by
the supporting fund for foreign experts grant wgxz2022021L. ZX is also
supported in part by NNSFC grant No. 12150610460.

##  Appendix A Correction to the kinetic energy fraction parameterisation

The kinetic energy fraction is often parameterised as

$K=\frac{\kappa\alpha}{1+\alpha}.$ (32)

This parameterisation introduces approximation to the fundamental definition
[22, 25, 9]

$K=\frac{\rho_{\text{kin}}(T_{*})}{\rho_{\text{tot}}(T_{*})},$ (33)

where $\rho_{\text{kin}}$ is the fluid kinetic energy. In the following we
assume $\rho$ and $p$ are renormalised such that the ground state energy
density vanishes. In this case, $\rho_{\text{tot}}=\rho_{f}$.

The inexact nature of eq. 32 was demonstrated in appendix B.2 of Ref. [22] and
implied in Ref. [25] (seen by comparing methods M2 and M3). A correction
$\delta$ can be applied such that [22]

$K=\frac{\kappa\alpha}{1+\alpha+\delta}.$ (34)

One can solve for $\delta$ by equating eq. 33 and eq. 34. If $\alpha$ is
calculated using the trace anomaly

$\theta=\frac{1}{4}\\!\left(\rho-3p\right)$ (35)

as in Ref. [22], one finds

$\delta=\frac{\theta_{t}}{3w_{f}}.$ (36)

If $\alpha$ is calculated using the pseudotrace [25]

$\bar{\theta}=\frac{1}{4}\\!\left(\rho-\frac{p}{c_{s,t}^{2}}\right),$ (37)

which reduces to the trace anomaly if $c_{s,t}^{2}=1/3$ (e.g. as in the bag
model), one instead finds

$\delta=\frac{4}{3w_{f}}\\!\left(\rho_{\text{tot}}-\Delta\bar{\theta}\right)-1.$
(38)

In our benchmark points we find $\delta\ll 1+\alpha$ such that the difference
between eq. 32 and eq. 34 is at most 1%. Thus, we do not include such
variations on the treatment of $K$ in our results.

## References

  * [1] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8 [2306.16213].
  * [2] H. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Res. Astron. Astrophys._ 23 (2023) 075024 [2306.16216].
  * [3] J. Antoniadis et al., _The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals_ , 2306.16214.
  * [4] D.J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6 [2306.16215].
  * [5] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Search for Signals from New Physics_ , _Astrophys. J. Lett._ 951 (2023) L11 [2306.16219].
  * [6] KAGRA, VIRGO, LIGO Scientific collaboration, _Open Data from the Third Observing Run of LIGO, Virgo, KAGRA, and GEO_ , _Astrophys. J. Suppl._ 267 (2023) 29 [2302.03676].
  * [7] P. Athron, C. Balázs, T.E. Gonzalo and M. Pearce, _Falsifying Pati-Salam models with LIGO_ , 2307.02544.
  * [8] F. Huang, V. Sanz, J. Shu and X. Xue, _LIGO as a probe of dark sectors_ , _Phys. Rev. D_ 104 (2021) 095001 [2102.03155].
  * [9] P. Athron, C. Balázs, A. Fowlie, L. Morris and L. Wu, _Cosmological phase transitions: from perturbative particle physics to gravitational waves_ , 2305.02357.
  * [10] P. Athron, C. Balazs, A. Fowlie, L. Morris, G. White and Y. Zhang, _How arbitrary are perturbative calculations of the electroweak phase transition?_ , _JHEP_ 01 (2023) 050 [2208.01319].
  * [11] D. Croon, O. Gould, P. Schicho, T.V.I. Tenkanen and G. White, _Theoretical uncertainties for cosmological first-order phase transitions_ , _JHEP_ 04 (2021) 055 [2009.10080].
  * [12] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [13] 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 [1304.2433].
  * [14] 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 [1809.08242].
  * [15] 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 [2003.07360].
  * [16] J. Ellis, M. Lewicki, J.M. No and V. Vaskonen, _Gravitational wave energy budget in strongly supercooled phase transitions_ , _JCAP_ 06 (2019) 024 [1903.09642].
  * [17] C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024 [1910.13125].
  * [18] H.-K. Guo, K. Sinha, D. Vagie and G. White, _Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories_ , _JCAP_ 01 (2021) 001 [2007.08537].
  * [19] 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) 103520 [1704.05871].
  * [20] D. Cutting, M. Hindmarsh and D.J. Weir, _Vorticity, kinetic energy, and suppressed gravitational wave production in strong first order phase transitions_ , _Phys. Rev. Lett._ 125 (2020) 021302 [1906.00480].
  * [21] M. Hindmarsh, _Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe_ , _Phys. Rev. Lett._ 120 (2018) 071301 [1608.04735].
  * [22] M. Hindmarsh and M. Hijazi, _Gravitational waves from first order cosmological phase transitions in the Sound Shell Model_ , _JCAP_ 12 (2019) 062 [1909.10040].
  * [23] X. Wang, F.P. Huang and Y. Li, _Sound velocity effects on the phase transition gravitational wave spectrum in the sound shell model_ , _Phys. Rev. D_ 105 (2022) 103513 [2112.14650].
  * [24] R.-G. Cai, S.-J. Wang and Z.-Y. Yuwen, _Hydrodynamic sound shell model_ , 2305.00074.
  * [25] F. Giese, T. Konstandin and J. van de Vis, _Model-independent energy budget of cosmological first-order phase transitions—A sound argument to go beyond the bag model_ , _JCAP_ 07 (2020) 057 [2004.06995].
  * [26] F. Giese, T. Konstandin, K. Schmitz and J. Van De Vis, _Model-independent energy budget for LISA_ , _JCAP_ 01 (2021) 072 [2010.09744].
  * [27] P. Athron, A. Fowlie, C.-T. Lu, L. Morris, L. Wu, Y. Wu et al., _Can supercooled phase transitions explain the gravitational wave background observed by pulsar timing arrays?_ , 2306.17239.
  * [28] P. Athron, C. Balázs and L. Morris, _Supercool subtleties of cosmological phase transitions_ , _JCAP_ 03 (2023) 006 [2212.07559].
  * [29] X. Wang, F.P. Huang and X. Zhang, _Phase transition dynamics and gravitational wave spectra of strong first-order phase transition in supercooled universe_ , _JCAP_ 05 (2020) 045 [2003.08892].
  * [30] H.-K. Guo, K. Sinha, D. Vagie and G. White, _The benefits of diligence: how precise are predicted gravitational wave spectra in models with phase transitions?_ , _JHEP_ 06 (2021) 164 [2103.06933].
  * [31] A.D. Linde, _Decay of the False Vacuum at Finite Temperature_ , _Nucl. Phys. B_ 216 (1983) 421.
  * [32] C.L. Wainwright, _CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields_ , _Comput. Phys. Commun._ 183 (2012) 2006 [1109.4189].
  * [33] K. Enqvist, J. Ignatius, K. Kajantie and K. Rummukainen, _Nucleation and bubble growth in a first order cosmological electroweak phase transition_ , _Phys. Rev. D_ 45 (1992) 3415.
  * [34] A. Mégevand and S. Ramírez, _Bubble nucleation and growth in very strong cosmological phase transitions_ , _Nucl. Phys. B_ 919 (2017) 74 [1611.05853].
  * [35] R.-G. Cai, M. Sasaki and S.-J. Wang, _The gravitational waves from the first-order phase transition with a dimension-six operator_ , _JCAP_ 08 (2017) 004 [1707.03001].
  * [36] P. Athron, C. Balázs, A. Fowlie and Y. Zhang, _PhaseTracer: tracing cosmological phases and calculating transition properties_ , _Eur. Phys. J. C_ 80 (2020) 567 [2003.02859].
  * [37] P. Athron, C. Balázs and L. Morris, _TransitionSolver: resolving cosmological phase histories_ , in preparation (2023).
  * [38] J.R. Espinosa, T. Konstandin, J.M. No and G. Servant, _Energy Budget of Cosmological First-order Phase Transitions_ , _JCAP_ 06 (2010) 028 [1004.4187].
  * [39] T. Alanne, T. Hugle, M. Platscher and K. Schmitz, _A fresh look at the gravitational-wave signal from cosmological phase transitions_ , _JHEP_ 03 (2020) 004 [1909.11356].
  * [40] A.F. Heckler, _The Effects of electroweak phase transition dynamics on baryogenesis and primordial nucleosynthesis_ , _Phys. Rev. D_ 51 (1995) 405 [astro-ph/9407064].
  * [41] A. Mégevand and S. Ramírez, _Bubble nucleation and growth in slow cosmological phase transitions_ , _Nucl. Phys. B_ 928 (2018) 38 [1710.06279].
  * [42] M.A. Ajmi and M. Hindmarsh, _Thermal suppression of bubble nucleation at first-order phase transitions in the early Universe_ , _Phys. Rev. D_ 106 (2022) 023505 [2205.04097].
  * [43] C. Caprini, R. Durrer and G. Servant, _The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition_ , _JCAP_ 12 (2009) 024 [0909.0622].
  * [44] LISA collaboration, _Laser Interferometer Space Antenna_ , 1702.00786.
  * [45] C. Caprini et al., _Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions_ , _JCAP_ 04 (2016) 001 [1512.06239].
  * [46] T. Robson, N.J. Cornish and C. Liu, _The construction and use of LISA sensitivity curves_ , _Class. Quant. Grav._ 36 (2019) 105011 [1803.01944].
  * [47] C. Caprini, D.G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni et al., _Reconstructing the spectral shape of a stochastic gravitational wave background with LISA_ , _JCAP_ 11 (2019) 017 [1906.09244].
  * [48] M.S. Turner, E.J. Weinberg and L.M. Widrow, _Bubble nucleation in first order inflation and other cosmological phase transitions_ , _Phys. Rev. D_ 46 (1992) 2384.