# A Precise Fitting Formula for Gravitational Wave Spectra from Phase
Transitions

Huai-ke Guo  Fazlollah Hajkarim  Kuver Sinha
Graham White  Yang Xiao

###### Abstract

Obtaining a precise form for the predicted gravitational wave (GW) spectrum
from a phase transition is a topic of great relevance for beyond Standard
Model (BSM) physicists. Currently, the most sophisticated semi-analytic
framework for estimating the dominant contribution to the spectrum is the
sound shell model; however, full calculations within this framework can be
computationally expensive, especially for large-scale scans. The community
therefore generally manages with fit functions to the GW spectrum, the most
widely used of which is a single broken power law. We provide a more precise
fit function based on the sound shell model: our fit function features a
double broken power law with two frequency breaks corresponding to the two
characteristic length scales of the problem – inter-bubble spacing and
thickness of sound shells, the second of which is neglected in the single
broken power law fit. Compared to previously proposed fits, we demonstrate
that our fit function more faithfully captures the GW spectrum coming from a
full calculation of the sound shell model, over most of the space of the
thermodynamic parameters governing the phase transition. The physical origins
of the fit parameters and their dependence on the thermodynamic parameters are
studied in the underlying sound shell model: in particular, we perform a
series of detailed scans for these quantities over the plane of the strength
of the phase transition ($\alpha$) and the bubble wall velocity ($v_{w}$).
Wherever possible, we comment on the physical interpretations of these scans.
The result of our study can be used to generate accurate GW spectra with our
fit function, given initial inputs of $\alpha$, $v_{w}$, $\beta/H$ (nucleation
rate parameter) and $T_{n}$ (nucleation temperature) for the relevant BSM
scenario.

##  1 Introduction

One of the most plausible cosmological mechanisms for producing a stochastic
gravitational wave (GW) background is a first order phase transition [1] . The
most common motivation for such an epoch is electroweak symmetry breaking,
which can become strongly first order in extensions to the Standard Model [2,
3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61,
62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
100, 101]. Other motivations for a cosmic first order phase transition include
dark matter [102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 112,
114, 115, 88, 116, 117, 118, 119, 120, 121, 35], symmetry breaking chains from
a grand unified group [122, 123, 124, 125, 126, 25, 127], or the QCD
transition [128, 129, 130, 131, 132, 133, 134, 135, 118, 136, 137, 117, 138,
119, 116, 139]. The GW spectrum resulting from a cosmological phase transition
has three contributions - a term arising from the collision of scalar shells,
a term arising from crashing sound waves in the plasma and the after party of
turbulence. For most of the parameter space, it is widely accepted that the
sound shell source produces the dominant contribution to the stochastic GW
background [1].

At the present time, the most widely used semi-analytical framework for
calculating the velocity spectrum sourced by sound waves is provided by the
sound shell model [140]. The corresponding GW spectrum was calculated for a
non-expanding Universe in [141], and subsequently extended to an expanding
Universe in [142]. Alongside the semi-analytical framework, simulations have
recently also made great progress in this area [143, 144, 145] and motivate
the use of a broken power law for the GW spectrum, which has been adopted by
the LISA cosmology working group [1]. The single broken power law comes from
the simplifying assumption of ignoring the the details of the velocity profile
of the sound waves, and instead taking the RMS 111Root mean square. of the
fluid velocity, $U_{f}$; the peak frequency is then set by the mean bubble
separation [146]. On the other hand, recent work with the sound shell model
[147, 148] motivates using a double broken power law, which comes from
including the full velocity profile. The use of the full velocity profile
results in qualitatively different physics: the spectrum is better described
by a double broken power law since the velocity profile now contains
information about the thickness of the profile, which implies that there is
more than one physically significant scale.

It should be noted that the dependence of the spectral features of a fit on
the underlying thermodynamic parameters of the phase transition (the
nucleation temperature $T_{n}$, nucleation rate parameter $\beta$, phase
transition strength $\alpha$, and wall speed $v_{w}$) is a difficult question
and progress can be made by a combination of analytical calculations and
numerical fits. In principle, one could imagine the possibility of
reconstructing the values of four thermal parameters from a hypothetical
spectrum [147], (though a more conservative study argued for three [149]).
Such a wealth of information, if it were to be available, would render GW
detectors almost competitive with colliders in the extent of details of the
underlying model they can reveal [79]. In a sense, the community is actually
already in the era where the degree of precision we use in calculating GW
spectra becomes important, since LIGO currently puts constraints on the Pati-
Salam unification scenario [150]. This urgency will grow as next generation
detectors cover many more decades in frequency at a much high strain
sensitivity (for a review see [151]). From a beyond Standard Model (BSM)
physicist’s point of view, the ideal workflow would be the following: for
every point on the parameter space of the model, obtain the thermodynamic
parameters $(T_{n},\,\beta,\,\alpha,\,v_{w})$, and compute the GW spectrum
using the sound shell model (or better still, use actual simulations); then
check how the parameter space is constrained by GW detectors. Needless to say,
full calculations of the sound shell model, not to speak of actual
simulations, are difficult and expensive; to perform large-scale scans over
the parameter space of particle physics models and obtain the predicted GW
spectrum, portable and simple fits to the spectrum, and the interpretation of
the fit parameters, become necessary.

The main purpose of this paper is to provide such a precise fit function.
Elaborating further, our goal in the current paper will be the following:
given a particle physics model with input values of
$(T_{n},\,\beta,\,\alpha,\,v_{w})$, an appropriate fit function for the GW
spectrum is provided, along with a public code. A historical progression of
fit functions can help contextualize our results: the cornerstone of GW
phenomenology has been a broken power law with a peak set by the mean bubble
separation; this was recently improved upon by the double broken power law fit
function introduced in [147, 148] with fit parameters corresponding to the
peak amplitude, two frequency breaks and the spectral slope between them; our
fit function is a further improvement, by introducing an extra parameter that
governs the IR behavior of the spectrum. Scanning over $(\alpha,v_{w})$, we
show that our fit function reproduces, much more faithfully, the GW spectrum
coming from the sound shell model (Fig. 2). We provide data files and scripts
in Python and Mathematica that can be directly used by a front-end user to
generate accurate GW spectra, given inputs of
$(T_{n},\,\beta,\,\alpha,\,v_{w})$ for a given beyond-SM scenario. This tool
is expected to allow the community to fully investigate the extent to which
each individual BSM model can be constrained by future detectors, without
actually performing calculations with the sound shell model.

The second purpose of our paper is to attempt to study the physical
interpretation of the fit parameters by way of detailed scans on the plane of
$(\alpha,v_{w})$. We provide scans over $(\alpha,v_{w})$ in the deflagration,
hybrid, and detonation regimes and study the velocity and enthalpy profiles,
as well as a host of quantities: the maximum velocity $v_{max}$, the thermal
efficiency factor $\kappa$, the peak frequency $f_{p}$, the maximum relic
density $\Omega_{p}$, the width of the velocity profile $\delta_{\xi}$, and
the fit parameters $\tilde{b}$ and $\tilde{a}$. The dependence of these
quantities on $(\alpha,v_{w})$ can be highly non-linear. Wherever possible, we
make qualitative statements about this dependence, as well as correlations
with the relative behaviour of other parameters.

This paper is organized as follows. In Section 2, we discuss the general
framework for the calculation of GWs in a first order phase transition and
show, in turn, the standard broken power law approximation (Section 2.1) and
the improved double broken power law fit (Section 2.2). In Section 3 we
introduce our new fit function and discuss its faithfulness in capturing the
GW spectrum compared to the older fits. The dependence of various physical and
fit parameters on $(\alpha,v_{w})$ is discussed in Section 4, where we also
display the results of our scans. Discussions and code for using the new fit
function 222It can be downloaded from this GitHub link. are relegated to the
Appendix A and B. Finally we summarize and conclude in Section 5.

##  2 General Framework for First Order Phase Transition

In this section we briefly review the framework for the calculation of GWs
from sound waves during a first order phase transition (FOPT) [152, 153, 154,
155, 156, 157, 141, 158, 142, 147, 140]. The details of a gravitational wave
spectrum produced by sound waves depend upon a number of macroscopic features
of the phase transition. For instance, it matters how far apart the bubbles
are during the phase transition, how fast the bubbles expand, the fraction of
energy released during the phase transition that is dumped into kinetic energy
modes, the details of the fluid veloicty and of course the temperature at
which all this occurs.

Since the background is expanding, the key parameter for describing the
nucleation of bubbles during the phase transition is not the nucleation rate -
as the phase transition always begins when the nucleation rate is fast enough
so that the nucleation of bubbles can keep pace with the expanding Universe.
Instead, it is customary to consider how fast the nucleation rate increases
from this bare minimum rate

$\Gamma\sim e^{-S(t)}\sim e^{-S_{0}+\beta(t-t_{n})}$ (2.1)

where $t_{n}$ is the time of nulceation at $S_{0}$ is the value needed to have
at least one bubble per hubble volume . Unless there is a significant amount
of supercooling [150], it is $\beta$ that actually controls the time scale of
the transition. It is straightforward to relate this parameter to the mean
bubble separation [141, 147] , $R_{\ast}$

$\displaystyle\frac{\beta}{H_{n}}=\left(8\pi\right)^{\frac{1}{3}}\frac{v_{w}}{r_{*}}\,.$
(2.2)

Here, the wall velocity is denoted by $v_{w}$ and

$\displaystyle r_{*}=H_{n}R_{*}\,,$ (2.3)

where $H_{n}$ is the Hubble rate at the nucleation time. The fraction of the
kinetic energy dumped into the plasma, as well as the details of the velocity
profile tend to depend upon how strong the transition is. This is governed by
the change in the trace anomaly, $\Delta\theta$ in the energy momentum tensor
between the phases, normalized to the enthalpy of the symmetric phase,
$w_{s}$,

$\alpha=\left.\frac{4}{3}\frac{\Delta\theta}{w_{\mathrm{s}}}\right|_{T=T_{\mathrm{n}}}\,.$
(2.4)

Making predictions in the sound shell model require solving the hydrodynamic
equations and using the macroscopic thermal parameters of a phase transition
as boundary conditions. We use a bag model to match solutions of thermodynamic
variables at the boundary between the two phases. We consider a non-vanishing
bag constant, which we denote as $\epsilon$, [159] which has the value

$\displaystyle\alpha=\frac{4}{3}\frac{\epsilon}{w_{n}}\,.$ (2.5)

There are three qualitatively distinct types of phase transitions which are
delineated by the value of the bubble wall velocity, compared to the Jouget
velocity, $v_{j}$ and the speed of sound, $c_{s}$ [146, 160, 141].
Specifically, the Jouget velocity is defined as

$c_{J}=c_{\text{s}}\frac{\left(1+\sqrt{\alpha(2+3\alpha)}\right)}{\left(1+\alpha\right)}\,,$
(2.6)

where the detonation regime occurs when $v_{w}>c_{J}>c_{s}$ , a deflagration
of subsonic regime is when $v_{w}<c_{J}<c_{s}$. A hybrid transition involves
both deflagration and detonation phases. Initially, $v_{w}<c_{s}$ and
$v_{w}<c_{J}$ indicating a deflagration. However, as the phase transition
evolves, energy accumulates in the bubble wall or in a shock front preceding
the wall, potentially allowing the bubble to accelerate to a velocity $v_{w}$
such that $c_{s}<v_{w}<c_{J}$ or even $v_{w}>c_{J}$, at which point a
detonation may occur [146, 160, 141].

The relative density of GW for each value of $z=kR_{*}$ corresponding to a
specific frequency $k=2\pi f$ is controlled by the kinetic energy fraction
$K$, as well as [141, 147] the lifetime of the soundshell source,
$H_{n}\tau_{v}$ and its characteristic scale, $H_{n}R_{\ast}$, for which we
use the shorthand

$J=H_{n}R_{*}H_{n}\tau_{v}=r_{*}\left(1-\frac{1}{\sqrt{1+2x}}\right)\equiv
r_{*}\Upsilon\,,$ (2.7)

where $x=H_{n}R_{*}/\sqrt{K}$ and $r_{*}$ defined in Eq. (2.3). We can then
write the gravitaitonal wave power spectrum at the time it is produced as,

$\Omega_{\mathrm{GW}}(z)=3K^{2}(v_{w},\alpha)J\frac{z^{3}}{2\pi^{2}}\tilde{P}_{\mathrm{GW}}\left(z\right).$
(2.8)

If one takes the RMS fluid velocity, rather than calculating the gravitational
wave spectrum from the full velocity profile, one can write the frequency
dependent analytical form of spectral density $z^{3}\tilde{P}_{\mathrm{GW}}$
for GW from Refs. [161, 146] as a broken power law

$\displaystyle S_{\rm SW}(f)=\left(\frac{f}{f_{\rm
SW}}\right)^{3}\left[\frac{7}{4+3(f/f_{\rm SW})^{2}}\right]^{7/2}.$ (2.9)

We will be comparing the predictions of the soundshell model to this single
broken power law throughout, as this simplified case is ubiquitous in the
literature.

In such a case we take the kinetic energy fraction as
$K=\Gamma\bar{U_{f}}^{2}$ where $\Gamma=\bar{w}/\bar{v}$ is the ratio of the
average enthalpy to the average fluid velocity. Further, the mean squared
velocity defines [162, 141]

$\overline{U}_{f}^{2}=\frac{3}{\bar{w}v_{w}^{3}}\int
w\xi^{2}\frac{v^{2}}{1-v^{2}}d\xi$ (2.10)

that can be approximated
$\overline{U}_{f}\approx\sqrt{\frac{3}{4}\frac{\kappa\alpha}{1+\alpha}}$ and
the quantity $\kappa$ is given by [141]

$\displaystyle\kappa=\frac{3}{\epsilon w^{3}}\int
d\xi\,\xi^{2}w^{2}\gamma^{2}\nu^{2}\ .$ (2.11)

where we have defined

$\xi\,\equiv\,r/t\,\,,$ (2.12)

with $r$ being the distance from the center of the bubble and $t$ being the
time from the onset of the phase transition and nucleation.

The sound shell model differs from the results of numerical simulations as
some sound shells collide before they reach a self similar solution, that is a
solution that depends upon the ratio $\xi=r/t$ only. The effect of this was
numerically approximated in Ref. [147] by introducing an error factor

$\Omega_{\rm
GW}(z)=\Omega_{\mathrm{GW}}^{\mathrm{SSM}}(z)\Sigma(v_{w},\alpha),$ (2.13)

where $\Omega_{\mathrm{GW}}^{\mathrm{SSM}}(z)$ is the relic from sound wave
model at redshift $z$ and $\Sigma(v_{w},\alpha)$ compensates from the sound
shell model’s overestimation of the gravitational waves due to energy lost
into vorticity modes as demonstrated in Ref. [163]. As the precise nature of
this factor for all possible thermal parameters is not known, we assume it as
$\Sigma(v_{w},\alpha)=1$ in order to focus on other aspects of sound shell
model that affect the GW signal. The reader should note that it is a simple
fix to include the results of Ref. [163] if one is working in a regime where
interpolations of their result can be used.

After formation the power spectra will redshift like radiation leading to its
dilution and a shift in the peak frequency such that,

$\Omega^{0}_{\mathrm{GW}}(f)=F_{\mathrm{GW}}^{0}\Omega_{\text{GW}}(z(f)),$
(2.14)

where the $F_{\mathrm{GW}}^{0}$ factor that redshifts the GW relic from its
moment of production in the early Universe to today is

$F_{\mathrm{GW}}^{0}=\Omega_{\gamma,0}\left(\frac{g_{*s,0}}{g_{*s}}\right)^{\frac{4}{3}}\frac{g_{*}}{g_{*,0}}\simeq(3.57\pm
0.05)\times 10^{-5}{\bigg{(}\frac{100}{g_{*}}\bigg{)}}^{\frac{1}{3}}\,\,,$
(2.15)

assuming $g_{*}\simeq g_{*s}$ at $T\gg T_{\nu}$.The peak frequency is set by
the mean bubble separation and a redshift factor

$f=\frac{z}{r_{*}}f_{*,0},$ (2.16)

where $z=kR_{*}$ and the redshift factor is given by

$f_{*,0}=2.6\times
10^{-6}\,\textrm{Hz}\left(\frac{T_{n}}{100\,\textrm{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,.$
(2.17)

Here we use the degrees of freedom of the thermal background of SM particles
for energy density $g_{*}$ and entropy density $g_{*s}$ from Ref. [164].

Putting everything together, the relic density of GW from sound waves is given
by [142, 146]

$\Omega_{\rm GW}h^{2}\simeq 1.2\times
10^{-6}\left(\frac{100}{g_{*}}\right)^{1/3}K^{2}\Upsilon\left(\frac{H_{n}}{\beta}\right)v_{w}S_{\rm
SW}(f)\,.$ (2.18)

The low frequency tail of the GW spectrum from sound waves can either have
$k^{3}$ or $k^{9}$ dependence based on recent studies [165, 166, 167]. We will
address this issue in the future.

###  2.1 Broken Power Law Approximation

If the velocity profile is replaced with the RMS fluid velocity, there is a
single scale left in the mean bubble separation. Thus under this
approximation, the gravitational wave power spectrum is a broken power law
(BPL) with a peak set by the mean bubble separation. This approximation has
been the cornerstone of gravitational wave phenomenology and we will briefly
review it here, following Ref. [155]. The peak amplitude is set by a thermal
efficiency factor, which ref [155] develop an analytic approximation for,

$\displaystyle\kappa(\alpha,v_{w})\equiv\begin{cases}v_{w}^{6/5}\frac{6.9\alpha}{1.36-0.037\sqrt{\alpha}+\alpha}~{},&~{}~{}0\lesssim
v_{w}\lesssim 0.2,\\\
\frac{\alpha^{2/5}}{0.017+(0.997+\alpha)^{2/5}}~{},&0.2\lesssim v_{w}\lesssim
0.8,\\\ \frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha}~{},&0.8\lesssim
v_{w}\lesssim 1.\end{cases}$ (2.19)

In the above we have given (very) rough limits of validity of each expression
with approximate values of $v_{w}$. The peak frequency of the spectrum is
given by

$\displaystyle f_{p}^{\rm BPL}=1.19\times
10^{-6}\frac{1}{v_{w}}\left(\frac{g_{*}}{100}\right)^{1/6}\left(\frac{T}{100}\right)\times
0.7\times\frac{\beta}{H_{n}}\,,$ (2.20)

with the mean bubble separation being the only scale in the problem

$\displaystyle R_{*}=(8\pi)^{1/3}\frac{\beta}{v_{w}}\,.$ (2.21)

The peak assuming adiabatic index $\Gamma\approx 4/3$ is then set in our
approximation mostly by the RMS fluid velocity [146]

$\displaystyle\Omega_{p}^{\rm BPL}\approx 1.2\times
10^{-6}\left(\frac{100}{g_{*}}\right)^{1/3}\left(\frac{4}{3}\right)^{2}\bar{U}_{f}(\alpha,v_{w})^{4}\frac{H_{n}}{\beta}v_{w}\Upsilon\,,$
(2.22)

which can be related to the thermal efficiency factor and trace anomaly,

$\displaystyle\bar{U}_{f}\approx\sqrt{\frac{3}{4}\kappa(\alpha,v_{w})\alpha}\,.$
(2.23)

The approximation can be rather poor. In Fig. 1, we show a comparison of the
GW relic density coming from the full sound shell model and that calculated
using the broken power law given by Eq. (2.22). The left (right) panel of the
first row shows the ratio of the peak relic densities (frequencies) as a
function of $v_{w}$ for several fixed values of $\alpha$. The second row
depicts the same quantities as a function of $\alpha$, for several fixed
values of $v_{w}$. There are two metrics of comparison to keep track of.
Firstly, it is clear that $\Omega_{p}/\Omega_{p}^{\rm BPL}\sim 1$ and
$f_{p}/f_{p}^{\rm BPL}\sim 1$ is only achieved for a narrow range of values of
$\\{v_{w},\alpha\\}$. For example, from the top and bottom left panels, it is
clear that $\Omega_{p}/\Omega_{p}^{\rm BPL}\sim 1$ is achieved for $\alpha\sim
0.2-0.4$ with $v_{w}\sim 0.4$. Similarly, from the top and bottom right
panels, it is clear that $0.1\lesssim f_{p}/f_{p}^{\rm BPL}\lesssim 1$ is
achieved for $v_{w}\lesssim 0.3$ and for $f_{p}/f_{p}^{\rm BPL}\sim 0.3$ for
$v_{w}\gtrsim 0.9$ while being relatively stable across values of $\alpha$.

Figure 1: Ratio of $\Omega_{p}/\Omega_{p}^{\rm BPL}$ and $f_{p}/f_{p}^{\rm
BPL}$ in different cases of either fixed $\alpha$’s or $v_{w}$’s. For the
definition of $\Omega_{p}^{\rm BPL}$ see Eq. (2.22).

###  2.2 Double Broken Power Law Fit

Using the full velocity profile in calculating the gravitational wave spectrum
(re)introduces a second scale, the thickness of the sound shell. A double
broken power law fit function to match the calculated GW spectrum was
therefore introduced in Refs. [147, 148]

$\displaystyle\Omega_{\rm GW,fit}=F^{0}_{\rm GW}\Omega_{\rm p}M(s,r_{b},b)\,,$
(2.24)

Here, $\Omega_{p}$ is the peak value of the GW relic density, and the
corresponding frequency is $f_{p}$. Also, $F^{0}_{\rm GW}$ is determined in
Eq. (2.15). The parameter $r_{b}=f_{\text{b}}/f_{p}$ gives the ratio between
the two breaks in the spectrum. The parameter $b$ gives the spectral slope
between the two breaks. The variable $s=f/f_{p}$ is introduced to normalize
the frequency with respect to the peak frequency. The function $M$ is given by
[141, 147, 148]

$\displaystyle
M(s,r_{b},b)=s^{9}{\left(\frac{1+r_{b}^{4}}{r_{b}^{4}+s^{4}}\right)}^{(9-b)/4}\left(\frac{b+4}{b+4-m+ms^{2}}\right)^{(b+4)/2}\,,$
(2.25)

where the choice of $m$ ensures that for $r_{b}<1$ one gets a peak at $s=1$
and $M(1,r_{b},b)=1$. This leads to

$\displaystyle m=\left(9{r_{b}}^{4}+b\right)/\left({r_{b}}^{4}+1\right)\,.$
(2.26)

In our work, when using the fit formula from [147] we use PySwarm optimization
package to determine $b$ and $r_{b}$. We consider a range of values for $b$
and $r_{b}$: $b\in[-200,200]$ and $r_{b}\in[0,200]$.

##  3 New Fit Formula

We will find that the fit function Eq. 2.25 described in the previous Section
gives a poor approximation when calculating the accoustic contribution to the
gravitational wave spectrum where the hydrodynamic equations are solved using
a non-vanishing bag constant $\epsilon$. However, physically we still expect a
doubly broken power law to work. We therefore, consider the following new fit
function that is different from [147]

$\Omega_{\rm
GW}h^{2}=\Omega_{p}\cdot\left(\frac{f}{\tilde{s}_{0}}\right)^{9}\cdot\frac{\left(2+\tilde{r}_{b}^{-12+\tilde{b}}\right)}{\left[\left(\frac{f}{\tilde{s}_{0}}\right)^{\tilde{a}}+\left(\frac{f}{\tilde{s}_{0}}\right)^{\tilde{b}}+\tilde{r}_{b}^{-12+\tilde{b}}\cdot\left(\frac{f}{\tilde{s}_{0}}\right)^{12}\right]}\,\,.$
(3.1)

Here the parameters $\Omega_{p}$, ${\tilde{s}}_{0}$, $\tilde{a}$, $\tilde{b}$
and ${\tilde{r}}_{b}=f_{b}/f_{p}$ are calculated by fitting to the numerical
solution to the GW spectrum. The parameter, $\tilde{r}_{b}$ controls the ratio
between the peaks, $\tilde{b}$ adjusts the spectral slope between the two
breaks in the frequency, while $\tilde{a}$ governs the IR behavior of the GW
spectrum. It should be noted that the two spectral breaks correspond to the
two characteristic length scales in the system: the peak frequency $f_{p}$
corresponds to the mean bubble separation $R_{*}$, while the frequency $f_{\rm
b}$ corresponds to the thickness of the sound shell $\Delta R_{*}$. Compared
to the fit formula in Eq. 2.25, an extra parameter $\tilde{a}$ has been
introduced. The values of all parameters in Eq. 3.1 are produced from our
calculations of the GW relic density based on the sound shell model [141, 147,
148]. The dependence of the fit parameters on the underlying physical
parameters such as $\alpha$ and $v_{w}$ governing the phase transition are
explored in Section 4.

Figure 2: Histogram for comparison of fits using the fit function in Eq.
(2.24) motivated by Ref. [141, 147, 148] and our proposed fit formula in Eq.
(3.1).

Before proceeding to compare the different fit formulas, we pause to give
details about our computation procedure. A Monte Carlo approach to fitting Eq.
3.1 to the GW spectrum of the sound shell model is adopted. Relevant variables
$\Omega_{p}$, ${\tilde{s}}_{0}$, $\tilde{a}$, $\tilde{b}$ and
${\tilde{r}}_{b}$ are first initialized. Our algorithm then executes a loop
100,000 times, each iteration randomly generating the parameters introduced in
Eq. 3.1 to match the GW spectrum from the sound shell model data. The
appropriateness of the fit is measured by calculating residuals - the squared
differences between the logarithms of actual data for $\Omega_{\rm GW}h^{2}$
and the predictions from Eq. 3.1. The program continually updates the
parameters that result in the lowest residuals, storing the best-fit
parameters. After completing all iterations, it uses these parameters to
generate the final predicted model values and calculates the total error
between this fit model and the original sound shell data. The objective is to
robustly determine the parameters that best explain how the GW spectrum
behaves across various frequencies, optimizing the fit to minimize error.

A comparison of the fit to GW spectrum from [141] as given by Eq. 2.24 to our
new fit in Eq. 3.1 is illuminating. We do a $\chi^{2}$ analysis based on the
fit equations and the GW spectrum produced by the sound shell model, using the
full bubble velocity profile based on [141, 142]. The two histograms in Fig. 2
show the results.

Figure 3: A benchmark point of GW from FOPT showing the relic density of GW
with respect to frequency for “old fit" ([147]), “new fit" (the current work),
and sound shell model ([141, 147, 148]) is plotted. Various current and future
observational constraints are also shown and labeled using the integrated
sensitivity curves from [168] for different experiments. Here we assume
$T_{n}=100$ GeV and $\beta/H=1$.

As can be seen from the left panel of Fig. 2, the $\chi^{2}$ values for Eq.
2.24 are much larger than one and most values are above a hundred. On the
other hand, the right panel of Fig. 2 shows that most data points have
$\chi^{2}$ values smaller than one for the new fit function in Eq. 3.1
introduced in this paper. While some data points do have $\chi^{2}$ large than
one, the number of such points decreases as $\chi^{2}$ increases. In Fig. 3,
we consider four benchmark points and show the $\chi^{2}$ values for the old
and new fits in each panel. These benchmarks have these set of FOPT
parameters: $[6.13\times 10^{-2},3.08\times 10^{-3},1,100]$, $[1.28\times
10^{-1},8.66\times 10^{-2},1,100]$, $[3.32\times 10^{-1},6.36\times
10^{-2},1,100]$ and $[5.99\times 10^{-1},3.47\times 10^{-1},1,100]$ for this
FOPT parameters $[v_{w},\alpha,\beta/H,T_{n}({\rm GeV})]$. For these
benchmarks it is apparent that the new fit in Eq. 3.1 matches better with the
sound shell model than the fit formula in Eq. 2.24. Overall, based on these
figures, we have validated the quality of the fit function for most choices of
$v_{w}$ and $\alpha$. This implies that our proposed fit in Eq. 3.1 more
accurately captures the features of GW spectrum from the sound shell model and
can be used by the community (see Appendix B for the fit function code).

##  4 Physical Behavior of Fit and FOPT Parameters

Having demonstrated the appropriateness of our fit function, we proceed
towards a deeper physical understanding of the different parameters that we
have introduced. We perform a scan over the following input parameters: the
wall velocity $v_{w}$ and the strength of transition $\alpha$, keeping the
nucleation temperature $T_{n}=100$ GeV and the bubble nucleation rate
$\beta/H=1$ fixed. Using these parameters, we calculate: $(i)$ the maximum
velocity $v_{\rm max}$, $(ii)$ the width of velocity profile $\delta_{\xi}$
(defined as the width of the velocity profile in a window $95~{}\%$ of $v_{\rm
max}$), $(iii)$ the corresponding value of $\xi$ as $\xi_{\rm max}$, and
$(iv)$ the difference between the values of the enthalpy before the wall and
after it i.e. $\delta_{w}$. Our results are displayed in Figs. 4 \- 8 (also
see Appendix A for more scanning plots).

We begin with a discussion of Fig. 4. Here, we depict values of
${\tilde{r}}_{b}$ (top panel), $\tilde{b}$ (middle panel), and $\tilde{a}$
(bottom panel) on the plane of $\\{v_{w},\alpha\\}$ for three regimes:
deflagration (red), detonation (green) and hybrid (blue). The value of
$\alpha$ is chosen to vary from $10^{-3}-1$, since we are interested in
studying the global behavior of the parameters. For the plot of
${\tilde{r}}_{b}$, it is clear that there is some structure: for example, in
the deflagration regime, the values of ${\tilde{r}}_{b}$ for a given $\alpha$
are small at small values of $v_{w}$, increase around $v_{w}\sim 0.2-0.4$,
before decreasing again. Similarly, in the detonation regime, the values of
${\tilde{r}}_{b}$ for a given value of $v_{w}$ decrease with increasing
$\alpha$. A somewhat opposite behavior can be discerned in the hybrid regime,
where increasing $\alpha$ corresponds to lower ${\tilde{r}}_{b}$. In contrast
to the results for ${\tilde{r}}_{b}$, it is difficult to discern any structure
in the scans of $\tilde{b}$ and $\tilde{a}$. It is therefore fruitful to study
the dependencies in a series of two dimensional plots, which we turn to next.

Figure 4: From top to bottom: scatter plots for $v_{w}$, $\alpha$ and
$\tilde{r}_{b}$; for $v_{w}$, $\alpha$ and $\tilde{b}$; and for $v_{w}$,
$\alpha$ and $\tilde{a}$, respectively. The type of FOPT is shown in the
legend. Here we assume $T_{n}=100$ GeV and $\beta/H=1$.

In Fig. 5, we show the velocity and enthalpy profiles with respect to $\xi$
for different choices of $\alpha$. The scan has been performed over different
values of wall velocity $0<v_{w}<1$. The velocity and enthalpy profiles
exhibit the behavior expected for the deflagration, hybrid, and detonation
regimes, respectively [79]. In the deflagration mode, the plasma in front of
the bubble wall flows outward while remaining static inside the bubble. The
velocity profiles to the left of the maximum in all panels of Fig. 5 exhibit
this behavior. As $v(\xi)$ increases, a discontinuity appears and
$v(\xi)\rightarrow 0$; this is the shock front, beyond which a supersonic
deflagration (hybrid) mode develops when the velocity exceeds the sound speed.
The value of $\alpha$ determines the maximum value of $v_{\rm max}$ and
maximum of $w_{\rm max}$ where the velocity profile maximises over different
values of $\xi$. Each value of $v_{\rm max}$ occurs when $\xi=v_{w}$.

Fig. 6 depicts the behavior of the various physical quantities, as well as the
fit parameters, as a function of $v_{w}$ for several benchmark values of
$\alpha=0.001,0.01,0.1,0.3,$ and $0.5$. The behaviour of $v_{max}$ w.r.t.
$v_{w}$ is shown in the bottom left panel of Fig. 6. It reaches a maximum that
depends on the value $\alpha$; a similar pattern is exhibited by the thermal
efficiency factor $\kappa$ defined in Eq. 2.11 and depicted in the top right
panel. The results of this panel, as well as the top right panel depicting
$\kappa$, are consistent with the left panels of Fig. 5. The peak frequency
$f_{p}$ is depicted in the bottom right panel of Fig. 6 and shows a local
maximum where $v_{\rm max}$ is maximized. The dependence of $\Omega_{p}$
w.r.t. $v_{w}$ is shown in the second row of the left panel of Fig. 6.
$\Omega_{p}$ is a decreasing function of $v_{w}$ in the deflagration regime;
on reaching $v_{\rm max}$, it rises somewhat and then continues falling in the
hybrid and detonation regimes. We have defined the width of velocity profile
$\delta_{\xi}$ as the difference between the values of $\xi$ when $\xi_{\rm
max}=v_{w}$ and $v(\xi_{\rm width})=~{}0.95~{}v_{\rm max}$. This quantity is
shown in the left panel of the first row. Its behavior matches with
expectations from the velocity profiles shown in Fig. 5. The value of
$\xi_{\rm max}$ increases monotonously with $v_{w}$, independently of
$\alpha$, as shown in the right of panel of the fourth row of Fig. 6.

We now turn to a discussion of the fit parameters $\tilde{b}$ and $\tilde{a}$.
Firstly, it should be noted that smaller values of $\tilde{b}$ correspond to a
situation where more energy resides in the contribution from the thickness of
the sound shell rather than from the bubble separation. Conversely, larger
values of $\tilde{b}$ correspond to more energy residing in the contribution
from the bubble separation. While it is difficult to discern structure, we do
find a linearly increasing trend in $\tilde{b}$ as $v_{w}$ increases, shown by
trend lines, albeit crude. Physically, it is reasonable that as $v_{w}$
increases, more energy should reside in the bubble separation contribution.
The behavior of the fit parameter $\tilde{a}$ in the left panel of the third
row is interesting. It should be noted that smaller values of $\tilde{a}$
correspond to sharper IR spectra, while larger values of $\tilde{a}$
correspond to shallower IR spectra. It is clear that the value of $\tilde{a}$
is insensitive to variations in $v_{w}$ in the deflagration mode, but swiftly
decreases after $v_{\rm max}$, for all values of $\alpha$.

We now turn to a discussion of Fig. 8, where a scan over $\alpha$ between $0$
and $1$ for fixed values of $v_{w}=0.1,0.4,0.7,0.95$ was performed. The plot
of the ${\tilde{a}}$ parameter versus $\alpha$ is shown in the fourth row,
left panel. ${\tilde{a}}$ appears to be constant with a value of around $4$,
which reduces when it reaches the maximum of $v_{\rm max}$. This is consistent
with the behavior in Fig. 6. The plot of ${\tilde{b}}$ versus $\alpha$ is
shown in the second row of the right panel of Fig. 7. It is difficult to
discern a specific pattern here; nevertheless, we have made a crude linear
fit. For $v_{w}=0.1,0.4,$ and $0.7$, we see that ${\tilde{b}}$ decreases very
mildly with increasing $\alpha$.

##  5 Conclusions

The sound shell model currently provides the most sophisticated semi-
analytical framework for calculating the GW spectrum in models of BSM physics.
In conjunction with numerical simulations, a precision GW spectrum frontier is
then already somewhat within reach for particle physicists; one can then ask
why fit functions of the type studied in this paper are even necessary, or why
the physical interpretation of fit parameters is an interesting question. The
response is that full calculations of the sound shell model, not to speak of
actual simulations, are difficult and expensive; to perform large-scale scans
over the parameter space of particle physics models and obtain the predicted
GW spectrum, portable and simple fits to the spectrum, and the interpretation
of the fit parameters, become necessary. Providing such a fit – currently the
one that most faithfully captures the full results coming from the sound shell
model – has been the goal of this paper. We have shared the results of our
work in a set of CSV, Python and Mathematica files333GitHub link.. These files
can be used by the BSM community to efficiently depict the GW spectrum at
various points in the parameter space of their models.

Studies in GW physics are rapidly entering a precision frontier. Indeed, as
pointed out by a subset of the present authors in [158], uncertainties of
several orders of magnitude can be introduced in the calculation of the GW
spectrum if one neglects careful treatments of the source lifetime; mean
bubble separation; employing a beyond the bag model approximation when solving
the hydrodynamics equations and explicitly calculating the fraction of energy
in the fluid from these equations rather than using a fit. The incorporation
of these effects constitutes, currently, the “highest level of diligence" (in
the parlance of [158]) in obtaining the GW spectrum. Calculating the full GW
spectrum at the highest level of diligence, and providing the relevant fit
functions, constitute concrete future directions.

Figure 5: Velocity and ethalpy profiles for fixed values of $\alpha$. Each row
of these panels corresponds to a color that are blue, orange, green, red and
violet equivalent to a specific value of
$\alpha=0.001,\,0.01,\,0.1,\,0.3,\,0.5$ shown in the legends.

Figure 6: Plots for different variables with respect to $\alpha$ for fixed
$\alpha$’s are shown. Each of colors blue, orange, green, red and violet
corresponds to a specific value of $\alpha=0.001,\,0.01,\,0.1,\,0.3,\,0.5$
shown in the legends. We describe them in the main text. Here we assume
$T_{n}=100$ GeV and $\beta/H=1$.

Figure 7: Plots for different variables with respect to $\alpha$ for fixed
$v_{w}$’s are shown. Each of colors blue, orange, green and red corresponds to
a specific value of $v_{w}=0.1,\,0.4,\,0.7,\,0.95$. We describe them in the
main text. Here we assume $T_{n}=100$ GeV and $\beta/H=1$.

Figure 8: Velocity and enthalpy profiles for fixed values of wall velocity
$v_{w}$. Each row of these panels corresponds to a color that are blue,
orange, green and red equivalent to a specific value of
$v_{w}=0.1,\,0.4,\,0.7,\,0.95$ shown in the legends.

The precision GW frontier has several other challenges that lie outside the
scope of this work. The transitions between the different power law regimes in
the GW spectrum, the existence of new regimes, the shape of the GW spectrum
near the peak, and the peak frequency depend on several factors: fluctuations
in the local temperature affect the distribution of nucleated bubbles [169];
energy lost to vorticity can suppress the spectrum [163]; and dissipative
effects encoded by the shear viscosity, bulk viscosity, and thermal conduction
can also suppress the spectrum [170]. Calculating these effects would make
further material progress on the precision frontier.

## Acknowledgments

We are grateful to Jose Miguel No for collaboration during the initial stages
of the project, and to Daniel Vagie for useful discussions. F.H. and K.S.
thank the organizers of the Mitchell Conference in May 2024 at Texas A & M
University for their hospitality and support during the final stages of this
project. They are also grateful to the organizers of workshop of Center for
Theoretical Underground Physics and Related Areas (CETUP* - 2024), The
Institute for Underground Science at Sanford Underground Research Facility
(SURF), Lead, South Dakota for their hospitality and financial support. K.S.
would like to thank the Aspen Center for Theoretical Physics, supported by the
National Science Foundation grant PHY-2210452, for hospitality during the
course of this work.

##  Appendix A Other Scanning Plots

We have shown other scanning plots over $\delta_{\xi}$ and $\kappa$ parameter
space in this section in Fig. 9. We also represented the variation of some of
our new fit parameters: $\tilde{r}_{b}$, $\tilde{b}$ and $\tilde{a}$ as
introduced in Eq (3.1). In first panel we have shown the plot for
$\tilde{r}_{b}$ i.e. the ratio of frequencies between two peaks. As Fig. 9
shows large values of $\tilde{r}_{b}$ are accumulating between
$0\lesssim\delta_{\xi}\lesssim 0.01$ and $10^{-2}\lesssim\kappa\lesssim 0.5$
for all regimes and $\tilde{r}_{b}$ decreases as $\delta_{\xi}$ and $\kappa$
are changing from these regions. Also, the second panel shows the changes
parameter $\tilde{b}$ w.r.t. $\delta_{\xi}$ and $\kappa$. For deflagration and
hybrid regimes increasing $\delta_{\xi}$ and $\kappa$ makes $\tilde{b}$
larger. However, for the detonation regime $\tilde{b}$ slightly reduces by
decreasing of $\kappa$. In the last panel we examine the variation of
parameter $\tilde{a}$ that deviates from $4$ depending on the values of
$\delta_{\xi}$ and $\kappa$. As it can be seen the points in deflagration and
hybrid regimes are almost independent of $\delta_{\xi}$ and $\kappa$. However,
as $\kappa$ increases $\tilde{a}$ also increases for detonation regime.

Figure 9: From top to bottom: scatter plots for $\delta_{\xi}$, $\kappa$ and
$\tilde{r}_{b}$; for $\delta_{\xi}$, $\kappa$ and $\tilde{b}$; and for
$\delta_{\xi}$, $\kappa$ and $\tilde{a}$. The type of FOPT is shown in the
legend. Here we assume $T_{n}=100$ GeV and $\beta/H=1$.

##  Appendix B GW Spectrum for a Specific Choice of FOPT Paramters

Here we briefly explain how one can use our script to calculate the spectrum
of GW from sound waves of FOPT. We provide two scripts written in Python and
Mathematica.

###  B.1 Mathematica script

The Mathematica script is displayed below:

⬇

1(*Load the CSV file*)

2data = Import[NotebookDirectory[] <> "data-fopt-fit.csv"];

3

4(*Extract and clean column names from the first row*)

5columnNames = StringTrim /@ First[data];

6

7(*Remove the header row to isolate the data*)

8data = Rest[data];

9

10(*Print column names to check*)

11Print[columnNames];

12

13(*Assign each column to a separate variable using safe access to \

14Position*)

15vw = If[Length[Position[columnNames, "vw"]] > 0,

16 data[[All, Position[columnNames, "vw"][[1, 1]]]], {}];

17alpha = If[Length[Position[columnNames, "alpha"]] > 0,

18 data[[All, Position[columnNames, "alpha"][[1, 1]]]], {}];

19betaH = If[Length[Position[columnNames, "betaH"]] > 0,

20 data[[All, Position[columnNames, "betaH"][[1, 1]]]], {}];

21Tn = If[Length[Position[columnNames, "Tn"]] > 0,

22 data[[All, Position[columnNames, "Tn"][[1, 1]]]], {}];

23fp = If[Length[Position[columnNames, "f_p"]] > 0,

24 data[[All, Position[columnNames, "f_p"][[1, 1]]]], {}];

25Omp = If[Length[Position[columnNames, "Om_p"]] > 0,

26 data[[All, Position[columnNames, "Om_p"][[1, 1]]]], {}];

27at = If[Length[Position[columnNames, "at"]] > 0,

28 data[[All, Position[columnNames, "at"][[1, 1]]]], {}];

29bt = If[Length[Position[columnNames, "bt"]] > 0,

30 data[[All, Position[columnNames, "bt"][[1, 1]]]], {}];

31rbt = If[Length[Position[columnNames, "rbt"]] > 0,

32 data[[All, Position[columnNames, "rbt"][[1, 1]]]], {}];

33s0t = If[Length[Position[columnNames, "s0t"]] > 0,

34 data[[All, Position[columnNames, "s0t"][[1, 1]]]], {}];

35omega0t =

36 If[Length[Position[columnNames, "omega_0t"]] > 0,

37 data[[All, Position[columnNames, "omega_0t"][[1, 1]]]], {}];

38K = If[Length[Position[columnNames, "K"]] > 0,

39 data[[All, Position[columnNames, "K"][[1, 1]]]], {}];

40

41(*Define the data points*)

42dataVWAlphaFP = Transpose[{vw, alpha, fp}];

43dataVWAlphaOmP = Transpose[{vw, alpha, Omp}];

44dataVWAlphaAt = Transpose[{vw, alpha, at}];

45dataVWAlphaBt = Transpose[{vw, alpha, bt}];

46dataVWAlphaRbt = Transpose[{vw, alpha, rbt}];

47dataVWAlphaS0t = Transpose[{vw, alpha, s0t}];

48dataVWAlphaOmega0t = Transpose[{vw, alpha, omega0t}];

49dataVWAlphaK = Transpose[{vw, alpha, K}];

50

51

52(*Create interpolating functions*)

53intFP = Interpolation[dataVWAlphaFP, InterpolationOrder -> 1];

54intOmP = Interpolation[dataVWAlphaOmP, InterpolationOrder -> 1];

55intAt = Interpolation[dataVWAlphaAt, InterpolationOrder -> 1];

56intBt = Interpolation[dataVWAlphaBt, InterpolationOrder -> 1];

57intRbt = Interpolation[dataVWAlphaRbt, InterpolationOrder -> 1];

58intS0t = Interpolation[dataVWAlphaS0t, InterpolationOrder -> 1];

59intOmega0t =

60 Interpolation[dataVWAlphaOmega0t, InterpolationOrder -> 1];

61intK = Interpolation[dataVWAlphaK, InterpolationOrder -> 1];

62

63(*Evaluate at a point*)

64vw0 = 0.5;

65alpha0 = 0.5;

66fp0 = intFP[vw0, alpha0];

67at0 = intAt[vw0, alpha0];

68bt0 = intBt[vw0, alpha0];

69rbt0 = intRbt[vw0, alpha0];

70s0t0 = intS0t[vw0, alpha0];

71Omegap0 = intOmega0t[vw0, alpha0];

72K0 = intK[vw0, alpha0];

73

74(*Degrees of freedom*)

75gstarData = Import[NotebookDirectory[] <> "gstar.txt", "Table"];

76gHighTemp = {{1.*10^14*gstarData[[1, 1]], gstarData[[1, 2]]}};

77dataG1 = Join[gHighTemp, gstarData];

78TData = dataG1[[All, 1]];

79gstarData = dataG1[[All, 2]];

80gstarFun = Interpolation[Transpose[{TData, gstarData}]];

81

82(*Frequency domain*)

83fdomainDefault = 10^Range[-24, 24, 48/999] // N;

84

85(*Generating fit parameters from the interpolated function*)

86fb0 = rbt0*fp0;

87(*This part can be modified based on different nucleation \

88tamperatures and bubble nucleation rates.*)

89betaHDefault = 1;

90betaHNew = 10;

91TnDefault = 100;

92TnNew = 1000;

93

94rstarDefault = (8 \\[Pi])^(1/3)*vw0/betaHDefault;

95rstarNew = (8 \\[Pi])^(1/3)*vw0/betaHNew;

96gstarDefault = gstarFun[TnDefault];

97gstarNew = gstarFun[TnNew];

98

99xDefault = rstarDefault/(K0^(1/2));

100JDefault = rstarDefault*(1 - 1/(1 + 2*xDefault)^(1/2));

101FgwDefault = 3.57*10^-5*(100/gstarDefault)^(1/3);

102xNew = rstarNew/(K0^(1/2));

103JNew = rstarNew*(1 - 1/(1 + 2*xNew)^(1/2));

104FgwNew = 3.57*10^-5*(100/gstarNew)^(1/3);

105

106f0Default = 2.6*10^-6*(TnDefault/100)*(gstarDefault/100)^(1/6);

107kR = fdomainDefault*(rstarDefault/f0Default);

108

109

110(*GW spectrum-default*)

111omegaFitDefault =

112 Omegap0*(fdomainDefault/

113 s0t0)^9*((2 +

114 rbt0^(-12 + bt0))/((fdomainDefault/s0t0)^

115 at0 + (fdomainDefault/s0t0)^bt0 +

116 rbt0^(-12 + bt0)*(fdomainDefault/s0t0)^12)); // Quiet

117

118

119

120f0New = 2.6*10^-6*(TnNew/100)*(gstarNew/100)^(1/6);

121fdomainNew = kR/(rstarNew/f0New);

122omegaFitNew = (omegaFitDefault/(FgwDefault*JDefault))*(FgwNew*JNew);

123

124

125(*Plotting*)

126ListLogLogPlot[{Transpose[{fdomainDefault, omegaFitDefault}],

127 Transpose[{fdomainNew, omegaFitNew}]},

128 PlotRange -> {{1.*10^-12, 1.*10^4}, {1.*10^-22, 1.*10^-2}},

129 Frame -> True,

130 FrameLabel -> {"f (Hz)",

131 "\\!\\(\\*SubscriptBox[\\(\\[CapitalOmega]\\), \

132\\(GW\\)]\\)\\!\\(\\*SuperscriptBox[\\(h\\), \\(2\\)]\\)"},

133 PlotLabel -> "FOPT - GW", GridLines -> Automatic, Joined -> True,

134 PlotLegends -> {"Tn = 100 GeV, beta/H = 1",

135 "Tn = 1000 GeV, beta/H = 10"}, FrameStyle -> 18, ImageSize -> 800]

Listing 1: Mathematica code to calculate the spectrum of GW from sound waves
of FOPT.

###  B.2 Python script

We have written the following code in Python to calculate GW from the phase
transitions:

⬇

1# Python script for this manuscript

2import numpy as np

3from scipy.interpolate import LinearNDInterpolator

4from scipy.interpolate import interp1d

5import matplotlib.pyplot as plt

6import pandas as pd

7

8# Load the CSV file into a DataFrame

9df = pd.read_csv(’data-fopt-fit.csv’)

10

11# Assign each column to a separate variable

12vw = df[’vw’].to_numpy()

13alpha = df[’alpha’].to_numpy()

14betaH = df[’betaH’].to_numpy()

15Tn = df[’Tn’].to_numpy()

16f_p = df[’f_p’].to_numpy()

17Om_p = df[’Om_p’].to_numpy()

18at = df[’at’].to_numpy()

19bt = df[’bt’].to_numpy()

20rbt = df[’rbt’].to_numpy()

21s0t = df[’s0t’].to_numpy()

22omega_0t = df[’omega_0t’].to_numpy()

23K = df[’K’].to_numpy()

24

25# Interpolatoin over data (t : tilde)

26int_f_p = LinearNDInterpolator(list(zip(vw, alpha)), f_p)

27int_Om_p = LinearNDInterpolator(list(zip(vw, alpha)), Om_p)

28int_at = LinearNDInterpolator(list(zip(vw, alpha)), at)

29int_bt = LinearNDInterpolator(list(zip(vw, alpha)), bt)

30int_rbt = LinearNDInterpolator(list(zip(vw, alpha)), rbt)

31int_s0t = LinearNDInterpolator(list(zip(vw, alpha)), s0t)

32int_omega_0t = LinearNDInterpolator(list(zip(vw, alpha)), omega_0t)

33int_K = LinearNDInterpolator(list(zip(vw, alpha)), K)

34

35vw0 = 0.5

36alpha0 = 0.5

37

38# Generatng fit parameters from the interpolated function

39f_p0 = int_f_p(vw0, alpha0)

40at0 = int_at(vw0, alpha0)

41bt0 = int_bt(vw0, alpha0)

42rbt0 = int_rbt(vw0, alpha0)

43s0t0 = int_s0t(vw0, alpha0)

44omega_p0 = int_omega_0t(vw0, alpha0)

45K0 = int_K(vw0, alpha0)

46

47# Degrees of freedom

48# arXiv: https://arxiv.org/pdf/1503.03513

49filename = "gstar.txt"

50data_g = np.loadtxt(filename, skiprows=0)

51g_high_temp = np.array([[1.e14 * data_g[0, 0], data_g[0, 1]]])

52data_g1 = np.concatenate((g_high_temp, data_g))

53T_data = data_g1[:, 0] # [GeV]

54gstar_data = data_g1[:, 1] # [GeV]

55gstar_fun = interp1d(T_data, gstar_data)

56#

57f_b0 = rbt0 * f_p0

58# This part can be modified based on different nucleation tamperatures and
bubble nucleation rates.

59# beta over H

60betaH_default = 1

61betaH_new = 10

62# Nucleation temperature

63Tn_default = 100 # GeV

64Tn_new = 1000 # GeV

65#

66rstar_default = (8*np.pi)**(1/3)*vw0/betaH_default

67rstar_new = (8*np.pi)**(1/3)*vw0/betaH_new

68# Degrees of Freedom

69gstar_default = gstar_fun(Tn_default)

70gstar_new = gstar_fun(Tn_new)

71#

72x_default = rstar_default/(K0**(1/2))

73J_default = rstar_default*(1-1/(1+2*x_default)**(1/2))

74Fgw_default = 3.57*1e-5*(100/gstar_default)**(1/3)

75x_new = rstar_new/(K0**(1/2))

76J_new = rstar_new*(1-1/(1+2*x_new)**(1/2))

77Fgw_new = 3.57*1e-5*(100/gstar_new)**(1/3)

78

79# Frequency

80fdomain_default = 10**np.linspace(-24,24,1000)

81

82f0_default = 2.6*1.e-6*(Tn_default/100)*(gstar_default/100)**(1/6)

83omega_p0_default = omega_p0

84kR = fdomain_default * (rstar_default/f0_default)

85

86# GW spectrum - default

87omega_fit_default = omega_p0_default * (fdomain_default/s0t0)**9 * (

88 (2 + rbt0**(-12 + bt0)) /

89 ((fdomain_default/s0t0)**at0 + (fdomain_default/s0t0)**bt0 +

90 rbt0**(-12 + bt0) * (fdomain_default/s0t0)**12))

91

92f0_new = 2.6*1.e-6*(Tn_new/100)*(gstar_new/100)**(1/6)

93fdomain_new = kR / (rstar_new/f0_new)

94omega_fit_new = (omega_fit_default / (Fgw_default*J_default) ) *
(Fgw_new*J_new)

95

96

97# Plotting

98plt.figure(figsize=(10, 6))

99plt.loglog(fdomain_default, omega_fit_default, label=’FOPT - Sound Waves -
GW - $T_n = 100$ GeV and $\\\beta /H = 1$’)

100plt.loglog(fdomain_new, omega_fit_new, label=’FOPT - Sound Waves - GW’)

101

102# Set limits for frequency and GW relic axes

103plt.xlim(1.e-12, 1.e4)

104plt.ylim(1.e-22, 1.e-2)

105plt.legend()

106plt.xlabel(’f (Hz)’)

107plt.ylabel(’$\Omega_{GW}h^2$(f)’)

108plt.grid(True)

109plt.show()

Listing 2: Python code to calculate the spectrum of GW from sound waves of
FOPT.

Using the above script and the files provided along with this paper one can
generate the predicted GW spectrum by inserting the input parameters from a
beyond Standard Model scenario and the cosmological scale that the phase
transition occurs at. These input parameters are $\alpha$, $v_{w}$, $\beta/H$
and $T_{n}$: with these, one can produce the GW spectrum from the above Python
and Mathematica scripts 444GitHub link..

## References

  * [1] C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024 [1910.13125].
  * [2] R. Caldwell et al., _Detection of early-universe gravitational-wave signatures and fundamental physics_ , _Gen. Rel. Grav._ 54 (2022) 156 [2203.07972].
  * [3] A. Alves, D. Gonçalves, T. Ghosh, H.-K. Guo and K. Sinha, _Di-Higgs Production in the $4b$ Channel and Gravitational Wave Complementarity_, _JHEP_ 03 (2020) 053 [1909.05268].
  * [4] G. C. Dorsch, S. J. Huber and J. M. No, _A strong electroweak phase transition in the 2HDM after LHC8_ , _JHEP_ 10 (2013) 029 [1305.6610].
  * [5] P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt and A. Wlotzka, _Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited_ , _JHEP_ 02 (2017) 121 [1612.04086].
  * [6] G. C. Dorsch, S. J. Huber, K. Mimasu and J. M. No, _Hierarchical versus degenerate 2HDM: The LHC run 1 legacy at the onset of run 2_ , _Phys. Rev. D_ 93 (2016) 115033 [1601.04545].
  * [7] G. C. Dorsch, S. J. Huber, T. Konstandin and J. M. No, _A Second Higgs Doublet in the Early Universe: Baryogenesis and Gravitational Waves_ , _JCAP_ 05 (2017) 052 [1611.05874].
  * [8] J. Bernon, L. Bian and Y. Jiang, _A new insight into the phase transition in the early Universe with two Higgs doublets_ , _JHEP_ 05 (2018) 151 [1712.08430].
  * [9] G. C. Dorsch, S. J. Huber, K. Mimasu and J. M. No, _The Higgs Vacuum Uplifted: Revisiting the Electroweak Phase Transition with a Second Higgs Doublet_ , _JHEP_ 12 (2017) 086 [1705.09186].
  * [10] J. O. Andersen, T. Gorda, A. Helset, L. Niemi, T. V. I. Tenkanen, A. Tranberg et al., _Nonperturbative Analysis of the Electroweak Phase Transition in the Two Higgs Doublet Model_ , _Phys. Rev. Lett._ 121 (2018) 191802 [1711.09849].
  * [11] K. Kainulainen, V. Keus, L. Niemi, K. Rummukainen, T. V. I. Tenkanen and V. Vaskonen, _On the validity of perturbative studies of the electroweak phase transition in the Two Higgs Doublet model_ , _JHEP_ 06 (2019) 075 [1904.01329].
  * [12] X. Wang, F. P. Huang and X. Zhang, _Gravitational wave and collider signals in complex two-Higgs doublet model with dynamical CP-violation at finite temperature_ , _Phys. Rev. D_ 101 (2020) 015015 [1909.02978].
  * [13] W. Su, A. G. Williams and M. Zhang, _Strong first order electroweak phase transition in 2HDM confronting future Z & Higgs factories_, _JHEP_ 04 (2021) 219 [2011.04540].
  * [14] H. Davoudiasl, I. M. Lewis and M. Sullivan, _Multi-TeV signals of baryogenesis in a Higgs troika model_ , _Phys. Rev. D_ 104 (2021) 015024 [2103.12089].
  * [15] T. Biekötter, S. Heinemeyer, J. M. No, M. O. Olea and G. Weiglein, _Fate of electroweak symmetry in the early Universe: Non-restoration and trapped vacua in the N2HDM_ , _JCAP_ 06 (2021) 018 [2103.12707].
  * [16] Z. Zhang, C. Cai, X.-M. Jiang, Y.-L. Tang, Z.-H. Yu and H.-H. Zhang, _Phase transition gravitational waves from pseudo-Nambu-Goldstone dark matter and two Higgs doublets_ , _JHEP_ 05 (2021) 160 [2102.01588].
  * [17] M. Aoki, T. Komatsu and H. Shibuya, _Possibility of a multi-step electroweak phase transition in the two-Higgs doublet models_ , _PTEP_ 2022 (2022) 063B05 [2106.03439].
  * [18] D. Gonçalves, A. Kaladharan and Y. Wu, _Electroweak phase transition in the 2HDM: Collider and gravitational wave complementarity_ , _Phys. Rev. D_ 105 (2022) 095041 [2108.05356].
  * [19] V. Q. Phong, N. M. Anh and H. N. Long, _Dual electroweak phase transition in the two-Higgs-doublet model with the S3 discrete symmetry_ , _Phys. Rev. D_ 107 (2023) 035020 [2209.14672].
  * [20] T. Biekötter, S. Heinemeyer, J. M. No, M. O. Olea-Romacho and G. Weiglein, _The trap in the early Universe: impact on the interplay between gravitational waves and LHC physics in the 2HDM_ , _JCAP_ 03 (2023) 031 [2208.14466].
  * [21] Anisha, L. Biermann, C. Englert and M. Mühlleitner, _Two Higgs doublets, effective interactions and a strong first-order electroweak phase transition_ , _JHEP_ 08 (2022) 091 [2204.06966].
  * [22] O. Atkinson, M. Black, C. Englert, A. Lenz, A. Rusov and J. Wynne, _The flavourful present and future of 2HDMs at the collider energy frontier_ , _JHEP_ 11 (2022) 139 [2202.08807].
  * [23] T. Biekötter, S. Heinemeyer, J. M. No, K. Radchenko, M. O. O. Romacho and G. Weiglein, _First shot of the smoking gun: probing the electroweak phase transition in the 2HDM with novel searches for A → ZH in ${\ell}^{+}{\ell}^{-}t\overline{t}$ and $\nu\nu b\overline{b}$ final states_, _JHEP_ 01 (2024) 107 [2309.17431].
  * [24] D. Gonçalves, A. Kaladharan and Y. Wu, _Gravitational waves, bubble profile, and baryon asymmetry in the complex 2HDM_ , _Phys. Rev. D_ 108 (2023) 075010 [2307.03224].
  * [25] L. Gráf, S. Jana, A. Kaladharan and S. Saad, _Gravitational wave imprints of left-right symmetric model with minimal Higgs sector_ , _JCAP_ 05 (2022) 003 [2112.12041].
  * [26] G. Arcadi, N. Benincasa, A. Djouadi and K. Kannike, _Two-Higgs-doublet-plus-pseudoscalar model: Collider, dark matter, and gravitational wave signals_ , _Phys. Rev. D_ 108 (2023) 055010 [2212.14788].
  * [27] N. Blinov, S. Profumo and T. Stefaniak, _The Electroweak Phase Transition in the Inert Doublet Model_ , _JCAP_ 07 (2015) 028 [1504.05949].
  * [28] F. P. Huang and J.-H. Yu, _Exploring inert dark matter blind spots with gravitational wave signatures_ , _Phys. Rev. D_ 98 (2018) 095022 [1704.04201].
  * [29] A. Paul, B. Banerjee and D. Majumdar, _Gravitational wave signatures from an extended inert doublet dark matter model_ , _JCAP_ 10 (2019) 062 [1908.00829].
  * [30] S. Fabian, F. Goertz and Y. Jiang, _Dark matter and nature of electroweak phase transition with an inert doublet_ , _JCAP_ 09 (2021) 011 [2012.12847].
  * [31] N. Benincasa, L. Delle Rose, K. Kannike and L. Marzola, _Multi-step phase transitions and gravitational waves in the inert doublet model_ , _JCAP_ 12 (2022) 025 [2205.06669].
  * [32] A. Paul, D. Majumdar and B. Banerjee, _Gravitational Wave Signatures from First-Order Phase Transitions in an Extended Inert Doublet Dark Matter Model_ , _Springer Proc. Phys._ 277 (2022) 73.
  * [33] S. Jiang, F. P. Huang and X. Wang, _Bubble wall velocity during electroweak phase transition in the inert doublet model_ , _Phys. Rev. D_ 107 (2023) 095005 [2211.13142].
  * [34] M. D. Astros, S. Fabian and F. Goertz, _Minimal Inert Doublet benchmark for dark matter and the baryon asymmetry_ , _JCAP_ 02 (2024) 052 [2307.01270].
  * [35] N. Benincasa, A. Hryczuk, K. Kannike and M. Laletin, _Phase transitions and gravitational waves in a model of $\mathbb{Z}$3 scalar dark matter_, _JHEP_ 02 (2024) 207 [2312.04627].
  * [36] M. Carena, M. Quiros and C. E. M. Wagner, _Opening the window for electroweak baryogenesis_ , _Phys. Lett. B_ 380 (1996) 81 [hep-ph/9603420].
  * [37] J. R. Espinosa, _Dominant two loop corrections to the MSSM finite temperature effective potential_ , _Nucl. Phys. B_ 475 (1996) 273 [hep-ph/9604320].
  * [38] D. Delepine, J. M. Gerard, R. Gonzalez Felipe and J. Weyers, _A Light stop and electroweak baryogenesis_ , _Phys. Lett. B_ 386 (1996) 183 [hep-ph/9604440].
  * [39] J. M. Cline and K. Kainulainen, _Supersymmetric electroweak phase transition: Beyond perturbation theory_ , _Nucl. Phys. B_ 482 (1996) 73 [hep-ph/9605235].
  * [40] M. Losada, _High temperature dimensional reduction of the MSSM and other multiscalar models_ , _Phys. Rev. D_ 56 (1997) 2893 [hep-ph/9605266].
  * [41] M. Laine, _Effective theories of MSSM at high temperature_ , _Nucl. Phys. B_ 481 (1996) 43 [hep-ph/9605283].
  * [42] D. Bodeker, P. John, M. Laine and M. G. Schmidt, _The Two loop MSSM finite temperature effective potential with stop condensation_ , _Nucl. Phys. B_ 497 (1997) 387 [hep-ph/9612364].
  * [43] B. de Carlos and J. R. Espinosa, _The Baryogenesis window in the MSSM_ , _Nucl. Phys. B_ 503 (1997) 24 [hep-ph/9703212].
  * [44] J. M. Cline and G. D. Moore, _Supersymmetric electroweak phase transition: Baryogenesis versus experimental constraints_ , _Phys. Rev. Lett._ 81 (1998) 3315 [hep-ph/9806354].
  * [45] M. Losada, _The Two loop finite temperature effective potential of the MSSM and baryogenesis_ , _Nucl. Phys. B_ 537 (1999) 3 [hep-ph/9806519].
  * [46] M. Laine and K. Rummukainen, _Two Higgs doublet dynamics at the electroweak phase transition: A Nonperturbative study_ , _Nucl. Phys. B_ 597 (2001) 23 [hep-lat/0009025].
  * [47] M. Carena, G. Nardini, M. Quiros and C. E. M. Wagner, _The Baryogenesis Window in the MSSM_ , _Nucl. Phys. B_ 812 (2009) 243 [0809.3760].
  * [48] A. Delgado, G. Nardini and M. Quiros, _The Light Stop Scenario from Gauge Mediation_ , _JHEP_ 04 (2012) 137 [1201.5164].
  * [49] M. Carena, G. Nardini, M. Quiros and C. E. M. Wagner, _MSSM Electroweak Baryogenesis and LHC Data_ , _JHEP_ 02 (2013) 001 [1207.6330].
  * [50] D. J. H. Chung, A. J. Long and L.-T. Wang, _125 GeV Higgs boson and electroweak phase transition model classes_ , _Phys. Rev. D_ 87 (2013) 023509 [1209.1819].
  * [51] W. Huang, J. Shu and Y. Zhang, _On the Higgs Fit and Electroweak Phase Transition_ , _JHEP_ 03 (2013) 164 [1210.0906].
  * [52] M. Laine, G. Nardini and K. Rummukainen, _Lattice study of an electroweak phase transition at $m_{h}\backsimeq$ 126 GeV_, _JCAP_ 01 (2013) 011 [1211.7344].
  * [53] S. Liebler, S. Profumo and T. Stefaniak, _Light Stop Mass Limits from Higgs Rate Measurements in the MSSM: Is MSSM Electroweak Baryogenesis Still Alive After All?_ , _JHEP_ 04 (2016) 143 [1512.09172].
  * [54] H. H. Patel and M. J. Ramsey-Musolf, _Stepping Into Electroweak Symmetry Breaking: Phase Transitions and Higgs Phenomenology_ , _Phys. Rev. D_ 88 (2013) 035013 [1212.5652].
  * [55] M. Chala, M. Ramos and M. Spannowsky, _Gravitational wave and collider probes of a triplet Higgs sector with a low cutoff_ , _Eur. Phys. J. C_ 79 (2019) 156 [1812.01901].
  * [56] S. Baum, M. Carena, N. R. Shah, C. E. M. Wagner and Y. Wang, _Nucleation is more than critical: A case study of the electroweak phase transition in the NMSSM_ , _JHEP_ 03 (2021) 055 [2009.10743].
  * [57] M. J. Kazemi and S. S. Abdussalam, _Electroweak Phase Transition in an Inert Complex Triplet Model_ , _Phys. Rev. D_ 103 (2021) 075012 [2103.00212].
  * [58] I. Baldes, S. Blasi, A. Mariotti, A. Sevrin and K. Turbang, _Baryogenesis via relativistic bubble expansion_ , _Phys. Rev. D_ 104 (2021) 115029 [2106.15602].
  * [59] A. Azatov, M. Vanvlasselaer and W. Yin, _Baryogenesis via relativistic bubble walls_ , _JHEP_ 10 (2021) 043 [2106.14913].
  * [60] M. Jiang, L. Bian, W. Huang and J. Shu, _Impact of a complex singlet: Electroweak baryogenesis and dark matter_ , _Phys. Rev. D_ 93 (2016) 065032 [1502.07574].
  * [61] C.-W. Chiang, M. J. Ramsey-Musolf and E. Senaha, _Standard Model with a Complex Scalar Singlet: Cosmological Implications and Theoretical Considerations_ , _Phys. Rev. D_ 97 (2018) 015005 [1707.09960].
  * [62] A. Ahriche, K. Hashino, S. Kanemura and S. Nasri, _Gravitational Waves from Phase Transitions in Models with Charged Singlets_ , _Phys. Lett. B_ 789 (2019) 119 [1809.09883].
  * [63] N. Chen, T. Li, Y. Wu and L. Bian, _Complementarity of the future $e^{+}e^{-}$ colliders and gravitational waves in the probe of complex singlet extension to the standard model_, _Phys. Rev. D_ 101 (2020) 075047 [1911.05579].
  * [64] G.-C. Cho, C. Idegawa and E. Senaha, _Electroweak phase transition in a complex singlet extension of the Standard Model with degenerate scalars_ , _Phys. Lett. B_ 823 (2021) 136787 [2105.11830].
  * [65] P. Schicho, T. V. I. Tenkanen and G. White, _Combining thermal resummation and gauge invariance for electroweak phase transition_ , _JHEP_ 11 (2022) 047 [2203.04284].
  * [66] P. Di Bari and M. H. Rahat, _The split majoron model confronts the NANOGrav signal_ , 2307.03184.
  * [67] J. R. Espinosa and M. Quiros, _The Electroweak phase transition with a singlet_ , _Phys. Lett. B_ 305 (1993) 98 [hep-ph/9301285].
  * [68] S. Profumo, M. J. Ramsey-Musolf and G. Shaughnessy, _Singlet Higgs phenomenology and the electroweak phase transition_ , _JHEP_ 08 (2007) 010 [0705.2425].
  * [69] J. M. Cline, G. Laporte, H. Yamashita and S. Kraml, _Electroweak Phase Transition and LHC Signatures in the Singlet Majoron Model_ , _JHEP_ 07 (2009) 040 [0905.2559].
  * [70] J. M. Cline and K. Kainulainen, _Electroweak baryogenesis and dark matter from a singlet Higgs_ , _JCAP_ 01 (2013) 012 [1210.4196].
  * [71] J. M. Cline, K. Kainulainen, P. Scott and C. Weniger, _Update on scalar singlet dark matter_ , _Phys. Rev. D_ 88 (2013) 055025 [1306.4710].
  * [72] A. Katz and M. Perelstein, _Higgs Couplings and Electroweak Phase Transition_ , _JHEP_ 07 (2014) 108 [1401.1827].
  * [73] S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright and P. Winslow, _Singlet-catalyzed electroweak phase transitions and precision Higgs boson studies_ , _Phys. Rev. D_ 91 (2015) 035018 [1407.5342].
  * [74] V. Vaskonen, _Electroweak baryogenesis and gravitational waves from a real scalar singlet_ , _Phys. Rev. D_ 95 (2017) 123515 [1611.02073].
  * [75] K. Hashino, M. Kakizaki, S. Kanemura, P. Ko and T. Matsui, _Gravitational waves and Higgs boson couplings for exploring first order phase transition in the model with a singlet scalar field_ , _Phys. Lett. B_ 766 (2017) 49 [1609.00297].
  * [76] W. Chao, H.-K. Guo and J. Shu, _Gravitational Wave Signals of Electroweak Phase Transition Triggered by Dark Matter_ , _JCAP_ 09 (2017) 009 [1702.02698].
  * [77] T. Matsui, _Gravitational waves from the first order electroweak phase transition in the $Z_{3}$ symmetric singlet scalar model_, _EPJ Web Conf._ 168 (2018) 05001 [1709.05900].
  * [78] Z. Kang, P. Ko and T. Matsui, _Strong first order EWPT $\&$ strong gravitational waves in Z3-symmetric singlet scalar extension_, _JHEP_ 02 (2018) 115 [1706.09721].
  * [79] A. Alves, T. Ghosh, H.-K. Guo, K. Sinha and D. Vagie, _Collider and Gravitational Wave Complementarity in Exploring the Singlet Extension of the Standard Model_ , _JHEP_ 04 (2019) 052 [1812.09333].
  * [80] V. R. Shajiee and A. Tofighi, _Electroweak Phase Transition, Gravitational Waves and Dark Matter in Two Scalar Singlet Extension of The Standard Model_ , _Eur. Phys. J. C_ 79 (2019) 360 [1811.09807].
  * [81] A. Beniwal, M. Lewicki, M. White and A. G. Williams, _Gravitational waves and electroweak baryogenesis in a global study of the extended scalar singlet model_ , _JHEP_ 02 (2019) 183 [1810.02380].
  * [82] T. Matsui, _Gravitational waves from the first order electroweak phase transition in the $Z_{3}$ symmetric singlet scalar model_, _PoS_ CORFU2017 (2018) 092.
  * [83] K. Kannike, K. Loos and M. Raidal, _Gravitational wave signals of pseudo-Goldstone dark matter in the $\mathbb{Z}_{3}$ complex singlet model_, _Phys. Rev. D_ 101 (2020) 035001 [1907.13136].
  * [84] S. Demidov, D. Gorbunov and E. Kriukova, _Gravitational waves from first-order electroweak phase transition in a model with light sgoldstinos_ , _JHEP_ 07 (2022) 061 [2112.06083].
  * [85] 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) 123529 [2102.12490].
  * [86] Q.-H. Cao, K. Hashino, X.-X. Li and J.-H. Yue, _Multi-step phase transition and gravitational wave from general $\mathbb{Z}_{2}$ scalar extensions_, 2212.07756.
  * [87] P. Athron, C. Balázs, A. Fowlie, L. Morris and L. Wu, _Cosmological phase transitions: From perturbative particle physics to gravitational waves_ , _Prog. Part. Nucl. Phys._ 135 (2024) 104094 [2305.02357].
  * [88] T. Alanne, N. Benincasa, M. Heikinheimo, K. Kannike, V. Keus, N. Koivunen et al., _Pseudo-Goldstone dark matter: gravitational waves and direct-detection blind spots_ , _JHEP_ 10 (2020) 080 [2008.09605].
  * [89] L. Bian, Y. Wu and K.-P. Xie, _Electroweak phase transition with composite Higgs models: calculability, gravitational waves and collider searches_ , _JHEP_ 12 (2019) 028 [1909.02014].
  * [90] K.-P. Xie, L. Bian and Y. Wu, _Electroweak baryogenesis and gravitational waves in a composite Higgs model with high dimensional fermion representations_ , _JHEP_ 12 (2020) 047 [2005.13552].
  * [91] C. Grojean, G. Servant and J. D. Wells, _First-order electroweak phase transition in the standard model with a low cutoff_ , _Phys. Rev. D_ 71 (2005) 036001 [hep-ph/0407019].
  * [92] D. Bodeker, L. Fromme, S. J. Huber and M. Seniuch, _The Baryon asymmetry in the standard model with a low cut-off_ , _JHEP_ 02 (2005) 026 [hep-ph/0412366].
  * [93] C. Delaunay, C. Grojean and J. D. Wells, _Dynamics of Non-renormalizable Electroweak Symmetry Breaking_ , _JHEP_ 04 (2008) 029 [0711.2511].
  * [94] C. Balazs, G. White and J. Yue, _Effective field theory, electric dipole moments and electroweak baryogenesis_ , _JHEP_ 03 (2017) 030 [1612.01270].
  * [95] J. de Vries, M. Postma, J. van de Vis and G. White, _Electroweak Baryogenesis and the Standard Model Effective Field Theory_ , _JHEP_ 01 (2018) 089 [1710.04061].
  * [96] J. De Vries, M. Postma and J. van de Vis, _The role of leptons in electroweak baryogenesis_ , _JHEP_ 04 (2019) 024 [1811.11104].
  * [97] M. Chala, C. Krause and G. Nardini, _Signals of the electroweak phase transition at colliders and gravitational wave observatories_ , _JHEP_ 07 (2018) 062 [1802.02168].
  * [98] S. A. R. Ellis, S. Ipek and G. White, _Electroweak Baryogenesis from Temperature-Varying Couplings_ , _JHEP_ 08 (2019) 002 [1905.11994].
  * [99] M. Postma and G. White, _Cosmological phase transitions: is effective field theory just a toy?_ , _JHEP_ 03 (2021) 280 [2012.03953].
  * [100] B. Batell, A. Ghalsasi, M. Low and M. Rai, _Gravitational Waves from Nnaturalness_ , _JHEP_ 01 (2024) 148 [2310.06905].
  * [101] T. Ghosh, A. Ghoshal, H.-K. Guo, F. Hajkarim, S. F. King, K. Sinha et al., _Did we hear the sound of the Universe boiling? Analysis using the full fluid velocity profiles and NANOGrav 15-year data_ , 2307.02259.
  * [102] P. Schwaller, _Gravitational Waves from a Dark Phase Transition_ , _Phys. Rev. Lett._ 115 (2015) 181101 [1504.07263].
  * [103] I. Baldes and C. Garcia-Cely, _Strong gravitational radiation from a simple dark matter model_ , _JHEP_ 05 (2019) 190 [1809.01198].
  * [104] 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 [1811.11175].
  * [105] D. Croon, V. Sanz and G. White, _Model Discrimination in Gravitational Wave spectra from Dark Phase Transitions_ , _JHEP_ 08 (2018) 203 [1806.02332].
  * [106] E. Hall, T. Konstandin, R. McGehee, H. Murayama and G. Servant, _Baryogenesis From a Dark First-Order Phase Transition_ , _JHEP_ 04 (2020) 042 [1910.08068].
  * [107] I. Baldes, _Gravitational waves from the asymmetric-dark-matter generating phase transition_ , _JCAP_ 05 (2017) 028 [1702.02117].
  * [108] D. Croon, A. Kusenko, A. Mazumdar and G. White, _Solitosynthesis and Gravitational Waves_ , _Phys. Rev. D_ 101 (2020) 085010 [1910.09562].
  * [109] E. Hall, T. Konstandin, R. McGehee and H. Murayama, _Asymmetric matter from a dark first-order phase transition_ , _Phys. Rev. D_ 107 (2023) 055011 [1911.12342].
  * [110] W. Chao, X.-F. Li and L. Wang, _Filtered pseudo-scalar dark matter and gravitational waves from first order phase transition_ , _JCAP_ 06 (2021) 038 [2012.15113].
  * [111] J. B. Dent, B. Dutta, S. Ghosh, J. Kumar and J. Runburg, _Sensitivity to dark sector scales from gravitational wave signatures_ , _JHEP_ 08 (2022) 300 [2203.11736].
  * [112] A. J. Helmboldt, J. Kubo and S. van der Woude, _Observational prospects for gravitational waves from hidden or dark chiral phase transitions_ , _Phys. Rev. D_ 100 (2019) 055025 [1904.07891].
  * [113] M. Aoki and J. Kubo, _Gravitational waves from chiral phase transition in a conformally extended standard model_ , _JCAP_ 04 (2020) 001 [1910.05025].
  * [114] D. Croon, J. N. Howard, S. Ipek and T. M. P. Tait, _QCD baryogenesis_ , _Phys. Rev. D_ 101 (2020) 055042 [1911.01432].
  * [115] D. Croon, R. Houtz and V. Sanz, _Dynamical Axions and Gravitational Waves_ , _JHEP_ 07 (2019) 146 [1904.10967].
  * [116] J. Garcia-Bellido, H. Murayama and G. White, _Exploring the early Universe with Gaia and Theia_ , _JCAP_ 12 (2021) 023 [2104.04778].
  * [117] W.-C. Huang, M. Reichert, F. Sannino and Z.-W. Wang, _Testing the dark SU(N) Yang-Mills theory confined landscape: From the lattice to gravitational waves_ , _Phys. Rev. D_ 104 (2021) 035005 [2012.11614].
  * [118] J. Halverson, C. Long, A. Maiti, B. Nelson and G. Salinas, _Gravitational waves from dark Yang-Mills sectors_ , _JHEP_ 05 (2021) 154 [2012.04071].
  * [119] Z. Kang, J. Zhu and S. Matsuzaki, _Dark confinement-deconfinement phase transition: a roadmap from Polyakov loop models to gravitational waves_ , _JHEP_ 09 (2021) 060 [2101.03795].
  * [120] B. Fornal and E. Pierre, _Asymmetric dark matter from gravitational waves_ , _Phys. Rev. D_ 106 (2022) 115040 [2209.04788].
  * [121] F. Costa, S. Khan and J. Kim, _A two-component dark matter model and its associated gravitational waves_ , _JHEP_ 06 (2022) 026 [2202.13126].
  * [122] K. Hashino, M. Kakizaki, S. Kanemura, P. Ko and T. Matsui, _Gravitational waves from first order electroweak phase transition in models with the U(1) X gauge symmetry_, _JHEP_ 06 (2018) 088 [1802.02947].
  * [123] F. P. Huang and X. Zhang, _Probing the gauge symmetry breaking of the early universe in 3-3-1 models and beyond by gravitational waves_ , _Phys. Lett. B_ 788 (2019) 288 [1701.04338].
  * [124] D. Croon, T. E. Gonzalo and G. White, _Gravitational Waves from a Pati-Salam Phase Transition_ , _JHEP_ 02 (2019) 083 [1812.02747].
  * [125] V. Brdar, L. Graf, A. J. Helmboldt and X.-J. Xu, _Gravitational Waves as a Probe of Left-Right Symmetry Breaking_ , _JCAP_ 12 (2019) 027 [1909.02018].
  * [126] W.-C. Huang, F. Sannino and Z.-W. Wang, _Gravitational Waves from Pati-Salam Dynamics_ , _Phys. Rev. D_ 102 (2020) 095025 [2004.02332].
  * [127] B. Fornal, K. Garcia and E. Pierre, _Testing unification and dark matter with gravitational waves_ , _Phys. Rev. D_ 108 (2023) 055022 [2305.12566].
  * [128] P. Costa, C. A. de Sousa, M. C. Ruivo and H. Hansen, _The QCD critical end point in the PNJL model_ , _EPL_ 86 (2009) 31001 [0801.3616].
  * [129] F. Marquez and R. Zamora, _Critical end point in a thermomagnetic nonlocal NJL model_ , _Int. J. Mod. Phys. A_ 32 (2017) 1750162 [1702.04161].
  * [130] V. Vovchenko, M. I. Gorenstein, C. Greiner and H. Stoecker, _Hagedorn bag-like model with a crossover transition meets lattice QCD_ , _Phys. Rev. C_ 99 (2019) 045204 [1811.05737].
  * [131] R. Câmara Pereira, J. a. Moreira and P. Costa, _The strange critical endpoint and isentropic trajectories in an extended PNJL model with eight Quark interactions_ , _Eur. Phys. J. A_ 56 (2020) 214 [2006.02385].
  * [132] F. Gao and J. M. Pawlowski, _QCD phase structure from functional methods_ , _Phys. Rev. D_ 102 (2020) 034027 [2002.07500].
  * [133] F. Gao and J. M. Pawlowski, _Chiral phase structure and critical end point in QCD_ , _Phys. Lett. B_ 820 (2021) 136584 [2010.13705].
  * [134] F. Gao and I. M. Oldengott, _Cosmology Meets Functional QCD: First-Order Cosmic QCD Transition Induced by Large Lepton Asymmetries_ , _Phys. Rev. Lett._ 128 (2022) 131301 [2106.11991].
  * [135] F. Gao, J. Harz, C. Hati, Y. Lu, I. M. Oldengott and G. White, _Sphaleron freeze-in baryogenesis with gravitational waves from the QCD transition_ , 2309.00672.
  * [136] F. Bigazzi, A. Caddeo, A. L. Cotrone and A. Paredes, _Fate of false vacua in holographic first-order phase transitions_ , _JHEP_ 12 (2020) 200 [2008.02579].
  * [137] F. Bigazzi, A. Caddeo, A. L. Cotrone and A. Paredes, _Dark Holograms and Gravitational Waves_ , _JHEP_ 04 (2021) 094 [2011.08757].
  * [138] X. Wang, F. P. Huang and X. Zhang, _Bubble wall velocity beyond leading-log approximation in electroweak phase transition_ , 2011.12903.
  * [139] E. Morgante, N. Ramberg and P. Schwaller, _Gravitational waves from dark SU(3) Yang-Mills theory_ , _Phys. Rev. D_ 107 (2023) 036010 [2210.11821].
  * [140] 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].
  * [141] M. Hindmarsh and M. Hijazi, _Gravitational waves from first order cosmological phase transitions in the Sound Shell Model_ , _JCAP_ 12 (2019) 062 [1909.10040].
  * [142] 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].
  * [143] 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) 123513 [1802.05712].
  * [144] D. Cutting, E. G. Escartin, M. Hindmarsh and D. J. Weir, _Gravitational waves from vacuum first order phase transitions II: from thin to thick walls_ , _Phys. Rev. D_ 103 (2021) 023531 [2005.13537].
  * [145] O. Gould, S. Sukuvaara and D. Weir, _Vacuum bubble collisions: From microphysics to gravitational waves_ , _Phys. Rev. D_ 104 (2021) 075039 [2107.05657].
  * [146] D. J. Weir, _Gravitational waves from a first order electroweak phase transition: a brief review_ , _Phil. Trans. Roy. Soc. Lond. A_ 376 (2018) 20170126 [1705.01783].
  * [147] C. Gowling and M. Hindmarsh, _Observational prospects for phase transitions at LISA: Fisher matrix analysis_ , _JCAP_ 10 (2021) 039 [2106.05984].
  * [148] C. Gowling, M. Hindmarsh, D. C. Hooper and J. Torrado, _Reconstructing physical parameters from template gravitational wave spectra at LISA: first order phase transitions_ , _JCAP_ 04 (2023) 061 [2209.13551].
  * [149] F. Giese, T. Konstandin and J. van de Vis, _Finding sound shells in LISA mock data using likelihood sampling_ , _JCAP_ 11 (2021) 002 [2107.06275].
  * [150] P. Athron, C. Balázs, T. E. Gonzalo and M. Pearce, _Falsifying Pati-Salam models with LIGO_ , _Phys. Rev. D_ 109 (2024) L061303 [2307.02544].
  * [151] R. Roshan and G. White, _Using gravitational waves to see the first second of the Universe_ , 2401.04388.
  * [152] A. Mazumdar and G. White, _Review of cosmic phase transitions: their significance and experimental signatures_ , _Rept. Prog. Phys._ 82 (2019) 076901 [1811.01948].
  * [153] 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].
  * [154] 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) 123009 [1504.03291].
  * [155] 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].
  * [156] 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].
  * [157] 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].
  * [158] 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].
  * [159] F. Giese, T. Konstandin, K. Schmitz and J. van de Vis, _Model-independent energy budget for LISA_ , _JCAP_ 01 (2021) 072 [2010.09744].
  * [160] T. Krajewski, M. Lewicki and M. Zych, _Bubble-wall velocity in local thermal equilibrium: hydrodynamical simulations vs analytical treatment_ , 2402.15408.
  * [161] S. J. Huber and T. Konstandin, _Gravitational Wave Production by Collisions: More Bubbles_ , _JCAP_ 09 (2008) 022 [0806.1828].
  * [162] 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].
  * [163] 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].
  * [164] M. Drees, F. Hajkarim and E. R. Schmitz, _The Effects of QCD Equation of State on the Relic Density of WIMP Dark Matter_ , _JCAP_ 06 (2015) 025 [1503.03513].
  * [165] A. Roper Pol, S. Procacci and C. Caprini, _Characterization of the gravitational wave spectrum from sound waves within the sound shell model_ , 2308.12943.
  * [166] A. Roper Pol, A. Neronov, C. Caprini, T. Boyer and D. Semikoz, _LISA and $\gamma$-ray telescopes as multi-messenger probes of a first-order cosmological phase transition_, 2307.10744.
  * [167] R. Sharma, J. Dahl, A. Brandenburg and M. Hindmarsh, _Shallow relic gravitational wave spectrum with acoustic peak_ , _JCAP_ 12 (2023) 042 [2308.12916].
  * [168] K. Schmitz, _New Sensitivity Curves for Gravitational-Wave Signals from Cosmological Phase Transitions_ , _JHEP_ 01 (2021) 097 [2002.04615].
  * [169] R. Jinno, T. Konstandin, H. Rubira and J. van de Vis, _Effect of density fluctuations on gravitational wave production in first-order phase transitions_ , _JCAP_ 12 (2021) 019 [2108.11947].
  * [170] H.-K. Guo, _Dissipative Effects as New Observables for Cosmological Phase Transitions_ , 2310.10927.