###### Abstract

An extension of the relativistic density functional approach to the equation
of state for strongly interacting matter is suggested that generalizes a
recently developed modified excluded-volume mechanism to the case of
temperature- and density-dependent available-volume fractions. A
parametrization of this dependence is presented for which, at low temperatures
and suprasaturation densities, a first-order phase transition is obtained. It
changes for increasing temperatures to a crossover transition via a critical
endpoint. This provides a benchmark case for studies of the role of such a
point in hydrodynamic simulations of ultrarelativistic heavy-ion collisions.
The approach is thermodynamically consistent and extendable to finite isospin
asymmetries that are relevant for simulations of neutron stars, their mergers,
and core-collapse supernova explosions.

###### keywords:

Equation of state; QCD matter; phase transition; critical point; modified
excluded-volume mechanism

x 10.3390/—— xx Academic Editor: name Received: date; Accepted: date;
Published: date Phenomenological Equation of State of Strongly Interacting
Matter with First-Order Phase Transitions and Critical Points Stefan Typel
1,2* and David Blaschke 3,4,5 Stefan Typel and David Blaschke Correspondence:
<EMAIL_ADDRESS>Tel.: +49-6159-71-2744

##  1 Introduction

The simulation of astrophysical phenomena, such as core-collapse supernovae
(CCSN) or neutron-star (NS) mergers, requires a careful modeling of strongly
interacting matter in a wide range of densities and temperatures. The same
applies to the theoretical description of heavy-ion collisions (HIC) that
study compressed baryonic matter in the laboratory from low to high beam
energies. The properties of such matter are represented by the equation of
state (EoS) that provides information on pressure, entropy, energies, and
other thermodynamic variables of interest.

A particular feature of QCD matter is the supposed phase transition (PT) from
hadronic matter to quark matter when density or temperature increase to
sufficiently high values. A strong first-order PT could allow for the
existence of a third branch of compact stars and the occurrence of the twin-
star phenomenon Alford et al. (2013); Benic et al. (2015); Alvarez-Castillo
and Blaschke (2017). Signals of the PT might also have direct consequences in
dynamical processes when matter in the quark phase expands and cools down,
e.g., the release of a second neutrino burst in CCSN Sagert et al. (2009);
Fischer et al. (2009, 2011). For a recent review on the role of the EoS in
CCSN simulations, see Fischer et al. (2017).

The theoretical description of the hadron–quark PT in strongly interacting
matter often relies on a construction employing different models for the two
phases. With such an approach, the coexistence line of the first-order PT will
usually connect a point on the zero-temperature axis at finite baryon chemical
potential $\mu_{B}$ with a point at finite temperature $T$ on the zero baryon
chemical potential axis. By this construction, the QCD hadron–quark PT is of
first order in the whole temperature–density plane, see for instance Blaschke
et al. (2010). However, from lattice QCD studies, it is known that there is a
smooth crossover at $\mu_{B}=0$ with increasing $T$ Bazavov et al. (2014);
Borsanyi et al. (2014), so at least one critical point at finite $\mu_{B}$ and
$T$ is expected to exist. Other possibilities are that the character of the
transition is crossover all over the QCD phase diagram Bratovic et al. (2013)
or, as is advocated in studies of the BEC-BCS crossover transition in low-
temperature QCD, that a second critical endpoint exists Baym et al. (2017);
Abuki et al. (2010); Hatsuda et al. (2006). Since lattice QCD studies are
presently incapable of exploring the EoS close to the presumed critical point
with much confidence, unified models are needed that can account for the
existence of these features, see, e.g., Klähn and Fischer (2015). There are
dedicated microscopic models available that incorporate the major expected
features in the QCD phase diagram, e.g., chiral mean-field models Dexheimer
and Schramm (2010) or parity-doublet quark–hadron models Mukherjee et al.
(2017). Simulations of CCSN or HIC that are based on a hydrodynamic
description of matter during dynamical evolution use the thermodynamic
properties of matter encoded in the EoS as an input. Such data can be provided
by phenomenological models that do not need to incorporate all of the details
of the underlying physics.

In this work, a novel approach is introduced to provide a phenomenological EoS
of baryonic matter that exhibits a first-order PT and a critical point at
densities and temperatures expected in QCD matter. The parameters of the model
can be adjusted to place the coexistence line at arbitrary positions in the
phase diagram. The description uses an extension of a relativistic energy
density functional for hadronic matter assuming a medium-dependent change in
the number of degrees of freedom. This approach employs a recently developed
version of a modified excluded-volume (EV) mechanism that gives a
thermodynamically consistent EoS with nuclear matter properties that are
consistent with present constraints. Here, we concentrate on the hadron–quark
transition but not on the liquid–gas PT, which is also contained in our model.
The model allows us to study the PT for arbitrary isospin asymmetries;
however, only isospin-symmetric matter is considered in this first exploratory
study for simplicity. In the present work, no attempt was made to reproduce
the EoS of QCD matter at vanishing baryon chemical potential obtained in
lattice QCD studies. With appropriately chosen EV parameters, the crossover
transition with increasing temperature can be well modeled, even for imaginary
chemical potentials, e.g., in a hadron resonance gas model Vovchenko et al.
(2017). With improved parametrizations, the structure of the phase diagram in
the full space of variables, i.e., temperature, baryon density/chemical
potential, and isospin asymmetry, can be investigated in the future.

The theoretical formalism of the model is presented in Section 2, which
includes the main equations that define the relevant thermodynamic quantities
in Section 2.1. In Section 2.2, details of the parametrization of the
interaction and of the effective degeneracy factors are given. They account
for the change in the number of degrees of freedom with density and
temperature. The phase transitions are explored in Section 3 for isospin-
symmetric matter. Conclusions follow in Section 4.

##  2 Theoretical Model

The theoretical description of strongly interacting matter in the present work
is adapted from the model introduced in Typel (2016). It combines a
relativistic mean-field (RMF) approach for hadronic matter with density-
dependent nucleon–meson couplings and a modified EV mechanism. Here, it is
sufficient to provide only the main equations without a detailed derivation.
The essential quantities that determine the position of the PT and the
critical point in the phase diagram are the effective degeneracy factors that
depend on the number densities of the particles and the temperature.

###  2.1 Relativistic Energy Density Functional with Modified Excluded-Volume
Mechanism

The present model assumes neutrons and protons as well as their antiparticles
as the basic degrees of freedom. These particles interact by the exchange of
mesons, and the model effectively describes the short-range repulsion
($\omega$ meson), the intermediate-range attraction ($\sigma$ meson), and the
isospin dependence of the nuclear interaction ($\rho$ and $\delta$ mesons), as
is common of RMF models. The contribution of leptons or other degrees of
freedom like nuclei, hyperons or photons, as required for multi-purpose EoS
for astrophysical applications Oertel et al. (2017), is not considered here.

The nucleons $i=n,p,\bar{n},\bar{p}$ with rest masses $m_{i}$ are treated as
quasi-particles of energy,

$E_{i}(k)=\sqrt{k^{2}+\left(m_{i}-S_{i}\right)^{2}}+V_{i}$ (1)

which depends on the particle momentum $k$ and the scalar ($S_{i}$) and vector
($V_{i}$) potentials. Denoting the particle chemical potentials with
$\mu_{i}$, the contribution of the quasi-particles to the total pressure

$p=\sum_{i}p_{i}+p_{\rm meson}-p^{(r)}$ (2)

of the system can be written as

$p_{i}=Tg_{i}^{\rm(eff)}\int\frac{d^{3}k}{(2\pi)^{3}}\>\ln\left[1+\exp\left(-\frac{E_{i}(k)-\mu_{i}}{T}\right)\right]$
(3)

where the medium-dependent effective degeneracy factors

$g_{i}^{\rm(eff)}=g_{i}\Phi_{i}$ (4)

are a product of the usual degeneracy factor $g_{i}=2$ for nucleons and the
available-volume fraction $\Phi_{i}$, which is defined in Section 2.2.

The meson contribution

$p_{\rm
meson}=\frac{1}{2}\left(C_{\omega}n_{\omega}^{2}+C_{\rho}n_{\rho}^{2}-C_{\sigma}n_{\sigma}^{2}-C_{\delta}n_{\delta}^{2}\right)$
(5)

to the total pressure in Equation (2) contains the coupling factors of the
mesons

$C_{j}=\frac{\Gamma_{j}^{2}}{m_{j}^{2}}$ (6)

given as a ratio of the density-dependent coupling functions $\Gamma_{j}$ and
the meson masses $m_{j}$. The source densities

$n_{j}=\sum_{i}g_{ij}n_{i}^{(v)}$ (7)

for vector mesons ($j=\omega,\rho$) and

$n_{j}=\sum_{i}g_{ij}n_{i}^{(s)}$ (8)

for scalar mesons ($j=\sigma,\delta$) in Equation (5) are obtained from the
quasi-particle vector densities

$n_{i}^{(v)}=g_{i}^{\rm(eff)}\int\frac{d^{3}k}{(2\pi)^{3}}\>f_{i}(k)$ (9)

and scalar densities

$n_{i}^{(s)}=g_{i}^{\rm(eff)}\int\frac{d^{3}k}{(2\pi)^{3}}\>f_{i}(k)\frac{m_{i}-S_{i}}{\sqrt{k^{2}+\left(m_{i}-S_{i}\right)^{2}}}$
(10)

with the Fermi-Dirac distribution function

$f_{i}(k)=\left[\exp\left(\frac{E_{i}(k)-\mu_{i}}{T}\right)+1\right]^{-1}\>.$
(11)

The scaling factors

$\displaystyle g_{n\omega}=g_{p\omega}=-g_{\bar{n}\omega}=-g_{\bar{p}\omega}$
$\displaystyle=$ $\displaystyle 1$ (12) $\displaystyle
g_{n\rho}=-g_{p\rho}=-g_{\bar{n}\rho}=g_{\bar{p}\rho}$ $\displaystyle=$
$\displaystyle 1$ (13) $\displaystyle
g_{n\sigma}=g_{p\sigma}=g_{\bar{n}\sigma}=g_{\bar{p}\sigma}$ $\displaystyle=$
$\displaystyle 1$ (14) $\displaystyle
g_{n\delta}=-g_{p\delta}=g_{\bar{n}\delta}=-g_{\bar{p}\delta}$
$\displaystyle=$ $\displaystyle 1$ (15)

in Equations (7) and (8) determine the coupling between mesons and nucleons.
They also appear in the vector potential

$V_{i}=C_{\omega}g_{i\omega}n_{\omega}+C_{\rho}g_{i\rho}n_{\rho}+B_{i}V_{\rm
meson}^{(r)}+V_{i}^{(r)}$ (16)

and the scalar potential

$S_{i}=C_{\sigma}g_{i\sigma}n_{\sigma}+C_{\delta}g_{i\delta}n_{\delta}+S_{i}^{(r)}$
(17)

in the quasi-particle energy (Equation (1)). The rearrangement potential

$V_{\rm
meson}^{(r)}=\frac{1}{2}\left(C_{\omega}^{\prime}n_{\omega}^{2}+C_{\rho}^{\prime}n_{\rho}^{2}-C_{\sigma}^{\prime}n_{\sigma}^{2}-C_{\delta}^{\prime}n_{\delta}^{2}\right)$
(18)

contributes to the vector potential (Equation (16)) because the couplings
$\Gamma_{j}$ in Equation (6) are assumed to depend on the baryon density
$n_{B}=\sum_{i}B_{i}n_{i}^{(v)}$, where $B_{i}=g_{i\omega}$ is the baryon
number of particle $i$, and the quantities $C_{j}^{\prime}=dC_{j}/dn_{B}$ are
the derivatives of the coupling factors.

The dependence of the available-volume fractions $\Phi_{i}$ in the effective
degeneracy factor (Equation (4)) on the vector or scalar quasi-particle
densities (9) and (10) also generates rearrangement contributions

$V_{i}^{(r)}=-\sum_{j}p_{j}\frac{\partial\ln\Phi_{j}}{\partial n_{i}^{(v)}}$
(19)

and

$S_{i}^{(r)}=\sum_{j}p_{j}\frac{\partial\ln\Phi_{j}}{\partial n_{i}^{(s)}}$
(20)

in the potentials (Equations (16) and (17)), respectively. Furthermore, there
is a rearrangement term

$p^{(r)}=p^{(r)}_{\rm meson}+p^{(r)}_{\Phi}$ (21)

in the total pressure (Equation (2)) with two contributions from the density
dependence of the couplings

$p^{(r)}_{\rm meson}=-V^{(r)}_{\rm meson}n_{B}$ (22)

and

$p^{(r)}_{\Phi}=\sum_{i}\left(n_{i}^{(s)}S_{i}^{(r)}-n_{i}^{(v)}V_{i}^{(r)}\right)$
(23)

from the EV effects.

The free energy density of the system

$f=\sum_{i}\mu_{i}n_{i}^{(v)}-p$ (24)

is obtained with the total pressure (Equation (2)) and the chemical potentials
of the particles $\mu_{i}$. They are not independent since, for nucleons with
baryon number $B_{i}$ and charge number $Q_{i}$, they are given by

$\mu_{i}=B_{i}\mu_{B}+Q_{i}\mu_{Q}$ (25)

with the baryon chemical potential $\mu_{B}$ and the charge chemical potential
$\mu_{Q}$. Only the latter two are independent quantities. For the internal
energy density

$\varepsilon=f+Ts$ (26)

the entropy density

$s=-\sum_{i}g_{i}^{\rm(eff)}\int\frac{d^{3}k}{(2\pi)^{3}}\>\left[f_{i}\ln
f_{i}+\left(1-f_{i}\right)\ln\left(1-f_{i}\right)\right]+\sum_{i}p_{i}\frac{\partial\ln\Phi_{i}}{\partial
T}$ (27)

is needed. Besides the standard contribution depending on the distribution
functions $f_{i}$, there is a term from the possible temperature dependence of
the available-volume fractions $\Phi_{i}$. In order to guarantee the third law
of thermodynamics, i.e., $\lim_{T\to 0}s=0$, the temperature derivative of the
available-volume fractions has to vanish for $T\to 0$. After solving the
equations above for a given $T$, $\mu_{B}$, and $\mu_{Q}$, a fully consistent
thermodynamic EoS is obtained. For practical purposes, however, the baryon
density $n_{B}$ and the hadronic charge fraction

$Y_{q}=\frac{\sum_{i}Q_{i}n_{i}^{(v)}}{n_{B}}$ (28)

are used as independent variables instead of $\mu_{B}$ and $\mu_{Q}$.

A possible shortcoming of models that consider EV effects is the potential
appearance of a superluminal speed of sound in certain regions of the space of
thermodynamic variables, see, e.g., Rischke et al. (1991); Satarov et al.
(2009). This causality constraint has to be checked case by case depending on
the specific implementation of the EV mechanism.

###  2.2 Available-Volume Fractions and Model Parameters

For a quantitative evaluation of the EoS in the present approach, the
functional forms of the meson–nucleon couplings $\Gamma_{j}$ and the
available-volume fractions $\Phi_{i}$ have to be specified as well as all
parameters. Here, we use the masses of nucleons and mesons and the coupling
functions of the DD2 parametrization presented in Typel et al. (2010). It only
considers $\omega$, $\rho$, and $\sigma$ mesons for the effective description
of the nuclear interaction but not the $\delta$ mesons. The parameters were
obtained by fitting observables (binding energies, radii, etc.) of selected
nuclei. With this set, the EoS of nuclear matter at zero temperature exhibits
characteristic nuclear matter parameters, e.g., the saturation density
($n_{\rm sat}=0.149065$ fm-3), binding energy at saturation ($B=16.02$ MeV),
incompressibility ($K=242.7$ MeV), symmetry energy ($J=32.73$ MeV), and slope
($L=57.94$ MeV), that are consistent with modern constraints from experiment
and theory.

EV effects are frequently employed to describe an effective repulsive
interaction between particles, in particular in calculations of the EoS in the
framework of hadron resonance gas models, see, e.g., Vovchenko et al. (2017).
For zero baryon density, a comparison of the resulting EoS with results from
lattice QCD calculations can be used to fix the volume parameters. If a finite
volume $v_{i}$ is attributed to each particle $i$, the available volume for
the motion of the particle is reduced from the total system volume $V$ to
$V\Phi^{(cl)}$ with the classical available-volume fraction

$\Phi^{(cl)}=1-\sum_{i}v_{i}n_{i}^{(v)}\>.$ (29)

Clearly, there is a limiting density above which a compression of the system
becomes impossible. In general, the volumes and available-volume fractions can
depend on the particle species, and the EV mechanism can be used to suppress
particles, e.g., nuclei, in a mixture, causing them to disappear above a
certain density, see, e.g., Hempel and Schaffner-Bielich (2010) for
applications to the low-temperature and low-density EoS in astrophysical
simulations.

In Typel (2016), the interpretation of the EV mechanism was changed by moving
the factor $\Phi_{i}$ from the system volume $V$ to the degeneracy factor
$g_{i}$ as in Equation (4) and allowing the available-volume fractions to be
arbitrary functions of the particle densities and temperature. The medium
dependence of the effective degeneracy factors is interpreted as a change in
the effective number of degrees of freedom. A decrease in $g_{i}^{\rm(eff)}$
has the effect of a repulsive interaction between the particles, whereas an
increase can be seen as the action of an attractive interaction, c.f., the
softening of the nuclear EoS when hyperons are included, see, e.g., Chatterjee
and Vidaña (2016); Schaffner-Bielich (2008). This freedom leaves the room to
modify the properties of an EoS in a favored way.

In the present application of the modified EV mechanism, the available-volume
fraction is defined to be identical for all particles $i$ as

$\Phi_{i}=1+sg_{1}(T)\theta(x)\exp\left(-\frac{1}{2x^{2}}\right)$ (30)

depending on the temperature $T$ and an auxiliary quantity

$x=v\left(\sum_{j}B_{j}n_{j}^{(v)}-g_{8}(T)n_{\rm cut}\right)$ (31)

depending on $T$ and the quasi-particle vector densities $n_{i}^{(v)}$ with
parameters $s$, $v$, and a cutoff density of $n_{\rm cut}$. The $\theta$
function in Equation (30) guarantees that $\Phi_{i}=1$ for $x\leq 0$ and the
EV mechanism has no effect on the EoS. The functions $g_{1}$ in Equation (30)
and $g_{8}$ in Equation (31) are defined as

$g_{t}(T)=\theta\left(T_{0}-T\right)\exp\left[-\frac{t}{2}\left(\frac{T}{T_{0}-T}\right)^{2}\right]$
(32)

with parameters $t$ and $T_{0}$. In the limit $T\to 0$, the function $g_{t}$
approaches one and it decreases with increasing temperature. Furthermore, the
derivatives $\partial\Phi_{i}/\partial T$ approach zero for $T\to 0$, as
required for the thermodynamic consistency, because of the choice of the
function $g_{t}(T)$. For $T\to T_{0}$, the function $g_{t}$ vanishes very
smoothly and there are no effects at higher temperatures because $\Phi_{i}=1$.
In order to reproduce the correct high-temperature limit, given by a
Stefan–Boltzmann-type behavior, a modification of the available-volume
fractions for temperatures well above $T_{0}$ is required. This is left to
future extensions of the model. According to Equation (31), the quantity $x$
is only positive for baryon densities larger than $g_{8}n_{\rm cut}$, a
quantity that decreases with increasing temperature. There are no artificial
singularities due to the presence of the $\theta$ functions in Equations (30)
and (32) because all derivatives of the exponential functions are zero when
the arguments of the $\theta$ functions vanish. The actual values of the
parameters for the modified EV mechanism used in the present study are $s=3$,
$v=2$ fm3, $T_{0}=270$ MeV, and $n_{\rm cut}=n_{\rm sat}$ of the DD2
parametrization.

##  3 Results

In order to illustrate the characteristic effects of the modified EV mechanism
on the EoS, we limit the presentation to the case of symmetric matter, i.e.,
$Y_{q}=0.5$. The pressure $p$ and baryon chemical potential $\mu_{B}$ are
calculated as a function of the baryon density $n_{B}$. Due to the increase in
the available-volume fractions $\Phi_{i}$ or the effective degeneracy factors
$g_{i}^{\rm(eff)}$, a considerable softening of the EoS is observed in a
certain range of densities. Below the critical temperature $T_{\rm crit}$, the
pressure is not a monotonous function of the baryon chemical potential for an
isotherm. A Maxwell construction is used to determine the two densities of the
coexisting phases where $p$ and $\mu_{B}$ are identical. The pressure $p$ as a
function of the baryon density $n_{B}$ is depicted in Figure 1 for selected
temperatures. Note that an integral part of the underlying RMF model is a
detailed description of the liquid–gas PT in nuclear matter and the formation
and dissociation of nuclear clusters in compact-star matter. For details, see,
e.g., Typel et al. (2014).

Figure 1:  Isotherms in isospin-symmetric strongly interacting matter in the
pressure–baryon density diagram at temperatures from $0$ to $160$ MeV in steps
of $20$ MeV (dashed colored lines) and at the critical temperature $T_{\rm
crit}$ of the pseudo hadron–quark phase transition (black dot-dashed line).
The binodals and critical points are denoted by full black lines and a full
(open) circle of the pseudo hadron–quark (liquid-gas) phase transition,
respectively.

In the coexistence region of the pseudo hadron–quark PT between the low and
high density phases, the pressure is constant as typical for a first-order PT.
The area of coexistence is enclosed by the binodal, and, above the critical
temperature $T_{\rm crit}\approx 155.5$ MeV, there is no PT anymore. The
peculiar shape of the binodal is a result of the specific form (30) of the
available-volume fractions. It can be adjusted with appropriate changes in the
functional form and parameters.

The binodals of the liquid–gas and pseudo hadron–quark PT in the
temperature–baryon density plane are shown in panel (a) of Figure 2. At
vanishing temperatures, the coexistence region of the pseudo hadron–quark PT
covers a density range from $0.270$ to $0.408$ fm-3, well above the nuclear
saturation density $n_{\rm sat}$. At higher temperatures, it moves to lower
densities with an almost constant extension in baryon density except for
temperatures close to $T_{\rm crit}$. Here, the critical density is found as
$0.201$ fm-3, still above $n_{\rm sat}$. The dashed line in panel (a) marks
the boundary between regions without (lower left) and with (upper right)
effects of the modified EV mechanism in the present parametrization. It
corresponds to the condition $x=0$. There is another region in the phase
diagram without modified EV effects at temperatures above $T_{0}$, outside the
figure.

Figure 2:  Binodals (full lines) and critical points (full and open circles)
of isospin-symmetric strongly interacting matter in (a) the temperature–baryon
density diagram and (b) the temperature–baryon chemical potential diagram. The
dashed line in panel (a) separates the region without effects of the modified
excluded-volume mechanism (lower left) from the region with effects (upper
right). Results for the liquid–gas phase transition are shown at subsaturation
densities.

Panel (b) of Figure 2 depicts the lines of the first-order PT in the
temperature–baryon chemical potential diagram ending in critical points. With
increasing temperature, the baryon chemical potential at the pseudo
hadron–quark PT reduces from $979.1$ to $591.8$ MeV at the critical point. By
crossing the transition line, an abrupt change in the density occurs that
becomes continuous at the critical point.

##  4 Conclusions

The extension of the modified EV approach to a density- and temperature-
dependent parametrization of the available-volume fractions as introduced in
this work was successful in achieving the main goal of this study: As a
generic structure of the QCD phase diagram, a first-order pseudo hadron–quark
phase transition at low temperatures and a crossover for low baryon densities
could be modeled that also includes a critical endpoint at $T_{\rm
crit}=155.5$ MeV and $\mu_{B,{\rm crit}}=591.8$ MeV. Other patterns of the QCD
phase diagram that have been theoretically motivated could also be modeled
within the present approach. Further extensions of the model are
straightforward. They include the extension to a larger number of components
of the hadron resonance gas in the underlying RMF model and an isospin
dependence of the EV model. It would be worthwhile to study further
thermodynamic quantities such as the speed of sound, heat capacities, or the
susceptibilities in such an enlarged model. It would also be interesting to
elaborate on a parametrization that would result in a second endpoint at low
temperatures, as suggested by Hatsuda et al. Hatsuda et al. (2006). The so-
generalized parametrization of the QCD EoS can be used in Bayesian analysis
studies for astrophysical applications pertaining to compact stars Steiner et
al. (2010); Alvarez-Castillo et al. (2016), their mergers, and core-collapse
supernova explosions, as well as heavy-ion collisions analogous to those
studied in Pratt et al. (2015).

###### Acknowledgements.

This work was supported by the Russian Science Foundation under contract
number 17-12-01427. S.T. was supported in part by the DFG through grant No.
SFB1245. D.B. and S.T. conceived and designed the concept of this work. S.T.
worked out the theoretical framework, performed the numerical calculations and
wrote the main body of the paper. D.B. contributed abstract, discussion,
conclusions and references. The authors declare no conflict of interest. The
following abbreviations are used in this manuscript:
CCSN | core-collapse supernova
---|---
EoS | equation of state
EV | excluded-volume
HIC | heavy-ion collision
NS | neutron star
PT | phase transition
QCD | quantum chromodynamics
RMF | relativistic mean-field
no multiple References yes

## References

  * Alford et al. (2013) Alford, M.G.; Han, S.; Prakash, M.  Generic conditions for stable hybrid stars.  Phys. Rev. 2013, D88, 083013, [arXiv:astro-ph.SR/1302.4732].
  * Benic et al. (2015) Benic, S.; Blaschke, D.; Alvarez-Castillo, D.E.; Fischer, T.; Typel, S.  A new quark-hadron hybrid equation of state for astrophysics - I. High-mass twin compact stars.  Astron. Astrophys. 2015, 577, A40, [arXiv:astro-ph.HE/1411.2856].
  * Alvarez-Castillo and Blaschke (2017) Alvarez-Castillo, D.E.; Blaschke, D.B.  High-mass twin stars with a multipolytrope equation of state.  Phys. Rev. 2017, C96, 045809, [arXiv:nucl-th/1703.02681].
  * Sagert et al. (2009) Sagert, I.; Fischer, T.; Hempel, M.; Pagliara, G.; Schaffner-Bielich, J.; Mezzacappa, A.; Thielemann, F.K.; Liebendörfer, M.  Signals of the QCD phase transition in core-collapse supernovae.  Phys. Rev. Lett. 2009, 102, 081101, [arXiv:astro-ph/0809.4225].
  * Fischer et al. (2009) Fischer, T.; Whitehouse, S.C.; Mezzacappa, A.; Thielemann, F.K.; Liebendörfer, M.  The neutrino signal from protoneutron star accretion and black hole formation.  Astron. Astrophys. 2009, 499, 1, [arXiv:astro-ph/0809.5129].
  * Fischer et al. (2011) Fischer, T.; Sagert, I.; Pagliara, G.; Hempel, M.; Schaffner-Bielich, J.; Rauscher, T.; Thielemann, F.K.; Käppeli, R.; Martinez-Pinedo, G.; Liebendörfer, M.  Core-collapse supernova explosions triggered by a quark-hadron phase transition during the early post-bounce phase.  Astrophys. J. Suppl. 2011, 194, 39, [arXiv:astro-ph.HE/1011.3409].
  * Fischer et al. (2017) Fischer, T.; Bastian, N.U.; Blaschke, D.; Cerniak, M.; Hempel, M.; Klähn, T.; Martínez-Pinedo, G.; Newton, W.G.; Röpke, G.; Typel, S.  The state of matter in simulations of core-collapse supernovae – Reflections and recent developments.  2017, [arXiv:astro-ph.HE/1711.07411].
  * Blaschke et al. (2010) Blaschke, D.B.; Sandin, F.; Skokov, V.V.; Typel, S.  Accessibility of Color Superconducting Quark Matter Phases in Heavy-ion Collisions.  Acta Phys. Polon. Supp. 2010, 3, 741–746, [arXiv:hep-ph/1004.4375].
  * Bazavov et al. (2014) Bazavov, A.; others.  Equation of state in ( 2+1 )-flavor QCD.  Phys. Rev. 2014, D90, 094503, [arXiv:hep-lat/1407.6387].
  * Borsanyi et al. (2014) Borsanyi, S.; Fodor, Z.; Hoelbling, C.; Katz, S.D.; Krieg, S.; Szabo, K.K.  Full result for the QCD equation of state with 2+1 flavors.  Phys. Lett. 2014, B730, 99–104, [arXiv:hep-lat/1309.5258].
  * Bratovic et al. (2013) Bratovic, N.M.; Hatsuda, T.; Weise, W.  Role of Vector Interaction and Axial Anomaly in the PNJL Modeling of the QCD Phase Diagram.  Phys. Lett. 2013, B719, 131–135, [arXiv:hep-ph/1204.3788].
  * Baym et al. (2017) Baym, G.; Hatsuda, T.; Kojo, T.; Powell, P.D.; Song, Y.; Takatsuka, T.  From hadrons to quarks in neutron stars 2017.  [arXiv:astro-ph.HE/1707.04966].
  * Abuki et al. (2010) Abuki, H.; Baym, G.; Hatsuda, T.; Yamamoto, N.  The NJL model of dense three-flavor matter with axial anomaly: the low temperature critical point and BEC-BCS diquark crossover.  Phys. Rev. 2010, D81, 125010, [arXiv:hep-ph/1003.0408].
  * Hatsuda et al. (2006) Hatsuda, T.; Tachibana, M.; Yamamoto, N.; Baym, G.  New critical point induced by the axial anomaly in dense QCD.  Phys. Rev. Lett. 2006, 97, 122001, [arXiv:hep-ph/hep-ph/0605018].
  * Klähn and Fischer (2015) Klähn, T.; Fischer, T.  Vector interaction enhanced bag model for astrophysical applications.  Astrophys. J. 2015, 810, 134, [arXiv:nucl-th/1503.07442].
  * Dexheimer and Schramm (2010) Dexheimer, V.A.; Schramm, S.  A Novel Approach to Model Hybrid Stars.  Phys. Rev. 2010, C81, 045201, [arXiv:astro-ph.SR/0901.1748].
  * Mukherjee et al. (2017) Mukherjee, A.; Steinheimer, J.; Schramm, S.  Higher-order baryon number susceptibilities: interplay between the chiral and the nuclear liquid-gas transitions.  Phys. Rev. 2017, C96, 025205, [arXiv:nucl-th/1611.10144].
  * Vovchenko et al. (2017) Vovchenko, V.; Pasztor, A.; Fodor, Z.; Katz, S.D.; Stoecker, H.  Repulsive baryonic interactions and lattice QCD observables at imaginary chemical potential.  Phys. Lett. 2017, B775, 71–78, [arXiv:hep-ph/1708.02852].
  * Typel (2016) Typel, S.  Variations on the excluded-volume mechanism.  Eur. Phys. J. 2016, A52, 16.
  * Oertel et al. (2017) Oertel, M.; Hempel, M.; Klähn, T.; Typel, S.  Equations of state for supernovae and compact stars.  Rev. Mod. Phys. 2017, 89, 015007, [arXiv:astro-ph.HE/1610.03361].
  * Rischke et al. (1991) Rischke, D.H.; Gorenstein, M.I.; Stoecker, H.; Greiner, W.  Excluded volume effect for the nuclear matter equation of state.  Z. Phys. 1991, C51, 485–490.
  * Satarov et al. (2009) Satarov, L.M.; Dmitriev, M.N.; Mishustin, I.N.  Equation of state of hadron resonance gas and the phase diagram of strongly interacting matter.  Phys. Atom. Nucl. 2009, 72, 1390–1415, [arXiv:hep-ph/0901.1430].
  * Typel et al. (2010) Typel, S.; Röpke, G.; Klähn, T.; Blaschke, D.; Wolter, H.H.  Composition and thermodynamics of nuclear matter with light clusters.  Phys. Rev. 2010, C81, 015803, [arXiv:nucl-th/0908.2344].
  * Vovchenko et al. (2017) Vovchenko, V.; Gorenstein, M.I.; Stoecker, H.  van der Waals Interactions in Hadron Resonance Gas: From Nuclear Matter to Lattice QCD.  Phys. Rev. Lett. 2017, 118, 182301, [arXiv:hep-ph/1609.03975].
  * Hempel and Schaffner-Bielich (2010) Hempel, M.; Schaffner-Bielich, J.  Statistical Model for a Complete Supernova Equation of State.  Nucl. Phys. 2010, A837, 210–254, [arXiv:nucl-th/0911.4073].
  * Chatterjee and Vidaña (2016) Chatterjee, D.; Vidaña, I.  Do hyperons exist in the interior of neutron stars?  Eur. Phys. J. 2016, A52, 29, [arXiv:nucl-th/1510.06306].
  * Schaffner-Bielich (2008) Schaffner-Bielich, J.  Hypernuclear Physics for Neutron Stars.  Nucl. Phys. 2008, A804, 309–321, [arXiv:astro-ph/0801.3791].
  * Typel et al. (2014) Typel, S.; Wolter, H.H.; Röpke, G.; Blaschke, D.  Effects of the liquid-gas phase transition and cluster formation on the symmetry energy.  Eur. Phys. J. 2014, A50, 17, [arXiv:nucl-th/1309.6934].
  * Steiner et al. (2010) Steiner, A.W.; Lattimer, J.M.; Brown, E.F.  The Equation of State from Observed Masses and Radii of Neutron Stars.  Astrophys. J. 2010, 722, 33–54, [arXiv:astro-ph.HE/1005.0811].
  * Alvarez-Castillo et al. (2016) Alvarez-Castillo, D.; Ayriyan, A.; Benic, S.; Blaschke, D.; Grigorian, H.; Typel, S.  New class of hybrid EoS and Bayesian M-R data analysis.  Eur. Phys. J. 2016, A52, 69, [arXiv:nucl-th/1603.03457].
  * Pratt et al. (2015) Pratt, S.; Sangaline, E.; Sorensen, P.; Wang, H.  Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions.  Phys. Rev. Lett. 2015, 114, 202301, [arXiv:nucl-th/1501.04042].


# Probing Exotic Phases Via Stochastic Gravitational Wave Spectra

Joshua Berger  Colorado State University, Fort Collins, Colorado, 80523, USA
Amit Bhoonah  Department of Physics and Astronomy and IQ Initiative,
University of Pittsburgh, Pittsburgh, PA 15260, USA  Biswajit Padhi  Colorado
State University, Fort Collins, Colorado, 80523, USA

###### Abstract

Stochastic backgrounds of gravitational waves (GWs) from the pre-BBN era offer
a unique opportunity to probe the universe beyond what has already been
achieved with the Cosmic Microwave Background (CMB). If the source is short in
duration, the low frequency tail of the resulting GW spectrum follows a
universal frequency scaling dependent on the equation of state of the universe
when modes enter the horizon. We demonstrate that the distortion of the
equation of state due to massive particles becoming non-relativistic can lead
to an observable dip in the GW spectrum. To illustrate this effect, we
consider a first order chiral symmetry breaking phase transition in the weak-
confined Standard Model (WCSM). The model features a large number of pions and
mostly elementary fermions with masses just below the critical temperature for
the phase transition. These states lead to a 20$\%$ dip in the GW power. We
find potential sensitivity to the distortions in the spectrum to future GW
detectors such as LISA, DECIGO, BBO, and $\mu$Ares.

##  I Introduction

The prospects for gravitational waves (GWs) to probe exotic particle phenomena
have gained increasing attention in light of their discovery by the LIGO
scientific collaboration [1]. The proof of principle that they provided has
led to the proposal for a large number of new GW telescopes covering a wide
range of frequencies, from (10-100) nHz with Pulsar Timing Arrays (PTAs) [2,
3, 4] to the mHz to 10 kHz window using ground and space based telescopes [5,
6, 7]. The range between PTAs and telescopes is challenging to cover, but
there exist proposals such as the $\mu$Ares detector [8] and through the use
of asteroids [9]. At much higher frequencies, in the MHz-GHz range, novel
laboratory techniques have already constrained gravitational wave backgrounds
[10, 11] and there are plans for future detectors capable of extending current
sensitivity further [12, 13]. These observatories will significantly broaden
the scope of new physics that can be explored through gravitational waves in
the coming decades.

This development has the potential to radically alter our understanding of the
universe in much the way detailed studies of the Cosmic Microwave Background
(CMB) [14] have over the past few decades. In studying the CMB [15, 16, 17,
18], we were able to infer with great precision the properties of the early
universe fluid by studying the spectrum of background photons, as well as the
anisotropies therein. It would be remarkable indeed if the same level of
detailed study could be performed by GW experiments. For the purposes of
understanding the early universe, it is particularly interesting to learn what
can be inferred from GWs emitted prior to big bang nucleosynthesis (BBN) [19].
The CMB can only directly probe the universe back to recombination at
$T\sim\text{eV}$, but observations of the primordial abundance of light
elements probes the universe back to $T\sim\text{MeV}$. Above this, we thus
far have no real probes and the best hope for accessing this era is in GWs;
the universe is essentially transparent to them since they only interact via
gravity and this makes them ideal messengers for probing the pre-BBN era. This
is, however, not without challenges because the universe was really tiny at
these times and effects like the Sachs-Wolfe effect [20, 21], that have proven
essential to determining our picture of the universe from CMB photons, are
suppressed.

Several early universe events, such as first order phase transitions [22, 23,
24], hybrid preheating after inflation [25], dynamics of topological defects
[26], and soliton collisions [27], can produce GWs in the pre-BBN era. These
are referred to as “stochastic” and would constitute a nearly isotropic
background of gravitational radiation just like the CMB is a nearly isotropic
background of electromagnetic radiation. In this work, we focus for
concreteness on GWs produced via a first order phase transition, though none
of the qualitative results of our analysis rely on much beyond assuming that
the GWs are produced nearly instantaneously and isotropically on cosmological
times.

A first order phase transition proceeds by nucleation of the new vacuum phase
as bubbles that expand. Eventually, those bubbles collide violently enough to
provide a first source of GWs. Colliding bubbles can also percolate and
interact with the surrounding thermal plasma, causing bulk motion in the fluid
which, in turn, provides two further sources of GWs: acoustic (sound wave) and
magnetohydrodynamic disturbances [28]. The spectrum of GW emission from these
processes have been simulated numerically and depends on multiple parameters
such as the duration of the phase transition, the amount of latent heat
released, and the fraction of energy which gets converted into the bulk motion
of the plasma (see, for example, Refs. [28, 29, 30] for more detailed
discussions). The precise spectrum of gravitational waves from these three
sources is an area of ongoing research.

Despite these model dependent features, a pair of related recent papers [31,
32] pointed out the interesting fact that the low frequency tail of GWs from
any “instantaneous” event would have a universal spectrum which only depends
on features of the universe when those modes enter the horizon. Intriguingly,
the spectrum is sensitive to the equation of state of the plasma. Provided the
universe is in a radiation dominated phase, as expected in a standard
cosmology, we expect $w\approx 1/3$. Nevertheless, there can be small
deviations from $w=1/3$ and a careful measurement of the spectrum of GWs would
be sensitive to such deviations. The equation of state has proven to be a
useful probe of early universe phases, as illustrated in e.g. [33, 34].

There are several potential sources such a distortion, but all require a
breaking of scaling invariance. Any scale-invariant theory will have $w=1/3$
exactly as

$\partial_{\mu}J^{\mu}_{\lambda}=T^{\mu}_{\mu}=\rho-3\,p,$ (1)

where $J_{\lambda}$ is the scale transformation current. One could imagine
several small sources of breaking in the radiation-dominated era. The recent
works metioned above, Refs. [31, 32], considered interactions inducing a scale
anomaly, leading to $w-1/3\propto\beta$, where $\beta$ is the beta function of
the theory, and free-streaming of particles that decouple from the plasma.

In this work, we focus primarily on scale-invariance breaking due to the
masses of particles. We demonstrate that when several degrees of freedom
simultaneously become non-relativistic, the resulting distortion of the
equation of state is sufficiently significant to be visible in future GW
observations. Since there are hundreds of bosonic degrees of freedom at high
temperatures, the resulting distortion goes like $g/g_{*}$, where $g$ is the
number of degrees of freedom that are becoming non-relativistic and $g_{*}$ is
the total number of relativistic degrees of freedom. Unlike the signals
considered previously, this distortion is transient. If the mass of the
particles is not too far from the temperature at which the GWs are sourced,
then a localized distortion of the long-wavelength tail of the GW spectrum
could be observed.

To demonstrate how this could work and could help determine qualitative
information about the non-standard early universe cosmology, we consider the
weak-confined standard model (WCSM) presented in Ref. [35, 36]. In that model,
a scalar field coupling to the $SU(2)_{L}$ gauge kinetic operator leads to an
effective shift in the scalar field coupling. This shift occurs when the
scalar field has a non-zero expectation value, a mechanism first presented in
Ref. [37]. If the shift is toward a stronger $SU(2)_{L}$ coupling, then it is
possible for the weak force to become strongly coupled at $T>T_{c}$, the
critical temperature for the electroweak crossover of standard cosmology. In a
strongly coupled $SU(2)_{L}$ phase, $B+L$-violating transitions would occur
readily, opening the door for a possible baryogenesis mechanism. For such a
mechanism to work, the standard cosmology phase would need to be restored at
$T<T_{c}$ so that the electroweak crossover does not washout the generated
asymmetry.

More to the point of this work, the transition to strong coupling could be a
first order phase transition, leading naturally to a source of stochastic
gravitational waves. What makes the WCSM a particularly interesting source is
that the nature of the phase of strongly-coupled weak force is in question. It
is unclear whether this would be a chiral symmetry-breaking confined phase, as
with standard model Quantum Chromodynamics (QCD), or an infrared conformal
phase. In the latter case, the only distortions away from $w=1/3$ would be due
to the scale anomaly of the remaining interactions (QCD, electromagnetic, and
Yukawa). On the other hand, the former would break scale invariance more
strongly. A number of species would have masses near and a bit below the
critical temperature, leading to modifications of the equation of state in the
region below the phase transition. We show that in a chiral symmetry-breaking
phase, pions becoming non-relativistic below the phase transition lead to more
than 10% deviation from $w=1/3$, correspondingly leading to an order 10% dip
in the characteristic strain or 20% dip in the GW power at a frequency
corresponding to modes that enter the horizon just as the pions are becoming
non-relativistic. For a WCSM phase transition above the would-be electroweak
crossover, such a dip could yield potentially observable gravitational wave
signatures.

In this work, we develop a calculation of the distortion of the gravitational
waves due to particle masses in the WCSM. This requires modifying the results
of Refs. [31, 32] to account for a time-dependent equation of state. We solve
the resulting gravitational wave equations numerically for benchmark models.
We find potential sensitivity to the distortions in the spectrum to future GW
detectors such as the Laser Interferometer Space Antenna (LISA) [38], Deci-
hertz Interferometer Gravitational wave Observatory (DECIGO) [5], Big Bang
Observer (BBO) [6], and the $\mu$Ares detector.

The remainder of this paper is structured as follows. In Section II, we
develop the formalism required to study time-dependent distortions of the
equation of state and demonstrate the effect of a small number of non-
relativistic species on the gravitational wave spectrum. We then review the
relevant details of the WCSM model in Section III. The results of our study
are presented in Section IV. Finally, in Section V, we discuss the
implications of our results for future measurements.

##  II General Formalism

###  II.1 Gravitational Wave Equation

Our starting point is to consider a Friedmann-Robertson-Lemaître-Walker
universe with metric

$ds^{2}=dt^{2}-a(t)^{2}\,d\mathbf{x}^{2}.$ (2)

Following Refs. [31, 32], we work in conformal time $\tau=\int dt/a(t)$ and
expand the FLRW metric in small perturbations to obtain the linearized
gravitational wave equation in comomving momentum space:

$h^{\prime\prime}+2\,\mathcal{H}\,h^{\prime}+k^{2}\,h=0,$ (3)

where $\prime$ denotes derivatives with respect to conformal time,
$\mathcal{H}=a^{\prime}/a$ is the conformal Hubble rate, $k$ is the comoving
wavenumber, and $h$ is the metric perturbation amplitude.

Suppose that at a time $\tau_{*}$ an instantaneous event (meaning that the
duration of the event $\Delta\tau_{*}\ll 1/\mathcal{H}_{*}$, where
$\mathcal{H}^{*}\equiv\mathcal{H}(\tau_{\star})$) generates gravitational
waves. Then one must add a source on the right-hand side of (3), and the
resulting equation can be solved using Green’s function methods, which is
equivalent to solving the source-free equation with initial conditions

$h(\tau_{*})=0,\qquad h^{\prime}(\tau_{*})=J_{*},$ (4)

where $J_{*}$ is related to the dimensionless anisotropic stress inducing the
gravitational waves,

$J_{*}=a^{2}\,\frac{32\,\pi\,G\,\rho}{3}\,\epsilon^{ij}\,\Pi_{ij},$ (5)

projecting onto the GW polarization tensor $\epsilon$.

When $k\ll\mathcal{H}_{*}$, the modes are initially frozen outside the
horizon. As they enter the horizon, they become sensitive the equation of
state determining $\mathcal{H}$. When $k\gg\mathcal{H}$, the modes lose
sensitivity to the equation of state and simply oscillate. Thus, if a number
of degrees of freedom are becoming non-relativistic at the time when a given
mode is entering the horizon, the spectrum for that mode gets distorted.

We demonstrate that future experiments should have sensitivity to this
distortion. To see this, we must first generalize the calculation of [31, 32].
In that work, the authors considered the equation of state to be approximately
constant. For the effect we are interested in, this approximation is not
accurate as the modification to the equation of state is happening on a Hubble
time, which is the same as the time scale for the modes to enter the horizon.

Defining the equation of state $w=p/\rho$ and plugging the first Friedmann
equation $\mathcal{H}^{2}=8\,\pi\,G\,\rho\,a^{2}/3$ into the acceleration
Friedmann equation, we get an equation for $a$ in terms of $w$:

$2\,a\,a^{\prime\prime}=(1-3\,w)\,(a^{\prime})^{2}.$ (6)

This equation can be solved to determine $\mathcal{H}$ as

$\mathcal{H}=\frac{2}{\tau+3\,\int_{0}^{\tau}w(\tau_{1})\,d\tau_{1}}.$ (7)

As $w(\tau)$ can be determined from the matter content of the universe, we can
solve the wave equation (3), at least numerically.

It will be helpful to change variables in equation (3), defining

$\tilde{h}=\tilde{\tau}\,h,\qquad\tilde{\tau}=k\,\tau,\qquad\tilde{\mathcal{H}}=\frac{\mathcal{H}}{k}.$
(8)

In terms of these variables, equation (3) becomes

$\tilde{h}^{\prime\prime}+2\,\left(\tilde{\mathcal{H}}-\frac{1}{\tilde{\tau}}\right)\,\tilde{h}^{\prime}+\left(1+\frac{2}{\tilde{\tau}}\,\left[\frac{1}{\tilde{\tau}}-\tilde{\mathcal{H}}\right]\right)\,\tilde{h}=0,$
(9)

with primes now denoting derivatives with respect to $\tilde{\tau}$. In terms
of these variables, for a purely conformal equation of state $w=1/3$,
$\tilde{H}=1/\tilde{\tau}$, and the equation reduces to that of a harmonic
oscillator with angular frequency 1.

At these late times, the solution goes like

$h=A(k)\frac{J_{*}\tau_{*}\sin(k(\tau-\tau_{*})+\phi)}{k\tau},$ (10)

where $\phi$ is a phase shift induced due to the change in equation of state
and $A$ quantifies the change in amplitude relative to the case of no change
in equation of state. The ensemble-averaged power is approximately a period-
averaged power at this time, so that we can write

$\frac{dP}{dk}\approx\langle(h^{\prime})^{2}\rangle\approx\frac{A^{2}(k)J_{*}^{2}\tau_{*}^{2}}{2k\tau^{2}},$
(11)

independent of the phase shift. We can therefore examine $A^{2}(k)$ as an
observable measure of the effect of the distortion of the equation of state.

###  II.2 Time-Dependent Equation of State

Each species in equilibrium with temperature $T$ in the universe contributes
an energy density

$\rho_{i}=g_{i}\,\frac{T^{4}}{2\,\pi^{2}}\,\int_{x}^{\infty}dy\,\frac{y^{2}\,\sqrt{y^{2}-x^{2}}}{e^{y}\pm
1},$ (12)

and a pressure

$p_{i}=g_{i}\,\frac{T^{4}}{6\,\pi^{2}}\,\int_{x}^{\infty}dy\,\frac{(y^{2}-x^{2})^{3/2}}{e^{y}\pm
1},$ (13)

where $+$ is for fermions and $-$ is for bosons, $g_{i}$ is the number of
degrees of freedom corresponding to the species and $x=m_{i}/T$ for a species
of mass $m_{i}$. The total energy density and pressure are obtained by simply
summing over the species,

$\rho_{\text{tot}}=\sum_{i}\rho_{i},\qquad p_{\text{tot}}=\sum_{i}p_{i}.$ (14)

From this, we can determine the equation of state parameter $w$ as a function
of temperature via

$w(T)=\frac{p_{\text{tot}}}{\rho_{\text{tot}}}.$ (15)

Then, we can calculate $\delta w$ as a function of temperature using

$\delta w(T)=w(T)-\frac{1}{3}$ (16)

To illustrate the change in the equation of state, we first consider a toy
model. We consider a scenario in which a thermal bath contains 2 different
bosonic species. One of them has 40 relativistic degrees of freedom and the
other species has 60 degees of freedom of mass $10^{4}$ GeV. The behavior of
$\delta w$ with respect to temperature in this toy scenario is shown in Figure
1. In this toy model, a peak change in $w$ is found at $T=0.38m$ with $\delta
w=0.03$.

Figure 1: $\delta w$ for a thermal bath containing 2 bosonic species. The
degrees of freedom with a mass of $10^{4}$ GeV decouple from the bath at
temperatures of the order of their mass causing $w$ to deviate from 1/3.

To write $w$ as a function of $\tau$, we need the relation between conformal
time and temperature, obtained by combining the first Friedmann equation with
the conservation of entropy,

$\frac{d}{dT}(s\,a^{3})=0.$ (17)

We can then obtain the solution

$\tau(T)=\int_{T}^{\infty}\frac{ds(T^{\prime})/dT^{\prime}}{2\,s(T^{\prime})\,\sqrt{6\,\pi
G\,\rho_{\text{tot}}(T^{\prime})}}\,dT^{\prime}.$ (18)

The entropy is given by the usual Euler relation

$s=\frac{\rho_{\text{tot}}+p_{\text{tot}}}{T}.$ (19)

Eqs. (18) and (19) allow us to express the conformal time as a function of
temperature. This relation can be inverted to get the temperature as a
function of conformal time, which allows us to express $\delta w$ as a
function of $\tau$.

###  II.3 Gravitational Wave Model

The phase for the gravitational wave is arbitrary and depends on the random
initial conditions of the universe. To estimate the sensitivity, we therefore
consider the ensemble-averaged power in the gravitational waves,

$\rho_{\text{gw}}=\frac{1}{32\,\pi\,G\,a^{2}}\,\sum_{\text{helicities}}\langle
h^{\prime}_{ij}(\tau,\mathbf{x})h^{\prime}_{ij}(\tau,\mathbf{x})\rangle$ (20)

In Fourier space, we define the power spectrum by

$\langle
h(\tau,k)^{2}\rangle=(2\,\pi)^{3}\,\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})\,P_{h}(k).$
(21)

From this, we can define the differential power by

$\frac{d\Omega_{\text{gw}}}{d\log
k}=\frac{1}{\rho_{c}}\,\frac{k^{5}\,P_{h}(\tau,k)}{2\,(2\,\pi)^{3}\,a^{2}\,G}.$
(22)

As discussed before, for a first order phase transition there are three
contributions to the gravitational wave spectrum, from collisions of the
bubbles themselves and from acoustic and magnetohydrodynamic distrubances
induced by their expansion in the thermal plasma. The three sources are
expected to combine linearly,

$h^{2}\Omega_{\text{gw}}=h^{2}\Omega_{\text{col}}+h^{2}\Omega_{\text{sw}}+h^{2}\Omega_{\text{MHD}}.$
(23)

To determine which one of these sources dominate we must consider the system
undergoing a first order phase transition. In a strongly interacting theory
such as the WCSM, one can reasonably expect that the bubble wall interacts
sufficiently with the thermal plasma to cause bulk motion of the latter. In
this scenario, numerous numerical studies (see for example [29, 39]) have
found that sound waves are the dominant source of GWs compared to the
collisions of bubbles themselves.111In the case of a supercooled electroweak
phase transition, the dominant source comes from bubble wall collisions, but
in our scenario, we expect the sound wave contribution to dominate. The peak
frequency is given by [29, 30]

$f_{\text{sw}}=1.9\times
10^{-5}\frac{1}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100\
\text{GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\text{Hz},$ (24)

where $\beta$ is the inverse duration of the phase transition, $H_{*}$ and
$g_{*}$ are respectively the Hubble parameter and number of relativistic
degrees of freedom at a temperature $T_{*}$. The latter is usually taken to be
the nucleation temperature $T_{n}$ when the bubble nucleation rate is one per
Hubble volume since that is when the production of GWs from all three sources
is most significant. The dimensionless gravitational wave spectrum today
obtained numerically in [39] is well fitted by [29]

$\begin{split}h^{2}\Omega_{sw}&=2.65\times
10^{-6}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}\
v_{w}\ S_{sw}(f)\\\
S_{sw}(f)&=\left(\frac{f}{f_{sw}}\right)^{3}\frac{7}{4+3\left(f/f_{sw}\right)^{2}},\end{split}$
(25)

where $\alpha$ is the ratio of the vacuum energy released during the phase
transition to the energy of the radiation bath, and $\kappa$ is the fraction
of the energy of the vacuum energy that gets converted into the bulk motion of
the fluid. Aside from $\alpha$ and $\kappa$, the density is parameterized in
terms of other properties of the phase transition such as the ratio of the
inverse duration of the phase transition to the conformal Hubble at the phase
transition $\mathcal{H}_{*}$, the bubble wall velocity $v_{w}$, and the number
of relativistic degrees of freedom at $T_{*}$.

We consider only long wavelength waves, with $f<f_{\text{sw}}$. The $f^{3}$
scaling is a universal property of short duration phase transition. Rather
than directly calculate this power spectrum in our model, we perform a simple
rescaling by

$h^{2}\Omega_{\text{gw}}=\frac{\langle\tilde{h}^{2}\rangle}{\langle\tilde{h}^{2}\rangle_{0}}h^{2}\Omega_{\text{sw}}\\\
,$ (26)

where we define the averaging procedure here by averaging over a period in
$\tilde{\tau}$ of $2\pi$ at $\tau\to\infty$, when the gravitational waves are
well within the horizon and lose their sensitivity to the equation of state.
We denote by $\langle\tilde{h}^{2}\rangle_{0}$ the averaged gravitational wave
in a purely radiation-dominated universe. Taking the now-arbitrary
normalization $J_{*}=1$, we find

$\langle\tilde{h}^{2}\rangle_{0}=\frac{\tilde{\tau}^{2}_{*}}{2}.$ (27)

The goal for much of the remainder of this paper is then to solve the equation
(9), with the initial conditions $\tilde{h}(k\,\tau_{*})=0$ and
$\tilde{h}^{\prime}(k\,\tau_{*})=\tilde{\tau}_{*}$.

Before going into that calculation, we note on the observable spectrum in the
present day universe. The gravitational waves are a bit colder than the
radiation, as they decouple immediately after the phase transition. The
spectrum gets scaled by a factor of $a(\tau_{*})^{-4}$ and we need to convert
from co-moving wave number to physical frequency. We work in a convention with
$a(\tau_{*})=1$, in which case $f=k/(2\,\pi\,a(\tau_{0}))$.

##  III Particle Physics Models

|  Particle |  Degrees of freedom |  Mass (GeV)
---|---|---|---
|  $\Pi^{0}_{3}$ |  3 |  $1.95\times 10^{5}$
|  |  6 |  $2.18\times 10^{5}$
|  $\Pi_{2}^{\pm}$ |  8 |  $6.38\times 10^{4}$
|  |  8 |  $1.59\times 10^{5}$
|  |  12 |  $1.55\times 10^{5}$
Pions |  |  4 |  $7.23\times 10^{4}$
|  $\Pi_{1}^{\pm}$ |  6 |  $8.90\times 10^{4}$
|  $\Pi_{1}^{0}$ |  2 |  $9.82\times 10^{4}$
|  |  4 |  $1.07\times 10^{3}$
|  |  1 |  5.32
|  |  2 |  0
|  $\Xi_{1}^{0}$ |  6 |  0
|  $\Xi_{1}^{\pm}$ |  4 |  $2.31\times 10^{-7}$
|  |  4 |  $2.81\times 10^{-2}$
Fermions |  |  4 |  $6.04\times 10$
|  $\Xi_{2}^{\pm}$ |  8 |  $2.11\times 10^{-7}$
|  |  8 |  $2.48\times 10^{-2}$
|  |  8 |  $1.42\times 10^{2}$
|  $A^{\prime}$ |  2 |  0
|  $G$ |  6 |  0
Gauge bosons |  $W^{\prime\pm}$ |  12 |  $1.19\times 10^{5}$
|  $Z^{\prime}$ |  3 |  $2.09\times 10^{5}$
Table 1: Numerical spectrum of the WCSM at $f_{\pi}=$ 80 TeV. Additionally,
there are 4 light tachyonic degrees of freedom as well.

The Weak Confined Standard Model (WCSM) is a model in which the
$\text{SU}(2)_{\text{L}}$ component of the electroweak force is strongly
coupled and weak isospin is confined. To achieve this, a scalar is coupled to
the kinetic term for the gauge bosons,

$\mathcal{L}=-\frac{\Phi}{\Lambda}W^{i}_{\mu\nu}W^{i\mu\nu},$ (28)

where $\Lambda$ is the scale of the operator. If the scalar $\Phi$ acquires an
expectation value in the early universe, the value of the weak coupling
$g^{\prime}$ is shifted. Since the weak coupling constant is asymptotically
free, it will run to strong coupling if the nominal confinement scale
$\Lambda_{\text{W}}$ is larger than the elecotroweak phase transition
temperature $T_{c}\lesssim v$. In order to have Big Bang Nucleosynthesis as
observed, the weak coupling should return to the SM phase at a temperature
$T_{\text{SM}}\gtrsim 10~{}\text{MeV}$. Further details of the WCSM can be
found in Ref. [35].

The SM has 12 $SU(2)_{\text{L}}$ doublet Weyl fermions and therefore has a
$SU(12)$ flavor symmetry with respect to this gauge group. This symmetry gets
broken to $Sp(12)$ due to chiral symmetry breaking assumed to occur at strong
$SU(2)_{L}$ coupling. As a result, there are 65 massless Goldstone bosons
(pions) in the spectrum. The flavor symmetry is explicitly broken by both the
$SU(3)_{C}\times U(1)_{Y}$ gauge couplings and the fermion Yukawa couplings.
The non-isospin gauge symmetry is spontaneously broken as $SU(3)_{C}\times
U(1)_{Y}\rightarrow SU(2)_{C}\times U(1)_{Q}$. Out of 9 gauge bosons, 5 gauge
bosons become massive and five pions are eaten. Out of the remaining 60 pions,
58 pions acquire masses through gauge interactions and Yukawa interactions.
Four of these pion masses are tachyonic, but we expect these pions to be light
after minimizing their potential.

In addition to the massive pions, there are fermionic states in the WCSM which
contribute to modifying the equation of state. The left-handed fermions of the
SM become mostly composite, while the right-handed states are mostly
elementary. The composite fermions have masses around the confinement scale
$\Lambda_{\text{W}}\sim 4\pi f_{\pi}$. Only the elementary fermions, with
masses $\ll\Lambda_{\text{W}}$, contribute significantly to modifying the
equation of state at relevant temperatures.

The mass spectrum of the pions and the fermions is sensitive to dimensionless
operator coefficients which encode loop corrections and are expected to be
$\mathcal{O}$(1) numbers. The values of these coefficients can be taken to be
$\pm 1$ for the purpose of calculating the numerical spectrum of the WCSM
given in Table 1.

With this numerical spectrum, we can repeat the procedure outlined in Sec.
II.2 to calculate $\delta w$ for the WCSM. The behaviour of $\delta w$ over a
range of values of $T$ for a benchmark value of $f_{\pi}=80$ TeV is shown in
Fig. 2.

Figure 2: Deviations from $w=\frac{1}{3}$ in the WCSM for a benchmark scenario
where $f_{\pi}\simeq$ 80 TeV. The peak near $T/f_{\pi}\simeq$ 1 corresponds
particles of mass around $f_{\pi}$ \- mainly 49 pions which get their mass
predominantly from the strong coupling or the top Yukawa - decoupling from the
thermal bath.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: The GW spectrum of an assumed first order WCSM phase transition
compared to the reach of future experiments like LISA (purple), DECIGO
(orange), BBO (green), and $\mu$Ares (brown). The stochastic background of GWs
from galactic binaries is shown in gray. The dotted lines in the spectrum
denote frequencies with model-dependent effects. The left panels show the
spectra for BI and the right panels the spectra for BII as described in the
text.

For the gravitational wave spectrum, we illustrate three benchmark values for
$f_{\pi}$ = $\left(\text{0.08, 8, 80}\right)$ TeV, which sets the critical
temperature of the chiral phase transition $T_{C}\simeq f_{\pi}$. For the
gravitational wave parameters in (25), the precise value of $\beta/H_{*}$ is
obtainable only through a numerical lattice study of the WCSM. That study is
beyond the scope of this work and we therefore adopt three benchmarks,
$\beta/H_{*}$ = 1, 10, 100, noting that since we are assuming that the phase
transition occurs quasi-instantaneously, $\beta/H_{*}$ = 1 is a lower limit.
Furthermore, we adopt two benchmark (BI and BII) values for $\alpha$ and
$\kappa$. BII assumes values for the latter that are fully consistent with the
numerical simulation of [39] while BII is an extrapolation of the latter
assuming that the GW spectrum scaling with respect $\alpha$ and $\kappa$ is
unchanged. For both benchmarks, we assume the wall velocity to be $v_{w}$ =
0.83 c, the upper limit considered [39]. With these, we show the GW spectrum
of an assumed first order WCSM phase transition and compare them to the reach
of future experiments like LISA, DECIGO, and BBO. We also show the expected GW
spectrum from galactic binaries, the dominant background for our model. The
dotted lines in the spectrum are used to describe modes entering the horizon
when the validity of the equation of state is uncertain (we describe this in
more detail in the next section).

##  IV Results

We have solved (9) numerically in the WCSM phase with initial conditions
described in Section (II). As discussed above, the particle physics effects of
the WCSM on this equation describing the propagation of GWs are encoded in the
small changes to the Hubble parameter due to deviations of the equation of
state from $w=1/3$ during the radiation domination era. The latter is
determined assuming the validity of chiral perturbation theory, but there is
some uncertainty, inherent to the strongly coupled nature of this theory, on
what the critical temperature for the phase is, what the properties and nature
of the phase transition are, and what distortions of the equation of state
would occur due to relatively strong pion interactions close to the
confinement scale. A robust answer to these questions can only come from a
detailed numerical lattice study, which is beyond the scope of this work.

As such, we trace the high frequency end of our spectrum within an order of
magnitude of the cutoff scale $\Lambda_{\chi}=4\pi f_{\pi}$ for chiral
perturbation theory, as a dashed line. Furthermore, any modes that would not
quite reach a point where they lose sensitivity to the equation of state
before the universe deconfines back to the SM phase are subject to
modification by that deconfining phase transition and changes to the equation
of state just below that phase transition. Therefore, modes for which
$10\pi/k<\tau_{\text{SM}}$, the conformal time of the transition to the SM
phase, are also drawn as dashed lines.

In Figure 4, we plot the ratio of the time-averaged gravitational wave
amplitude today for $\delta w\neq$ 0 compared to the $\delta w$ = 0 case for
$f_{\pi}\simeq$ 80 TeV. Since this value of $f_{\pi}$ corresponds to the
variation of the equation of state with respect to temperature (upper axis)
shown in Figure 2, we can compare the two figures and see that the largest
change occurs at $k\lesssim$ 10-9 GeV and corresponds to the mode that enters
the horizon when particles of mass around $f_{\pi}$ (mainly 49 pions which get
their mass predominantly from the strong coupling or the top Yukawa) become
non-relativistic.

Figure 4: Average gravitational wave power relative to that in the absence of
changes to the equation of state as a function of comoving wavenumber in
(GeV). The vertical axis is the temperature of the thermal bath at the moment
of horizon entry for each mode. The largest change - the dip at $k\lesssim$
10-7 GeV corresponds to the mode that enters the horizon when particles of
mass $\mathcal{O}(f_{\pi})$ \- mainly 49 pions - decouple from the thermal
bath.

We also show a zoomed in version of the $f_{\pi}$ = 80 TeV, $\alpha$ = 0.5,
$\kappa$ = 0.1, $\beta/H_{*}$ = 1 benchmark (Fig. 3f) in Fig. 5. To make the
$f^{3}$ scaling at low frequencies apparent, we plot
$h^{2}\Omega_{\text{gw}}/f^{3}$ instead of $h^{2}\Omega_{\text{gw}}$. We show
the $\delta w=0$ case, corresponding to no conformal symmetry breaking, and
the WCSM with chiral symmetry breaking, where conformal symmetry is strongly
broken. We see that LISA has some sensitivity to the two curves in a region
where the conformal and non-conformal scenarios differ which lies at the upper
end of the region where the causality limited $f^{3}$ scaling is valid. There
are also pronounced dips well into the region where the $f^{3}$ scaling is
valid, but this lies in the $\mu\text{Hz}$ gap described in Ref. [9].
Unfortunately, the sensitivity of the proposed detection scheme in this work
still does not quite reach the sensitivity required to see the deviation from
the expected $f^{3}$ scaling for this benchmark, but we hope that our model
adds further motivation to thinking about future detectors capable of probing
this frequency range, where there is already significant interest for
astrophysical sources of GWs.

Figure 5: A zoomed in plot Figure 3f with only the $\beta/H_{*}$ = 1 benchmark
case shown. To make the $f^{3}$ scaling at low frequencies apparent, we plot
$h^{2}\Omega_{\text{gw}}/f^{3}$ instead of $h^{2}\Omega_{\text{gw}}$. The blue
curve represents the $\delta w=0$ case, corresponding to no conformal symmetry
breaking, the dark blue curve represents WCSM with chiral symmetry breaking.
We see that LISA and, especially, $\mu$Ares have some sensitivity to the two
curves in a region where they defer and the universal $f^{3}$ scaling is
valid, making such an effect potentially observable.

##  V Discussion and Outlook

The signal we describe here is quite general and applies to any scenario with
a long wavelength gravitational wave background and a significant number of
particles becoming non-relativistic. For example, one could examine
supersymmetric scenarios with a large spectrum of superpartners at similar
masses. The signal also bears some resemblance to that studied by Ref. [41].
One could imagine slight tweaks of that scenario to get similar dips.

One further intriguing possibility is to modify the QCD phase transition such
that it is first order. One scenario in which this would occur is a
supercooled electroweak phase transition [42], as the QCD transition would
then occur with six light flavors of quarks. After the QCD phase transition,
the amount of radiation in the universe is rather low. As such, the process of
the muons becoming non-relativistic is sufficient to induce a roughly 5%
change in $w$, which is more challenging to observe.

While the prediction of a $f^{3}$ scaling of the power, up to modifications
due to the equation of state, is universal, the normalization of the spectrum
and peak frequency are dependent on the properties of the phase transition.
Further calculations are needed to determine these properties of the phase
transition. Lattice gauge theory is rather challenging for this model due to
the addition of the Higgs field, but this would be the most robust way to
study the phase transition. A lattice study would also elucidate the nature of
the strongly interacting phase, which we assume be chiral symmetry breaking as
in QCD. It is possible that the gravitational wave properties can be also
estimated from a Nambu-Jona-Lasinio model [43, 44].

Note that any dimensionful parameters in the model could in principle lead to
such a deviation in the spectrum. In particular a super-renormalizable scalar
cubic coupling or non-renormalizable operator will lead to a small shift in
the power of the low frequency spectrum.

The promising signal type presented in this work will ultimately require
further study by gravitational wave telescopes.

Acknowledgments. We thank Akshay Ghalsasi for useful discussions and Andrew
Long for feedback on the manuscript. JB would like to thank the Mainz
Institute for Theoretical Physics and the Aspen Center for Physics for their
hospitality during the completion of parts of this work. JB and BP are
supported by the National Science Foundation under Grant No. 2112789.

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Manchester [2006] R. N. Manchester, Chin. J. Astron. Astrophys. Suppl. 6, 02139 (2006), arXiv:astro-ph/0604288 .
  * Kramer and Champion [2013] M. Kramer and D. J. Champion, Class. Quant. Grav. 30, 224009 (2013).
  * Arzoumanian  _et al._ [2020] Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Kawamura  _et al._ [2006] S. Kawamura _et al._ , Class. Quant. Grav. 23, S125 (2006).
  * Harry  _et al._ [2006] G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner, and E. S. Phinney, Class. Quant. Grav. 23, 4887 (2006), [Erratum: Class.Quant.Grav. 23, 7361 (2006)].
  * Mei  _et al._ [2021] J. Mei _et al._ (TianQin), PTEP 2021, 05A107 (2021), arXiv:2008.10332 [gr-qc] .
  * Sesana  _et al._ [2021] A. Sesana _et al._ , Exper. Astron. 51, 1333 (2021), arXiv:1908.11391 [astro-ph.IM] .
  * Fedderke  _et al._ [2022] M. A. Fedderke, P. W. Graham, and S. Rajendran, Phys. Rev. D 105, 103018 (2022), arXiv:2112.11431 [gr-qc] .
  * Chou  _et al._ [2017] A. S. Chou _et al._ (Holometer), Phys. Rev. D 95, 063002 (2017), arXiv:1611.05560 [astro-ph.IM] .
  * Akutsu  _et al._ [2008] T. Akutsu _et al._ , Phys. Rev. Lett. 101, 101101 (2008), arXiv:0803.4094 [gr-qc] .
  * Ito and Soda [2022] A. Ito and J. Soda, (2022), arXiv:2212.04094 [gr-qc] .
  * Berlin  _et al._ [2023] A. Berlin, D. Blas, R. Tito D’Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn, J. Schütte-Engel, and M. Wentzel, (2023), arXiv:2303.01518 [hep-ph] .
  * Penzias and Wilson [1965] A. A. Penzias and R. W. Wilson, Astrophys. J. 142, 419 (1965).
  * Smoot  _et al._ [1992] G. F. Smoot _et al._ (COBE), Astrophys. J. Lett. 396, L1 (1992).
  * Bennett  _et al._ [2013] C. L. Bennett _et al._ (WMAP), Astrophys. J. Suppl. 208, 20 (2013), arXiv:1212.5225 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A1 (2020), arXiv:1807.06205 [astro-ph.CO] .
  * Swetz  _et al._ [2011] D. S. Swetz _et al._ , Astrophys. J. Suppl. 194, 41 (2011), arXiv:1007.0290 [astro-ph.IM] .
  * Wagoner  _et al._ [1967] R. V. Wagoner, W. A. Fowler, and F. Hoyle, Astrophys. J. 148, 3 (1967).
  * Sachs and Wolfe [1967] R. K. Sachs and A. M. Wolfe, Astrophys. J. 147, 73 (1967).
  * Rees and Sciama [1968] M. J. Rees and D. W. Sciama, Nature (London) 217, 511 (1968).
  * Turner and Wilczek [1990] M. S. Turner and F. Wilczek, Phys. Rev. Lett. 65, 3080 (1990).
  * Hogan [1986] C. J. Hogan, MNRAS 218, 629 (1986).
  * Kamionkowski  _et al._ [1994] M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D 49, 2837 (1994), arXiv:astro-ph/9310044 .
  * Garcia-Bellido and Figueroa [2007] J. Garcia-Bellido and D. G. Figueroa, Phys. Rev. Lett. 98, 061302 (2007), arXiv:astro-ph/0701014 .
  * Vilenkin [1981] A. Vilenkin, Phys. Lett. B 107, 47 (1981).
  * Helfer  _et al._ [2019] T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, Phys. Rev. D 99, 044046 (2019), arXiv:1802.06733 [gr-qc] .
  * Caprini  _et al._ [2009] C. Caprini, R. Durrer, T. Konstandin, and G. Servant, Phys. Rev. D 79, 083519 (2009), arXiv:0901.1661 [astro-ph.CO] .
  * Caprini  _et al._ [2016] C. Caprini _et al._ , JCAP 04, 001 (2016), arXiv:1512.06239 [astro-ph.CO] .
  * Caprini  _et al._ [2020] C. Caprini _et al._ , JCAP 03, 024 (2020), arXiv:1910.13125 [astro-ph.CO] .
  * Hook  _et al._ [2021] A. Hook, G. Marques-Tavares, and D. Racco, JHEP 02, 117 (2021), arXiv:2010.03568 [hep-ph] .
  * Brzeminski  _et al._ [2022] D. Brzeminski, A. Hook, and G. Marques-Tavares, (2022), arXiv:2203.13842 [hep-ph] .
  * Chung  _et al._ [2011] D. Chung, A. Long, and L.-T. Wang, Phys. Rev. D 84, 043523 (2011), arXiv:1104.5034 [astro-ph.CO] .
  * Chung and Long [2011] D. J. H. Chung and A. J. Long, Phys. Rev. D 84, 103513 (2011), arXiv:1108.5193 [astro-ph.CO] .
  * Berger  _et al._ [2019] J. Berger, A. J. Long, and J. Turner, Phys. Rev. D 100, 055005 (2019), arXiv:1906.05157 [hep-ph] .
  * Lohitsiri and Tong [2019] N. Lohitsiri and D. Tong, SciPost Phys. 7, 059 (2019), arXiv:1907.08221 [hep-th] .
  * Ipek and Tait [2019] S. Ipek and T. M. P. Tait, Phys. Rev. Lett. 122, 112001 (2019), arXiv:1811.00559 [hep-ph] .
  * Amaro-Seoane  _et al._ [2012] P. Amaro-Seoane _et al._ , Class. Quant. Grav. 29, 124016 (2012), arXiv:1202.0839 [gr-qc] .
  * Hindmarsh  _et al._ [2015] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 92, 123009 (2015), arXiv:1504.03291 [astro-ph.CO] .
  * Note [1] In the case of a supercooled electroweak phase transition, the dominant source comes from bubble wall collisions, but in our scenario, we expect the sound wave contribution to dominate.
  * Dienes  _et al._ [2022] K. R. Dienes, L. Heurtier, F. Huang, D. Kim, T. M. P. Tait, and B. Thomas, Phys. Rev. D 105, 023530 (2022), arXiv:2111.04753 [astro-ph.CO] .
  * Iso  _et al._ [2017] S. Iso, P. D. Serpico, and K. Shimada, Phys. Rev. Lett. 119, 141301 (2017), arXiv:1704.04955 [hep-ph] .
  * Aoki  _et al._ [2017] M. Aoki, H. Goto, and J. Kubo, Phys. Rev. D 96, 075045 (2017), arXiv:1709.07572 [hep-ph] .
  * Wang  _et al._ [2022] Y. Wang, M. Kawaguchi, S. Matsuzaki, and A. Tomiya, Phys. Rev. D 106, 095010 (2022), arXiv:2208.03975 [hep-ph] .
  * Weir [2017] D. J. Weir, Phil. Trans. Roy. Soc. Lond. A 376, 20170126 (2017), arXiv:1705.01783 [hep-ph] .
  * von Harling and Servant [2018] B. von Harling and G. Servant, JHEP 01, 159 (2018), arXiv:1711.11554 [hep-ph] .
  * Hayakawa  _et al._ [2013] M. Hayakawa, K. I. Ishikawa, S. Takeda, M. Tomii, and N. Yamada, Phys. Rev. D 88, 094506 (2013), arXiv:1307.6696 [hep-lat] .
  * Leino  _et al._ [2018] V. Leino, K. Rummukainen, J. M. Suorsa, K. Tuominen, and S. Tähtinen, Phys. Rev. D 97, 114501 (2018), arXiv:1707.04722 [hep-lat] .