NORDITA 2022-031

DESY-22-091

Speed of sound in cosmological phase transitions and
effect on gravitational waves

Tuomas V. I<EMAIL_ADDRESS>and Jorinde van de
<EMAIL_ADDRESS>
${}^{\mbox{\scriptsize a}}$Nordita, KTH Royal Institute of Technology and
Stockholm University,
Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

${}^{\mbox{\scriptsize b}}$Tsung-Dao Lee Institute & School of Physics and
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China

${}^{\mbox{\scriptsize c}}$Shanghai Key Laboratory for Particle Physics and
Cosmology, Key Laboratory for Particle Astrophysics and Cosmology (MOE),
Shanghai Jiao Tong University, Shanghai 200240, China

${}^{\mbox{\scriptsize d}}$Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany

Abstract

The energy budget for gravitational waves of a cosmological first order phase
transitions depends on the speed of sound in the thermal plasma in both phases
around the bubble wall. Working in the real-singlet augmented Standard Model,
which admits a strong two-step electroweak phase transition, we compute higher
order corrections to the pressure and sound speed. We compare our result to
lower-order approximations to the sound speed and the energy budget and
investigate the impact on the gravitational wave signal. We find that
deviations in the speed of sound from $c_{s}^{2}=1/3$ are enhanced up to
$\mathcal{O}(5\%)$ in our higher-order computation. This results in a
suppression in the energy budget of up to $\mathcal{O}(50\%)$ compared to
approximations assuming $c_{s}^{2}=1/3$. The effect is most significant for
hybrid and detonation solutions. We generalise our discussion to the case of
multiple inert scalars and the case of a reduced number of fermion families in
order to mimic hypothetical dark sector phase transitions. In this sector with
modified field content, the sound speed can receive significant suppression,
with potential order-of-magnitude impact on the gravitational wave amplitude.

###### Contents

  1. 1 Introduction
  2. 2 Model dependence of the gravitational wave signal: pressure and speed of sound
  3. 3 Computation of the pressure and sound speed
  4. 4 Impact on gravitational wave predictions
  5. 5 Speed of sound at a dark sector
  6. 6 Discussion
  7. A Dimensional reduction and pressure with $O(N)$ singlet
    1. A.1 Renormalisation
    2. A.2 Dimensional reduction
    3. A.3 Coefficient of the unit operator, or pressure, in the symmetric phase
    4. A.4 The effective potential, or pressure, in the broken phase
    5. A.5 Order parameters in perturbation theory and gauge invariance
    6. A.6 Relation to lower order computations

##  1 Introduction

#### Gravitational waves from first order phase transitions

A stochastic gravitational wave (GW) background generated by a cosmological
phase transition can offer a window to the early universe predating the
recombination era. In a first order phase transition, a barrier separates the
two degenerate minima of the free energy of the thermal plasma at the critical
temperature $T_{c}$. When the temperature ($T$) of the plasma drops below the
critical temperature, the system transitions from the high temperature to the
lower temperature phase via thermal or quantum fluctuations. Bubbles of the
low-$T$ phase nucleate, expand due to the pressure difference between the
phases inside and outside the bubble, and collide eventually, filling the
whole universe with the low-$T$ phase. A fraction of the latent heat of the
transition gets converted into kinetic energy in the form of sound waves
within the plasma, which act as a source of shear stress producing
gravitational waves, resulting in a stochastic background [1, 2, 3, 4, 5, 6].
A detection of such a hypothetical background is one of the science goals of
the space-based GW observatory LISA [7] and other similar next-generation GW
experiments.

The Standard Model (SM) of particle physics does not predict any cosmological
first order phase transitions, as both the electroweak (EW) and QCD phase
transitions occur via a smooth crossover [8, 9, 10, 11]. Hence, the existence
of a cosmic GW background would point to physics beyond the Standard Model
(BSM). The GW signal would offer a probe complementary to high energy
colliders [12, 13, 14]. Furthermore, a first order EW phase transition is a
crucial ingredient in EW baryogenesis, a dynamical mechanism aiming to explain
the matter/antimatter imbalance of the present day universe [15, 16, 17, 18,
19]. Multiple BSM theories have been proposed that are believed to incorporate
cosmic phase transitions, see e.g. the references listed in [14, 20].

#### Modelling the signal – importance of the sound speed

The gravitational wave signal of a cosmological phase transition receives
contributions from the kinetic energy stored in the bubble walls [21, 22, 23,
24], sound waves in the plasma, and turbulence [25, 26, 27, 28, 29]. In the
case of weak to moderate phase transitions, the sound wave contribution is
expected to be dominant [2, 3, 6], and it will be the main focus of this work.

The gravitational wave signal from sound waves has been determined in a range
of numerical simulations [2, 4, 5, 30]. The results of [4, 5] were used to
obtain a fitting function [6] which approximates the gravitational wave signal
as a function of the nucleation temperature ($T_{n}$), the inverse transition
duration ($\beta$), the bubble wall speed ($v_{w}$), the sound speed ($c_{s}$)
– which is obtained from the pressure, and usually assumed to be $c_{s}\sim
1/\sqrt{3}$ – and the kinetic energy fraction, or energy budget, $K$. The
value of $K$ can be determined from the fluid velocity profile of an isolated
bubble, but is usually obtained from a fit provided by [31], which gives the
value of $K$ as a function of the phase transition strength $\alpha$ and
$v_{w}$.

The results of [31] were obtained in the ‘bag model’, which assumes that both
the broken and the symmetric phase are composed of radiation degrees of
freedom only. The same assumption underlies the choice of $c_{s}=1/\sqrt{3}$.
In a more realistic model, these assumptions do not hold, as some degrees of
freedom are usually massive around the phase transition temperature. The
procedure sketched above thus neglects some of the model-dependence of the
signal, by describing it fully in terms of $T_{n},\alpha,\beta$ and $v_{w}$.
This leads to an inaccuracy in the estimated gravitational wave signal, but
also possibly to an unnecessary degeneracy between different models.

Recently, it was demonstrated in [32, 33] that further model-dependence of $K$
can, to a good approximation, be captured by a modification of the phase
transition strength (see Eq.(2.12) below) and the values of the speed of sound
of both phases of the phase transition. A code snippet provided by [33] serves
as an alternative to the fit based on the bag model by [31]. It was
demonstrated for a set of toy models, that the code snippet approximates the
numerically obtained energy budget with an error of $\lesssim 5\%$, compared
to an error of $\lesssim 50\%$ for a mapping onto the bag model.

Although some studies of the effect of the sound speed on the shape and
strength of the gravitational wave signal now exist [32, 33, 34, 35], an
accurate computation of the _value_ of the speed of sound in a first order
phase transition is still lacking. Careful studies of the sound speed thus far
have been limited to QCD [36, 37] and the SM electroweak crossover [38]. The
goal of this work is to improve the accuracy of the computation of the
pressure – for a representative BSM model that accommodates a strong
electroweak or dark sector phase transition [39] – and in particularly to
determine what are realistic values of the sound speed. To understand the
importance of our analysis, we will compare the resulting amplitude of the
gravitational wave spectrum with the result obtained in the approximation
$c_{s}=1/\sqrt{3}$. In particular, we will study the SM extended by $N$
singlets with $O(N)$ symmetry as a proxy for first order electroweak phase
transitions with an extended Higgs sector. We will also investigate the effect
of the removal of some of the light SM fermions, mimicking a dark sector with
a smaller number of degrees of freedom.

The reason that we include multiple inert singlets is that one can expect
deviations in the speed of sound when additional massive particles are
present. In particular we expect that the deviation from $c_{s}^{2}=1/3$
increases with increasing $N$. In addition, we also manually remove some of
the light degrees of freedom by reducing number of fermion families ($N_{f}$)
to one instead of three. In this case, one can expect the sound speed to drop,
because the many light degrees of freedom all push the speed of sound
(squared) towards 1/3. We indeed observe such behaviour, and can conclude that
deviations in the speed of sound can be relevant for dark sector phase
transitions.

#### BSM model

While the features that we discuss are generic for thermal phase transitions
in BSM scalar extensions, practically we work in the model where the SM is
augmented with an $O(N)$ symmetric scalar singlet [40]. The Langrangian reads

$\displaystyle\mathcal{L}$
$\displaystyle=\mathcal{L}_{\text{SM}}+\frac{1}{2}\mu^{2}_{S}S^{\text{T}}S+\frac{1}{4}\lambda_{S}(S^{\text{T}}S)^{2}+\frac{1}{2}\lambda_{m}S^{\text{T}}SH^{\dagger}H,$
(1.1)

where the singlet decomposes as $S=(S_{1},...,S_{N})^{\text{T}}$ such that
each component $S_{i}$ is real. We will focus on a two-step phase transition.
We assume that in the first step of the phase transition only one of the
components of the singlet gets a non-zero vev and participates dynamically in
the phase transition, while the other components remain inert. Therefore, for
$N=1$ this setup matches the $Z_{2}$-symmetric xSM studied in [41, 42, 20, 43]
in the context of dimensionally reduced effective field theories. We
generalise the results of these works by including the inert singlet
contributions for general $N$. As another novelty, we also compute the
coefficient of the unit operator that describes the hard and soft scale
contributions to the symmetric phase pressure. These contributions were not
discussed in earlier singlet dimensional reduction literature as therein the
intrest was purely in the pressure difference between different phases, and
corrections to the speed of sound were not determined.

#### Outline of this work

This article is organised as follows. In Sec. 2 we review the setup for
gravitational wave production, paying attention to the role of the sound
speed. In Sec. 3 we discuss the computation of the pressure and inclusion of
several higher order corrections. In Sec. 4 we present our results for
quantities important for the GW wave power spectrum. In Sec. 5 we discuss how
our findings could generalise to phase transitions in a dark sector. In Sec. 6
we discuss our conclusions and their impact. Appendix A collects multiple
technical details of our analysis.

##  2 Model dependence of the gravitational wave signal: pressure and speed of
sound

For a BSM model with a first order phase transition, the gravitational wave
signal from sound waves can be estimated as [6]

$\displaystyle\frac{d\Omega_{{\rm gw},0}}{d\ln{(f)}}=0.687F_{{\rm
gw},0}K^{2}\left(H_{*}\tau_{s}\right)H_{*}R_{*}/c_{s}\tilde{\Omega}_{\rm
gw}C\left(\frac{f}{f_{{\rm p},0}}\right),$ (2.1)

where $F_{{\rm gw},0}$ accounts for the redshift from the generation of the
signal until now, $H_{*}$ is the Hubble parameter at the phase transition
temperature, $R_{*}$ is the mean bubble separation at the moment of collision,
given by $R_{*}\sim(8\pi)^{1/3}v_{w}/\beta$. $c_{s}$ is the speed of sound in
the plasma, $\tilde{\Omega}_{\rm gw}$ is a numerical factor, $C$ is a function
that determines the shape of the spectrum and $f_{{\rm p},0}$ is the peak
frequency. $\tau_{s}$, the lifetime of the source, is estimated as [44]

$\displaystyle\tau_{s}=\begin{cases}1/H_{*}&\qquad H_{*}\tau_{\rm sh}<1\,,\\\
R_{*}/K^{1/2}&\qquad H_{*}\tau_{\rm sh}>1\,,\end{cases}$ (2.2)

where $\tau_{\rm sh}$ is the time of shock formation.111The shock formation
time can be estimated as [45, 5] $\tau_{\rm
sh}\sim\frac{R_{*}}{\bar{U}_{f}}\,,\qquad\bar{U}_{f}=\left(\frac{K}{\Gamma}\right)^{1/2}\,,$
where $\bar{U}_{f}$ is the enthalpy-weighted root-mean-square fluid velocity
and $\Gamma\sim 4/3$ is the adiabatic index.  When shock formation is
relatively slow, the gravitational wave spectrum scales with $K^{2}$, and for
fast-developing shocks the spectrum scales as $K^{3/2}$. The kinetic energy
budget $K$ is the fraction of energy that is available for the production of
gravitational waves. It can be determined by solving the hydrodynamic
equations for the plasma around a single expanding bubble. We will now give a
summary of this computation, for a more complete description, see e.g. [46,
47, 48, 31].

The hydrodynamics follow from the energy-momentum tensor of the plasma:

$\displaystyle T_{\mu\nu}=wu_{\mu}u_{\nu}-g_{\mu\nu}p,$ (2.3)

where $u_{\mu}$ is the four-velocity of the plasma and $p$ and $w$ are the
pressure and enthalpy respectively. The pressure is equal to the negative of
the free energy of the plasma, and the enthalpy and energy density $e$ can be
obtained from the pressure via

$\displaystyle w=T\frac{\partial p}{\partial T},\qquad e=T\frac{\partial
p}{\partial T}-p.$ (2.4)

The fluid equations are simply the continuity equations
$\partial^{\mu}T_{\mu\nu}$. By projecting these along the directions parallel
and perpendicular to the fluid flow, and by introducing the radial coordinate
$\xi=r/t$, with $r$ the distance from the center of the bubble and $t$ the
time since nucleation, the following hydrodynamic equations are found

$\displaystyle\frac{dv}{d\xi}$
$\displaystyle=\frac{2v(1-v^{2})}{\xi(1-v\xi)}\left(\frac{\mu(\xi,v)^{2}}{c_{s}^{2}}-1\right)^{-1},$
(2.5) $\displaystyle\frac{dw}{d\xi}$
$\displaystyle=w\left(1+\frac{1}{c_{s}^{2}}\gamma^{2}\mu(\xi,v)\frac{dv}{d\xi}\right),$
(2.6)

where $v(\xi)$ is the fluid velocity, $\gamma$ the Lorentz factor and
$\mu(\xi,v)$ the boosted velocity

$\displaystyle\mu(\xi,v)=\frac{\xi-v}{1-\xi v}.$ (2.7)

The speed of sound can be derived from the pressure via

$\displaystyle c_{s}^{2}=\frac{dp/dT}{de/dT}.$ (2.8)

We compute the sound speed separately for each metastable phase.222
Consequently, $c_{s}^{2}$ does not vanish at the transition point, in contrary
to the prescription in [49].  The boundary conditions at the bubble wall are
given by

$\displaystyle\frac{v_{+}}{v_{-}}$
$\displaystyle=\frac{e_{b}(T_{-})+p_{s}(T_{+})}{e_{s}(T_{+})+p_{b}(T_{-})},$
(2.9) $\displaystyle v_{+}v_{-}$
$\displaystyle=\frac{p_{s}(T_{+})-p_{b}(T_{-})}{e_{s}(T_{+})-e_{b}(T_{-})},$
(2.10)

where the $+(-)$ sign denotes a quantity right in front of (behind) the bubble
wall. The subscript $s(b)$ refers to the equation of state of the symmetric
(broken) phase. The boundary condition for $w$ is that $w=w_{n}$ in the region
in front of the bubble where the fluid is at rest. We use the subscript $n$ to
denote quantities evaluated at the nucleation temperature $T_{n}$. The kinetic
energy fraction $K$ is obtained from the solution of the hydrodynamic
equations via

$\displaystyle K=\frac{3}{e_{n}v_{w}^{3}}\int d\xi\xi^{2}v^{2}\gamma^{2}w,$
(2.11)

where $v_{w}$ is the bubble wall velocity.

Although both the fluid equations (2.5) and (2.6) and the boundary conditions
(2.9) and (2.10) depend on the equation of state, Refs. [32, 33] show that
this dependence can, to a good approximation, be captured by the following
parameters only: the speed of sound in the broken phase $c_{s,b}$, the speed
of sound in the symmetric phase $c_{s,s}$ (both evaluated at the nucleation
temperature) and the phase transition strength

$\displaystyle\alpha_{\bar{\theta}}=\frac{D\bar{\theta}(T_{n})}{3w_{n}},\qquad{\rm
with}\quad\bar{\theta}=e-\frac{p}{c_{s,b}^{2}},\qquad{\rm and}\quad
D\theta(T_{n})=\theta_{s}(T_{n})-\theta_{b}(T_{n}).$ (2.12)

Consequently, instead of solving the full hydrodynamic equations for each BSM
model, $K$ can be obtained as a function of $c_{s,b},c_{s,s}$, $v_{w}$ and
$\alpha_{\bar{\theta}}$ in a simplified model using the code snippet provided
by [33]. This procedure is a generalisation of the fit in terms of the phase
transition strength and $v_{w}$ provided by Ref. [31], using the bag equation
of state. The new approximation works particularly well when the temperature-
dependence of the sound speed is weak, but approximates the full numerical
solution to an accuracy of $<5\%$ for all toy models considered in [32, 33].

To quantify the importance of an accurate computation of the sound speed, in
Sec. 4 we will determine the kinetic energy fraction for our benchmark point
using three different methods (these methods are identical to methods 1-3 in
[33]):

  * •

Method 1: full solution of Eqs. (2.5) and (2.6), using the complete,
temperature-dependent speed of sound.

  * •

Method 2: solution of Eqs. (2.5) and (2.6), assuming constant sound speeds in
the broken and symmetric phases.

  * •

Method 3: assuming the speed of sound is $c_{s}=1/\sqrt{3}$. This corresponds
to a mapping onto the bag equation of state. Note that this is the most common
method in the literature, and it was also used in Ref.[50], where the
uncertainty in the GW spectrum associated to the other thermodynamic
quantities was estimated.

The difference between Method 1 and Method 3 gives a measure of the importance
of a careful computation of the speed of sound. The difference between Method
1 and Method 2 gives a measure of the temperature-dependence of the speed of
sound.

In this work we will assume that the gravitational wave spectrum is fully
described by the broken power law of Eq. (2.1), and only investigate the
effect of our improved sound speed computation on the overall amplitude.
According to Eq. (2.1), the position of the peak is determined by the size of
the bubbles at the moment of collision. However, results of the sound shell
model [51, 52] and hybrid simulations [53] suggest that the width of the sound
shells also gets imprinted onto the spectrum, resulting in a doubly broken
power law shape. The results of [51, 52, 53] were all obtained under the
assumption that the sound speed is $c_{s}=1/\sqrt{3}$. Since deviations on the
sound speed affect the shape of the fluid profile [54, 33] it is expected that
the shape of the gravitational wave spectrum is also affected by variations in
the sound speed. This effect was demonstrated for the sound shell model in
[35], but will not be considered further in our analysis.

##  3 Computation of the pressure and sound speed

In this section we review how to compute the pressure perturbatively in high
temperature field theory. We present this discussion for the SM augmented with
$N$ singlets with O($N$) symmetry, but this discussion readily generalises to
other BSM theories. In fact, many of the computations presented in this work
can be obtained with the automated DRalgo package for generic, user defined
models [55]. Many details of our computation in the model in question are
relegated to Appendix A.

#### Using dimensional reduction to determine the phase transition parameters

In order to improve the accuracy of the computation of the sound speed, we
will use the method of high temperature dimensional reduction [56, 57].
Therein, the perturbative computation is organised in the effective field
theory language [58, 59], based on a chain of scale hierarchies at high
temperature. In the rest of this article, we will refer to the different
scales as hard ($\pi T$), soft ($gT$) and ultrasoft ($g^{2}T$), where $g$ is a
weak coupling, often identified with the weak gauge coupling in EW
theories.333 In QCD the soft and ultrasoft scales are often referred to as the
electric and magnetic scales, respectively.  This method was recently applied
in [50, 60] to obtain higher order corrections of $T_{n}$, $\alpha$ and
$\beta$. Ref. [50] reported an alarming leftover renormalisation scale
dependence plaguing the conventional analysis based on a mere one-loop
determination of the effective potential – that describes the free energy and
pressure of the plasma – which in general signals slow convergence of
perturbation theory. In particular, there could be a multiple order-of-
magnitude uncertainty in the predicted amplitude of the GW power spectrum,
location of the peak frequency, and the consequent signal-to-noise (SNR) ratio
for experiments such as LISA. In [60], it was reviewed and discussed in detail
why conventional one-loop analyses fail to provide better accuracy: due to
slower convergence of the perturbative expansion at high temperatures, several
two-loop contributions have to be included to enable leading renormalisation
group (RG) improvement. After such RG improvement, order-by-order in
perturbation theory, the renormalisation scale-dependence cancels between the
running of lower order contributions and explicit logarithms at higher orders.
The inclusion of these two-loop level contributions is standard in a three-
dimensional effective field theory (3d EFT) description of the phase
transition thermodynamics, from which equilibrium properties of the
transition, such as the pressure, and subsequent quantities of interest such
as the critical temperature, latent heat, and the speed of sound can be
derived.

In this work at hand, we focus on the computation of the pressure and the
speed of sound, and investigate the importance of its higher order corrections
to the energy budget $K$, and hence to the GW power spectrum. In practice, we
compute the pressure by following [61, 62, 38], and extend the SM computation
of these works to a BSM setup. For a similar, related discussion of QCD and
the speed of sound in a quark-gluon plasma, see e.g [36, 37] (c.f. also the
review [63]).

Although a necessary improvement in the accuracy of equilibrium quantities can
be made by including higher order corrections to the pressure, or the
effective potential, it was concluded in [60] that the most limiting source of
theoretical uncertainty in the phase transition parameters is the bubble
nucleation rate. This rate is used to derive $T_{n}$ and all other parameters
have to be evaluated at this temperature. In [60], the bubble nucleation rate
was determined only at leading order within the 3d EFT. However, recent work
has described how to extend such a computation to higher orders [64, 65, 66,
67] (also c.f. [68, 69, 70]). In the work at hand, we will not focus on such
an improvement for the bubble nucleation rate, but instead treat $T_{n}$ as an
input parameter that we vary; once its computation at higher orders is
available, our results for the pressure, sound speed and the kinetic energy
fraction can be used to determine the actual effect on GW power spectra.

#### Pressure in perturbation theory

In perturbation theory, the pressure admits the following form and subsequent
formal expansion444 The analogous expansion in QCD has been worked out, order
by order: $g^{2}$ [71, 72], $g^{3}$ [73], $g^{4}\ln g$ [74], $g^{4}$ [75, 76],
$g^{5}$ [77, 78] and $g^{6}\ln g$ [79, 80]. Here, logarithmic contributions
are highlighted separately, since they appear in a result expanded in $g$,
when the 3d EFT formalism is not used. For the case of non-zero quark chemical
potentials, see [81].  in the high temperature limit:

$\displaystyle p(T)=p_{\text{sym}}-V_{\text{eff}}(v,s)\simeq
T^{4}\Big{(}a+b_{1}g^{2}+cg^{3}+dg^{4}+eg^{5}+\mathcal{O}(g^{6})\Big{)}+m^{2}T^{2}\Big{(}b_{2}g^{2}\Big{)},$
(3.1)

where $p_{\text{sym}}$ describes the pressure in the symmetric phase where all
scalar background fields vanish and $V_{\text{eff}}$ is the effective
potential as a function of the background fields, and $a$-$f$ represent
constants. The $\overline{\rm MS}$ mass parameters are assumed to scale as
$m^{2}\sim(gT)^{2}$.

We comment below on our conventions for the normalisations of both terms. A
formal expansion in $g$ follows when a power counting rule is associated to
each coupling. In particular we assume that all quartic couplings $\lambda$
scale as $g^{2}$, and we assume that there is a barrier present in the leading
order effective potential; for radiatively generated transitions, the
corresponding expansion is different [82].555 For example, for a radiatively
generated transition, the 3d scalar mass parameter scales as $\mu^{2}_{3}\sim
g^{3}T^{2}$, which leads to fractional powers $g^{4.5}$ and $g^{5.5}$ coming
from the one-loop ultrasoft contribution [62].  Compared to the perturbative
expansion of the effective potential at zero temperature, the structure of the
perturbative expansion of the pressure is inflicted by several finite-$T$
peculiarities: odd powers of $g$ appear, logarithmic terms include ratios of
various thermal mass scales, and crucially the $\mathcal{O}(g^{6})$ term is
not even attainable in perturbation theory, due to the Linde problem [83]. In
addition, the coupling expansion misaligns with the loop expansion due to the
enhancement of IR contributions.

We emphasise, however, that in practice the pressure is not computed directly
order-by-order in the above expansion in $g$. The reason is, that the pressure
is not directly computable in such an expansion, since the direct use of
perturbation theory is inhibited by various IR singularities. These have to be
tackled by thermal resummations, based on a scale hierarchy of the different
mass scales (hard/soft/ultrasoft) at high temperature. The perturbative result
will include logarithms of all three mass scales in the problem, and it is not
possible to choose a UV cutoff that would simultaneously remove all large
logarithms. In particular, these logarithms are of the form $\ln(\mu/T)$,
$\ln(T/m_{\mbox{\tiny\rm{E}}})$ and
$\ln(m_{\mbox{\tiny\rm{E}}}/m_{\mbox{\tiny\rm{M}}})$, where
$m_{\mbox{\tiny\rm{E}}}\sim gT$ and $m_{\mbox{\tiny\rm{M}}}\sim g^{2}T$ are
the soft and ultrasoft – or electric (E) and magnetic (M) – mass scales. These
logarithms lead formally to orders $g^{4}\ln(g)$ and $g^{6}\ln(g^{6})$, where
logarithms are large for $g\rightarrow 0$.

Dimensional reduction to a chain of EFTs resolves this problem systematically:
the thermal contributions of the hard scale are resummed as contributions to
the parameters of the soft and ultrasoft scales. The running within the EFT
introduces new RG scales ($\mu_{3}$ and $\bar{\mu}_{3}$), which can be chosen
independently, allowing to avoid large logarithms. In our computation, we
associate $p_{\text{sym}}$ with the coefficient of the unit operator [59] in
dimensional reduction, and this contribution describes the hard and soft
contributions to the symmetric phase pressure. Contributions of the ultrasoft
scale are, by our convention, encoded in $V_{\text{eff}}$, for both symmetric
and broken phases.

Let us take a closer look at how different order terms in the pressure arise,
and how we label them:

  * LO ($ag^{0}$): one-loop hard contributions to $p_{\text{sym}}$.

  * NLO ($bg^{2}$): tree-level terms in $V_{\text{eff}}$ and two-loop hard pieces in $p_{\text{sym}}$ .

  * NNLO ($cg^{3}$): one-loop soft (ultrasoft) terms in $p_{\text{sym}}$ ($V_{\text{eff}}$).

  * $\text{N}^{3}$LO ($dg^{4}$): two-loop (three-loop) soft (hard) pieces in $p_{\text{sym}}$ and two-loop ultrasoft pieces in $V_{\text{eff}}$. Tree-level $V_{\text{eff}}$ includes contributions at this order, via the resummed parameters in the EFT.

  * $\text{N}^{4}$LO ($eg^{5}$): three-loop soft (ultrasoft) contributions to $p_{\text{sym}}$ ($V_{\text{eff}}$).

Note that terms that are formally of higher order are resummed inside the 3d
EFT parameters in different contributions and, in accord with the EFT
construction, we will keep them untruncated in the final result [84, 36]. We
remind that the $\mathcal{O}(g^{6})$-contribution is non-perturbative,
although the $\mathcal{O}(g^{6}\ln g)$ logarithmic term is still computable at
4-loop level [79, 80].

The LO result for the sound speed, $c^{2}_{s}=1/3$, follows from the LO result
for the pressure. Typical phase transition analyses, based on the one-loop
daisy resummed thermal effective potential (c.f. Appendix A.6), produce a
pressure that is correct to NNLO at $\mathcal{O}(g^{3})$ for the broken phase
contributions, but are lacking two-loop $\mathcal{O}(g^{2})$ terms in the
symmetric phase. Also, such analyses include some $\mathcal{O}(g^{4})$
corrections, but are not complete at $\text{N}^{3}$LO. In our computation in
this article at hand, we provide complete accuracy at $\text{N}^{3}$LO: this
requires a three-loop computation for the hard contributions, and a two-loop
computation of the soft contributions and the broken phase effective
potential; technical details can be found in Appendix A. At $\text{N}^{4}$LO
the $\mathcal{O}(g^{5})$ correction would require a three-loop computation in
the 3d EFT, which we do not include here, but leave to future work. It is
tempting to pursue this correction, as it is the last term that can still be
obtained in perturbation theory, c.f. [79].

Figure 1:  Minima of the real part of the effective potential as a function of
the temperature, determined in a one-loop approximation (left) and at two-loop
within the 3d EFT with NLO dimensional reduction (right). Blue (red) denotes
the Higgs (singlet) direction, and dots (crosses) present global (local)
minima. Local minima describe the metastable phases. The different dark and
light colours present two choices for the 4d RG scale $\mu$: $0.5\pi T$ and
$2\pi T$. The vertical lines depict the critical temperatures for the second
transition ($T_{c,\phi}$), and the band between them the theoretical
uncertainty in their determination. We note that the critical temperatures for
the first transition ($T_{c,S}$) is very sensitive to two-loop corrections,
while $T_{c,\phi}$ is not.  Figure 2:  Pressure of each phase as function of
temperature, normalised by the leading order pressure $p_{0}$. Bands
correspond to varying the 4d RG scale $\mu$ between $0.5\pi T$ and $2\pi T$.
Darker colours correspond to $\text{N}^{3}$LO computation, and lighter colours
to a mere one-loop approximation. The result of the $\text{N}^{3}$LO
computation is less sensitive to the RG scale, and corrections to $p_{0}$ in
all phases are mildly smaller than in the mere one-loop approximation.  Figure
3:  Percentual deviation in the speed of sound from $c_{s}^{2}=1/3$ in the
singlet phase (red) and Higgs phase (blue) as a function of the temperature.
The horizontal line indicates the value of a relativistic gas of particles.
Colour codes are the same as in Fig. 2

#### Numerical study

For our numerical analysis, we select a single, representative benchmark point
with a two-step phase transition [85, 86, 87, 88, 89] with a strong second
transition to the Higgs phase.

$\displaystyle\text{BM:}\quad(m_{S},\lambda_{S},\lambda_{m})$
$\displaystyle=(160\,\text{GeV},1.0,1.6),$ (3.2)

where $m_{S}$ is the singlet mass.

The order parameters of the different phases are depicted in Fig. 1. The Higgs
(singlet) background field is denoted by $v$ ($s$). The pressure and sound
speed – computed from Eq. (2.8) – are shown in Figs. 2 and 3, respectively. In
all these plots, we present a comparison of our full $\text{N}^{3}$LO
computation with the frequently used, sole one-loop approximation [90] (cf.
Appendix A.6). In addition, we vary the RG scale $\mu$ from $0.5\,\pi T$ to
$2\,\pi T$, to monitor convergence, and we find in all the aforementioned
figures that in the full $\text{N}^{3}$LO computation the sensitivity to the
4d theory RG scale is milder, as expected. For the 3d RG scales we have used
$\bar{\mu}_{3}=\mu_{3}=T$ in Fig. 1 and
$\bar{\mu}_{3}=\mu_{3}=1.25\bar{g}^{2}_{3}$ in Figs. 2 and 3. We note that our
use of minima of the effective potential as order parameters leads to an
undesired gauge-dependence that signals an inconsistent perturbative
expansion. We carefully comment on this issue in Appendix A.5 and argue why
this does not compromise our discussion regarding the sound speed.

We find that the deviations in the speed of sound squared from $c_{s}^{2}=1/3$
are small in the singlet phase, where they remain subpercent over the entire
temperature range for the one-loop result. Deviations are somewhat more
apparent in the $\text{N}^{3}$LO analysis, of the order of $\gtrsim 1\%$. Note
that the regime where the deviations are largest is the regime with most
supercooling, where the high-temperature expansion might not be reliable. In
the Higgs phase, the one-loop result deviates from the relativistic value by
$\sim 2\%$ over a large range of temperatures. Again, the deviations in the
$\text{N}^{3}$LO analysis are larger, $\sim 3-4\%$. In both cases, the speed
of sound exceeds the relativistic value in the less trustworthy regime of
strong supercooling.

##  4 Impact on gravitational wave predictions

Figure 4:  Ratio of the kinetic energy fraction computed in Method 2 and
Method 3 to the result of Method 1. The left plot shows a hybrid with
$v_{w}=0.65$ and the right plot a detonation with wall velocity $v_{w}=0.99$.

We have seen in Sec. 2 that the gravitational wave spectrum depends on the
sound speed through the kinetic energy fraction $K$ and through an explicit
appearance in Eq. (2.1). Fig. 3 demonstrates that, for $N=1$, the deviations
from $c_{s}=1/\sqrt{3}$ are moderate, so we do not expect major deviations in
the gravitational wave spectrum from the explicit dependence on $c_{s}$. As
demonstrated in [33] however, a small deviation in the sound speed can already
lead to a significant effect on the kinetic energy budget. Bear in mind that
the gravitational wave amplitude scales with $K^{2}$ or $K^{3/2}$, which
amplifies the dependence. We will quantify the effect of the accurate
computation of the sound speed on the gravitational wave spectrum by comparing
the kinetic energy budget computed in the bag model (Method 3) with the
computations accounting for the deviations in the sound speed (Method 1 and
2). In all following plots, we demonstrate the results obtained in the 3d EFT
with fixed RG scales $\mu=1.25\;\pi T$ and
$\mu_{3}=\bar{\mu}_{3}=1.25\;\bar{g}^{2}_{3}$.

It should be noted that an accurate prediction of the gravitational wave
spectrum depends on the nucleation temperature, which enters in
$H_{*},R_{*},f_{\rm p,0}$ and $K$. A self-consistent determination of the
nucleation temperature requires a computation of the bubble nucleation rate
within the 3d nucleation EFT [65] (for applications, see [64, 69]). This goes
beyond the scope of this work. Instead we will present our results for $K$ as
a function of an undetermined nucleation temperature.

One distinguishes three types of hydrodynamic solutions, characterised by the
wall velocity: (subsonic) deflagrations, in which the fluid is at rest inside
the bubble, and forms a shock wave in front of the bubble wall; (supersonic)
detonations, for which the fluid in front of the bubble is at rest and a
rarefaction wave forms inside; hybrids, or supersonic deflagrations, which
consist of a rarefaction wave and a shock. We treat the wall velocity as an
external parameter, but in reality, the computation of its value is
challenging. Recently, it has been argued that solutions likely fall into one
of two categories: deflagrations and ultrarelativistic detonations, but that
no solutions exist between the Jouguet velocity and wall velocities with
$\gamma_{w}\lesssim 10$ [91, 92].

Figure 5:  Absolute value of $K$ as a function of the nucleation temperature
for a hybrid ($v_{w}=0.65$). The black lines indicate the results of the
$\text{N}^{3}$LO analysis and the gray lines show the result of the one-loop
analysis. The dashed lines show the results based on a mapping onto the bag
equation of state.

We have determined the ratio of the kinetic energy fraction $K$ computed in
Method 2, $K_{2}$ and Method 3, $K_{3}$, to the result in Method 1, $K_{1}$
for our benchmark point. The results are demonstrated in Fig. 4 for a hybrid
and detonation solution. The temperature has been rescaled by the critical
temperature of the transition to the singlet to the Higgs phase. We found that
for the deflagration solution, the differences between the three methods are
very small ($<\mathcal{O}(1\%)$) for the entire temperature range. This
implies that the gravitational wave spectrum can be estimated accurately with
$c_{s}=1/\sqrt{3}$ and a mapping onto the bag equation of state. For a hybrid,
the story is different. Especially when the nucleation temperature is
relatively close to the critical temperature,666Note that the upper limit of
$T/T_{c,\phi}$ for the hybrid is smaller than unity, since the solutions with
$T/T_{c,\phi}\gtrsim 0.9$ and $v_{w}=0.65$ correspond to detonations. the
difference between $K_{3}$ and the other methods is more pronounced,
approaching $\mathcal{O}(40\%)$. The estimate based on the bag equation of
state typically overestimates the kinetic energy fraction. Over the entire
temperature range, $K_{2}$ reproduces Method 1 very well, meaning that the
temperature dependence of the sound speed plays an insignificant role, and a
full numerical solution of the hydrodynamics is not required. For detonations,
the effect of the sound speed is again smaller, with a maximum of
$\mathcal{O}(8\%)$ for $T_{n}\sim T_{c,\phi}$. For significant supercooling,
we see that $K_{2}$ starts to deviate slightly from $K_{1}$, implying that the
temperature-dependence of the sound speed mildly affects the result.

Fig. 5 showcases the absolute value of $K$ for a hybrid solution as a function
of the nucleation temperature, for the $\text{N}^{3}$LO (black) and sole one-
loop result (grey). The solid line displays $K_{1}$, and the dashed line
$K_{3}$. We see that the difference between $K_{1}$ and $K_{3}$ is larger in
the $\text{N}^{3}$LO analysis, which could be expected from Fig. 3. We see
that, for a given nucleation temperature, the $\text{N}^{3}$LO computation
yields a smaller value of $K$ than the one-loop analysis (but of course the
two approaches do not yield the same nucleation temperature). In addition, our
accurate treatment of the speed of sound suppresses the kinetic energy
fraction compared to the treatment with $c_{s}^{2}=1/3$. For a given
temperature, the determination of the pressure at one loop and a determination
of $K$ via the bag model together overestimate $K$ by $3\%$ in the strong
supercooled regime up to $50\%$ in the regime around $0.9T_{c,\phi}$.

##  5 Speed of sound at a dark sector

Figure 6:  Sound speed as a function of temperature in the singlet (left) and
the Higgs phases (right). The solid lines correspond to the case with all
fermionic families ($N_{f}=3$) for $N=1,3,5$, and dashed lines show result
with only one fermionic family ($N_{f}=1$). The endpoints of the graphs are
set by the critical temperature $T_{c,\phi}$.

We have seen in the previous section that the sound speed only deviates from
$c_{s}^{2}=1/3$ by a few percent. The main reason for this is that the light
SM fermions give the dominant contribution to the pressure. In this section,
we study how the the sound speed behaves when the particle content deviates
further from the SM. We investigate two different effects:

  * •

Increasing the number of singlets $N$.

  * •

Reducing the number of fermions by setting $N_{f}=1$. This modification mimics
a dark sector with a smaller amount of fermions than the SM.

Fig. 6 demonstrates the result for the sound speed. The solid lines show the
result for the SM + $N$ singlets, with the usual SM fermionic content. As
expected, increasing the number of singlets suppresses the sound speed. Again,
the deviations are largest in the broken phase, and approach $5\%$. The dotted
lines demonstrate the sound speed in the dark sector with only 1/3 of the
fermionic content of the SM. We see that this leads to an even stronger
suppression in the sound speed, approaching $8\%$ in the singlet phase and
$15\%$ in the Higgs phase. This is a significant result, as deviations of that
size can suppress the gravitational wave signal by an order of magnitude [33,
35] and also affect its shape. An even stronger suppression is expected if the
particle content gets reduced further. From this simple exercise we thus
conclude that especially in dark sectors with a particle content that strongly
deviates from the SM, a computation of the gravitational wave spectrum
requires a careful computation of the sound speed.

##  6 Discussion

In this work, we have computed the pressure and the corresponding sound speed
at $\text{N}^{3}$LO in formal expansion in $g$ (c.f. Eq. (3.1) and discussion
below it) in the Standard Model augmented by $N$ singlets with
$\mathcal{O}(N)$ symmetry. The speed of sound enters in the computation of the
energy budget of gravitational waves. Whenever the sound speed deviates from
the often-assumed value $c_{s}^{2}=1/3$, the amplitude and shape of the
gravitational wave spectrum are modified. In order to obtain $\text{N}^{3}$LO
accuracy, we used a dimensionally reduced EFT. We worked with a benchmark
point defined in Eq. (3.2) giving rise to a two-step phase transition. We
expect our qualitative results to hold for other parameter choices, and other
BSM models that accommodate similar strong two-step phase transitions,
relevant for observable GW wave signatures.

As known from earlier comparisons between the mere one-loop analysis and
higher order analysis in the 3d EFT, thermodynamic quantities such as the
critical temperature and latent heat differ significantly between the two
computations, potentially leading to multiple order of magnitude differences
in the gravitational wave signal. The one-loop result suffers from a
disturbing leftover RG-scale dependence, which gets significantly reduced in
the result obtained from the 3d EFT. We observe the same reduction in the RG-
scale dependence in the pressure and the sound speed. In the $\text{N}^{3}$LO
analysis, suppressions of the sound speed compared to the LO value
$c_{s}^{2}=1/3$ are more significant than in the mere one-loop computation.
This suppression affects the energy budget, which sets the amplitude of the
gravitational wave signal: in a phase transition with $c_{s}^{2}\neq 1/3$ the
energy budget should be computed with the methods of [32, 33] instead of the
commonly-used bag model. The error associated with using the bag model or an
inaccurately computed $c_{s}^{2}$ is not as large as the error associated to
some of the other parameters of the phase transition[50, 60]. Nevertheless,
especially in the case of hybrid and detonation solutions, the error is non-
negligible, and can be as high as $\mathcal{O}(50\%)$ for the energy budget.

In Sec. 5 we further modified the particle content, by increasing the number
of singlets, and by removing 2/3 of the fermionic content. These modifications
allowed us to envision the possible behaviour of the sound speed in phase
transition in dark sectors. We found that both an increase in $N$, as well as
the removal of the 2/3 of the fermions cause larger deviations in the pressure
compared to its LO value. As a result, the sound speed decreases with $N$, and
especially small values are obtained in the fewer fermion case. The effect was
strongest in the Higgs phase, where the sound speed could decrease by almost
$15\%$. Such strong deviations from $c_{s}^{2}=1/3$ lead to a significant
suppression in the gravitational wave signal, by possibly an order of
magnitude, and also significantly modify its shape. We thus conclude that
computations of the gravitational wave signal from phase transitions in a dark
sector require an accurate computation of the sound speed in order to even
obtain the right order of magnitude of the signal.

## Acknowledgements

We thank Sebastian Bruggisser, Andreas Ekstedt, Oliver Gould, Lauri Niemi,
Philipp Schicho, Aleksi Vuorinen and Juuso Österman for illuminating
discussions and Mikko Laine for reading and commenting on the manuscript. TT
is supported by National Science Foundation of China grant no. 19Z103010239.
JvdV is supported by the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306, and thanks
Nordita for the hospitality during the completion of this work.

##  Appendix A Dimensional reduction and pressure with $O(N)$ singlet

In this appendix we collect multiple details of our computation, in particular
the dimensional reduction matching relations for the SM + $O(N)$ singlet, as
well as the computation of the pressure. For the SM parameters, we use the
same power counting in terms of $g$ as in [41, 20], and in analogy to these
references we use

$\displaystyle\mu^{2}_{S}\sim(gT)^{2},\quad\lambda_{S},\lambda_{m}\sim g^{2}.$
(A.1)

The detailed structure and exact definitions of the 3d EFTs at the soft and
ultrasoft scales are discussed in [41, 20, 43]. We do not repeat the
discussion here, as the only new aspect is the singlet having $N$ components
instead of one. Parameters of the 3d EFTs will be denoted with an additional
subscript 3, and parameters of the ultrasoft theory with an additional bar. We
emphasise, that compared to [62] in which the SM pressure is computed, we use
a different power counting for the 3d Higgs mass parameter. In this reference,
in the context of pure SM, the proper counting is $\mu^{2}_{h,3}\sim
g^{3}_{3}$, which corresponds to a radiatively generated barrier (see for a
similar discussion the recent work [69, 82]), whereas here we use
$\mu^{2}_{h,3}\sim g^{2}_{3}$ which corresponds to the case with a tree-level
barrier in a two-step phase transition of the Higgs and singlet.

###  A.1 Renormalisation

We start with the zero $T$ renormalisation. All the relevant counterterms are
listed in [41] and [20], where the case of SM with a single additional singlet
was considered. Here, we consider the case of $N$ singlets, and list the
counterterms that get modified by the presence of these additional singlets.
We use the modified minimal subtraction scheme ($\overline{\rm MS}$).

$\displaystyle\delta\lambda_{h}$
$\displaystyle=\frac{1}{16\pi^{2}\epsilon}\left(\frac{3}{16}(3g^{4}+2g^{2}{g^{\prime}}^{2}+{g^{\prime}}^{4})-3g_{Y}^{4}+12\lambda_{h}^{2}+\frac{N}{4}\lambda_{m}^{2}-\frac{1}{2}\lambda_{h}(3g^{2}\xi_{2}+{g^{\prime}}^{2}\xi_{1})\right)\,,$
(A.2) $\displaystyle\delta\lambda_{S}$
$\displaystyle=\frac{1}{16\pi^{2}\epsilon}\left(\lambda_{m}^{2}+\lambda_{S}^{2}(N+8)\right)\,,$
(A.3) $\displaystyle\delta\lambda_{m}$
$\displaystyle=\frac{1}{16\pi^{2}\epsilon}\lambda_{m}\left(6\lambda_{h}+2\lambda_{m}+(N+2)\lambda_{S}\right)\,,$
(A.4) $\displaystyle\delta\mu^{2}_{h}$
$\displaystyle=\frac{1}{16\pi^{2}\epsilon}\left(6\lambda_{h}\mu_{h}^{2}+\frac{N}{2}\lambda_{m}\mu_{S}^{2}\right)\,,$
(A.5) $\displaystyle\delta\mu_{S}^{2}$
$\displaystyle=\frac{1}{16\pi^{2}\epsilon}\left((N+2)\lambda_{S}\mu_{S}^{2}+2\lambda_{m}\mu_{h}^{2}\right)\,,$
(A.6)

with $\lambda_{h}$ the Higgs quartic coupling and $\mu_{h}$ its mass
parameter, $g$ and $g^{\prime}$ the SU(2) and U(1) gauge couplings, $g_{Y}$
the Yukawa coupling and $\xi_{2}$ and $\xi_{1}$ the SU(2) and U(1) gauge
fixing parameters, in general covariant gauge, respectively.

The one-loop RG-equations, which encode the running of parameters as a
function of the RG-scale $\mu$ to order $\mathcal{O}(g^{4})$, read

$\displaystyle\mu\frac{d}{d\mu}\lambda_{h}$
$\displaystyle=\Big{(}\mu\frac{d}{d\mu}\lambda_{h}\Big{)}_{\text{SM}}+\frac{1}{(4\pi)^{2}}\Big{(}\frac{1}{2}N\lambda^{2}_{m}\Big{)},$
(A.7) $\displaystyle\mu\frac{d}{d\mu}\mu^{2}_{h}$
$\displaystyle=\Big{(}\mu\frac{d}{d\mu}\mu^{2}_{h}\Big{)}_{\text{SM}}+\frac{1}{(4\pi)^{2}}\Big{(}N\lambda_{m}\mu^{2}_{S}\Big{)},$
(A.8) $\displaystyle\mu\frac{d}{d\mu}\lambda_{S}$
$\displaystyle=\frac{1}{(4\pi)^{2}}\Big{(}2\lambda^{2}_{m}+2(N+8)\lambda^{2}_{S}\Big{)},$
(A.9) $\displaystyle\mu\frac{d}{d\mu}\lambda_{m}$
$\displaystyle=\frac{2}{(4\pi)^{2}}\lambda_{m}\Big{(}-\frac{3}{4}(3g^{2}+{g^{\prime}}^{2})+3g^{2}_{Y}+2\lambda_{m}+6\lambda_{h}+(N+2)\lambda_{S}\Big{)},$
(A.10) $\displaystyle\mu\frac{d}{d\mu}\mu^{2}_{S}$
$\displaystyle=\frac{1}{(4\pi)^{2}}\Big{(}2(N+2)\lambda_{S}\mu^{2}_{S}+4\lambda_{m}\mu^{2}_{h}\Big{)}.$
(A.11)

For $N=1$, they agree with the RG-equations of [41], but be aware of the
different convention for the sign of $\mu_{h}^{2}$ and $\mu_{S}^{2}$.

All $\overline{\rm MS}$ parameters in the Lagrangian can be related to pole
masses and other physical parameters as in [58, 93, 43], but in this work at
hand, we merely use leading, tree-level relations without one-loop
corrections.

###  A.2 Dimensional reduction

For the dimensional reduction procedure, the Feynman rules of the symmetric
phase are listed in [41]. Here we only list the Feynman rules for the singlet.
The singlet propagator reads

$\displaystyle\langle
S_{a}(P)S_{b}(Q)\rangle=\delta_{ab}\frac{\raise-0.20004pt\hbox{$\bar{}$}{\delta}(P+Q)}{P^{2}+\mu^{2}_{S}},$
(A.12)

where $P,Q$ denote Euclidean four-momentum and
$\,\raise-0.20004pt\hbox{$\bar{}$}{\delta}(K)\equiv
T^{-1}\delta_{K_{0},0}(2\pi)^{d}\delta^{(d)}({\bf k})$ where
$\delta_{K_{0},0}$ is the Kronecker delta for vanishing zero component. The
vertex Feynman rules read

$\displaystyle V(S_{a}S_{b}S_{c}S_{d})$
$\displaystyle=-2\lambda_{S}(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}),$
(A.13) $\displaystyle V(S_{a}S_{b}\phi^{\dagger}_{i}\phi_{j})$
$\displaystyle=-\lambda_{m}\delta_{ij}\delta_{ab},$ (A.14)

where the indices $a,b,c,d$ label singlet components, for which
$\delta_{aa}=N$, while $i,j$ are fundamental SU(2) indices for the Higgs. The
3d matching relations, or short distance coefficients, can be obtained
following [58, 59, 41, 42, 20, 43].777 At the time of our computation,
awakening of DRalgo [55] in its lair was still in the future.  We list the
matching relations that involve the singlet fields here at NLO, i.e. at
$\mathcal{O}(g^{4})$:

$\displaystyle\lambda_{h,3}=$
$\displaystyle\Big{(}\lambda_{h,3}\Big{)}_{\text{SM}}+\frac{1}{(4\pi)^{2}}L_{b}(\mu)\Big{(}-\frac{1}{4}N\lambda_{m}^{2}\Big{)},$
(A.15) $\displaystyle\mu^{2}_{h,3}(\mu_{3})=$
$\displaystyle\Big{(}\mu^{2}_{h,3}(\mu_{3})\Big{)}_{\text{SM}}+\frac{T^{2}}{24}N\lambda_{m}(\mu)+\epsilon\;T^{2}\beta_{\lambda_{m}2}\lambda_{m}-\frac{L_{b}(\mu)}{(4\pi)^{2}}\Big{(}\frac{1}{2}N\lambda_{m}\mu_{S}^{2}\Big{)}$
$\displaystyle+\frac{1}{(4\pi)^{2}}\bigg{(}\frac{3}{4}(3g^{2}+g^{\prime
2})L_{b}(\mu)-3g_{\mbox{\tiny\rm{$Y$}}}^{2}L_{f}\bigg{)}\bigg{(}\frac{T^{2}}{24}N\lambda_{m}\bigg{)}-\frac{1}{(4\pi)^{2}}\frac{1}{2}N\lambda^{2}_{m,3}\Big{(}c+\ln\Big{(}\frac{3T}{\mu_{3}}\Big{)}\Big{)}$
$\displaystyle-\frac{T^{2}}{(4\pi)^{2}}L_{b}(\mu)N\lambda_{m}\Big{(}\frac{1}{4}\lambda_{h}+\frac{5}{24}\lambda_{m}+\frac{(N+2)}{24}\lambda_{S}\Big{)},$
(A.16) $\displaystyle\lambda_{S,3}=$ $\displaystyle
T\bigg{[}\lambda_{S}(\mu)-\frac{1}{(4\pi)^{2}}L_{b}(\mu)\Big{(}\lambda_{m}^{2}+(N+8)\lambda_{\sigma}^{2}\Big{)}\bigg{]},$
(A.17) $\displaystyle\lambda_{m,3}=$ $\displaystyle
T\bigg{[}\lambda_{m}(\mu)+\frac{\lambda_{m}}{(4\pi)^{2}}\bigg{(}L_{b}(\mu)\Big{(}\frac{3}{4}(3g^{2}+g^{\prime
2})-6\lambda_{h}-2\lambda_{m}-(N+2)\lambda_{S}\Big{)}-3L_{f}(\mu)g_{\mbox{\tiny\rm{$Y$}}}^{2}\bigg{)}\bigg{]},$
(A.18) $\displaystyle\mu^{2}_{S,3}(\mu_{3})=$
$\displaystyle\mu_{S}^{2}(\mu)+T^{2}\Big{(}\frac{1}{6}\lambda_{m}(\mu)+\frac{1}{12}(N+2)\lambda_{S}(\mu)\Big{)}+\epsilon\;T^{2}\Big{(}\gamma_{\lambda_{m}2}\lambda_{m}+\gamma_{\lambda_{s}2}\lambda_{s}\Big{)}$
$\displaystyle-\frac{L_{b}(\mu)}{(4\pi)^{2}}\Big{(}2\lambda_{m}\mu_{h}^{2}+(N+2)\lambda_{S}\mu_{S}^{2}\Big{)}$
$\displaystyle+\frac{1}{(4\pi)^{2}}\Big{(}(3g^{2}_{3}+g^{\prime
2}_{3})\lambda_{m,3}-2\lambda_{m,3}^{2}-2(N+2)\lambda_{S,3}^{2}\Big{)}\Big{(}c+\ln\Big{(}\frac{3T}{\Lambda_{{\mbox{\tiny\rm{3d}}}}}\Big{)}\Big{)}$
$\displaystyle+\frac{T^{2}}{(4\pi)^{2}}\bigg{[}\frac{(2+3L_{b}(\mu))}{24}(3g^{2}+g^{\prime
2})\lambda_{m}-\frac{1}{4}(3L_{b}(\mu)-L_{f}(\mu))g_{\mbox{\tiny\rm{$Y$}}}^{2}\lambda_{m}$
$\displaystyle-
L_{b}(\mu)\bigg{(}\Big{(}\lambda_{h}+\frac{1}{12}(N+6)\lambda_{m}+\frac{1}{6}(N+2)\lambda_{S}\Big{)}\lambda_{m}+\frac{1}{12}(N+2)(N+8)\lambda_{S}^{2}\bigg{)}\bigg{]},$
(A.19)

where $\mu_{3}$ denotes the RG-scale of the 3d EFT. Here we have used the
shorthand notation

$\displaystyle c$
$\displaystyle\equiv\frac{1}{2}\left(\ln{\left(\frac{8\pi}{9}\right)}+\frac{\zeta^{\prime}(2)}{\zeta(2)}-2\gamma\right),$
(A.20) $\displaystyle L_{b}$ $\displaystyle\equiv 2\left(\ln{\frac{\mu}{4\pi
T}}+\gamma\right),\qquad\qquad\qquad L_{f}\equiv L_{b}+4\ln 2,$ (A.21)
$\displaystyle I_{b}$ $\displaystyle\equiv\bigg{(}1+\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\frac{\zeta^{\prime}(-1)}{\zeta(-1)}\bigg{)}.$ (A.22)

For the mass parameter we have included the
$\mathcal{O}(\epsilon\;T^{2}g^{2})$ level pieces, as these contribute to the
ultrasoft pressure $p_{\text{M2}}$ when multiplied with $1/\epsilon$ poles
therein. The coefficients $\beta_{\lambda_{m}2},\gamma_{\lambda_{m}2}$ and
$\gamma_{\lambda_{s}2}$ read

$\displaystyle\beta_{\lambda_{m}2}$
$\displaystyle=\frac{N}{12}I_{b},\qquad\gamma_{\lambda_{m}2}=\frac{1}{3}I_{b},\qquad\gamma_{\lambda_{s}2}=\frac{N+2}{6}I_{b}.$
(A.23)

Similar pure SM contributions to the Higgs mass parameter can be found in [62]
in Appendix B, there denoted by $\beta_{A2},\beta_{B2},\beta_{\lambda
2},\beta_{Y2}$. The 3d gauge parameters do not receive any singlet
contributions at NLO and can be read off from [58, 41]. The Debye masses for
the temporal scalars read

$\displaystyle m^{2}_{D}$
$\displaystyle=\Big{(}m^{2}_{D}\Big{)}_{\text{SM}}+\frac{N}{(4\pi)^{2}}\frac{T^{2}}{12}g^{2}\lambda_{m},$
(A.24) $\displaystyle{m}^{\prime 2}_{D}$ $\displaystyle=\Big{(}{m}^{\prime
2}_{D}\Big{)}_{\text{SM}}+\frac{N}{(4\pi)^{2}}\frac{T^{2}}{12}{g^{\prime}}^{2}\lambda_{m}.$
(A.25)

The SM parts at two-loop order can be found in Refs. [61, 62]. The singlet
interacts with the temporal scalar fields through [20]

$\displaystyle y_{3}$
$\displaystyle=T\frac{1}{(4\pi)^{2}}\frac{1}{2}g^{2}\lambda_{m},$ (A.26)
$\displaystyle y_{3}^{\prime}$
$\displaystyle=T\frac{1}{(4\pi)^{2}}\frac{1}{2}{g^{\prime}}^{2}\lambda_{m}.$
(A.27)

However, since these interactions only start at $\mathcal{O}(g^{4})$, they do
not contribute to any of our further expressions.

###  A.3 Coefficient of the unit operator, or pressure, in the symmetric phase

In addition to the matching of the couplings and masses, we also need the
matching of the unit operator (c.f. [59] Sec. 3A), which amounts to computing
vacuum contributions to the pressure from the hard and soft modes, as in [61].
Note that by convention, we include the similar ultrasoft pieces of the
symmetric phase in the effective potential.

#### Leading order pressure

At leading order, the pressure is given in terms of the sum-integrals [77]

$\displaystyle\mathcal{I}^{\prime}_{0}$
$\displaystyle\equiv-\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\ln(P^{2})=\frac{\pi^{2}}{45}T^{4}+\mathcal{O}(\epsilon),$
(A.28) $\displaystyle\widetilde{\mathcal{I}}^{\prime}_{0}$
$\displaystyle\equiv-\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{\\{P\\}}\ln(P^{2})=\frac{7}{8}\frac{\pi^{2}}{45}T^{4}+\mathcal{O}(\epsilon),$
(A.29)

arising from hard mode contributions of one-loop vacuum bubble diagrams. For
the definition of the sum-integration measure, see [77]. Following [61],
different fields contribute to the total pressure as

(fermions) $\displaystyle=-\Big{(}\underbrace{1}_{\text{lepton
singlets}}+\underbrace{2N_{c}}_{\text{quark singlets}}$
$\displaystyle+d_{F}(\underbrace{1}_{\text{lepton
doublets}}+\underbrace{N_{c}}_{\text{quark
doublets}})\Big{)}N_{f}\;\widetilde{\mathcal{I}}^{\prime}_{0},$ (A.30) (gauge
fields) $\displaystyle=\frac{1}{2}D\Big{(}\underbrace{d_{A}}_{\text{weak
bosons}}+\underbrace{1}_{\text{photon}}+\underbrace{(N^{2}_{c}-1)}_{\text{gluons}}\Big{)}\mathcal{I}^{\prime}_{0},$
(A.31) (ghosts) $\displaystyle=-\Big{(}\underbrace{d_{A}}_{\text{weak
bosons}}+\underbrace{1}_{\text{photon}}+\underbrace{(N^{2}_{c}-1)}_{\text{gluons}}\Big{)}\mathcal{I}^{\prime}_{0},$
(A.32) (Higgs) $\displaystyle=d_{F}\mathcal{I}^{\prime}_{0},$ (A.33) (singlet)
$\displaystyle=\frac{N}{2}\mathcal{I}^{\prime}_{0}.$ (A.34)

To label different contributions, we have used the following notation of [61]
for the constants: $N_{f}=3$ for the number of fermion families, $N_{c}=3$ for
the SU(3) colour group, $d_{F}=2$ for the dimensionality of the fundamental
representation of SU(2) and $d_{A}=3$ the dimensionality of the adjoint
representation of SU(2). In total,

$\displaystyle p^{\text{LO}}(T)$
$\displaystyle=\frac{\pi^{2}}{90}T^{4}\bigg{(}1+d_{A}+(N^{2}_{c}-1)+d_{F}+N+\frac{7}{8}N_{f}\Big{(}1+d_{F}+N_{c}(d_{F}+2)\Big{)}\bigg{)}$
$\displaystyle=(28+N+26.25N_{f})\frac{\pi^{2}}{90}T^{4}.$ (A.35)

This is the leading order pressure in the symmetric phase.

#### Higher order corrections

The full pressure in the symmetric phase has the form

$\displaystyle
p_{\text{sym}}(T)=p_{\text{E}}(T)+p_{\text{M1}}(T)+p_{\text{QCD}}(T)+\mathcal{O}(g^{5}T^{4}).$
(A.36)

Here, different contributions are organised in a manner similar to [61, 62]:
$p_{\text{E}}$ collects the hard mode contributions, and $p_{\text{M1}}$ the
soft mode contributions from temporal scalars. The subscripts denote electric
and magnetic contributions, as in hot EQCD [77]. We note that by convention,
the leading order QCD contributions are included in $p_{\text{E}}$, while
higher order QCD corrections are collected in $p_{\text{QCD}}$ [77, 75, 76,
78, 79]. We do not include these QCD corrections in our computation, as
consecutive orders are known to fluctuate around the ideal gas pressure – due
to the large values of the strong coupling $g_{s}$ and $g_{Y}$ – unless the
temperature is asymptotically large, and here our focus is solely on the EW
corrections.

In addition, for $p_{\text{sym}}$ there would also be an ultrasoft, or
additional magnetic contribution $p_{\text{M2}}(T)$ [61, 62] that arises from
one- and two-loop diagrams with ultrasoft scalars. However, by convention we
include these contributions in the effective potential part of the pressure,
i.e. $p_{\text{M2}}(T)=-TV^{\text{3d}}_{\text{eff}}(0,0)$ (c.f. Sec. A.4),
where both Higgs and singlet background fields vanish in the symmetric phase.

To reach the $\text{N}^{3}$LO, or $\mathcal{O}(g^{4})$ precision for the
vacuum diagrams requires a computation up to – and including – three-loop
topologies, for the hard scale contributions. These can be split into

$\displaystyle
p_{\text{E}}(T)=p^{\text{SM}}_{\text{E}}(T)+p^{\text{singlet}}_{\text{E}}(T),$
(A.37)

where the SM contribution can be read off from Eq. (20) and Appendix A in
[61].888 We note that Ref. [38] presents minor corrections to the results of
[61, 62], and for the parts relevant for our computation, we have confirmed
and used the results of [38].  In this work, we include the singlet
contributions

$\displaystyle p^{\text{singlet}}_{\text{E}}(T)$
$\displaystyle=T^{4}\bigg{[}\alpha^{\text{singlet}}_{E1}+\lambda_{m}(\mu)\alpha_{E\lambda_{m}}+\lambda_{S}(\mu)\alpha_{E\lambda_{S}}$
$\displaystyle+\frac{1}{(4\pi)^{2}}\bigg{(}g^{2}\lambda_{m}\alpha_{EA\lambda_{m}}+{g^{\prime}}^{2}\lambda_{m}\alpha_{EB\lambda_{m}}+g^{2}_{Y}\lambda_{m}\alpha_{EY\lambda_{m}}$
$\displaystyle+\lambda^{2}_{S}\alpha_{E\lambda_{S}\lambda_{S}}+\lambda^{2}_{m}\alpha_{E\lambda_{m}\lambda_{m}}+\lambda_{S}\lambda_{m}\alpha_{E\lambda_{S}\lambda_{m}}+\lambda_{h}\lambda_{m}\alpha_{E\lambda_{h}\lambda_{m}}\bigg{)}\bigg{]}$
$\displaystyle+T^{2}\bigg{[}\mu^{2}_{S}(\mu)\alpha_{E\mu^{2}_{S}}+\frac{1}{(4\pi)^{2}}\bigg{(}\lambda_{m}\mu^{2}_{S}\alpha_{E\lambda_{m}\mu^{2}_{S}}+\lambda_{m}\mu^{2}_{h}\alpha_{E\lambda_{m}\mu^{2}_{h}}+\lambda_{S}\mu^{2}_{S}\alpha_{E\lambda_{S}\mu^{2}_{S}}\bigg{)}\bigg{]}$
$\displaystyle+\frac{\mu^{4}_{S}}{(4\pi)^{2}}\alpha_{E\mu^{4}_{S}}+\mathcal{O}(g^{6}),$
(A.38)

where we used notation similar to [61, 62]. We list the results for all
coefficients $\alpha$ in Eqs. (A.57)-(A.71).

The soft contribution from the temporal scalars reads

$\displaystyle p_{\text{M1}}(T)=p^{\text{SM}}_{\text{M1}}(T),$ (A.39)

i.e. it is of exactly the same functional form as in the SM, see Eq. (12) in
[62]. Note, however, that the singlet contributes to the two-loop Debye masses
$m_{D}$ and $m_{D}^{\prime}$ as in Eqs. (A.24) and (A.25), and hence to
$p_{\text{M1}}$. Note that three-loop diagrams contribute to $p_{\text{M1}}$
at $\text{N}^{4}$LO, or $\mathcal{O}(g^{5})$, but we include them nonetheless,
since they provide the complete soft contribution at $\mathcal{O}(g^{5})$.
There are three-loop diagrams involving the singlet and temporal scalars, but
since their respective couplings are loop-induced, they are relatively
suppressed and do not contribute at $\mathcal{O}(g^{5})$.

The ultrasoft contribution to the pressure reads

$\displaystyle
p_{\text{M2}}(T)=p^{\text{SM}}_{\text{M2}}(T)+p^{\text{singlet}}_{\text{M2}}(T),$
(A.40)

where the SM piece can be read off from Eq. (15) in [62] and

$\displaystyle p^{\text{singlet}}_{\text{M2}}(T)$
$\displaystyle=T\bigg{(}\frac{N}{12\pi}\bar{\mu}^{3}_{s,3}-\frac{N}{(4\pi)^{2}}\Big{(}\frac{1}{4}(N+2)\bar{\lambda}_{s,3}\bar{\mu}^{2}_{s,3}+\bar{\lambda}_{m,3}\bar{\mu}_{s,3}\bar{\mu}_{h,3}\Big{)}\bigg{)}.$
(A.41)

As in [62], we compute $p_{\text{M2}}(T)$ at two-loop order, providing
$\mathcal{O}(g^{4})$ accuracy.

Finally, we comment that the divergent $1/\epsilon$ poles in $p_{\text{E}}$,
$p_{\text{M1}}$ and $p_{\text{M2}}$ cancel when all these terms are summed
together, i.e. the cancellation happens between contributions from different
scales.999 Here lies a technical subtlety. All $1/\epsilon$ poles indeed
cancel, but this cancellation happens exactly only when the ultrasoft pieces
are expanded to $\mathcal{O}(g^{4})$. As we have explained in Sec. 3, we do
not expand the ultrasoft pieces, but in practice we remove $1/\epsilon$ terms
manually.  Note that $p^{\text{singlet}}_{\text{M2}}$ is finite and
$p_{\text{M1}}$ does not have divergences related to the singlet
contributions. Therefore all leftover divergent singlet pieces from the hard
modes are cancelled by the ultrasoft Higgs-gauge field contribution that is
proportional to $\bar{\mu}^{2}_{h,3}$ and that encodes the contribution from
the singlet portal coupling.

#### Vacuum Feynman diagrams at the hard scale with the singlet

Figure 7:  Diagrams contributing to $p^{\text{singlet}}_{\text{E}}$. Dashed
(solid) lines denote the Higgs (singlet), wiggly line gauge fields and the
solid arrowed line the top quark. Vertices in these diagrams are bare
vertices, i.e. in renormalised perturbation theory there are counterterm
contributions in (a)-(e).

The computation of the singlet contributions depicted in Fig. 7 is relatively
straightforward – even at three-loop level – compared to the computation of
the SM contributions: gauge, ghost and fermion topology diagrams do not appear
at all at two-loop, and even at three-loop there are only three of them in
total. Furthermore, all these reduce trivially to lower order master
integrals, since the momentum in the singlet one-loop bubble separates from
the other momenta. All other three-loop topologies with the singlet are pure
scalar diagrams of the singlets and Higgs, and therefore straightforward to
compute. As in Appendix D of Ref. [61] for the SM contributions, we list the
explicit midstage results for all different new diagrams containing singlets.
Note that here we show also pure Higgs diagrams that include counterterms with
singlet couplings.

$\displaystyle(a)$ $\displaystyle=-d_{F}\delta\mu^{2}_{h}\mathcal{I}_{1}\;,$
(A.42) $\displaystyle(b)$
$\displaystyle=\frac{1}{2}N\Big{(}\mathcal{I}_{0}^{\prime}-(\mu^{2}_{S}+\delta\mu^{2}_{S})\mathcal{I}_{1}+\frac{1}{2}\mu^{4}_{S}\mathcal{I}_{2}\Big{)}\;,$
(A.43) $\displaystyle(c)$
$\displaystyle=-d_{F}(d_{F}+1)\delta\lambda_{h}\mathcal{I}^{2}_{1}\;,$ (A.44)
$\displaystyle(d)$
$\displaystyle=-\frac{1}{4}N(N+2)(\lambda_{S}+\delta\lambda_{S})\Big{(}\mathcal{I}^{2}_{1}-2\mu^{2}_{S}\mathcal{I}_{2}\mathcal{I}_{1}\Big{)}\;,$
(A.45) $\displaystyle(e)$
$\displaystyle=-\frac{1}{2}d_{F}N(\lambda_{m}+\delta\lambda_{m})\Big{(}\mathcal{I}^{2}_{1}-(\mu^{2}_{h}+\mu^{2}_{S})\mathcal{I}_{2}\mathcal{I}_{1}\Big{)}+\frac{1}{2}d_{F}N\lambda_{m}\delta
Z_{\phi}\mathcal{I}^{2}_{1}\;,$ (A.46) $\displaystyle(f)$
$\displaystyle=\frac{1}{4}N(N+2)^{2}\lambda^{2}_{S}\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.47) $\displaystyle(g)$
$\displaystyle=\frac{1}{2}d_{F}N(N+2)\lambda_{m}\lambda_{S}\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.48) $\displaystyle(h)$
$\displaystyle=\frac{1}{8}d_{F}N^{2}\lambda^{2}_{m}\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.49) $\displaystyle(i)$
$\displaystyle=Nd_{F}(d_{F}+1)\lambda_{h}\lambda_{m}\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.50) $\displaystyle(j)$
$\displaystyle=\frac{1}{4}Nd^{2}_{F}\lambda^{2}_{m}\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.51) $\displaystyle(k)$
$\displaystyle=\frac{1}{4}N(N+2)\lambda^{2}_{S}\mathcal{M}_{0,0}\;,$ (A.52)
$\displaystyle(l)$
$\displaystyle=\frac{1}{4}d_{F}N\lambda^{2}_{m}\mathcal{M}_{0,0}\;,$ (A.53)
$\displaystyle(m)$
$\displaystyle=\frac{1}{8}(d_{A}g^{2}+{g^{\prime}}^{2})d_{F}\lambda_{m}N(d+\xi)\mathcal{I}^{2}_{1}\mathcal{I}_{2}\;,$
(A.54) $\displaystyle(n)$
$\displaystyle=\;-\frac{1}{8}(d_{A}g^{2}+{g^{\prime}}^{2})d_{F}\lambda_{m}N\xi\mathcal{I}^{2}_{1}\mathcal{I}_{2},,$
(A.55) $\displaystyle(o)$
$\displaystyle=2\lambda_{m}NN_{c}g^{2}_{Y}\mathcal{I}_{1}\widetilde{\mathcal{I}}_{1}\mathcal{I}_{2}.\;.$
(A.56)

The results and definitions for the master integrals $\mathcal{I}_{1,2}$,
$\widetilde{\mathcal{I}}_{1}$ and $\mathcal{M}_{0,0}$ can be found in Appendix
D of Ref. [61] and Appendix A of Ref. [77]. For diagrams with gauge fields, we
have used the Fermi gauge, or general covariant gauge, with the identical
choice for the SU(2) and U(1) sector gauge fixing parameters
$\xi_{2}=\xi_{1}=\xi$. It can be seen how the gauge fixing parameter $\xi$
cancels for the sum of diagrams (m)$+$(n). Note that the SM parts in [61] were
computed in Feynman gauge, and the authors did not explicitly check the gauge
invariance of their final result for the pressure. We note that it would be a
valuable crosscheck of the correctness of the computation to confirm the gauge
invariance of the $p^{\text{SM}}_{E}(T)$, and for this the automated
computation developed and used in Refs. [94, 50, 20] could be utilised.

From our computation in Eqs. (A.42) - (A.56), we find the following results
for the different contributions from the hard modes to Eq. (A.38)

$\displaystyle\alpha^{\text{singlet}}_{E1}$
$\displaystyle=\frac{\pi^{2}}{90}N,$ (A.57)
$\displaystyle\alpha_{E\lambda_{m}}$ $\displaystyle=-\frac{d_{F}}{288}N,$
(A.58) $\displaystyle\alpha_{E\lambda_{S}}$
$\displaystyle=-\frac{1}{576}N(N+2),$ (A.59)
$\displaystyle\alpha_{EA\lambda_{m}}$
$\displaystyle=\frac{C_{F}d_{F}N}{144}\bigg{(}\frac{3}{\epsilon}+11+3\gamma+12\frac{\zeta^{\prime}(-1)}{\zeta(-1)}+15\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}\bigg{)},$ (A.60) $\displaystyle\alpha_{EB\lambda_{m}}$
$\displaystyle=\frac{d_{F}N}{144}\frac{1}{4}\bigg{(}\frac{3}{\epsilon}+11+3\gamma+12\frac{\zeta^{\prime}(-1)}{\zeta(-1)}+15\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}\bigg{)},$ (A.61) $\displaystyle\alpha_{EY\lambda_{m}}$
$\displaystyle=\frac{N_{c}N}{72}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}-\ln(2)+\gamma\bigg{)},$ (A.62)
$\displaystyle\alpha_{E\lambda_{S}\lambda_{S}}$
$\displaystyle=\frac{N(N+2)}{144}\bigg{(}\frac{31}{10}+6\frac{\zeta^{\prime}(-1)}{\zeta(-1)}-3\frac{\zeta^{\prime}(-3)}{\zeta(-3)}+\frac{1}{2}(N+2)\gamma+\frac{1}{2}(N+8)\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}\bigg{)},$ (A.63) $\displaystyle\alpha_{E\lambda_{m}\lambda_{m}}$
$\displaystyle=\frac{d_{F}N}{288}\bigg{(}\frac{31}{5}+12\frac{\zeta^{\prime}(-1)}{\zeta(-1)}-6\frac{\zeta^{\prime}(-3)}{\zeta(-3)}+\frac{1}{2}(N+2d_{F})\gamma+\frac{1}{2}(12+N+2d_{F})\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}\bigg{)},$ (A.64) $\displaystyle\alpha_{E\lambda_{S}\lambda_{m}}$
$\displaystyle=\frac{d_{F}N(N+2)}{144}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\gamma\bigg{)},$ (A.65)
$\displaystyle\alpha_{E\lambda_{h}\lambda_{m}}$
$\displaystyle=\frac{d_{F}(d_{F}+1)N}{72}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\gamma\bigg{)},$ (A.66) $\displaystyle\alpha_{E\mu^{2}_{S}}$
$\displaystyle=-\frac{1}{24}N,$ (A.67)
$\displaystyle\alpha_{E\lambda_{m}\mu^{2}_{S}}$
$\displaystyle=\frac{d_{F}N}{12}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\gamma\bigg{)},$ (A.68)
$\displaystyle\alpha_{E\lambda_{m}\mu^{2}_{h}}$
$\displaystyle=\frac{d_{F}N}{12}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\gamma\bigg{)},$ (A.69)
$\displaystyle\alpha_{E\lambda_{S}\mu^{2}_{S}}$
$\displaystyle=\frac{N(N+2)}{12}\bigg{(}\ln\Big{(}\frac{\mu}{4\pi
T}\Big{)}+\gamma\bigg{)},$ (A.70) $\displaystyle\alpha_{E\mu^{4}_{S}}$
$\displaystyle=-\frac{N}{2}\bigg{(}\ln\Big{(}\frac{\mu_{S}}{4\pi
T}\Big{)}-\frac{3}{4}+\gamma\bigg{)}.$ (A.71)

We have used similar notation to that of [61], and $L_{b}$ was defined in Eq.
(A.22). In analogy to [61], the normalisation $p(T=0)=0$ in the symmetric
phase is accounted for by $\alpha_{E\mu^{4}_{S}}$. In practice, we have
subracted the zero temperature Coleman-Weinberg part from the hard mode
contribution

$\displaystyle\alpha_{E\mu^{4}_{S}}=\frac{1}{(4\pi)^{2}}\frac{1}{4}N\mu^{4}_{S}L_{b}-\frac{1}{(4\pi)^{2}}\frac{1}{4}N\mu^{4}_{S}\bigg{(}\ln\Big{(}\frac{\mu^{2}_{S}}{\mu^{2}}\Big{)}-\frac{3}{2}\bigg{)}.$
(A.72)

Here, the first term is the contribution from the hard modes and the second
term is the $T=0$ contribution, c.f. e.g. Eqs. (A.22) and (2.6) in [60]. Note
in particular that the dependence on the RG-scale vanishes in
$\alpha_{E\mu^{4}_{S}}$. The divergence related to $\mu^{4}_{S}$ is removed by
imposing a vacuum counterterm, that reads, together with the similar
contribution of the Higgs

$\displaystyle\delta
V=-\frac{1}{(4\pi)^{2}}\frac{1}{4\epsilon}(N\mu^{4}_{S}+2\mu^{4}_{h}).$ (A.73)

This completes our computation of $p^{\text{singlet}}_{\text{E}}$.

###  A.4 The effective potential, or pressure, in the broken phase

The pressure in the broken phases, with non-zero vacuum expectation values for
the Higgs or singlet, can be decomposed as

$\displaystyle
p(T)=p_{\text{sym}}(T)-V_{\text{eff}}(v_{\text{min}},s_{\text{min}}),$ (A.74)

where the effective potential is evaluated at its minima
$(v_{\text{min}},s_{\text{min}})$. Note that according to our EFT
construction, the effective potential is computed in the final ultrasoft scale
EFT. The parameters, background fields and the ultrasoft RG scale
($\bar{\mu}_{3}$) of this EFT are denoted by bars.

For the xSM the effective potential has been computed in [43], in Landau
gauge, and here we generalise this computation for $N$ singlets used in this
work at hand. We compute the effective potential assuming only one singlet
component gets a vev at high temperature, i.e.
$S\rightarrow(S_{1}+\bar{x},...,S_{N})^{\text{T}}$. $S_{1}$ mixes with the
Higgs field, and they constitute two mass eigenstates $h_{1}$ and $h_{2}$ as
in the xSM case, where $N=1$. The mixing angle is defined as

$\displaystyle\text{tan}2\bar{\theta}\equiv\frac{4\bar{v}\bar{x}\bar{\lambda}_{m,3}}{2(\bar{\mu}^{2}_{S,3}-\bar{\mu}^{2}_{h,3})+(\bar{\lambda}_{m,3}-6\bar{\lambda}_{h,3})\bar{v}^{2}-(\bar{\lambda}_{m,3}-6\bar{\lambda}_{h,3})\bar{x}^{2}}.$
(A.75)

The other $N-1$ inert singlets do not couple to gauge fields, and they merely
have interaction with $h_{1}$, $h_{2}$, Goldstones and among themselves. All
inert singlets have mass (squared) eigenvalue

$\displaystyle\bar{m}^{2}_{S}$
$\displaystyle\equiv\bar{m}^{2}_{S_{i}}=\bar{\mu}^{2}_{S,3}+\bar{\lambda}_{S,3}\bar{x}^{2}+\frac{1}{2}\bar{\lambda}_{m,3}\bar{v}^{2},$
(A.76)

for $2<i\leq N$.

In total, we compose the effective potential as

$\displaystyle
V^{{\mbox{\tiny\rm{3d}}}}_{\text{eff}}(\bar{v},\bar{x})=V_{0}+V_{1}+V_{2},$
(A.77)

where the indices denote loop-order. The tree-level potential

$\displaystyle V_{0}=V^{\text{xSM}}_{0},$ (A.78)

is the same as in the xSM. At one-loop, the inert singlets add a contribution
to the effective potential

$\displaystyle V_{1}=V^{\text{xSM}}_{1}+(N-1)J_{3}(\bar{m}_{S}),$ (A.79)

where

$\displaystyle
J_{3}(m_{3})\equiv\frac{1}{2}\int_{p}\ln(p^{2}+m^{2}_{3})=-\frac{1}{12\pi}(m^{2}_{3})^{\frac{3}{2}}+\mathcal{O}(\epsilon).$
(A.80)

At two-loop, the inert singlets add

$\displaystyle
V_{2}=V^{\text{xSM}}_{2}-\bigg{(}(\text{SSS})+(\text{SS})\bigg{)},$ (A.81)

where

$\displaystyle(\text{SSS})$
$\displaystyle=(N-1)\bigg{(}\frac{1}{4}C^{2}_{s_{i}s_{i}h_{1}}\mathcal{D}_{SSS}(\bar{m}_{S},\bar{m}_{S},\bar{m}_{h,1})+\frac{1}{4}C^{2}_{s_{i}s_{i}h_{2}}\mathcal{D}_{SSS}(\bar{m}_{S},\bar{m}_{S},\bar{m}_{h,2})\bigg{)},$
(A.82) $\displaystyle(\text{SS})$
$\displaystyle=(N-1)\bigg{(}\frac{1}{8}C_{s_{i}s_{i}s_{i}s_{i}}\Big{(}I^{3}_{1}(\bar{m}_{S})\Big{)}^{2}+\frac{1}{2}C_{s_{i}s_{i}G^{+}G^{-}}I^{3}_{1}(\bar{m}_{S})I^{3}_{1}(\bar{m}_{G})$
$\displaystyle+\frac{1}{4}C_{s_{i}s_{i}GG}I^{3}_{1}(\bar{m}_{S})I^{3}_{1}(\bar{m}_{G})+\frac{1}{4}C_{s_{i}s_{i}h_{1}h_{1}}I^{3}_{1}(\bar{m}_{S})I^{3}_{1}(\bar{m}_{h,1})$
$\displaystyle+\frac{1}{4}C_{s_{i}s_{i}h_{2}h_{2}}I^{3}_{1}(\bar{m}_{S})I^{3}_{1}(\bar{m}_{h,2})\bigg{)}+\frac{1}{2}(N-2)(N-1)\bigg{(}\frac{1}{4}C_{s_{i}s_{i}s_{j}s_{j}}\Big{(}I^{3}_{1}(\bar{m}_{S})\Big{)}^{2}\bigg{)},$
(A.83)

where the master integrals can be found in the supplementary material of [88]
and the vertex coefficients read

$\displaystyle C_{s_{i}s_{i}h_{1}}$
$\displaystyle=-\bar{v}\bar{\lambda}_{m,3}c_{\bar{\theta}}+2\bar{\lambda}_{S,3}\bar{x}s_{\bar{\theta}},$
(A.84) $\displaystyle C_{s_{i}s_{i}h_{2}}$
$\displaystyle=-\bar{v}\bar{\lambda}_{m,3}s_{\bar{\theta}}-2\bar{\lambda}_{S,3}\bar{x}c_{\bar{\theta}},$
(A.85) $\displaystyle C_{s_{i}s_{i}s_{i}s_{i}}$
$\displaystyle=-6\bar{\lambda}_{S,3},$ (A.86) $\displaystyle
C_{s_{i}s_{i}G^{+}G^{-}}$
$\displaystyle=C_{s_{i}s_{i}GG}=-\bar{\lambda}_{m,3},$ (A.87) $\displaystyle
C_{h_{1}h_{1}s_{i}s_{i}}$
$\displaystyle=-\bar{\lambda}_{m,3}c_{\bar{\theta}}^{2}-2\bar{\lambda}_{S,3}s_{\bar{\theta}}^{2},$
(A.88) $\displaystyle C_{h_{2}h_{2}s_{i}s_{i}}$
$\displaystyle=-\bar{\lambda}_{m,3}s_{\bar{\theta}}^{2}-2\bar{\lambda}_{S,3}c_{\bar{\theta}}^{2},$
(A.89) $\displaystyle C_{s_{i}s_{i}s_{j}s_{j}}$
$\displaystyle=-2\bar{\lambda}_{S,3},\quad(i\neq j),$ (A.90)

where we denote $c_{\bar{\theta}}\equiv\cos\bar{\theta}$ and
$s_{\bar{\theta}}\equiv\sin\bar{\theta}$. The 3d EFT counterterms read

$\displaystyle\delta\bar{\mu}^{2}_{h,3}$
$\displaystyle=\delta\bar{\mu}^{2}_{h,3,\text{SM}}+\frac{1}{(4\pi)^{2}}\frac{1}{8\epsilon}N\bar{\lambda}^{2}_{m,3},$
(A.91) $\displaystyle\delta\bar{\mu}^{2}_{S,3}$
$\displaystyle=\frac{1}{(4\pi)^{2}}\frac{1}{4\epsilon}\Big{(}-(3\bar{g}^{2}_{3}+\bar{g}^{\prime
2}_{3})\bar{\lambda}_{m,3}+2\bar{\lambda}^{2}_{m,3}+2(N+2)\bar{\lambda}^{2}_{S,3}\Big{)},$
(A.92) $\displaystyle\delta V_{3}$
$\displaystyle=-\frac{1}{(4\pi)^{2}}\frac{1}{4\epsilon}(3\bar{g}^{2}_{3}+\bar{g}^{\prime
2}_{3})\bar{\mu}^{2}_{h,3}.$ (A.93)

Note that the vacuum counterterm $\delta V_{3}$ does not include contributions
from the singlet, but since it was not listed in [43], we add it here. The
magnetic, or ultrasoft, contribution to the symmetric phase pressure is
encoded in the effective potential as
$p_{\text{M2}}(T)=TV^{{\mbox{\tiny\rm{3d}}}}_{\text{eff}}(0,0)$.

###  A.5 Order parameters in perturbation theory and gauge invariance

Figure 8:  Schematic illustration of order parameters in perturbation theory,
as a function of temperature. Left: minimum of the effective potential $(v)$.
At high $T$, $v$ is zero, but non-zero at low $T$. If there is a
discontinuity, the transition between the two phases is of first order (A). If
$v$ is continuous, but there is a discontinuity in the first derivative, the
transition is of second order (B). A smooth crossover (C), where all order
derivatives are continuous, is not possible. Furthermore, $v$ itself is not a
gauge invariant quantity. Right: the scalar condensate
$\langle\phi^{\dagger}\phi\rangle$ has a non-vanishing value even in the high
$T$ phase, and hence a smooth interpolation between the two phases is
possible, leading to a crossover between the phases. Condensates can be
computed in a gauge invariant manner, even in perturbation theory.

Determining the character of a phase transition – whether it is of first or
second order, or crossover – is a non-trivial endeavour. In a crossover type
transition the system transitions smoothly from one phase to another, and
there is no discontinuity in the order parameter at any order in its
derivatives. In conventional analyses in which the minima of the potential are
treated as order parameters, such a transition can not be recovered: in the
high temperature symmetric phase such minima vanish, but they are non-zero and
vary smoothly at lower temperatures, so there cannot be a smooth interpolation
between the phases. Furthermore, these minima are notoriously gauge
dependent,101010 Only the value of the effective potential at the minima is
gauge invariant, and corresponds to the pressure.  and hence cannot serve as
realistic order-parameters. More realistic order parameters can be defined in
terms of condensates, such as $\langle\phi^{\dagger}\phi\rangle$ [95]. They
have a non-zero value even in the symmetric phase and

can be computed in a gauge-invariant manner in perturbation theory as

$\displaystyle\langle\phi^{\dagger}\phi\rangle=\frac{dV^{{\mbox{\tiny\rm{3d}}}}_{\text{eff}}}{d\mu^{2}_{3}},$
(A.94)

where the effective potential is evaluated in an expansion around its leading
order minimum. Even better, the condensates can be determined in a lattice
Monte Carlo simulation, which does not require gauge fixing. We point out,
however, that the computation of these condensates is dependent on their UV
renormalisation, and hence their direct numerical values do not have a
physical meaning.111111 Physical quantities such as the latent heat released
in the transition, can be related to these condensates though, see e.g. [50].
These order parameter-like quantities are schematically illustrated in Fig. 8.

In this article at hand, we are interested in the sound speed for a strong
phase transition, and we have focused solely on the second transition of a
two-stepper. In such a case, the barrier separating the phases exists already
at tree-level, and the phase transition can hence safely be assumed to be of
first order, even without a dedicated non-perturbative study. Regardless, the
computation in perturbation theory should be arranged in a gauge-invariant
manner, in order to study the physical thermodynamic properties of the plasma.
Despite this, however, in this article at hand we have resorted to a naive
analysis in terms of a direct minimisation of the real part of the potential,
and treat the location of the minima as a first estimate of the order
parameters, based on the arguments presented in [43]. With this decision, none
of our results are properly gauge invariant. Furthermore, in [96] it has been
demonstrated that the use of either minima in Landau gauge or condensates does
not make a qualitative difference. Hence, we expect our results in Landau
gauge to match the gauge invariant analysis, qualitatively. We argue, that
since our main interest is to get a handle on higher order corrections to the
speed of sound and estimate their effect on the GW power spectrum, the
limitation of a gauge dependent analysis does not compromise the conclusion of
this work. For future work, that would also include a determination of the
bubble nucleation rate and the nucleation temperature, and would hence provide
a more complete analysis of the predicted GW signature, we envision an upgrade
of our current computation of the pressure,121212 In our computation, the
dimensional reduction step is gauge invariant [20], as well as the symmetric
phase pressure [62, 61], but the use of the effective potential is not.  in
terms of a gauge invariant analysis in perturbation theory, c.f. recent [38,
60, 50, 96] but also [97, 98]. In addition, in such an analysis, a consistent
perturbative expansion around the leading order minima, imaginary parts do not
show up. The effective potential itself has a spurious imaginary part, that we
have simply disregarded in our direct minimisation of the potential.

###  A.6 Relation to lower order computations

For readers not familiar with the 3d EFT approach, we describe here the
commonly adopted one-loop computation, in order to provide a comparison. For
the remainder of this section, we set $N=1$, i.e. we discuss the xSM. In a
one-loop computation, the EFT technology is not as apparent as in computations
required for higher order corrections. Regardless, even the one-loop effective
potential based on daisy resummation by Arnold and Espinosa [99] utilises the
EFT picture for the thermal scale hierarchy, as only the soft modes are
resummed, corresponding to screening of the hard scale. However, the
underlying EFT picture is not apparent in [99] and the computation is
performed directly in the 4d parent theory; a formal discussion in the EFT
language is formulated in [58, 59]. At leading order, the parameters of the
EFT have trivial relations to the 4d parameters, i.e. $\lambda_{3}\sim
T\lambda$ and $v^{2}_{3}\sim v^{2}/T$, and only the masses receive thermal
corrections, at one-loop level. Let us start by taking a closer look at the
all-order resummation of the leading daisy diagrams.

Figure 9:  Schematic illustration of leading daisy resummation. (a): division
of the one-loop vacuum bubble into soft (dashed line) and hard contributions
(double lines). (b): all-order resummation of the hard mode contributions to
the soft term. (c): resummation of the soft scalar mass parameter.

#### Leading order daisy resummation

In Fig. 9 we illustrate schematically the premise of daisy resummation. Fig. 9
(a) showcases the division of the one-loop vacuum bubble into soft (dashed
line) and hard contributions (double lines). While the hard modes are
regulated in the IR by non-zero Matsubara frequencies, the soft part is IR-
sensitive. Physically, this is related to macroscopic collective phenomena,
i.e. screening of different scales in the plasma. Technically, a problem
appears in loops including the zero-mode, since the propagator has the same
form as at $T=0$, while the integration measure has one less power of
momentum: this leads to worse behaviour in the IR compared to the $T=0$ case.
This is problem is demonstrated in Fig. 9 (b): consider a $k$-loop daisy
diagram where the inner loop has a soft loop-momentum with $k$ propagators
with unresummed mass $m^{2}\sim(gT)^{2}$, and $k$ loops with hard loop momenta
$P$. When each vertex contributes with $g^{2}$ this diagram has the form [100]

$\displaystyle\sim\bigg{[}T\int_{p}\frac{1}{(p^{2}+m^{2})^{k}}\bigg{]}_{\text{soft}}\bigg{[}g^{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}^{\prime}\frac{1}{P^{2}}\bigg{]}^{k}_{\text{hard}}$
$\displaystyle=\bigg{[}T\frac{1}{(-1)_{k-1}}\frac{d^{k-1}}{d(m^{2})^{k-1}}\int_{p}\frac{1}{(p^{2}+m^{2})}\bigg{]}_{\text{soft}}\bigg{[}\frac{1}{12^{k}}g^{2k}T^{2k}\bigg{]}_{\text{hard}}$
$\displaystyle=\bigg{[}-\frac{T}{4\pi}\frac{(\frac{1}{2})_{k-1}}{(-1)_{k-1}}(m^{2})^{\frac{3}{2}-k}\bigg{]}_{\text{soft}}\bigg{[}\frac{1}{12^{k}}g^{2k}T^{2k}\bigg{]}_{\text{hard}}$
$\displaystyle=\bigg{(}-\frac{1}{4\pi\times
12^{k}}\frac{(\frac{1}{2})_{k-1}}{(-1)_{k-1}}\bigg{)}\bigg{(}m^{3}T\Big{(}\frac{gT}{m}\Big{)}^{2k}\bigg{)}\sim
g^{3}\quad\text{for all}\;k.$ (A.95)

Results for the used one-loop integrals can be found in e.g. [77]. Here
$(x)_{n}$ denotes the falling factorial, but the overall numerical constant is
irrelevant for the observation, that regardless of $k$, any such diagram
contributes at $\mathcal{O}(g^{3})$, since $m\sim gT$. Even worse, for $k\geq
2$ and $m\rightarrow 0$ all these contributions are IR-divergent. However, for
massive fields – i.e. scalars and temporal scalar fields, but not for
massless, spatial gauge fields – there is a salvation by means of a
resummation. In Fig. 9 (c) the soft scalar propagator is resummed as a
geometric Dyson series. Formally, the scalar 2-point Green’s function $G$ can
be witten in terms of the tree-level propagator $G_{0}=1/(p^{2}+m^{2})$ and
the one-loop thermal correction $\Pi\sim g^{2}T^{2}$ as

$\displaystyle G$ $\displaystyle=G_{0}+G_{0}\Pi G_{0}+G_{0}(\Pi
G_{0}\Pi)G_{0}+\ldots=G_{0}\Big{(}\sum_{i=0}^{\infty}(\Pi G_{0})^{i}\Big{)}$
$\displaystyle=G_{0}\Big{(}\frac{1}{1-\Pi
G_{0}}\Big{)}=\frac{1}{p^{2}+m^{2}+\Pi},$ (A.96)

which leads to a thermally corrected mass for the three-dimensional zero-mode
$m^{2}_{3}\equiv m^{2}+\Pi$, which is resummed to all-orders, by the one-loop
hard correction. When all different daisy diagrams in Fig. 9 (b) are summed
together, the result for the one-loop soft vacuum bubble reads

$\displaystyle(\text{LO soft})=-\frac{m^{3}_{3}}{4\pi}\sim\mathcal{O}(g^{3}).$
(A.97)

When this contribution is combined with the tree-level terms and thermally
corrected mass at $\mathcal{O}(g^{2})$, the effective potential is consistent
at $\mathcal{O}(g^{3})$.

Figure 10:  Illustration of higher order resummations in NLO dimensional
reduction, going beyond leading one-loop daisy resummation of the mass. (a):
two-loop resummation of the mass parameter. (b): one-loop resummation of the
couplings. (c): one-loop resummation of the fields, from the part dependent on
the external momentum. Here a dashed line represents a soft field, and a
double line and wiggly line hard scalar and gauge field modes, respectively.
All NLO effects in dimensional reduction contibute at $\mathcal{O}(g^{4})$.

In Fig. 9 (a) in addition to the $\mathcal{O}(g^{3})$ soft contribution, there
are one-loop hard mode contributions – that include the $T=0$ Coleman-Weinberg
potential (see end of this section) – which appear at order
$\mathcal{O}(g^{4})$. Alas, for a complete computation at this order, also
higher order resummations are required. Such higher order resummations are
those depicted in Fig. 10: (a) two-loop corrections to the mass, (b) one-loop
corrections to couplings and (c) one-loop resummation of the fields, from
momentum dependent contributions to the 2-point function. The inclusion of all
these contributions is systematic in the EFT language, and they are resummed
to all-orders in dimensional reduction to 3d EFT, in analogy to daisy
resummation that resums one-loop thermal mass to all orders.

#### Comparison of 3d and 4d effective potentials

Schematically, the relation between the lower-order 4d effective potential and
the higher-order 3d effective potential is131313Here the labels 4d and 3d
merely denote the approach used to compute the free energy, that describes
physical properties of the plasma.

$\displaystyle V^{\text{4d}}_{\text{eff}}(v,s)\simeq
T\;V^{\text{3d}}_{\text{eff}}(v_{3},s_{3}).$ (A.98)

This relation is only an approximate equality, since
$V^{\text{3d}}_{\text{eff}}(v_{3},s_{3})$ contains higher-order contributions.
Technically, the one-loop $V^{\text{4d}}_{\text{eff}}$ with daisy resummation
is correct at $\mathcal{O}(g^{3})$, but incomplete at $\mathcal{O}(g^{4})$. By
comparing the above expressions, we have for the soft terms

$\displaystyle V^{\text{4d}}_{\text{soft}}=T\;V^{\text{3d}}_{\text{1-loop}},$
(A.99)

which holds exactly if only the LO matching relations of Eqs. (A.121) and
(A.122) are used, at $\mathcal{O}(g^{2})$. In practice one uses
$\mathcal{O}(g^{4})$ accurate relations in the 3d EFT, providing higher
accuracy. We emphasise, that this is exactly the term that corresponds to
leading daisy resummation. For the hard terms we have

$\displaystyle V^{\text{4d}}_{\text{tree}}+V^{\text{4d}}_{\text{hard}}\simeq
T\;V^{\text{3d}}_{\text{tree}}.$ (A.100)

Again, equality holds exactly if the expressions are computed to the same
order in the coupling expansion. In the 3d picture, hard contributions are
encoded in the parameters of the EFT. Here we emphasise, that in the 4d
approach, the leading contributions that depend on the external momenta –
corresponding to a screening of the soft field by the heat bath – are not
included. In the 3d computation these terms are included, and necessary for a
consistent computation at $\mathcal{O}(g^{4})$ and RG-improvement, see Ref.
[99, 101, 60]. As discussed in these references, at $\mathcal{O}(g^{4})$ also
two-loop contributions to the thermal masses are required, as well as a two-
loop effective potential for the soft or 3d terms.

For concreteness, below we construct again the one-loop effective potential
using the background field method. We need the mass squared eigenvalues for
fluctuating quantum fields

$\displaystyle m^{2}_{\chi}$
$\displaystyle=\mu^{2}_{h}+\lambda_{h}v^{2}+\frac{1}{2}\lambda_{m}s^{2},$
(A.101) $\displaystyle m^{2}_{\pm}$
$\displaystyle=\frac{1}{4}\bigg{\\{}2\mu^{2}_{h}+2\mu^{2}_{S}+(6\lambda_{h}+\lambda_{m})v^{2}+(\lambda_{m}+6\lambda_{S})s^{2}\pm\sqrt{A+B}\bigg{\\}},$
(A.102) $\displaystyle A$
$\displaystyle\equiv(-6\lambda_{h}+\lambda_{m})^{2}v^{4}+\Big{(}2\mu^{2}_{h}-2\mu^{2}_{S}+(\lambda_{m}-6\lambda_{S})s^{2}\Big{)}^{2},$
(A.103) $\displaystyle B$ $\displaystyle\equiv
2v^{2}\bigg{(}2(6\lambda_{h}-\lambda_{m})(\mu^{2}_{h}-\mu^{2}_{S})+\Big{(}6\lambda_{h}(\lambda_{m}-6\lambda_{S})+\lambda_{m}(7\lambda_{m}+6\lambda_{S})\Big{)}s^{2}\bigg{)},$
(A.104)

where $m_{\pm}$ are the two neutral scalar mass eigenvalues, and the Goldstone
mass eigenvalue $m^{2}_{\chi}$ is triply-degenerate. For the gauge field and
top quark pieces, we need the mass eigenvalues

$\displaystyle m^{2}_{W}=\frac{1}{4}g^{2}v^{2},\quad\quad
m^{2}_{Z}=\frac{1}{4}(g^{2}+{g^{\prime}}^{2})v^{2},\quad\quad
m^{2}_{t}=\frac{1}{2}g^{2}_{Y}v^{2}.$ (A.105)

The mass squared eigenvalues related to the temporal scalar fields read [101,
58]

$\displaystyle M^{2}_{3}$ $\displaystyle=m^{2}_{D}+h_{3}v^{2}_{3},$ (A.106)
$\displaystyle M^{2}_{\pm,3}$
$\displaystyle=\frac{1}{2}\bigg{\\{}m^{2}_{D}+m_{D}^{\prime
2}+(h_{3}+h_{3}^{\prime})v^{2}_{3}+\pm\sqrt{X+Y}\bigg{\\}},$ (A.107)
$\displaystyle X$ $\displaystyle\equiv(m_{D}^{2}-m_{D}^{\prime
2})^{2}-2(h_{3}-h_{3}^{\prime})(m_{D}^{\prime 2}-m_{D}^{2})v^{2}_{3},$ (A.108)
$\displaystyle Y$
$\displaystyle\equiv\Big{(}(h_{3}-h_{3}^{\prime})^{2}+h^{\prime\prime
2}_{3}\Big{)}v^{4}_{3},$ (A.109)

where $M^{2}_{3}$ is doubly-degenerate. Here all quantities are those of the
3d EFT, since temporal scalars live at the soft scale.141414 Temporal scalars
$A^{a}_{0}$ and $B_{0}$ are mixing through a $h_{3}^{\prime\prime}$ term, but
in practise one can linearise this mixing.  Coefficients
$h_{3},h_{3}^{\prime},h_{3}^{\prime\prime}$ are couplings between Higgs and
temporal scalars, see [58, 20].

Above, the mass eigenvalues are input for the one-loop master integral

$\displaystyle V^{\text{4d}}_{\text{1-loop}}$
$\displaystyle\sim\frac{1}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\ln(P^{2}+m^{2})=J_{\text{soft}}(m)+J_{\text{hard}}(m).$
(A.110)

We divide it into the soft zero-mode contribution

$\displaystyle J_{\text{soft}}(m_{3})$ $\displaystyle=TJ_{3}(m_{3}),$ (A.111)

with all quantities resummed, i.e. those of the 3d EFT, and denoted herein
with subscript 3, and contributions from the hard non-zero modes

$\displaystyle J_{\text{hard}}(m)$
$\displaystyle\equiv\frac{1}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}^{\prime}\ln(P^{2}+m^{2})=\frac{1}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}^{\prime}\ln
P^{2}+\frac{1}{2}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n}(m^{2})^{n}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}^{\prime}\frac{1}{(P^{2})^{n}},$
(A.112)

where we have used the high-$T$ expansion in $m^{2}/T^{2}$ for the last line.
Concretely

$\displaystyle J^{b}_{\text{hard}}$
$\displaystyle=-\frac{\pi^{2}}{90}T^{4}+\frac{1}{24}T^{2}m^{2}-\frac{m^{4}}{4(4\pi)^{2}}\Big{(}\frac{1}{\epsilon}+L_{b}\Big{)}+\mathcal{O}\Big{(}\frac{m^{6}}{T^{2}}\Big{)},$
(A.113) $\displaystyle J^{f}_{\text{hard}}$
$\displaystyle=-\frac{7}{8}\frac{\pi^{2}}{90}T^{4}-\frac{1}{48}T^{2}m^{2}-\frac{m^{4}}{4(4\pi)^{2}}\Big{(}\frac{1}{\epsilon}+L_{f}\Big{)}+\mathcal{O}\Big{(}\frac{m^{6}}{T^{2}}\Big{)}\;,$
(A.114)

where we give the result for both bosonic and fermionic fields. Typically only
the first three terms are kept. Higher order terms would result in higher
dimensional, marginal, operators in the 3d EFT [58].

In total, the one-loop effective potential reads (in Landau gauge)

$\displaystyle
V^{\text{4d}}_{\text{eff}}=V^{\text{4d}}_{\text{tree}}+V^{\text{4d}}_{\text{CT}}+V^{\text{4d}}_{\text{1-loop}},$
(A.115)

where

$\displaystyle V^{\text{4d}}_{\text{tree}}$
$\displaystyle=\frac{1}{2}\mu^{2}_{h}v^{2}+\frac{1}{4}\lambda_{h}v^{4}+\frac{1}{2}\mu^{2}_{S}s^{2}+\frac{1}{4}\lambda_{S}s^{4}+\frac{1}{4}\lambda_{m}v^{2}s^{2},$
(A.116) $\displaystyle V^{\text{4d}}_{\text{CT}}$
$\displaystyle=\frac{1}{2}\delta\mu^{2}_{h}v^{2}+\frac{1}{4}\delta\lambda_{h}v^{4}+\frac{1}{2}\delta\mu^{2}_{S}s^{2}+\frac{1}{4}\delta\lambda_{S}s^{4}+\frac{1}{4}\delta\lambda_{m}v^{2}s^{2},$
(A.117)

and

$\displaystyle
V^{\text{4d}}_{\text{1-loop}}=V^{\text{4d}}_{\text{soft}}+V^{\text{4d}}_{\text{hard}},$
(A.118)

with ($d=3-2\epsilon$)

$\displaystyle V^{\text{4d}}_{\text{hard}}$
$\displaystyle=3J^{b}_{\text{hard}}(m_{\chi})+J^{b}_{\text{hard}}(m_{+})+J^{b}_{\text{hard}}(m_{-})$
$\displaystyle+d\Big{(}2J^{b}_{\text{hard}}(m_{W})+J^{b}_{\text{hard}}(m_{Z})\Big{)}-4N_{c}J^{f}_{\text{hard}}(m_{t}),$
(A.119) $\displaystyle V^{\text{4d}}_{\text{soft}}$
$\displaystyle=3J_{\text{soft}}(m_{\chi,3})+J_{\text{soft}}(m_{+,3})+J_{\text{soft}}(m_{-,3})$
$\displaystyle+(d-1)\Big{(}J_{\text{soft}}(m_{W,3})+J_{\text{soft}}(m_{Z,3})\Big{)}$
$\displaystyle+2J_{\text{soft}}(M_{3})+J_{\text{soft}}(M_{+,3})+J_{\text{soft}}(M_{-,3}).$
(A.120)

The counterterms cancel the $1/\epsilon$ poles in $J_{\text{hard}}$. In the
soft pieces, the conventional daisy-resummed approach corresponds to utilising
the dimensional reduction matching relations at leading order, i.e. one-loop
for the masses and tree-level for the couplings

$\displaystyle\mu^{2}_{h,3}$
$\displaystyle=\mu^{2}_{h}+\Pi_{\phi},\quad\quad\mu^{2}_{S,3}=\mu^{2}_{S}+\Pi_{S},$
(A.121) $\displaystyle\lambda_{h,3}$
$\displaystyle=T\lambda_{h},\quad\quad\lambda_{S,3}=T\lambda_{S}\quad\quad\lambda_{m,3}=T\lambda_{m}.$
(A.122)

The background fields are related as $v_{3}=v/\sqrt{T}$ and
$s_{3}=s/\sqrt{T}$. The one-loop thermal mass contributions read

$\displaystyle\Pi_{\phi}$
$\displaystyle=\frac{T^{2}}{12}\Big{(}\frac{3}{4}(3g^{2}+{g^{\prime}}^{2})+3g^{2}_{Y}+6\lambda_{h}+\frac{1}{2}\lambda_{m}\Big{)},$
(A.123) $\displaystyle\Pi_{S}$
$\displaystyle=\frac{T^{2}}{12}\Big{(}2\lambda_{m}+3\lambda_{S}\Big{)},$
(A.124) $\displaystyle m^{2}_{D}$
$\displaystyle=g^{2}T^{2}\Big{(}\frac{5}{6}+\frac{1}{3}N_{f}\Big{)},$ (A.125)
$\displaystyle m_{D}^{\prime 2}$
$\displaystyle={g^{\prime}}^{2}T^{2}\Big{(}\frac{1}{6}+\frac{5}{3}N_{f}\Big{)}.$
(A.126)

This completes the construction of the one-loop effective potential. Note that
the zero temperature Coleman-Weinberg potential is implicitly included above,
as

$\displaystyle V^{\text{4d}}_{\text{1-loop}}$
$\displaystyle\sim\frac{1}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\ln(P^{2}+m^{2})=J_{\text{soft}}(m)+J_{\text{hard}}(m)=J_{\text{CW}}(m)+J_{T}(m),$
(A.127)

where

$\displaystyle J_{\text{CW}}(m)$
$\displaystyle=\frac{1}{2}\Big{(}\frac{\mu^{2}e^{\gamma}}{4\pi}\Big{)}^{\epsilon}\int\frac{d^{D}p}{(2\pi)^{D}}\ln(p^{2}+m^{2})=-\frac{1}{2}\Big{(}\frac{\mu^{2}e^{\gamma}}{4\pi}\Big{)}^{\epsilon}\frac{\Big{(}m^{2}\Big{)}^{\frac{D}{2}}}{(4\pi)^{\frac{D}{2}}}\frac{\Gamma(-\frac{D}{2})}{\Gamma(1)}\;,$
(A.128) $\displaystyle J_{T}(m)$ $\displaystyle=-\int_{p}T\ln\Big{(}1\pm
n_{{\mbox{\tiny\rm{B/F}}}}(E_{p},T)\Big{)}.$ (A.129)

Here $J_{\text{CW}}$ is the result at $T=0$ in $D=d+1=4-2\epsilon$ dimensions
and all temperature dependence is captured in the thermal function
$J_{T}$151515This separation can be obtained by a contour trick [102, 103],
where the sum over Matsubara frequencies is turned into a contour integration
in the complex plane. in which $n_{{\mbox{\tiny\rm{B/F}}}}$ is the Bose or
Fermi distribution. Note, however, that in these expressions the daisy
resummation is not yet implemented. To re-express Eq. (A.127) into the more
popular form with daisy resummation, we write

$\displaystyle J_{\text{CW}}(m)+J_{T}(m)$
$\displaystyle=J_{\text{CW}}(m)+\underbrace{\Big{(}J^{\prime}_{T}(m)+J_{\text{soft}}(m)\Big{)}}_{J_{T}}$
$\displaystyle\rightarrow
J_{\text{CW}}(m)+J^{\prime}_{T}(m)+J_{\text{soft}}(m_{3})$
$\displaystyle=J_{\text{CW}}(m)+J_{T}(m)+\underbrace{\Big{(}J_{\text{soft}}(m_{3})-J_{\text{soft}}(m)\Big{)}}_{\equiv
J_{\text{daisy}}}$
$\displaystyle=J_{\text{CW}}(m)+J_{T}(m)+J_{\text{daisy}}(m_{3},m),$ (A.130)

a form often encountered in the literature [90]. Here the prime denotes the
sole hard contributions, and the resummation is implemented by
$J_{\text{soft}}(m)\rightarrow J_{\text{soft}}(m_{3})$, i.e. resumming the
soft mode. The arguably counterintuitive form $J_{\text{daisy}}(m_{3},m)\equiv
J_{\text{soft}}(m_{3})-J_{\text{soft}}(m)$ simply results from the fact that
the unresummed zero mode needs to be subtracted, as it is still present in
$J_{T}(m)$. Despite this counterintuitive formulation, by nature this
resummation is still equal to LO dimensional reduction and use of the 3d EFT.
The unresummed thermal function $J_{T}(m)$ can be evaluated numerically
without using the high-$T$ expansion, while the resummation of the zero mode
is encoded in $J_{\text{daisy}}(m_{3},m)$. However, the zero mode is resummed
due to screening by hard non-zero modes, and still assumes the high
temperature scale hierarchy. Therefore, in this formulation it is not
straightforward to move to a low-temperature limit even if the thermal
function is handled numerically, without high temperature expansion.

For completeness, we list high-$T$ expansions of bosonic and fermionic thermal
functions in a form often encountered in the literature

$\displaystyle
J_{T,b}(m)=T^{4}\bigg{(}-\frac{\pi^{2}}{90}+\frac{1}{24}\frac{m^{2}}{T^{2}}-\frac{1}{4(4\pi)^{2}}\frac{m^{4}}{T^{4}}\ln\Big{(}\frac{m^{2}}{a_{b}T^{2}}\Big{)}+\mathcal{O}\Big{(}\frac{m^{6}}{T^{6}}\Big{)}\bigg{)},$
(A.131) $\displaystyle
J_{T,f}(m)=T^{4}\bigg{(}-\frac{7}{8}\frac{\pi^{2}}{90}-\frac{1}{48}\frac{m^{2}}{T^{2}}-\frac{1}{4(4\pi)^{2}}\frac{m^{4}}{T^{4}}\ln\Big{(}\frac{m^{2}}{a_{f}T^{2}}\Big{)}+\mathcal{O}\Big{(}\frac{m^{6}}{T^{6}}\Big{)}\bigg{)},$
(A.132)

where $a_{b}=(4\pi)^{2}\text{Exp}(\frac{3}{2}-2\gamma)$ and $a_{f}=a_{b}/16$.
In the high-$T$ expansion it is straightforward to check that indeed
$J_{\text{CW}}+J_{T}^{\prime}=J_{\text{hard}}$, i.e. the zero-$T$ Coleman-
Weinberg potential is included in $J_{\text{hard}}$.

Finally, we comment that above we did not yet include $\\#T^{4}$ contributions
from gluons and light fermions in the effective potential, but this can be
done as in Sec. A.3 for the leading order pressure.

## References

  * [1] C. J. Hogan, _Gravitational radiation from cosmological phase transitions_ , _Mon. Not. Roy. Astron. Soc._ 218 (1986) 629\.
  * [2] 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].
  * [3] C. Caprini et al., _Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions_ , _JCAP_ 04 (2016) 001 [1512.06239].
  * [4] 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].
  * [5] 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].
  * [6] C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024 [1910.13125].
  * [7] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender et al., _Laser Interferometer Space Antenna_ , _arXiv e-prints_ (2017) arXiv:1702.00786 [1702.00786].
  * [8] K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, _Is there a hot electroweak phase transition at $m_{H}\gtrsim m_{W}$?_, _Phys. Rev. Lett._ 77 (1996) 2887 [hep-ph/9605288].
  * [9] K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, _The Electroweak phase transition: A Nonperturbative analysis_ , _Nucl. Phys. B_ 466 (1996) 189 [hep-lat/9510020].
  * [10] F. Csikor, Z. Fodor and J. Heitger, _Endpoint of the hot electroweak phase transition_ , _Phys. Rev. Lett._ 82 (1999) 21 [hep-ph/9809291].
  * [11] M. A. Stephanov, _QCD phase diagram: An Overview_ , _PoS_ LAT2006 (2006) 024 [hep-lat/0701002].
  * [12] 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].
  * [13] 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].
  * [14] M. J. Ramsey-Musolf, _The electroweak phase transition: a collider target_ , _JHEP_ 09 (2020) 179 [1912.07189].
  * [15] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, _On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe_ , _Phys. Lett. B_ 155 (1985) 36.
  * [16] D. E. Morrissey and M. J. Ramsey-Musolf, _Electroweak baryogenesis_ , _New J. Phys._ 14 (2012) 125003 [1206.2942].
  * [17] T. Konstandin, _Quantum Transport and Electroweak Baryogenesis_ , _Phys. Usp._ 56 (2013) 747 [1302.6713].
  * [18] G. A. White, _A Pedagogical Introduction to Electroweak Baryogenesis_ , .
  * [19] D. Bodeker and W. Buchmuller, _Baryogenesis from the weak scale to the grand unification scale_ , _Rev. Mod. Phys._ 93 (2021) 035004 [2009.07294].
  * [20] P. M. Schicho, T. V. I. Tenkanen and J. Österman, _Robust approach to thermal resummation: Standard Model meets a singlet_ , _JHEP_ 06 (2021) 130 [2102.11145].
  * [21] R. Jinno and M. Takimoto, _Gravitational waves from bubble dynamics: Beyond the Envelope_ , _JCAP_ 01 (2019) 060 [1707.03111].
  * [22] T. Konstandin, _Gravitational radiation from a bulk flow model_ , _JCAP_ 03 (2018) 047 [1712.06869].
  * [23] 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].
  * [24] 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].
  * [25] 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].
  * [26] P. Niksa, M. Schlederer and G. Sigl, _Gravitational Waves produced by Compressible MHD Turbulence from Cosmological Phase Transitions_ , _Class. Quant. Grav._ 35 (2018) 144001 [1803.02271].
  * [27] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili and A. Kosowsky, _Numerical simulations of gravitational waves from early-universe turbulence_ , _Phys. Rev. D_ 102 (2020) 083512 [1903.08585].
  * [28] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, _Circular polarization of gravitational waves from early-Universe helical turbulence_ , _Phys. Rev. Res._ 3 (2021) 013193 [2011.05556].
  * [29] P. Auclair, C. Caprini, D. Cutting, M. Hindmarsh, K. Rummukainen, D. A. Steer et al., _Generation of gravitational waves from freely decaying turbulence_ , 2205.02588.
  * [30] 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].
  * [31] 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].
  * [32] F. Giese, T. Konstandin and J. van de Vis, _Model-independent energy budget of cosmological first-order phase transitions—A sound argument to go beyond the bag model_ , _JCAP_ 07 (2020) 057 [2004.06995].
  * [33] F. Giese, T. Konstandin, K. Schmitz and J. Van De Vis, _Model-independent energy budget for LISA_ , _JCAP_ 01 (2021) 072 [2010.09744].
  * [34] X. Wang, F. P. Huang and X. Zhang, _Energy budget and the gravitational wave spectra beyond the bag model_ , _Phys. Rev. D_ 103 (2021) 103520 [2010.13770].
  * [35] X. Wang, F. P. Huang and Y. Li, _Sound velocity effects on the phase transition gravitational wave spectrum in the Sound Shell Model_ , 2112.14650.
  * [36] M. Laine and Y. Schroder, _Quark mass thresholds in QCD thermodynamics_ , _Phys. Rev. D_ 73 (2006) 085009 [hep-ph/0603048].
  * [37] HotQCD collaboration, _Equation of state in ( 2+1 )-flavor QCD_ , _Phys. Rev. D_ 90 (2014) 094503 [1407.6387].
  * [38] M. Laine and M. Meyer, _Standard Model thermodynamics across the electroweak crossover_ , _JCAP_ 07 (2015) 035 [1503.04935].
  * [39] D. Croon, V. Sanz and G. White, _Model Discrimination in Gravitational Wave spectra from Dark Phase Transitions_ , _JHEP_ 08 (2018) 203 [1806.02332].
  * [40] A. Drozd, B. Grzadkowski and J. Wudka, _Multi-Scalar-Singlet Extension of the Standard Model - the Case for Dark Matter and an Invisible Higgs Boson_ , _JHEP_ 04 (2012) 006 [1112.2582].
  * [41] T. Brauner, T. V. I. Tenkanen, A. Tranberg, A. Vuorinen and D. J. Weir, _Dimensional reduction of the Standard Model coupled to a new singlet scalar field_ , _JHEP_ 03 (2017) 007 [1609.06230].
  * [42] O. Gould, _Real scalar phase transitions: a nonperturbative analysis_ , _JHEP_ 04 (2021) 057 [2101.05528].
  * [43] L. Niemi, P. Schicho and T. V. I. Tenkanen, _Singlet-assisted electroweak phase transition at two loops_ , _Phys. Rev. D_ 103 (2021) 115035 [2103.07467].
  * [44] 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].
  * [45] U.-L. Pen and N. Turok, _Shocks in the Early Universe_ , _Phys. Rev. Lett._ 117 (2016) 131301 [1510.02985].
  * [46] L. Landau and E. Lifshitz, _Fluid Mechanics_. Pergamon Press, New York, 1989\.
  * [47] M. Kamionkowski, A. Kosowsky and M. S. Turner, _Gravitational radiation from first order phase transitions_ , _Phys. Rev._ D49 (1994) 2837 [astro-ph/9310044].
  * [48] H. Kurki-Suonio and M. Laine, _Supersonic deflagrations in cosmological phase transitions_ , _Phys. Rev._ D51 (1995) 5431 [hep-ph/9501216].
  * [49] C. Schmid, D. J. Schwarz and P. Widerin, _Amplification of cosmological inhomogeneities from the QCD transition_ , _Phys. Rev. D_ 59 (1999) 043517 [astro-ph/9807257].
  * [50] D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen and G. White, _Theoretical uncertainties for cosmological first-order phase transitions_ , _JHEP_ 04 (2021) 055 [2009.10080].
  * [51] 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].
  * [52] M. Hindmarsh and M. Hijazi, _Gravitational waves from first order cosmological phase transitions in the Sound Shell Model_ , _JCAP_ 12 (2019) 062 [1909.10040].
  * [53] R. Jinno, T. Konstandin and H. Rubira, _A hybrid simulation of gravitational wave production in first-order phase transitions_ , _JCAP_ 04 (2021) 014 [2010.00971].
  * [54] L. Leitao and A. Megevand, _Hydrodynamics of phase transition fronts and the speed of sound in the plasma_ , _Nucl. Phys. B_ 891 (2015) 159 [1410.3875].
  * [55] A. Ekstedt, P. Schicho and T. V. I. Tenkanen, _DRalgo: a package for effective field theory approach for thermal phase transitions_ , 2205.08815.
  * [56] P. H. Ginsparg, _First Order and Second Order Phase Transitions in Gauge Theories at Finite Temperature_ , _Nucl. Phys._ B170 (1980) 388.
  * [57] T. Appelquist and R. D. Pisarski, _High-Temperature Yang-Mills Theories and Three-Dimensional Quantum Chromodynamics_ , _Phys. Rev. D_ 23 (1981) 2305.
  * [58] K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, _Generic rules for high temperature dimensional reduction and their application to the standard model_ , _Nucl. Phys. B_ 458 (1996) 90 [hep-ph/9508379].
  * [59] E. Braaten and A. Nieto, _Effective field theory approach to high temperature thermodynamics_ , _Phys. Rev. D_ 51 (1995) 6990 [hep-ph/9501375].
  * [60] O. Gould and T. V. I. Tenkanen, _On the perturbative expansion at high temperature and implications for cosmological phase transitions_ , _JHEP_ 06 (2021) 069 [2104.04399].
  * [61] A. Gynther and M. Vepsalainen, _Pressure of the standard model at high temperatures_ , _JHEP_ 01 (2006) 060 [hep-ph/0510375].
  * [62] A. Gynther and M. Vepsalainen, _Pressure of the standard model near the electroweak phase transition_ , _JHEP_ 03 (2006) 011 [hep-ph/0512177].
  * [63] J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, _Perturbative Thermal QCD: Formalism and Applications_ , _Phys. Rept._ 880 (2020) 1 [2002.10188].
  * [64] A. Ekstedt, _Higher-order corrections to the bubble-nucleation rate at finite temperature_ , _Eur. Phys. J. C_ 82 (2022) 173 [2104.11804].
  * [65] O. Gould and J. Hirvonen, _Effective field theory approach to thermal bubble nucleation_ , _Phys. Rev. D_ 104 (2021) 096015 [2108.04377].
  * [66] A. Ekstedt, _Bubble Nucleation to All Orders_ , 2201.07331.
  * [67] A. Ekstedt, _Convergence of the nucleation rate for first-order phase transitions_ , 2205.05145.
  * [68] J. Löfgren, M. J. Ramsey-Musolf, P. Schicho and T. V. I. Tenkanen, _Nucleation at finite temperature: a gauge-invariant, perturbative framework_ , 2112.05472.
  * [69] J. Hirvonen, J. Löfgren, M. J. Ramsey-Musolf, P. Schicho and T. V. I. Tenkanen, _Computing the gauge-invariant bubble nucleation rate in finite temperature effective field theory_ , 2112.08912.
  * [70] J. Hirvonen, _Intuitive method for constructing effective field theories_ , 2205.02687.
  * [71] E. V. Shuryak, _Theory of Hadronic Plasma_ , _Sov. Phys. JETP_ 47 (1978) 212.
  * [72] S. A. Chin, _Transition to Hot Quark Matter in Relativistic Heavy Ion Collision_ , _Phys. Lett. B_ 78 (1978) 552.
  * [73] J. I. Kapusta, _Quantum Chromodynamics at High Temperature_ , _Nucl. Phys. B_ 148 (1979) 461.
  * [74] T. Toimela, _The Next Term in the Thermodynamic Potential of QCD_ , _Phys. Lett. B_ 124 (1983) 407.
  * [75] P. B. Arnold and C.-X. Zhai, _The Three loop free energy for pure gauge QCD_ , _Phys. Rev. D_ 50 (1994) 7603 [hep-ph/9408276].
  * [76] P. B. Arnold and C.-x. Zhai, _The Three loop free energy for high temperature QED and QCD with fermions_ , _Phys. Rev. D_ 51 (1995) 1906 [hep-ph/9410360].
  * [77] E. Braaten and A. Nieto, _Free energy of QCD at high temperature_ , _Phys. Rev. D_ 53 (1996) 3421 [hep-ph/9510408].
  * [78] C.-x. Zhai and B. M. Kastening, _The Free energy of hot gauge theories with fermions through g**5_ , _Phys. Rev. D_ 52 (1995) 7232 [hep-ph/9507380].
  * [79] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, _The Pressure of hot QCD up to g6 ln(1/g)_ , _Phys. Rev. D_ 67 (2003) 105008 [hep-ph/0211321].
  * [80] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, _Four loop vacuum energy density of the SU(N(c)) + adjoint Higgs theory_ , _JHEP_ 04 (2003) 036 [hep-ph/0304048].
  * [81] A. Vuorinen, _The Pressure of QCD at finite temperatures and chemical potentials_ , _Phys. Rev. D_ 68 (2003) 054017 [hep-ph/0305183].
  * [82] A. Ekstedt, O. Gould and J. Löfgren, _Radiative first-order phase transitions to next-to-next-to-leading-order_ , 2205.07241.
  * [83] A. D. Linde, _Infrared Problem in Thermodynamics of the Yang-Mills Gas_ , _Phys. Lett._ 96B (1980) 289.
  * [84] J. P. Blaizot, E. Iancu and A. Rebhan, _On the apparent convergence of perturbative QCD at high temperature_ , _Phys. Rev. D_ 68 (2003) 025011 [hep-ph/0303045].
  * [85] 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].
  * [86] S. Inoue, G. Ovanesyan and M. J. Ramsey-Musolf, _Two-Step Electroweak Baryogenesis_ , _Phys. Rev. D_ 93 (2016) 015013 [1508.05404].
  * [87] N. Blinov, J. Kozaczuk, D. E. Morrissey and C. Tamarit, _Electroweak Baryogenesis from Exotic Electroweak Symmetry Breaking_ , _Phys. Rev. D_ 92 (2015) 035012 [1504.05195].
  * [88] L. Niemi, M. J. Ramsey-Musolf, T. V. I. Tenkanen and D. J. Weir, _Thermodynamics of a Two-Step Electroweak Phase Transition_ , _Phys. Rev. Lett._ 126 (2021) 171802 [2005.11332].
  * [89] N. F. Bell, M. J. Dolan, L. S. Friedrich, M. J. Ramsey-Musolf and R. R. Volkas, _Two-Step Electroweak Symmetry-Breaking: Theory Meets Experiment_ , _JHEP_ 05 (2020) 050 [2001.05335].
  * [90] M. Quiros, _Finite temperature field theory and phase transitions_ , in _ICTP Summer School in High-Energy Physics and Cosmology_ , pp. 187–259, 1, 1999, hep-ph/9901312.
  * [91] G. C. Dorsch, S. J. Huber and T. Konstandin, _A sonic boom in bubble wall friction_ , _JCAP_ 04 (2022) 010 [2112.12548].
  * [92] B. Laurent and J. M. Cline, _First principles determination of bubble wall velocity_ , 2204.13120.
  * [93] O. Gould, J. Kozaczuk, L. Niemi, M. J. Ramsey-Musolf, T. V. I. Tenkanen and D. J. Weir, _Nonperturbative analysis of the gravitational waves from a first-order electroweak phase transition_ , _Phys. Rev. D_ 100 (2019) 115024 [1903.11604].
  * [94] P. Schicho, _Multi-loop investigations of strong interactions at high temperatures_ , Ph.D. thesis, U. Bern, 2020.  10.24442/BORISTHESES.1988.
  * [95] K. Farakos, K. Kajantie, K. Rummukainen and M. E. Shaposhnikov, _3-d physics and the electroweak phase transition: A Framework for lattice Monte Carlo analysis_ , _Nucl. Phys. B_ 442 (1995) 317 [hep-lat/9412091].
  * [96] P. Schicho, T. V. I. Tenkanen and G. White, _Combining thermal resummation and gauge invariance for electroweak phase transition_ , 2203.04284.
  * [97] H. H. Patel and M. J. Ramsey-Musolf, _Baryon Washout, Electroweak Phase Transition, and Perturbation Theory_ , _JHEP_ 07 (2011) 029 [1101.4665].
  * [98] A. Ekstedt and J. Löfgren, _A Critical Look at the Electroweak Phase Transition_ , _JHEP_ 12 (2020) 136 [2006.12614].
  * [99] P. B. Arnold and O. Espinosa, _The Effective potential and first order phase transitions: Beyond leading-order_ , _Phys. Rev. D_ 47 (1993) 3546 [hep-ph/9212235].
  * [100] M. Laine, _Thermodynamics and kinetics of cosmological phase transitions: academic dissertation_. Kirjapaino Oy Like, 1994.
  * [101] K. Farakos, K. Kajantie, K. Rummukainen and M. E. Shaposhnikov, _3-D physics and the electroweak phase transition: Perturbation theory_ , _Nucl. Phys. B_ 425 (1994) 67 [hep-ph/9404201].
  * [102] J. I. Kapusta and C. Gale, _Finite-temperature field theory: Principles and applications_ , Cambridge Monographs on Mathematical Physics. Cambridge University Press, 2011, 10.1017/CBO9780511535130.
  * [103] M. Laine and A. Vuorinen, _Basics of Thermal Field Theory_ , vol. 925. Springer, 2016, 10.1007/978-3-319-31933-9, [1701.01554].