# Phase transitions in the early universe

MPP-2020-163, HIP-2020-27/TH

Phase transitions in the early universe

Mark Hindmarsh1,2* , Marvin Lüben3** , Johannes Lumma4$\dagger$, Martin
Pauly4$\ddagger$

1 Department of Physics and Helsinki Institute of Physics,

PL 64, FI-00014 University of Helsinki, Finland

2 Department of Physics and Astronomy, University of Sussex,

Brighton BN1 9QH, United Kingdom

3 Max-Planck-Institut für Physik (Werner-Heisenberg-Institut),

Föhringer Ring 6, 80805 Munich, Germany

4 Institut für Theoretische Physik, Ruprecht-Karls-Universität Heidelberg,

Philosophenweg 16, 69120 Heidelberg, Germany

*<EMAIL_ADDRESS>**<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

These lecture notes are based on a course given by Mark Hindmarsh at the 24th
Saalburg Summer School 2018 and written up by Marvin Lüben, Johannes Lumma and
Martin Pauly. The aim is to provide the necessary basics to understand first-
order phase transitions in the early universe, to outline how they leave
imprints in gravitational waves, and advertise how those gravitational waves
could be detected in the future. A first-order phase transition at the
electroweak scale is a prediction of many theories beyond the Standard Model,
and is also motivated as an ingredient of some theories attempting to provide
an explanation for the matter-antimatter asymmetry in our Universe.

Starting from bosonic and fermionic statistics, we derive Boltzmann’s equation
and generalise to a fluid of particles with field dependent mass. We introduce
the thermal effective potential for the field in its lowest order
approximation, discuss the transition to the Higgs phase in the Standard Model
and beyond, and compute the probability for the field to cross a potential
barrier. After these preliminaries, we provide a hydrodynamical description of
first-order phase transitions as it is appropriate for describing the early
Universe. We thereby discuss the key quantities characterising a phase
transition, and how they are imprinted in the gravitational wave power
spectrum that might be detectable by the space-based gravitational wave
detector LISA in the 2030s.

##  1 Introduction

These lecture notes are intended to provide an introduction to the topic of
phase transitions in the early universe, focusing on a possible first-order
phase transition at temperatures around the scale of electroweak symmetry-
breaking, which the universe reached at an age of around $10^{-11}$ s.

Phase transitions are a generic, but not universal, feature of gauge field
theories, like the Standard Model, which are based on elementary particle mass
generation by spontaneous symmetry-breaking [1, 2]. When there is a phase
transition in a gauge theory it is (except for special parameter choices)
first-order, which means that just below the critical temperature, the
universe transitions from a metastable quasi-equilibrium state into a stable
equilibrium state, through a process of bubble nucleation, growth, and merger
[3, 4, 5, 6]. Such a first-order phase transition in the early universe
naturally leads to the production of gravitational waves [7, 8]. If it took
place around the electroweak scale, by which we mean temperatures in the range
100 – 1000 GeV, the gravitational wave signal could lie in the frequency range
of the upcoming space-based gravitational wave detector LISA (Laser
Interferometer Space Antenna) [9]. The approval of the mission, and the
detection of gravitational waves [10], has generated enormous interest in
phase transitions in the early universe.

While the Standard Model has a crossover rather than a true phase transition
[11], many extensions of the Standard Model, e.g. with extra scalar fields,
lead to first-order phase transitions at the electroweak scale. Gravitational
wave signatures are therefore a fascinating new window towards new physics,
complementary to that provided by the Large Hadron Collider (see e.g. Ref.
[12] for a recent review).

A further motivation for studying electroweak phase transitions is that one of
the requirements to explain the matter-antimatter asymmetry in the universe
[13] is a departure from thermal equilibrium, which is inevitable in a first-
order phase transition. The asymmetry is quantified in terms of the net baryon
number of the universe, leading to the name baryogenesis. We will
unfortunately not have time to study electroweak baryogenesis in these
lectures, and refer the interested reader to e.g. Refs. [14, 15, 16, 17].

Let us also briefly mention that non-thermal phenomena might lead to a sizable
gravitational wave background. One example of this is the production of
gravitational waves associated with preheating; at the end of inflation the
inflaton decays to Standard Model particles. The resulting distributions might
strongly deviate from local thermal equilibrium. Their violent dynamics could
lead to a background of gravitational waves [18, 19, 20]. We will not touch on
this topic and refer the interested reader to Refs. [21, 17, 22].

A thorough study of early universe phase transitions, gravitational wave
production and detection, requires quite a lot of theoretical apparatus from
particle physics and cosmology, which could not be covered in a short lecture
course. It is assumed that the student has done advanced undergraduate courses
on statistical physics and general relativity, and has been introduced to
particle physics and cosmology. Wherever possible, the full machinery of
thermal quantum field theory is avoided. The aim is to provide a direct route
to some important results, and motivate further study and, we hope, research.

The main points we would like the reader to take away are: that the
gravitational wave power spectrum from a first-order phase transition is
calculable from a few thermodynamic properties of matter at very high
temperatures; that these parameters are computable from an underlying quantum
field theory; and that these parameters are measurable by LISA. The final
point we would like to make is that we can see in outline how the
computations, calculations and measurements could be done, but they are far
from concrete methods. There is therefore a lot of exciting work to be done in
the years leading up to LISA’s launch in 2034 to realise the mission’s
potential scientific reward. These notes are organised as follows. In Sec. 2
we review basic thermodynamics of non-interacting fields and we discuss the
different relevant thermodynamic quantities for both fermions and bosons. In
Sec. 3 we introduce weak interactions among the fields and derive the thermal
Higgs potential. Further, we summarize phase transitions in the Standard Model
as well as in models beyond the Standard Model. In Sec. 4 we consider the
distribution function of a relativistic fluid and derive the relativistic
Boltzmann equation. We generalise the preceding results and study the
hydrodynamics of a fluid with a field-dependent mass in Sec. 5 This set-up is
analogous to the hydrodynamics with electromagnetic forces, which is governed
by the Vlasov equation. In Sec. 6 we study the transition of the Higgs from
the false, symmetric phase to the new symmetry-breaking phase and apply the
process to the early universe. After these preliminaries, in Sec. 7 we provide
a hydrodynamical description of the phase transition in the early unverse. We
then discuss the different sources of gravitational waves during a first-order
phase transition and the expected power spectra in Sec. 8. Finally, we provide
a summary of these lectures and comment on open issues in Sec. 9.

Conventions. Throughout these notes we set $\hbar=k_{\text{B}}=c=1$ and just
re-introduce these constants occasionally. We try to stick to a $(-,+,+,+)$
metric signature. 4-vectors are denoted by roman letters, e.g, $x$, $p$, and
$F$ with greek letters as space-time indices, e.g., $\mu,\nu=0,1,2,3$. Spatial
indices are latin letters, e.g., $i,j=1,2,3$ and we denote 3-vectors with an
arrow, e.g., $\vec{x}$ and $\vec{p}$.

##  2 Thermodynamics of free fields

We start by studying thermodynamic properties of free bosonic and fermionic
fields. For both cases, we derive the partition function, from which all
thermodynamic quantities can be derived. We will be particularly interested in
the free energy, as it can be used to find the equilibrium states of a theory.

###  2.1 Basic thermodynamics - the bosonic harmonic oscillator

The basic object of thermodynamics is the partition function

$Z(T)=\Tr\left[e^{-\beta\hat{H}}\right]$ (2.1)

where $\hat{H}$ is the Hamiltonian operator and $\beta={1}/{T}$ the inverse
temperature with $T$ the temperature. The free energy, entropy and energy of
the system are given by

$\displaystyle F$ $\displaystyle=$ $\displaystyle-T\ln Z\,,$ (2.2)
$\displaystyle S$ $\displaystyle=$ $\displaystyle-\partialderivative{F}{T}\,,$
(2.3) $\displaystyle E$ $\displaystyle=$ $\displaystyle-\partialderivative{\ln
Z}{\beta}\,,$ (2.4)

respectively. First, we will study a single bosonic harmonic oscillator, the
simplest case. To compute its partition function we consider

$Z_{\text{bho}}=\sum_{n=0}^{\infty}\matrixelement{n}{e^{-\beta\hat{H}}}{n}$
(2.5)

where the $\ket{n}$ are the eigenstates of the Hamiltonian $\hat{H}$ of the
harmonic oscillator with angular frequency $\omega$, together satisfying

$\hat{H}\ket{n}=\omega\left(n+\frac{1}{2}\right)\ket{n}\,.$ (2.6)

For the partition function and the free energy this yields

$\displaystyle Z_{\text{bho}}(T,\omega)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\exp[-\beta\omega\left(n+\frac{1}{2}\right)\Bigg{]}=\frac{e^{-\beta{\omega}/{2}}}{1-e^{-\beta\omega}},$
(2.7) $\displaystyle F_{\text{bho}}(T,\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\omega+T\,\ln\left(1-e^{-\beta\omega}\right),$ (2.8)

where the first term in the free energy describes the ground state energy, and
the second term is the thermal contribution.

Next we turn to the partition function for a field, or equivalently for a
collection of harmonic oscillators. We consider the field operator
$\hat{\phi}(\vec{x},t)$ and decompose it into its Fourier modes

$\hat{\phi}(\vec{x},t)=\int\frac{\differential[3]{k}}{(2\pi)^{3}}\frac{1}{2\omega_{\vec{k}}}\left(\hat{a}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}}+\hat{a}^{\dagger}_{\vec{k}}e^{i\vec{k}\cdot\vec{x}}\right)\,,$
(2.9)

where we postulate that the operators $\hat{a},\hat{a}^{\dagger}$ satisfy the
commutation relation

$\displaystyle\commutator{\hat{a}_{\vec{k}}}{\hat{a}^{\dagger}_{\vec{k^{\prime}}}}$
$\displaystyle=$ $\displaystyle
2\omega_{\vec{k}}(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})\,$ (2.10)
$\displaystyle\commutator{\hat{a}_{\vec{k}}}{\hat{a}_{\vec{k^{\prime}}}}$
$\displaystyle=$
$\displaystyle\commutator{\hat{a}^{\dagger}_{\vec{k}}}{\hat{a}^{\dagger}_{\vec{k^{\prime}}}}=0\,.$
(2.11)

The equation of motion for the free field is the Klein-Gordon equation,

$\left(\Box+m^{2}\right)\hat{\phi}(\vec{x},t)=0$ (2.12)

which in terms of the Fourier modes reads

$(k^{0})^{2}=\omega_{\vec{k}}^{2}=\vec{k}^{2}+m^{2}\,.$ (2.13)

This dispersion relation does not involve different momenta and hence the
different modes are not coupled. The free scalar field is a collection of
independent harmonic oscillators, one for each momentum mode $|\vec{k}|$. The
partition function of a bosonic field (indicated by the subscript $B$) thus
factorizes into

$Z_{B}=\prod_{\vec{k}}Z_{\text{bho}}(T,\omega_{\vec{k}})\,,$ (2.14)

where the multiplication here is a symbolic notation for a product over all
wavenumbers. It can be given meaning by working in finite volume, with the
infinite volume limit taken at the end of the calculation.

The free energy of a free bosonic field is given by

$\displaystyle F_{B}$ $\displaystyle=$ $\displaystyle-T\,\ln
Z_{B}=-T\,\sum_{\vec{k}}\ln Z_{\text{bho}}(T,\omega_{\vec{k}})$ (2.15)
$\displaystyle=$
$\displaystyle\sum_{\vec{k}}\left[\frac{1}{2}\omega_{\vec{k}}+T\,\ln\left(1-e^{-\beta\omega_{\vec{k}}}\right)\right]\,,$
(2.16)

Again, the sum over all momenta $\vec{k}$ is defined over a finite volume
$\mathcal{V}$. In the infinite volume limit $\mathcal{V}\rightarrow\infty$,
the sum is replaced by an integration as

$\sum_{\vec{k}}\rightarrow{\mathcal{V}}\int\frac{\differential[3]{k}}{(2\pi)^{3}}\,.$
(2.17)

The free energy density, i.e., the free energy normalized to the volume
$\mathcal{V}$, then becomes

$f_{B}=\frac{F_{\text{B}}}{\mathcal{V}}=V_{0,B}+T\int\frac{\differential[3]{k}}{(2\pi)^{3}}\ln\left(1-e^{-\beta\omega_{\vec{k}}}\right)\,.$
(2.18)

The first term $V_{0,B}$ is the energy density of the zero-temperature ground
state, which is divergent, cf. Eq. (2.16). The same divergence is encountered
in quantum field theories at zero temperature. In the following, we assume
that it is regularized with an appropriate counter-term. The standard
renormalisation convention takes the zero-temperature ground state free energy
to be zero.

Figure 1: This figure shows the dimensionless function $J_{B}$ that is
proportional to the free energy of bosons as defined in Eq. (2.19), as a
function of mass-to-temperature ratio (thick line). Also the expansions for
large $T$ (dashed), Eq. (2.1) and small $T$ (dotted), Eq. (2.20) are shown.
The large-T expansion is performed up to order four in $m/T$, being a good
approximation up to $m/T\sim 1.1$.

Due to the integration over all momenta, $f_{B}$ can only depend on $T$ and
$m$, where $m$ only appears as ${m}/{T}$. From dimensional analysis we infer
that the free energy density hence takes the form

$f_{B}(T,m)=T^{4}J_{B}\left(\frac{m}{T}\right)\,,$ (2.19)

where $J_{B}({m}/{T})$ is a dimensionless function. While the integral in Eq.
(2.18) cannot be solved exactly for all values of ${m}/{T}$, analytic
approximations exist in the low and the high temperature regimes.

In the low temperature regime, ${m}/{T}\gg 1$, one expands in ${T}/{m}$ and
gets

$J_{B}\left(\frac{m}{T}\right)=-\left(\frac{m}{2\pi
T}\right)^{\frac{3}{2}}e^{-{m}/{T}}\left(1+\mathcal{O}\left(\frac{T}{m}\right)\right)\,,$
(2.20)

recovering the familiar distribution function of Maxwell-Boltzmann statistics.

In the high temperature case ${m}/{T}\ll 1$ one expands in ${m}/{T}$ to obtain

$\displaystyle J_{B}\left(\frac{m}{T}\right)$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}}{90}+\frac{1}{24}\left(\frac{m}{T}\right)^{2}-\frac{1}{12\pi}\left(\left(\frac{m}{T}\right)^{2}\right)^{\frac{3}{2}}$
$\displaystyle-\frac{1}{2(4\pi)^{2}}\left(\frac{m}{T}\right)^{4}\left[\ln\left(\frac{1}{4\pi}\frac{m}{T}e^{\gamma_{\rm
E}}\right)-\frac{3}{4}\right]$
$\displaystyle+\mathcal{O}\left(\left(\frac{m}{T}\right)^{6}\right)\,.$

Here $\gamma_{\rm E}\approx 0.57721$ is the Euler–Mascheroni constant. While
the first two terms follow in a relatively simple way using the
$\zeta$-function, the third and fourth terms are more complicated in nature:
they are non-analytic in the fundamental expansion parameter
${m^{2}}/{T^{2}}$, and can only be derived using more advanced methods. For
details, we refer to Ref. [23].

A numerical evaluation of the function $J_{B}(m/T)$ is depicted in Fig. 1,
along with its high and low temperature approximations. It can be seen that
the high temperature approximation is good even up to $m/T\simeq 2$.

###  2.2 The fermionic harmonic oscillator

For fermions the Pauli exclusion principle holds: a quantum state can only be
occupied by a single fermion at once. To compute the fermionic partition
function, we therefore sum over only the occupation numbers $0$ and $1$ in
order to respect the Pauli exclusion principle, arriving at

$Z_{\text{fho}}(T,\omega)=\sum_{n=0}^{1}\matrixelement{n}{e^{-\beta\hat{H}}}{n}=e^{\beta{\omega}/{2}}\left(1+e^{-\beta{\omega}/{2}}\right)\,,$
(2.22)

after using Eq. (2.6). Consequently, the free energy for fermions is given by

$F_{\text{fho}}(T,\omega)=-\frac{1}{2}\omega-T\,\ln\left(1+e^{-\beta\omega}\right)\,.$
(2.23)

In order to generalize the above expression to fermionic fields we introduce a
Dirac spinor field $\hat{\Psi}_{\alpha}(\vec{x},t)$, which creates and
destroys massive fermions. Every such spinor has four components, denoted by
the index $\alpha$. The components describe particles and antiparticles, each
of which have two spin degrees of freedom.

The spinor can be decomposed into Fourier modes, where the Fourier
coefficients are operators subject to a set of anticommutation relations,
which need to respect the fermionic nature of $\Psi_{\alpha}$. Analogously to
the bosonic case, the free fermionic field is a collection of independent
harmonic oscillators, four for each momentum mode $|\vec{k}|$. The partition
function for a fermionic field $F$ hence factorizes into

$Z_{F}=4\prod_{\vec{k}}Z_{\text{fho}}(T,\omega_{\vec{k}})\,,$ (2.24)

where $Z_{\text{fho}}$ is the partition function of a single fermionic
harmonic oscillator, cf. Eq. (2.23). This leads to the fermionic free energy

$F_{F}=-4\sum_{\vec{k}}\Bigg{[}\frac{1}{2}\omega_{\vec{k}}+T\,\ln\left(1+e^{-\beta\omega_{\vec{k}}}\right)\Bigg{]}\,.$
(2.25)

The factor of $4$ arises due to the four components of an individual spinor.
In the continuum limit (infinite volume) one needs to replace the sum by an
integration, as in Eq. (2.17). The free energy density for each fermionic
degree of freedom is thus

$f_{F}=-V_{0,{F}}-T\int\frac{\differential[3]{k}}{(2\pi)^{3}}\ln\left(1+e^{-\beta\omega_{\vec{k}}}\right)\,.$
(2.26)

The vacuum energy density $V_{0,\text{F}}$ is again divergent, but comes with
the opposite sign compared to the bosonic case. We assume that the vacuum
energy is regularised by an appropriate counter-term such that we can take it
to be zero in the following.

Figure 2: This figure shows the dimensionless function $J_{F}$ that is
proportional to the free energy of fermions as defined in Eq. (2.27) as a
function of mass-to-temperature ratio (thick line). Additionally, the
expansion for large $T$ (dashed), Eq. (2.2) and small $T$ (dotted) , in
analogy to Eq. (2.20) are shown. Note that in the small $T$ limit, both the
fermionic and the bosonic expansions agree. The large $T$ expansion is
performed up to order four in $m/T$, working well up to $m/T\sim 0.5$, hence
being sligthly worse than the bosonic high-T expansion, depicted in Fig. 1 .

In analogy to the bosonic case the free energy density can be written as

$f_{F}=T^{4}J_{F}\left(\frac{m}{T}\right)\,.$ (2.27)

In Fig. 2, a numerical evaluation of $J_{F}(m/T)$ is shown together with its
small and large temperature expansions. The integral in Eq. (2.26) cannot be
solved analytically for all values of $m/T$, but in the limits of small and
large temperatures. In the low temperature limit the expansion agrees with Eq.
(2.20). As expected, at low energies the quantum nature of the field becomes
irrelevant and one recovers the Maxwell-Boltzmann statistics for both,
fermions and bosons. In the high temperature limit the free energy can be
approximated as

$\displaystyle J_{F}\left(\frac{m}{T}\right)$ $\displaystyle=$
$\displaystyle-\frac{7}{8}\frac{\pi^{2}}{90}+\frac{1}{48}\left(\frac{m}{T}\right)^{2}$
$\displaystyle-\frac{1}{2(4\pi)^{2}}\left(\frac{m}{T}\right)^{4}\left[\ln\left(\frac{1}{\pi}\frac{m}{T}e^{\gamma_{\rm
E}}\right)-\frac{3}{4}\right]$
$\displaystyle+\mathcal{O}\left(\left(\frac{m}{T}\right)^{6}\right)\,.$

Compared to Eq. (2.1), the term that appeared with a power of $3/2$
disappeared, and the first term has a prefactor of $7/8$. Eqs. (2.27) and
(2.2) together give the free energy of a single Dirac fermion field.

##  3 Phase transition in field theory

After having discussed the free energy of free bosons and fermions, let us
introduce interactions among these fields. In a weakly interacting field
theory one can compute the free energy as a perturbation around the free
energy of a free field. In this section we will present a setup, which is
tailored to discuss phase transitions in the Standard Model. Phase transitions
in weakly coupled gauge theories were first discussed in Refs. [1, 2].

particle | mass [GeV] | $c_{i}$
---|---|---
$t$ | 172.76 | $y_{t}$
$H$ | 125.18 | $\sqrt{2\lambda}$
$Z$ | 91.19 | $\sqrt{g^{2}+{g^{\prime}}^{2}}/\sqrt{2}$
$W^{\pm}$ | 80.38 | $g/\sqrt{2}$
Table 1: The zero-temperature masses and the mass proportionality constants
$c_{i}$, defined in Eq. (3.1), for the most massive fields in the Standard
Model. Here, $y_{t}$ is the Higgs-Yukawa coupling of the top quark, $\lambda$
the Higgs self-coupling, and $g$ and $g^{\prime}$ the gauge couplings of the
$W^{a}_{\mu}$ and $B_{\mu}$ bosons.

###  3.1 The Standard Model at weak coupling

In the Standard Model the masses of fermions and gauge bosons $M_{i}$ depend
linearly on the magnitude of the Higgs field $\phi$,

$M_{i}(\phi)=c_{i}\phi\,,$ (3.1)

with the $c_{i}$ proportional to the dimensionless coupling constants. The
index $i$ labels the Standard Model fields that couple to the Higgs. Today,
the Higgs is in its broken phase and the Higgs field assumes its vacuum
expectation value, $\phi=v_{\text{EW}}\simeq 246\,\text{GeV}$, which
determines the particle masses we observe in experiments. This value for the
field is dynamically determined by the minimisation of the zero-temperature
potential for the Higgs field,

$V_{0}(\phi)=\frac{\lambda}{4}\left(\phi^{2}-v_{\text{EW}}^{2}\right)^{2}\,,$
(3.2)

The relation of the $c_{i}$ to the standard coupling constants of the Standard
Model are given for the most massive fields in Tab. 1. One can see that for
these fields, the values of $c_{i}$ are all $\mathcal{O}(1)$. For the other
fields of the standard model, $c_{i}\ll 1$.

The Higgs particle is a quantised fluctuation around the ground state, with
mass

$M_{H}=\sqrt{V^{\prime\prime}_{0}(v_{\text{EW}})}=\sqrt{2\lambda}v_{\text{EW}}.$

The Higgs field is unique in that its mass does not in general depend linearly
on $\phi$. At this level of treatment, we will not need to know that in the
Standard Model, the Higgs field is a two-component vector of complex scalar
fields $\Phi$, but for completeness we mention that
$\phi^{2}=\Phi^{\dagger}\Phi/{2}$.

energy scale | event
---|---
$100$ GeV |  $t$ non-relativistic
$1$ GeV |  $b$ non-relativistic
$500$ MeV |  $c$, $\tau$ non-relativistic
$200$ MeV | QCD phase transition
$30$ MeV |  $\mu$ non-relativistic
$2$ MeV |  $\nu$ freeze-out
$0.2$ MeV |  $e$ non-relativistic
$1$ eV | matter-radiation equality
$0.1$ eV | photon decoupling
Table 2: An overview over events happening at different energy scales in the
early universe. These determine the effective number of degrees of freedom in
the Standard Model at a certain energy scale. Figure 3: This figure shows the
effective number of relativistic degrees of freedom $g_{\text{eff}}$ of a
Standard Model plasma as a function of temperature, taking into account
interactions between particles, with both perturbative and lattice methods
[24].

The free energy density $f$ of a gas of Standard Model particles is given by
the zero temperature result (3.2), plus terms that arise due to the
interaction with the Higgs according to Eq. (3.1). This gives

$\displaystyle f$ $\displaystyle=$ $\displaystyle
V_{0}(\phi)+\sum_{B}f_{B}+\sum_{F}f_{F}$ (3.3) $\displaystyle=$ $\displaystyle
V_{0}(\phi)+T^{4}\Bigg{[}\sum_{B}J_{B}\left(\frac{M_{B}}{T}\right)+\sum_{F}J_{F}\left(\frac{M_{F}}{T}\right)\Bigg{]}\,.$
(3.4)

Here, we sum over all fermions $F$ and bosons $B$ that are relativistic at
temperature $T$. For large temperatures, we can write the free energy density
as

$f=-g_{\text{eff}}\frac{\pi^{2}}{90}T^{4}+V_{T}(\phi)\,.$ (3.5)

Here, $g_{\text{eff}}$ is the effective number of relativistic degrees of
freedom, given by

$g_{\text{eff}}=\frac{7}{8}4N_{\rm F}+3N_{\rm V}+2N_{{\rm V}0}+N_{\rm S}$
(3.6)

where $N_{\rm F}$ is the number of Dirac fermions, $N_{\rm V}$ is the number
of massive vectors, $N_{{\rm V}0}$ is the number of massless vectors and
$N_{\rm S}$ is the number of scalars. The prefactors account for the degrees
of freedom of each of the particles. In the case of only bosons or only
fermions this expression reduces to the first term in Eq. (2.1) or Eq. (2.2),
resp. For the Standard Model at high energies this value is
$g_{\text{eff}}=106.75$.111It is a good exercise to verify this. Note that the
neutrinos of the Standard Model count as $N_{F}=1/2$, as they are two-
component spinors.  As the temperature decreases, so does the effective number
of relativistic degrees of freedom, as more and more particles become non-
relativistic, or are bound together into hadrons. The function
$g_{\text{eff}}(T)$ for the Standard Model is shown in Fig. 3, using data
taken from Ref. [24], where interactions between particles (and not just the
mass generation effect of the Higgs) are also taken into account. Tab. 2
summarizes key temperatures that affect $g_{\rm eff}$.

The second, field-dependent term in Eq. (3.5), $V_{T}(\phi)$, is called the
thermal effective Higgs potential. According to Eqs. (2.1) and (2.2) it is
given by [23, 25, 26, 27]

$\displaystyle V_{T}(\phi)$ $\displaystyle=$ $\displaystyle
V_{0}(\phi)+\frac{T^{2}}{24}\left(\sum_{S}M_{{S}}^{2}(\phi)+3\sum_{V}M_{{V}}^{2}(\phi)+2\sum_{F}M_{{F}}^{2}(\phi)\right)$
(3.7)
$\displaystyle-\frac{T}{12\pi}\left(\sum_{S}\left(M_{{S}}^{2}(\phi)\right)^{\frac{3}{2}}+\sum_{V}\left(M_{{V}}^{2}(\phi)\right)^{\frac{3}{2}}\right)$
$\displaystyle+\text{higher order terms}\,.$

Here $M_{S}$, $M_{V}$, $M_{F}$ are the masses of the scalar fields $S$, vector
fields $V$ and fermionic fields $F$, which are related to the expectation
value of the Higgs as in Eq. (3.1).

For small field values the expression for the scalar masses in Eqn. (3.7) is
negative, leading to imaginary contributions to the thermal Higgs potential.
This problem does not arise if instead of the bare mass one takes into account
the thermal mass in (3.7). However, the thermal mass itself depends on the
thermal potential. A self-consistent solution then requires resummation
techniques [28, 29]. In addition, technical subtleties associated with gauge
choice in the Higgs sector arise [30] \- the resulting effective potential is
gauge dependent, while observable quantities remain gauge independent, as they
should. For further details we refer the reader to Refs. [31, 27].

For high temperatures, the thermal effective potential can be approximated by
an expansion in ${\phi}/{T}$ yielding

$V_{T}(\phi)=\frac{D}{2}\left(T^{2}-T_{0}^{2}\right)\phi^{2}-\frac{A}{3}T\phi^{3}+\frac{\lambda_{T}}{4!}\phi^{4}+\dots\,,$
(3.8)

where $A$, $D$ are constants and $\lambda_{T}$ depends only logarithmically on
the temperature. In the Standard Model we have

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{1}{12\pi\phi^{3}}\left(M_{H}^{3}+6M_{W}^{3}+3M_{Z}^{3}\right)$
(3.9) $\displaystyle D$ $\displaystyle=$
$\displaystyle\frac{1}{12\phi^{2}}\left(M_{H}^{2}+6M_{W}^{2}+3M_{Z}^{2}+6M_{t}^{2}\right)$
(3.10) $\displaystyle\lambda_{T}$ $\displaystyle\simeq$ $\displaystyle\lambda$
(3.11) $\displaystyle T_{0}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{1}{2D}}M_{H},$ (3.12)

where we have dropped the logarithmic dependence of $\lambda_{T}$ on $T$.
Here, the subscripts $H$, $W$, $Z$, and $t$ denote the Higgs-boson, $W$\- and
$Z$-bosons, and the top quark of the Standard Model. Notice that only bosons
contribute to the cubic term in the potential.

The form of this quartic potential is sketched for different values of the
temperature in Fig. 4. For large temperatures, $T\gg T_{c}$, the potential has
a minimum at $\phi=0$, which is the only ground state or equilibrium state of
the system. We will refer to the ground state where $\phi=0$ also as symmetric
phase. As the temperature drops, a second, but higher lying minimum develops
as represented by the dark green line. Both minima are degenerate at the
critical temperature

$T_{\rm c}=T_{0}\left(1-\frac{2}{9}\frac{A^{2}}{\lambda
D}\right)^{-\frac{1}{2}}.$ (3.13)

which is well-defined only if $2A^{2}<9\lambda D$. This case will be of
particular interest to us, as in this case the two minima at $\phi=0$ and
$\phi(T_{\rm c})=2AT_{\rm c}/3\lambda$ are separated by a free energy barrier,
signaling a first-order phase transition. Below the critical temperature, the
system can supercool, staying in the false ground state at $\phi=0$ for some
time, before transitioning to the global minimum. In the bosonic case, this is
reflected by the cubic term. We will refer to the ground state where the Higgs
field is non-zero also as Higgs phase. For $T=0$ the thermal corrections are
absent and the minimum is at $\phi=v_{\text{EW}}\simeq 246\,\text{GeV}$.

Figure 4:  The figure shows the thermal effective Higgs potential
$V_{T}(\phi)$ at different temperatures. For large temperatures $T\gg T_{\rm
c}$ (red) the potential has a minimum at $\phi=0$ and the ground state is
symmetric. Below the temperature $T_{1}>T_{\rm c}$ (dark green) a second, but
higher lying minimum develops. At the critical temperature $T_{c}$ (green)
both minima are degenerate. Below the critical temperature, the new minimum at
non-zero field value is the global minimum representing the true (stable)
ground state.

###  3.2 Breakdown of weak coupling

So far we have assumed that the free energy is only slightly altered by
interactions, an assumption of weak (i.e. small) couplings between particles.
Moreover, we only included interactions with the Higgs field, in their
simplest form of a mass generation effect. To properly include interactions,
one really has to study and apply thermal field theory [27]. In this section
we merely sketch where the assumption of weak coupling breaks down, with a
qualitative argument using the statistical mechanics of a field in 3
dimensions.

In an interacting theory, one splits the Hamiltonian
$\hat{H}=\hat{H}_{0}+\hat{H}_{I}$ into a free and an interacting Hamiltonian.
As an example inspired by the Standard Model consider a scalar field (not
necessarily the Higgs) with an interaction Hamiltonian

$\hat{H}_{I}=g^{2}\int\differential[3]{x}\hat{\phi}^{4}\,,$ (3.14)

with $g^{2}$ an arbitrary dimensionless coupling constant. We leave the mass
of this scalar field free. Weak coupling means that $g^{2}\ll 1$.

We can then try to compute the partition function

$Z=\Tr\left[e^{-\beta(\hat{H}_{0}+\hat{H}_{I})}\right],$ (3.15)

by expanding in powers of the coupling constant. This is a non-trivial
exercise, but it turns out that we are in fact expanding in the parameter

$\varepsilon=g^{2}f(\vec{k})$ (3.16)

with $f$ the phase space density. For a boson,

$f(\vec{k})=\frac{1}{e^{\beta\omega_{\vec{k}}}-1}$ (3.17)

which approaches ${T}/{\omega_{\vec{k}}}$ for frequencies low compared with
the temperature, $\omega_{\vec{k}}\ll T$. In this limit, the expansion
parameter reads

$\varepsilon=\frac{g^{2}T}{\omega_{\vec{k}}}$ (3.18)

which is greater than unity for $k\lesssim g^{2}T$. The expansion parameter
diverges as $|\vec{k}|\to 0$ (in the “infrared”) for massless bosons, cf. Eq.
(2.13). We therefore learn that in the case of massless bosons at zero
chemical potential a perturbative expansion in powers of $g$ breaks down in a
thermal state, at any temperature [32], for momenta $k\lesssim g^{2}T$.

However, thermal corrections contribute to the mass of a thermal state which
have to be taken into account. One can apply the above argument to the $W$,
$Z$ and gluons of the Standard Model, which have an interaction term of a
similar form.222Indeed, this infrared problem was first pointed out for gauge
bosons [32].

For gauge fields, one distinguishes between the electric mass (a mass
parameter appearing in the wave equation for the timelike component of the
gauge field $A_{0}$) and the corresponding magnetic mass for the spacelike
components $A_{i}$. The timelike component of the gauge field behaves like a
scalar field, both of which have a mass-temperature relation as

$m^{2}(T)=m_{0}^{2}+cg^{2}T^{2}\,,$ (3.19)

where $c$ is a theory-dependent constant, and $m_{0}$ is the mass of the field
at zero temperature. Therefore, the expansion parameter $\epsilon$ for
massless gauge bosons (such as the photon) with $m_{0}=0$ is of the order of
the coupling constant, $\epsilon\sim g\ll 1$. This means that for fields with
electric mass perturbation theory is trustworthy at any temperature for small
coupling constants. Physically, the electric mass makes the electric field at
a distance $r$ from a static charge behave as $r^{-1}\exp[-m(T)r]$, that is,
the electric field is screened. The electric mass is none other than the
inverse Debye length: the free charges in the plasma become polarised around a
source.

The magnetic mass, on the other hand, vanishes in perturbation theory.
Therefore, the expansion parameter $\epsilon$ is divergent in the IR and one
should be suspicious of perturbation theory. The vanishing of the perturbative
magnetic mass turns out not to matter for the photon, as it has no self-
interaction terms in its Hamiltonian, but for the other gauge bosons of the
Standard Model our naive perturbation theory definitely breaks down. To study
phase transitions, more involved methods such as a combination of advanced
resummation techniques and lattice simulations are required.

###  3.3 Beyond weak coupling and the Standard Model

Figure 5: The phase diagram of the Standard model. For Higgs masses of
$m_{H}\lesssim 75$ GeV the Standard Model undergoes a first-order phase
transition. For larger Higgs masses, there is no phase transition between the
symmetric phase $\phi=0$ and the Higgs phase $\phi=v_{\text{EW}}$, but a
cross-over. Including higher-order interactions changes the picture
significantly.

In more advanced calculations based on numerical computations of the partition
function, the following picture emerges [11, 33, 34, 35, 36, 37]. One can
study the Standard Model (in fact any gauge theory with spontaneous symmetry-
breaking) in the 2-dimensional space spanned by the temperature and the ratio
of the Higgs mass to the gauge boson mass. To simplify matters when discussing
the Standard Model, one can take the gauge boson mass to be the one of the
$W$-boson, $m_{W}\simeq 80$ GeV, and take one of the parameters to be the
Higgs mass. This leads to the picture presented in Fig. 5.

If the ratio of the Higgs mass to the gauge boson mass is small, the simple
picture outlined in the previous sections is correct: the perturbative
evaluation of the thermal potential is reasonably accurate, and there is
indeed a first-order phase transition. This would correspond to the case of
the Standard Model if the Higgs mass were much less than $80$ GeV. However, as
the ratio increases, the strength of the transition, as measured for example
by the latent heat, decreases. At a critical value for this ratio, the latent
heat goes to zero. Above this critical value the transition is a cross-over.

In case of a cross-over the system smoothly changes from the symmetric phase
to the Higgs phase. The situation is similar to water at high pressure, whose
density smoothly decreases with temperature, rather than making a sharp
transition from the vapour to the liquid phase.

Precisely at the critical ratio, the transition is second-order, meaning that
the effective thermal mass of the Higgs $M_{H}(T)$ goes to zero, and its
correlation length diverges. The same phenomenon happens with water at its
critical point (374∘ C and 218 atm), where the divergence of the correlation
length can be observed as the phenomenon of critical opalescence. In terms of
the zero-temperature Higgs mass, the critical point is at around 80 GeV. Given
the measured value for the Higgs mass of 125 GeV, the Standard Model is well
into the cross-over region.

However, in theories beyond the Standard Model, the electroweak phase
transition can be a first-order phase transition. Indeed already the inclusion
of a $\phi^{6}$ operator in the Higgs potential could lead to a first-order
phase transition [38, 39, 40]. Such a term is not allowed by the Standard
Model, but it could be part of an effective field theory describing new
physics.

The motivation to study models with new physics includes providing an
explanation for the matter-antimatter asymmetry in the Universe [41] to
explaining dark matter [42, 43]. There are further shortcomings of the
Standard Model that need to be addressed [43].

Many such extensions of the Standard Model include extra scalars, which can
give first-order phase transitions. Examples include coupling the Standard
Model to an extra Higgs SU(2) singlet, doublet (“2HDM”) or triplet. Further,
extensions of the Standard Model with supersymmetry automatically include
extra scalars, although it seems that the simplest such extensions do not have
a first-order phase transition. Possibilities exist beyond the framework of
weakly-coupled field theory. A review of Standard Model extensions with first-
order phase transitions is given in a recent LISA Cosmology Working Group
report [12].

In the following chapters we will study the details of the dynamics of such a
first-order phase transition.

##  4 Relativistic hydrodynamics

In order to understand the dynamics of a first-order phase transition, we need
an appropriate hydrodynamic description of the early universe. In this section
we study interacting bosonic and fermionic particles in the thermodynamic
limit. For the resulting distribution function of a relativistic fluid, we
will derive the relativistic Boltzmann equation. For now we focus on the case
where the particle masses are constant throughout spacetime. We mostly follow
Ref. [44] in this section.

###  4.1 Distribution function

In the presence of several interacting harmonic oscillators, we need to also
include a number operator $\hat{N}$ in the definition of the partition
function

$Z=\Tr\left[e^{-\beta(\hat{H}-\mu\hat{N})}\right]\,,$ (4.1)

where $\mu$ is the chemical potential. The number operator is obtained as

$\hat{N}=T\partialderivative{\mu}\ln Z\,.$ (4.2)

With the same procedure as described for the $1$-particle partition function
in Sec. 2, we arrive at the partition functions and number operators for free
bosonic, indicated by $B$, and fermionic, indicated by $F$, fields with
3-momentum $\vec{p}$:

bosons: $\displaystyle Z_{B}=\frac{e^{-\beta
E_{\vec{p}}/2}}{1-e^{-\beta(E_{\vec{p}}-\mu)}}\ \Longrightarrow\
N_{B}=\frac{1}{e^{\beta(E_{\vec{p}}-\mu)}-1}\,,$ (4.3) fermions:
$\displaystyle Z_{F}=\frac{e^{\beta
E_{\vec{p}}/2}}{1+e^{-\beta(E_{\vec{p}}-\mu)}}\ \Longrightarrow\
N_{B}=\frac{1}{e^{\beta(E_{\vec{p}}-\mu)}+1}\,.$ (4.4)

From now on we switch to the notation $\omega_{\vec{p}}\rightarrow
E_{\vec{p}}$ such that the dispersion relation reads
$E^{2}_{\vec{p}}=\vec{p}^{\,2}+m^{2}$. The particle number densities are

$n_{B}=\frac{N_{B}}{\mathcal{V}}=\int\frac{\differential^{3}p}{(2\pi)^{3}}\frac{1}{e^{\beta(E_{\vec{p}}-\mu)}-1},\
\ \
n_{F}=\int\frac{\differential^{3}p}{(2\pi)^{3}}\frac{1}{e^{\beta(E_{\vec{p}}-\mu)}+1}\,.$
(4.5)

Hence, let us introduce the $1$-particle distribution function as

$f_{\eta}(\vec{p})=\frac{1}{e^{\beta(E_{\vec{p}}-\mu)}-\eta},$ (4.6)

where $\eta=+1$ for bosonic fields and $\eta=-1$ for fermionic fields. The
$\mu=0$ case will be the relevant case for us; although the total particle
density is high in the early universe, the net particle number densities are
very small.

Our aim is to allow small departures from equilibrium, which can be described
by local changes in the distribution function, such that it becomes space and
time dependent. We will then obtain the evolution equation for the
distribution function $f(p,x)\,\differential[3]p\,\differential[3]x$ where
$f(p,x)$ describes the average number of particles of momentum $\vec{p}$ in a
$3$-phase space volume element
$(\vec{x},\vec{x}+\differential\vec{x})\times(\vec{p},\vec{p}+\differential\vec{p})$
at time $t=x^{0}$. From this function, we can define various quantities like
the number density $n(x)$ and the particle flux $j^{i}(x)$ as

$\displaystyle n(x)$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{3}p}{(2\pi)^{3}}\,f(p,x)\,,$ (4.7)
$\displaystyle j^{i}(x)$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{3}p}{(2\pi)^{3}}\,\frac{p^{i}}{p^{0}}f(p,x)\,.$
(4.8)

We can combine both quantities and define the particle current

$j^{\mu}(x)=\int\frac{\differential^{3}p}{(2\pi)^{3}}\,\frac{p^{\mu}}{p^{0}}f(p,x)\,.$
(4.9)

Furthermore, we define the energy density $e(x)$, the $3$-momentum density
$\Pi^{i}(x)$ and the $3$-momentum flux $\Pi^{ij}(x)$ as

$\displaystyle e(x)$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{3}p}{(2\pi)^{3}}\,p^{0}f(p,x)\,,$
(4.10) $\displaystyle\Pi^{i}(x)$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{3}p}{(2\pi)^{3}}\,p^{i}f(p,x)\,,$
(4.11) $\displaystyle\Pi^{ij}(x)$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{3}p}{(2\pi)^{3}}\,p^{i}\frac{p^{j}}{p^{0}}f(p,x)\,.$
(4.12)

We can combine the last three quantities to form the energy-momentum tensor,

$T^{\mu\nu}(x)=\int\frac{\differential^{3}p}{(2\pi)^{3}}\,\,\frac{p^{\mu}p^{\nu}}{p^{0}}f(p,x),$
(4.13)

where $T^{00}=e$, $T^{0i}=\Pi^{i}=T^{i0}$, and $T^{ij}=\Pi^{ij}$. The integral
measure transforms as a scalar under Lorentz transformations because we can
rewrite it as

$\int\frac{\differential[3]p}{(2\pi)^{3}}\frac{1}{p^{0}}=\int\frac{\differential[4]p}{(2\pi)^{4}}\delta(p^{2}+m^{2})\theta(p^{0})\,,$
(4.14)

where $\theta(p^{0})$ ensures positivity of the energy. Hence, for the energy-
momentum tensor $T^{\mu\nu}$ to transform as a $2$-tensor under Lorentz-
transformation (and the particle current $j^{\mu}$ as a $1$-tensor), the
distribution function $f(p,x)$ has to be a Lorentz-scalar.

###  4.2 Relativistic Boltzmann equation

Let us study how the particle distribution function changes in time. First,
let us assume that there are no collisions between the individual particles,
and follow the trajectory of each particle in phase space
$(x^{\mu}(\tau),p^{\mu}(\tau))$, which is parametrised by proper time $\tau$.
The position and momentum after a small proper time interval
$\differential\tau$ hence change as

$\displaystyle x^{\mu}(\tau)$ $\displaystyle\longrightarrow$ $\displaystyle
x^{\mu}(\tau)+\derivative{x^{\mu}}{\tau}\differential\tau=x^{\mu}(\tau)+\frac{p^{\mu}}{m}\differential\tau\,,$
(4.15) $\displaystyle p^{\mu}(\tau)$ $\displaystyle\longrightarrow$
$\displaystyle p^{\mu}(\tau)+F^{\mu}\differential\tau\,,$ (4.16)

where $F^{\mu}$ describes an external $4$-force. This is sketched in Fig. 6,
which shows the phase space at two time steps $\tau$ and
$\tau+\differential\tau$. The particle distribution functions at time $\tau$
and $\tau+\differential\tau$ must be the same because the particle number is
conserved in an infinitesimal phase space volume. This leads to

$f\left(p+F\differential\tau,x+\frac{p}{m}\differential\tau\right)=f(p,x).$
(4.17)

A Taylor expansion around $\differential\tau=0$ yields

$\left(p^{\mu}\partial_{\mu}+mF^{\mu}\partialderivative{p^{\mu}}\right)f(p,x)=0$
(4.18)

A consistent differential equation for the distribution function needs to
maintain the on-shell condition $p^{2}+m^{2}=0$. This is the case when the
force satisfies either (a) $F^{\mu}p_{\mu}=0$ or (b)
$F^{\mu}=-\partial^{\mu}m$. We introduce the on-shell condition by hand to
arrive at the collisionless relativistic Boltzmann equation

$\left(p^{\mu}\partial_{\mu}+mF^{\mu}\partialderivative{p^{\mu}}\right)\delta(p^{2}+m^{2})f(p,x)=0.$
(4.19)

$x$$p$$\differential x$$\differential x^{\prime}$$\differential
p$$\differential p^{\prime}$ Figure 6: This figure depicts how the phase space
volume occupied by particles within $(x,x+\differential x)$ and
$(p,p+\differential p)$ changes in a time step $\differential\tau$ to
$(x^{\prime},x^{\prime}+\differential x^{\prime})$ and
$(p^{\prime},p^{\prime}+\differential p^{\prime})$. The particle number in
both elements is the same.

Next, we want to incorporate collisions between the particles. We focus on
$2$-body collisions between classical particles only, as these dominate the
scattering in the energy regime that we are interested in. Furthermore, we
assume that there are no external forces, $F^{\mu}=0$. Let us denote initial
values, i.e., before the collision, without a prime and final values with a
prime (see Figure 7). Momentum conservation requires

$p_{1}+p_{2}=p_{1}^{\prime}+p_{2}^{\prime}\,,$ (4.20)

where the subscript refers to the particle. In the presence of collisions, the
particle distribution function is not necessarily the same after a time step
$\differential\tau$ because scattering can remove and add particles to the
phase space volume element $(\differential p,\differential x)$ around $(p,x)$.
Let $R$ ($\bar{R}$) describe a scattering in time $\differential t$ which
removes (adds) an initial (final) particle with momentum $p$ at position $x$.
The Boltzmann equation including collisions reads

$p^{\mu}\partial_{\mu}f(p,x)=\bar{R}(p,x)-R(p,x)\,.$ (4.21)

Let us quantify $R$ and $\bar{R}$ a bit further. An incoming particle with
momentum $p_{1}$ can scatter with any particle that is at the same position,
but with arbitrary momentum $p_{2}$ (and likewise for outgoing). That amounts
to writing

$\displaystyle R(p_{1},x)$ $\displaystyle=$
$\displaystyle\int\frac{\bar{\differential}^{3}p_{2}}{2E_{2}}\frac{\bar{\differential}^{3}p_{1}^{\prime}}{2E_{1}^{\prime}}\frac{\bar{\differential}^{3}p_{2}^{\prime}}{2E_{2}^{\prime}}f(p_{1},x)f(p_{2},x)W(p_{1},p_{2}|p_{1}^{\prime},p_{2}^{\prime})\,$
(4.22) $\displaystyle\bar{R}(p_{1},x)$ $\displaystyle=$
$\displaystyle\int\frac{\bar{\differential}^{3}p_{2}}{2E_{2}}\frac{\bar{\differential}^{3}p_{1}^{\prime}}{2E_{1}^{\prime}}\frac{\bar{\differential}^{3}p_{2}^{\prime}}{2E_{2}^{\prime}}f(p_{1}^{\prime},x)f(p_{2}^{\prime},x)W(p_{1}^{\prime},p_{2}^{\prime}|p_{1},p_{2})\,,$
(4.23)

where $W$ is called the scattering function. Here, we introduced the short-
hand notation $\bar{\differential}^{3}p=\differential^{3}p/(2\pi)^{3}$. We can
write $W$ in terms of the cross section $\sigma$ as

$W(p_{1},p_{2}|p_{1}^{\prime},p_{2}^{\prime})=s\
\sigma(s,\theta)\delta^{(4)}(p_{1}^{\prime}+p_{2}^{\prime}-p_{1}-p_{2})$
(4.24)

in terms of the Mandelstam variables $s=(p_{1}+p_{2})^{2}$,
$t=(p_{1}-p_{1}^{\prime})^{2}$, and the scattering angle $\theta$ is obtained
from $\cos\theta=1-2t/(s+4m^{2})$.

Let us introduce the collision function $C[f]=\bar{R}-R$, where the notation
reminds us that $R$ and $\bar{R}$ depend on the distribution function. We now
have derived the relativistic Boltzmann equation

$p^{\mu}\partial_{\mu}f=C[f]\,.$ (4.25)

$p_{2}^{\prime}$$p_{1}$$p_{1}^{\prime}$$p_{2}$ Figure 7: Visual depiction of a
$2$-body scattering event. The unprimed quantities describe the values before
the scattering, the primed ones after the scattering.

###  4.3 Conservation laws and collision invariants

It is natural to expect that conservation laws obeyed during the collisions
constrain the evolution of the distribution function. In our above system with
2-body collisions, we expect that the particle number is conserved, as well as
the momentum. One can show that for any function
$\psi(p,x)=a(x)+b_{\mu}(x)\,p^{\mu}$ with $a$ and $b_{\mu}$ arbitrary
functions of $x$, the following integral vanishes identically,

$\int\frac{\bar{\differential}^{3}p}{2E}\psi(p,x)C[f]=0\,,$ (4.26)

as a consequence of particle number and momentum conservation. The function
$\psi$ is a collision invariant.

Now we replace the collision function using the Boltzmann equation (4.25). For
the case $b_{\mu}=0$, the above integral then implies

$\displaystyle
0=\int\frac{\bar{\differential}^{3}p}{2E_{\vec{p}}}p^{\mu}\partial_{\mu}f=\partial_{\mu}j^{\mu}\,,$
(4.27)

where we took the partial derivative out of the integral and used Eq. (4.9).
Therefore, the particle current is conserved.

Instead of setting $b_{\mu}=0$, we can set $a=0$ to arrive at

$\displaystyle
0=\int\frac{\bar{\differential}^{3}p}{2E_{\vec{p}}}p^{\nu}p^{\mu}\partial_{\mu}f=\frac{1}{2}\partial_{\mu}T^{\mu\nu}$
(4.28)

where we took the partial derivative out of the integral and used the
definition in Eq. (4.13). This establishes conservation of the tensor,

$T^{\mu\nu}=\int\frac{\bar{\differential}^{3}p}{2E_{\vec{p}}}p^{\nu}p^{\mu}f,$
(4.29)

which is the energy-momentum tensor of the system of particles with
distribution function $f$.

###  4.4 Local equilibrium and perfect fluid

The fluid is in local equilibrium in the presence of collisions when
$R(p_{1},x)=\bar{R}(p_{1},x)$, i.e. $C[f]=0$. From Eqs. (4.22) and (4.23) we
find that the fluid is in local equilibrium when,
$\forall\,p_{1},p_{2},p_{1}^{\prime},p_{2}^{\prime}$,

$\displaystyle f_{1}f_{2}=f^{\prime}_{1}f^{\prime}_{2}$ (4.30)
$\displaystyle\Rightarrow$ $\displaystyle\ln f_{1}+\ln f_{2}=\ln
f^{\prime}_{1}+\ln f^{\prime}_{2}$ (4.31)

where we introduced the short-hand notation
$f_{i}^{(\prime)}=f(p_{i}^{(\prime)},x)$. This implies that the quantity $\ln
f_{1}+\ln f_{2}$ is conserved. Therefore it must be expressible in terms of
the conserved quantities of the system. According to Eq. (4.26), it can
therefore be written as

$\ln f^{\rm eq}(p,x)=a(x)+b_{\mu}(x)p^{\mu}\,,$ (4.32)

in local equilibrium. We rewrite $a(x)$ and $b_{\mu}(x)$ in a suggestive
notation, $a(x)=\beta(x)\mu(x)$ and $b_{\mu}(x)=\beta(x)u_{\mu}(x)$, with
$u^{2}=-1$, and we recover the classical equilibrium distribution function

$f^{\rm eq}(p,x)=\exp[\beta(x)(p\cdot u(x)+\mu(x))\Big{]}.$ (4.33)

We then see that $\beta(x)=1/T$ is the inverse temperature, $\mu(x)$ the
chemical potential, and $u^{\mu}(x)$ is the local $4$-velocity of the system
of particles which we can now start calling a fluid. In the fluid local rest
frame, $u^{\mu}=(1,\vec{0})^{\mu}$, $p\cdot u=-p^{0}$, where
$p^{0}=E_{\vec{p}}$ is the particle energy.

To extend this analysis to quantum scattering, we need to account for the fact
that two fermions cannot occupy the same quantum state (Fermi blocking), and
that bosonic wave functions add coherently (Bose enhancement). This can be
done by writing the particle production and destruction rates as

$\displaystyle R(p_{1},x)$ $\displaystyle=$
$\displaystyle\int\frac{\bar{\differential}^{3}p_{2}}{2E_{2}}\frac{\bar{\differential}^{3}p_{1}^{\prime}}{2E_{1}^{\prime}}\frac{\bar{\differential}^{3}p_{2}^{\prime}}{2E_{2}^{\prime}}f_{1}f_{2}(1\pm
f_{1}^{\prime})(1\pm f_{2}^{\prime})\overrightarrow{W}\,,$ (4.34)
$\displaystyle\bar{R}(p_{1},x)$ $\displaystyle=$
$\displaystyle\int\frac{\bar{\differential}^{3}p_{2}}{2E_{2}}\frac{\bar{\differential}^{3}p_{1}^{\prime}}{2E_{1}^{\prime}}\frac{\bar{\differential}^{3}p_{2}^{\prime}}{2E_{2}^{\prime}}f_{1}^{\prime}f_{2}^{\prime}(1\pm
f_{1})(1\pm f_{2})\overleftarrow{W}\,.$ (4.35)

The additional factors of $1\pm f_{i}^{(\prime)}$ implement the Bose
enhancement and Fermi blocking in $C[f]$.

In local equilibrium, the particle production and destruction rates are equal,
i.e. $C[f]=0$. Hence, the distribution function has to satisfy

$f_{1}f_{2}(1\pm f_{1}^{\prime})(1\pm
f_{2}^{\prime})=f_{1}^{\prime}f_{2}^{\prime}(1\pm f_{1})(1\pm f_{2})\,.$
(4.36)

Separating primed and unprimed variables reveals that $\ln f_{1}/(1\pm
f_{1})+\ln f_{2}/(1\pm f_{2})$ is conserved. As for the classical case, this
implies that the distribution function in equilibrium can be written as

$\ln\frac{f^{\rm eq}}{1\pm f^{\rm eq}}=a+b_{\mu}p^{\mu}\,,$ (4.37)

where it is understood that as before $a$ and $b$ depend on $x$. Rearranging
and rewriting $a$ and $b$ in terms of the more physical quantities $\mu$,
$\beta$ and $u$ yields the following expression for the distribution function,

$f^{\rm eq}=\frac{1}{e^{\beta(\mu+u\cdot p)}\pm 1}.$ (4.38)

In the early universe the chemical potential is negligible, $\mu\simeq 0$.
Then, the energy-momentum tensor in local equilibrium reads

$T^{\mu\nu}=(e+p)u^{\mu}u^{\nu}+pg^{\mu\nu},$ (4.39)

where all quantities depend on $x$. This is the energy-momentum tensor for a
perfect fluid with energy density $e$ and pressure333The pressure is defined
as
$p=\int\frac{\differential[3]p}{(2\pi)^{3}}\frac{\vec{p}^{\,2}}{2p^{0}}f(p,x)\,.$
(4.40) $p$. Despite an ambiguity in the notation, the pressure is not to be
confused with 4-momentum.

##  5 Hydrodynamics with field-dependent mass

In the early universe the value of the Higgs changes in time during the
electroweak phase transition. The Higgs transitions from the symmetric to the
broken phase. As discussed in Sec. 3, the particle masses depend on the value
of the Higgs. Since the transition does not happen at the same time in the
entire universe, the particle masses are space-time dependent due to their
field dependency. In this section, we generalize the results of the previous
section to the case of a fluid of particles with space-time dependent mass.

The theory of hydrodynamics with a field dependent mass is structurally
similar to hydrodynamics with electromagnetic forces and hence to the wide
field of plasma physics. This is because the field dependency of the mass
leads to a non-zero external force $F^{\mu}$. In the case of plasma physics
the motion of particles under Lorentz forces is described by the Vlasov
equation. In close analogy, we will derive Boltzmann’s equation for a fluid
with field dependent (and hence space-time dependent) mass in the following.

We start from the action for a single particle (note this is not the action
for the field $\phi$)

$S=-\int\differential{\tau}m\sqrt{-g_{\mu\nu}\derivative{x^{\mu}}{\tau}\derivative{x^{\nu}}{\tau}}\,,$
(5.1)

where $m=m(\phi(x(\tau)))$ is the field-dependent mass,
$x^{\mu}=x^{\mu}(\tau)$ is the space-time coordinate of a particle,
parameterised by the particle’s proper time $\tau$, and $g_{\mu\nu}$ is the
space-time metric. Varying this action we obtain the equation of motion

$\derivative{\tau}\left(m\derivative{x^{\mu}}{\tau}\right)+\partial^{\mu}\phi\derivative{m}{\phi}=0\,.$
(5.2)

From this equation we identify the 4-momentum
$p^{\mu}=m\differential{x^{\mu}}/\differential{\tau}$, and the force

$F^{\mu}=-\partial^{\mu}\phi\derivative{m}{\phi}\,,$ (5.3)

where we used the equation of motion (5.2). Hence, the field-dependence of the
mass yields a force acting on the particle. This has to be taken into account
when deriving conservation laws.

In the previous section, we derived the conservation of the particle current
and energy-momentum assuming there are no external sources. Now, consider
instead the Boltzmann equation with collisions and external forces,

$\left(p^{\mu}\partial_{\mu}+mF^{\mu}\partialderivative{p^{\mu}}\right)\Theta(p^{0})\delta(p^{2}+m^{2})f=C\left[f\right]\,,$
(5.4)

where we have introduced the on-shell condition again. We multiply both sides
with $p^{\nu}$, integrate over momenta, and use Eq. (4.26) to find

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int\frac{\differential^{4}p}{(2\pi)^{4}}p^{\nu}C[f]$ (5.5)
$\displaystyle=$
$\displaystyle\int\frac{\differential^{4}p}{(2\pi)^{4}}p^{\nu}\left(p^{\mu}\partial_{\mu}+mF^{\mu}\partialderivative{p^{\mu}}\right)\Theta(p^{0})\delta(p^{2}+m^{2})f$
(5.6) $\displaystyle=$
$\displaystyle\partial_{\mu}\int\frac{\differential^{3}p}{(2\pi)^{3}}\frac{p^{\mu}p^{\nu}}{2E_{\vec{p}}}f-mF^{\nu}\int\frac{\differential^{4}p}{(2\pi)^{4}}\Theta(p^{0})\delta(p^{2}+m^{2})f$
(5.7) $\displaystyle=$
$\displaystyle\frac{1}{2}\partial_{\mu}T^{\mu\nu}-mF^{\nu}\int\frac{\differential^{3}p}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}f\Bigg{|}_{p^{0}=E_{\vec{p}}}\,,$
(5.8)

where we assumed that collisions are local and preserve momenta. For the last
step, we used the definition of the energy-momentum tensor as in Eq. (4.13),
integrated the second term by parts and used that $\partial p^{\nu}/\partial
p^{\mu}=\delta^{\nu}_{\mu}$. Using Eq. (5.3), this yields

$\partial_{\mu}T^{\mu\nu}=-\partial^{\nu}\phi\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}f(p,x)\Bigg{|}_{p^{0}=E_{\vec{p}}}$
(5.9)

where the 4-momentum is constrained such that $p^{0}=E_{\vec{p}}$. The energy-
momentum tensor of the fluid is not conserved when the mass is field-
dependent, but sourced by a term proportional to the change of the mass
throughout spacetime. Total energy-momentum of the fluid-field system is
conserved, so there is a corresponding term in the energy-momentum
conservation equation for the field.

Note that in the derivation of Eq. (5.9) we assumed momentum conservation in
collisions. However, in the presence of gradients in the scalar field it is
not clear if momentum is conserved; there could be exchange of momentum with
the field during the collision. However, as long as the gradients in the
scalar field are small, by which we mean

$\lambda_{\text{mfp}}\frac{\partial\phi}{\phi}\lesssim\mathcal{O}(1)\,,$
(5.10)

where $\lambda_{\text{mfp}}$ is the mean free path in the fluid, conservation
of momentum in collisions should be a good approximation. In the opposite
limit, the particles are unlikely to interact in the wall, and conservation of
momentum in collisions is again recovered.

###  5.1 Complete model of scalar field fluid system

A more detailed description of a system consisting of a scalar field and a
fluid is given in Ref. [45]. As a starting point, split the full energy-
momentum tensor into two components, one for the fluid $T^{\mu\nu}_{\text{f}}$
and one for the Higgs $T^{\mu\nu}_{\phi}$, as

$\displaystyle T_{\text{f}}^{\mu\nu}$ $\displaystyle=$
$\displaystyle(e+p_{1})u^{\mu}u^{\nu}+p_{1}g^{\mu\nu}\,,$ (5.11)
$\displaystyle T_{\phi}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\partial^{\mu}\phi\partial^{\nu}\phi-g^{\mu\nu}\left(\frac{1}{2}(\partial\phi)^{2}+V_{0}(\phi)\right)\,,$
(5.12)

where we label the fluid pressure as $p_{1}$ for now. As an example we use the
symmetry breaking potential

$V_{0}(\phi)=\frac{\lambda}{4}\left(\phi^{2}-v^{2}\right)^{2}\,,$ (5.13)

for the scalar field. Note that there are some field-dependent terms in the
fluid pressure, which we define as

$p_{1}(\phi,T)=\frac{\pi^{2}}{90}g_{\text{eff}}T^{4}-V_{1}(\phi,T)\,,$ (5.14)

where the fluid pressure is the total contribution from all particles, equal
to the negative of the free energy densities (2.19) and (2.26)

$p_{1}(\phi,T)=-\sum_{B}f_{B}(m(\phi),T)-\sum_{F}f_{F}(m(\phi),T).$ (5.15)

The full thermally corrected potential then reads

$V_{T}(\phi)=V_{0}(\phi)+V_{1}(\phi,T)\,.$ (5.16)

In the high-temperature expansion, $V_{1}(\phi,T)$ contains all $T$-dependent
terms of the potential in Eq. (3.7).

We demand that the full energy-momentum tensor is conserved,

$\partial_{\mu}\left(T_{\text{f}}^{\mu\nu}+T_{\phi}^{\mu\nu}\right)=0\,,$
(5.17)

and hence non-conservation of $T^{\mu\nu}_{\text{f}}$ has to appear in
$T^{\mu\nu}_{\phi}$. Using Eq. (5.9) this implies

$\partial_{\mu}T_{\phi}^{\mu\nu}=+\partial^{\nu}\phi\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}f(p,x)\,.$
(5.18)

The energy-momentum tensor for the field is hence not conserved individually,
fluid contributions can source changes in its energy and momentum.

Computing the equation of motion for the scalar field one obtains

$\Box\phi-
V_{0}^{\prime}(\phi)=-\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}f(p,x)$
(5.19)

where prime denotes the derivative with respect to the argument $\phi$. Here
again the right-hand side stems from the fluid contributions and in a more
complicated theory should contain a sum over all massive degrees of freedom.

Let us study the situation, when the system is close to equilibrium and
parametrize the distribution function as a small perturbation around
equilibrium as $f(p,x)=f^{\text{eq}}(p,x)+\delta f(p,x)$. If the fluid were in
local equilibrium everywhere, Eq. (5.19) becomes [45]

$\Box\phi-
V_{0}^{\prime}(\phi)=-\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}f^{\text{eq}}(p,x)=V^{\prime}_{1}(\phi,T)$
(5.20)

according to the definition of the free energies in Eqs. (2.18) and (2.26).
Therefore, the equation of motion of the Higgs in exact equilibrium reads,

$\Box\phi-V_{T}^{\prime}(\phi)=0\,.$ (5.21)

Consequently the equation for small departures from fluid equilibrium (5.19)
becomes

$\Box\phi-
V_{T}^{\prime}(\phi)=-\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}\delta
f(p,x)$ (5.22)

to first order in $\delta f$.

Let us repackage the contribution from the effective potential into the fluid
by defining the fluid pressure as $p(\phi,T)=p_{1}(\phi,T)-V_{0}(\phi)$. The
new fluid energy-momentum tensor then reads

$T_{\text{f}}^{\mu\nu}=(e+p)u^{\mu}u^{\nu}+pg^{\mu\nu}\,$ (5.23)

which yields

$\partial_{\mu}T_{\text{f}}^{\mu\nu}+V^{\prime}_{T}(\phi)\partial^{\nu}\phi=-\partial^{\nu}\phi\derivative{m^{2}}{\phi}\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}\delta
f(p,x)$ (5.24)

for its conservation equation. The right hand side is generated by departures
from the equilibrium phase-space distribution $f^{\text{eq}}$, which in turn
have to be induced by gradients of the scalar field. We assume that these are
small, in the sense outlined at the end of the last section.

To express the dependence on the gradient in analogy to linear response
theory, we write the integral of the perturbed distribution function on the
right hand side of Eq. (5.24) as a term proportional to the gradient of the
scalar field,

$\int\frac{\differential[3]{p}}{(2\pi)^{3}}\frac{1}{2E_{\vec{p}}}\delta
f(p,x)=\tilde{\eta}u^{\mu}\partial_{\mu}\phi\,,$ (5.25)

where we use the fluid $4$-velocity $u^{\mu}$ to contract the open index in
order to respect isotropy. The proportionality factor $\tilde{\eta}$ might in
general depend on the field value, the temperature, and the fluid
$\gamma$-factor, or $\tilde{\eta}=\tilde{\eta}(T,\phi,\gamma)$. The overall
mass dimension of $\tilde{\eta}$ is zero, so it must depend on the ratio
$\phi/T$.

This leaves us with the equation of motion

$\Box\phi-
V_{T}^{\prime}(\phi)=-\tilde{\eta}\derivative{m^{2}}{\phi}u^{\mu}\partial_{\mu}\phi$
(5.26)

Note that in using the ansatz (5.25) we assumed full Lorentz symmetry.
However, invariance under boosts is broken in a plasma. The following argument
from entropic considerations supports the Lorentz-symmetric form.

Consider the scalar product of $u$ with both sides of Eq. (5.24),

$u_{\nu}\partial_{\mu}(wu^{\mu}u^{\nu}+pg^{\mu\nu})+V^{\prime}_{T}(\phi)u\cdot\partial\phi=\tilde{\eta}(u\cdot\partial\phi)^{2},$
(5.27)

where $w=e+p=Ts$ is the enthalpy density, and $s=\differential p/\differential
T$ is the entropy density. From this we can derive a conservation equation for
the entropy current $S^{\mu}=su^{\mu}$,

$\partial\cdot S=\frac{\tilde{\eta}}{T}(u\cdot\partial\phi)^{2}.$ (5.28)

The fact that this is always positive supports the claim that the term on the
right hand side of Eq. (5.25) is indeed proportional to $u\cdot\partial\phi$.
In this model entropy is generated by the interaction of the field and the
fluid.

##  6 The bubble nucleation rate

We now to study how, and how quickly, a phase transition happens. We focus on
first-order phase transitions, which are characterised by a critical
temperature, a latent heat, and mixed phases separated by phase boundaries
with a characteristic surface energy. An order parameter distinguishes between
the phases.

A simple example is the phase transition between water in its liquid and
vapour form, at a critical temperature of 373 K at 1 atmosphere pressure. The
order parameter is the density, which changes by a factor of about 1000 at the
phase transition.

In the laboratory, systems cooled below their critical temperature, proceed by
the nucleation of small bubbles or droplets of the new phase around
impurities. Especially pure systems can be cooled well below their critical
temperature before the transition occurs.

The early Universe is thought to be perfectly pure, so bubbles of the new
phase can appear only by thermal or quantum fluctuations.444This picture may
not be accurate if there are primordial black holes, which can act as
nucleation sites [46, 47].  The order parameter is a scalar field $\phi$, and
we will take the values of the field in the two phases to be $\phi=0$ in the
high temperature phase and $\phi=\phi_{\rm b}$ in the low temperature phase.

The transition proceeds via bubble nucleation: individual bubbles of the new
ground state $\phi=\phi_{\rm b}$ appear due to thermal fluctuations. If the
temperature of the universe is below the critical value $T_{\rm c}$, these
bubbles grow if they exceed a critical radius $R_{\rm c}$. Eventually they
merge, and the whole universe is in the new phase. In this section we will
compute the rate at which these bubbles appear and how the system changes from
the old to the new ground state. This process was first studied in the context
of statistical mechanics in Ref. [48], and in quantum field theory in Refs.
[49] and [50]. In turns out that at high temperatures higher than the masses
of the particles involved, the quantum field theory is well approximated by
statistical mechanical description.

We start with the partition function of a classical field $\phi$ at
temperature $\beta=1/T$, which is the functional integral

$Z_{\beta}=\int\mathcal{D}{\pi}\mathcal{D}{\phi}e^{-\beta H}\,,$ (6.1)

with the Hamiltonian

$H=\int\differential[3]{x}\Bigg{[}\frac{1}{2}\pi^{2}+\frac{1}{2}(\vec{\nabla}\phi)^{2}+V_{T}(\phi)\Bigg{]}\,,$
(6.2)

where $\pi$ is the conjugate momentum to the field $\phi$. We will use the
partition function to calculate the rate per unit volume of fluctuations over
a thermal barrier in the potential $V_{T}$.

###  6.1 Transition rate in $0$ dimensions

To get started, we discuss a simple toy model without spatial dimensions [51].
Therefore, spatial gradients $\vec{\nabla}\phi$ do not exist, and the
integrals over functions $\pi$ and $\phi$ become ordinary integrals. The
theory describes a particle with coordinate $\phi$ in a potential $V_{T}$ as
depicted in Fig. 8. The partition function for this simple system is

$Z_{\beta}=\int\differential{\pi}\differential{\phi}e^{-\beta\big{[}\pi^{2}/2+V_{T}(\phi)\big{]}}\,.$
(6.3)

We can immediately perform the integral over $\pi$ yielding

$Z_{\beta}=\sqrt{\frac{2\pi}{\beta}}\int\differential{\phi}e^{-\beta
V_{T}(\phi)}\,.$ (6.4)

To carry out the remaining integration over $\phi$ we use the saddle point
approximation. In principle the potential $V_{T}(\phi)$ has three stationary
points, two minima, a local one at $\phi=0$ and a global one at
$\phi=\phi_{\text{b}}$ as well as a maximum at $\phi=\phi_{\text{m}}$
separating the two minima (see Fig. 8). However, in our set-up we assume that
initially there is a thermal bath of particles only in the metastable ground
state at $\phi=0$. To infer the transition rate of particles fluctuating to
the stable ground state at $\phi=\phi_{\rm b}$, only two saddle-points are of
relevance: the false minimum and the maximum. Intuitively, this can be
understood from the fact that once a particle in the metastable ground state
makes it over the top of the potential barrier, due to thermal fluctuations,
it just rolls down into the true ground state. Hence, the transition rate is
dominated by the two aforementioned saddle-points. In particular, even though
with time particles settle in the global minimum, the flux of particles from
the stable to the metastable ground state will always be negligible.

Let us evaluate the partition function at the two saddle points. The first
point sits at $\phi=0$. We hence expand as $\phi=0+\delta\phi$ to obtain

$\displaystyle Z^{0}_{\beta}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2\pi}{\beta}}\int\differential{\phi}\exp[-\beta\left(V_{T}(0)+\frac{1}{2}V_{T}^{\prime\prime}(0)\delta\phi^{2}\right)\Bigg{]}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{2\pi}{\beta}}\sqrt{\frac{2\pi}{\beta
V_{T}^{\prime\prime}(0)}}e^{-\beta V_{T}(0)}\,.$

The second stationary point is at $\phi=\phi_{\text{m}}$ with
$V_{T}^{\prime\prime}(\phi_{\text{m}})<0$. As the second derivative at the
stationary point is negative here we need to deform the integration contour
into the complex plane, which results in the integral acquiring an imaginary
part. At this second stationary point we obtain, for the imaginary part,

$Z_{\beta}^{1}=\sqrt{\frac{2\pi}{\beta}}\sqrt{\frac{2\pi}{V_{T}^{\prime\prime}(\phi_{\text{m}})}}\frac{1}{2}e^{-\beta
V_{T}(\phi_{\text{m}})}\,.$ (6.6)

This expression is imaginary because of the negative sign for
$V_{T}^{\prime\prime}$. The factor of $1/2$ is due to the integration contour
transversing only one half of the complex plane. The leading term in the full
partition function is the sum of these two terms,
$Z_{\beta}=Z_{\beta}^{0}+Z_{\beta}^{1}$.

Computing the free energy we obtain

$F=-\frac{1}{\beta}\ln Z\approx
F^{0}_{\beta}\pm\frac{i}{2\beta}e^{-\beta\Delta
V_{T}}\sqrt{\frac{V_{T}^{\prime\prime}(0)}{\absolutevalue{V_{T}^{\prime\prime}(\phi_{\text{m}})}}}\,,$
(6.7)

where $\Delta V_{T}\equiv V_{T}(\phi_{\text{m}})-V_{T}(0)$ is the barrier
height in the potential, and we expanded the logarithm for $\beta\Delta
V_{T}\gg 1$. The two signs arise due to the choice of the integration contour.

Figure 8: The effective thermal potential $V_{T}(\phi)$ for a temperature
below the critical temperature, $T<T_{\rm c}$. The potential has two minima
$\phi=0$ an $\phi=\phi_{\rm b}$, and a maximum at $\phi=\phi_{\rm m}$, which
we use for the saddle point approximation. Due to thermal fluctuations, the
field can transition from the symmetric phase to the Higgs phase.

Next, we compute the the probability flux $\Gamma$ across the potential
barrier given by

$\Gamma=\frac{1}{Z_{\beta}}\int\differential{\pi}\differential{\phi}e^{-\beta
H}\delta(\phi-\phi_{\text{m}})\pi\Theta(\pi)\,.$ (6.8)

The $\delta$-distribution picks the saddle point at the maximum for us. The
Heaviside function only allows for positive velocities (increasing field
values) and hence we only consider particles that flow to the new minimum. The
whole expression describes a flux of phase space volume from the old to the
new ground state. Evaluating Eq. (6.8) in a similar manner as before yields

$\Gamma=\frac{1}{2\pi}\sqrt{V_{T}^{\prime\prime}(0)}e^{-\beta\Delta V_{T}}\,.$
(6.9)

Comparing this equation to the expression for the free energy (6.7)
establishes the relation

$\Gamma=\frac{\beta}{\pi}\sqrt{V_{T}^{\prime\prime}(\phi_{\rm
m})}\,\absolutevalue{\Im{F}}\,.$ (6.10)

The probability flux density is thus proportional to the imaginary part of the
free energy. To be explicit, it is given by

$\Im{F}=\frac{T}{2}\sqrt{\frac{V_{T}^{{}^{\prime\prime}}(0)}{\absolutevalue{V_{T}^{\prime\prime}(\phi_{\text{m}})}}}e^{-\beta\Delta
V_{T}}\,.$ (6.11)

Note that in Eq. (6.8) we implicitly assumed that the particles were moving
freely. This is equivalent to assuming that the particles do not interact with
the heat bath which maintains thermal equilibrium as they cross the barrier.
The timescale for crossing the barrier can be estimated as $\tau_{\rm
m}=1/\sqrt{\absolutevalue{V_{T}^{\prime\prime}(\phi_{\text{m}})}}$. If the
mean free time between interactions is much less than $\tau_{\rm m}$, then it
is more appropriate to study the diffusion of particles across the barrier,
which can be modelled by the Langevin equation

$\dot{\pi}=-\gamma\pi-V^{\prime}_{T}(\phi)+\xi(t).$ (6.12)

Here, $\gamma$ is the diffusion constant, and $\xi(t)$ represents the forces
exerted by the heat bath. It is a Gaussian random variable, obeying
$\langle\xi(t)\xi(t^{\prime})\rangle=2\gamma T\delta(t-t^{\prime})$, meaning
that it is uncorrelated from one instant to the next, but has a mean amplitude
$\sqrt{2\gamma T}$. In the case $\gamma\tau_{\rm m}\gg 1$, the rate of
particles crossing over the barrier is

$\Gamma=\frac{1}{2\pi}\frac{\sqrt{V_{T}^{\prime\prime}(0)\absolutevalue{V_{T}^{\prime\prime}(\phi_{\text{m}})}}}{\gamma}e^{-\beta\Delta
V_{T}}\,.$ (6.13)

This is a factor $1/\gamma\tau_{\rm m}$ smaller than the rate of free
particles crossing the barrier.

The problem of how classical particles escape over an energy barrier was
originally solved by Kramers [52] in the context of chemical reactions. The
formulation of the Kramers escape problem in quantum field theory was recently
studied in Ref. [53]. For an early reference relating the imaginary part of
the free energy to the tunneling rate, see Ref. [54].

###  6.2 Bubble nucleation in $3$ dimensions

After having derived the probability flux without spatial dimensions, in this
section we will repeat this procedure for $3$ dimensions.

####  6.2.1 The critical bubble

The starting point is the path integral which for a thermal field theory is
given by

$\mathcal{Z}_{\beta}=\int\mathcal{D}\pi\int\mathcal{D}\phi\,e^{-\beta
H[\pi,\phi]}\,,$ (6.14)

where the Hamiltonian $H$ is given by Eq. (6.2) now including spatial
gradients of the scalar field. We can perform the integration over $\pi$, as
it is Gaussian, leading to

$\mathcal{Z}_{\beta}=\mathcal{N}\int\mathcal{D}\phi\,e^{-\beta
E_{T}[\phi]}\,,$ (6.15)

where $\mathcal{N}$ is the result of the integration. It will play no role in
the following.

We will again perform a saddle point evaluation of this integral, which means
first identifying the stationary point of the function in the exponential.

At the stationary or saddle point, we have $\delta E_{T}[\phi]/\delta\phi=0$,
leading to the following differential equation

$-\vec{\nabla}^{2}\phi+\frac{\partial V_{T}}{\partial\phi}=0\,.$ (6.16)

It turns out that the energy functional $E_{T}$ is minimized for a radially
symmetric field configuration. In that case, the above equation of motion can
be rewritten as

$-\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}\left(r^{2}\frac{\differential\phi}{\differential
r}\right)+V^{{}^{\prime}}_{T}(\phi)=0\,.$ (6.17)

We will suppose a thermal effective potential of the form sketched in Fig. 8,
and assume a high temperature approximation of quartic form as in Eq. (3.8).
We are looking for a solution which represents a finite-energy fluctuation
away from the metastable phase $\phi=0$ at a temperature just below $T_{\rm
c}$. The asymptotics of the solution are therefore: as $r\rightarrow\infty$,
$\phi\rightarrow 0$, and as $r\rightarrow 0$,
$\phi^{{}^{\prime}}(r)\rightarrow 0$, where the prime denotes derivative with
respect to $r$.

Figure 9:  Left panel: the scaled quartic potential as a function of the
parameter $\bar{\lambda}$, defined in Eq. (6.18). Right panel: the critical
bubble solution to Eq. (6.17). Scripts by kind permission of D. Cutting, see
also Ref. [55].

A trivial solution $\phi(r)=0$ always exists. A non-trivial solutions can be
found by numerical integration, using a shooting method (see e.g. Ref. [55]).
Although there are three parameters, it turns out that the differential
equation can be rewritten so that it depends only on the combination

$\bar{\lambda}=\frac{9\lambda D}{2A^{2}},$ (6.18)

provided one scales the potential with the vacuum energy difference, the field
with its stable phase value $\phi_{\rm b}$, and radial distance with the
effective mass $M=D\sqrt{T^{2}-T_{0}^{2}}$. The results for selected values of
$\bar{\lambda}$ are shown in Fig. 9. We refer to a non-trivial solution of Eq.
(6.17) with these boundary conditions as a critical bubble and denote it by
$\bar{\phi}(r)$.

At the critical temperature, $\bar{\lambda}\to 1$, and the potential is
approximately

$V(\phi)=\frac{\lambda}{4}\phi^{2}\left(\phi-\phi_{\rm b}\right)^{2}$ (6.19)

This leads to the following approximate solution to the equation of motion,

$\bar{\phi}=\frac{\phi_{\rm b}}{2}\left(1-\tanh(\frac{r-R_{\rm c}}{\ell_{\rm
w}})\right)\,,$ (6.20)

where $R_{\rm c}$ is the radius of the critical bubble, $\ell_{\rm
w}=1/M(T_{\rm c})$ the thickness of the bubble wall, and $R_{\rm
c}\gg\ell_{\rm w}$. This solution is plotted in the right panel of Fig. 9 for
different values of $\bar{\lambda}$. We refer to this approximate solution as
a thin-wall bubble.

The energy of a critical bubble can be computed from the energy functional

$E_{T}[\phi]=\int\text{d}^{3}x\left(\frac{1}{2}\left(\vec{\nabla}\phi\right)^{2}+V_{T}(\phi)\right).$
(6.21)

This is potentially divergent in an infinite volume system. We are actually
interested in the difference between the energy of the field configuration of
the critical bubble $\bar{\phi}$ and the energy of the metastable state
$\phi=0$, or

$E_{\text{c}}=E_{T}[\bar{\phi}]-E_{T}[0]\,.$ (6.22)

This quantity is finite, and it is the energy a thermal fluctuation needs to
have in order to lift the field over the barrier in the spherical region of
radius $R$.

It is simplest to analyse this expression for a thin-wall. In that case, the
interior of the bubble represents the ’new’ minimum (for instance the Higgs
phase in the case of electroweak symmetry breaking), the exterior the ’old’
minimum (i.e. the metastable phase) and the wall of the bubble defines the
region where the field interpolates between the two minima.

If $\epsilon=\Delta V_{T}/V_{T}(\phi_{\rm m})\ll 1$, where $\Delta
V_{T}=V_{T}(0)-V_{T}(\phi_{\rm b})$ is the free-energy difference between the
two phases, there exists an approximate solution for $E_{\rm c}$ as a function
of radius $R$, given by

$E_{\rm c}(R)=4\pi\sigma R^{2}-\frac{4}{3}\pi\Delta V_{T}R^{3},$ (6.23)

Here, $\sigma$ denotes the surface tension given by

$\sigma=\frac{2^{3/2}}{3^{4}}\frac{A^{3}}{\lambda^{5/2}}T_{\rm c}^{3}$ (6.24)

with $A$ and $\lambda$ the constants appearing in the thermal Higgs potential
(3.8). It is possible to introduce a critical radius $R_{\rm c}$ at which the
bubble does neither shrink nor grow. This is the radius where the force, due
to the release of latent heat, pushing the bubble to grow, is balanced with
the interaction of the bubble with the plasma of particles creating friction.
Above the critical radius it is favorable for the bubble to grow. The critical
radius is defined as the stationary point of the energy $E$ as a function of
the radius $R$, i.e.

$\frac{\partial E}{\partial R}=8\pi\sigma R-4\pi\Delta
V_{T}R^{2}\Bigg{|}_{R=R_{\rm c}}\overset{!}{=}0\,,$ (6.25)

such that the critical radius is given by

$R_{\rm c}=\frac{2\sigma}{\Delta
V_{T}}=\frac{2\sigma}{V_{T}(\phi_{\text{m}})\epsilon}.$ (6.26)

Plugging this expression into Eq. (6.23), we find the energy of the critical
bubble as

$E_{\rm c}=E(R_{\rm c})=\frac{16\pi}{3}\frac{\sigma^{3}}{\Delta V_{T}^{2}}\,.$
(6.27)

In the thin-wall approximation, the energy of the critical bubble can be
computed as a function of temperature $T$ as [56]

$E_{\rm c}(T)=\frac{16\pi}{3}\frac{\sigma T_{\rm c}^{2}}{[L(T_{\rm
c}-T)]^{2}}\,,$ (6.28)

where $L$ is the latent heat. It is given by $L=T_{\rm c}(p_{\rm s}-p_{\rm
b})$ with $p_{\rm s}$ and $p_{\rm b}$ the fluid pressure in the symmetric and
in the broken phase, resp. Hence, the energy $E_{\rm c}$ diverges at the
critical temperature.

Note the different uses of the word “critical” in critical bubble, with radius
$R_{\rm c}$ and energy $E_{\rm c}$, and the critical temperature $T_{\rm c}$.

####  6.2.2 Saddle point evaluation

We now return to completing the saddle-point approximation of the integral
(6.15). We recall that the saddle point method consists in expanding around
stationary values of the energy functional $E_{T}[\phi]$. In the case at hand
the two previously mentioned classical solutions to the equations of motion,
the trivial solution $\phi(r)=0$ and the critical bubble solution
$\phi(r)=\bar{\phi}(r)$, extremise $E_{T}[\phi]$. In contrast to the
simplified case without spatial dimensions where we evaluated the partition
function at the extrema of the thermal potential, i.e., at the points $\phi=0$
and $\phi=\phi_{\rm m}$, we now evaluate the partition function at two
classical solutions out of the field configuration space of $\phi$. So the
partition function will be a sum of two terms $\mathcal{Z}^{0}$ and
$\mathcal{Z}^{1}$ evaluated around the extrema $\phi=0$ and the critical
bubble $\phi=\bar{\phi}$, respectively.

Let us first expand around the classical solution $\bar{\phi}$ as
$\phi=\bar{\phi}+\delta\phi$. Performing the actual saddle-point approximation
we find

$\mathcal{Z}\simeq\mathcal{N}\mathcal{D}(\delta\phi)e^{-\beta
E[\bar{\phi}]-\beta\int\delta\phi\left(-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(\bar{\phi})\right)\delta\phi/2}.$
(6.29)

To evaluate the functional integral it is useful to diagonalise the operator
$-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(\bar{\phi})$. In order to do so we
consider the following eigenvalue equation

$(-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(\bar{\phi}))\psi_{\alpha}=\lambda_{\alpha}\psi_{\alpha},$
(6.30)

with $\delta\phi=\sum_{\alpha}c_{\alpha}\psi_{\alpha}$, where the
$\psi_{\alpha}$ are normalized so that
$\beta\int\differential^{3}x|\psi_{\alpha}|^{2}=1.$ The measure then takes the
following form

$\mathcal{D}(\delta\phi)=\int\prod_{\alpha}\frac{\text{d}c_{\alpha}}{\sqrt{2\pi}}\,.$
(6.31)

For the non-zero eigenvalues, we have to perform Gaussian integrals, resulting
in

$\mathcal{Z}^{1}\simeq\mathcal{N}^{\prime}\int\prod_{i}\frac{\text{d}c_{i}}{\sqrt{2\pi}}\prod^{\prime}_{\alpha}\frac{1}{{\lambda_{\alpha}}^{\frac{1}{2}}}e^{-\beta
E[\bar{\phi}]}\,,$ (6.32)

where the index $i$ labels the zero eigenvalues, and the prime indicates that
the zero modes have been omitted from the product.

There are in fact three eigenfunctions with zero eigenvalues. Let us recall
the equation of motion for a solution $\bar{\phi}$, given by

$-\vec{\nabla}^{2}\bar{\phi}+V_{T}^{\prime}(\bar{\phi})=0\,.$ (6.33)

Taking the spatial gradient yields

$\displaystyle\left(-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(\bar{\phi})\right)\partial_{i}\bar{\phi}=0\,.$
(6.34)

These are the translational zero modes, as they correspond to the change in
the field when performing a small translation $x^{i}\to x^{i}+\delta x^{i}$,
for which

$\delta\bar{\phi}=\delta x^{i}\partial_{i}\bar{\phi}\,.$ (6.35)

For these modes, the normalised integration variables are

$\differential c^{i}=\differential
x^{i}\left[\frac{\beta}{3}\int\differential^{3}x\left(\vec{\nabla}\bar{\phi}\right)^{2}\right]^{\frac{1}{2}}\,,$
(6.36)

where we have used the spherical symmetry of the solution. One can show from
the stationarity of the solution under a scale transformation
$\delta\bar{\phi}=x^{i}\partial_{i}\bar{\phi}$ that the integral in the square
brackets in Eq. (6.36) is equal to the energy of the critical bubble. Hence

$\differential c^{i}=\differential x^{i}(\beta E_{\text{c}})^{\frac{1}{2}}\,,$
(6.37)

from which it follows that

$\mathcal{Z}^{1}\simeq\mathcal{N}^{\prime}\mathcal{V}\left(\frac{\beta E_{\rm
c}}{2\pi}\right)^{3/2}\prod^{\prime}_{\alpha}\frac{1}{{\lambda_{\alpha}}^{\frac{1}{2}}}e^{-\beta
E[\bar{\phi}]}\,.$ (6.38)

where $\mathcal{V}$ is the volume of the system.

Furthermore, also the trivial solution $\phi=0$ is an extremum. Therefore, if
we solve Eq. (6.15) at $\phi=0$ we obtain

$\mathcal{Z}^{0}\simeq\mathcal{N}^{\prime}\prod_{\alpha}\frac{1}{\left(\lambda_{\alpha}^{0}\right)^{1/2}}e^{-\beta
E[0]}\,.$ (6.39)

The solution $\phi=0$ has no zero modes, as it is translation invariant.

The partition function evaluated through a saddle-point approximation is now
given as a sum of only two terms,
$\mathcal{Z}\simeq\mathcal{Z}^{0}+\mathcal{Z}^{1}$.

In order to compute the transition rate from the trivial solution to the
critical bubble we have to consider the imaginary part of the free energy.
This comes from the expansion around the critical bubble, which has one
negative eigenvalue $\lambda_{-}$, linked to an unstable mode (i.e. expansion
or contraction of the bubble).

We can already see that the ratio $\mathcal{Z}^{1}/\mathcal{Z}^{0}$ contains
an exponentially small factor $\exp(-\beta E_{\text{c}})$. Therefore, we can
expand the logarithm in the free energy and approximate its imaginary part as

$\Im{F}\approx-T\frac{|\mathcal{Z}^{1}|}{\mathcal{Z}^{0}}$ (6.40)

If the interaction rate of the field with the thermal bath is small, we use
the example of the particle in a potential well to infer that the probability
flux per volume $\mathcal{V}$ is then given by

$\displaystyle\frac{\Gamma}{\mathcal{V}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{V_{T}^{\prime\prime}(0)}}{\pi}\left[\frac{\det(-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(0))}{|\det^{\prime}(-\vec{\nabla}^{2}+V_{T}^{\prime\prime}(\bar{\phi}))|}\right]^{1/2}\left(\frac{\beta
E_{\rm c}}{2\pi}\right)^{3/2}e^{-\beta E_{\rm c}}$ (6.41)

where we have introduced a standard notation for the eigenvalue products. The
symbol $\det^{\prime}$ means that the three zero modes have been omitted,
which means that the overall mass dimension of the square root of the ratio of
determinants is that of the square root of the product of three eigenvalues,
or 3. Hence the mass dimension of the right hand side is 4, as is appropriate
for a rate per volume.

We can estimate the order of magnitude of the dimension-4 prefactor as
$\left[\sqrt{V_{T}^{\prime\prime}(0)}\right]^{4}$ in a Standard-Model-like
plasma. Recalling the thermal potential in the high-temperature approximation
(3.8), and using Eq. (3.13), we find that $\sqrt{V_{T}^{\prime\prime}(0)}\sim
g^{2}T$ at $T_{\text{c}}$, where we have used the fact that the most important
couplings in the Standard Model are all of the same order of magnitude, or
$g^{2}\simeq\lambda\simeq y_{t}$. This seems to indicate that the prefactor
should go as $g^{8}T^{4}$. It turns out that for the Higgs in a Standard Model
plasma, the interactions with the $W$ and $Z$ fields are rather rapid, and it
evolves diffusively [57, 58], as in the Kramers problem. The diffusion
constant in a plasma of $W$ and $Z$ particles is $\gamma^{-1}\sim
g^{4}\ln(1/g)$ [59], bringing an extra factor of order $g^{2}\ln(1/g)$.

##  7 Dynamics of expanding bubbles

In this section we will consider what happens after the bubble has nucleated
and started to expand. As the bubble gets larger, it is appropriate to move to
the hydrodynamic description of the system.

We would like to know how the fluid reacts to the expanding bubble, and in
particular the amplitude of the shear stresses which are generated, as they
are the source of the gravitational waves.

A particular quantity of interest is the speed at which the phase boundary is
expanding, the so-called wall speed $v_{\text{w}}$. Having good knowledge of
the wall speed is important, as the flow set-up around the expanding bubble
depends strongly on it, and therefore the gravitational wave power. The flow
also depends on the strength of the transition parametrised by
$\alpha_{\text{n}}$, introduced in Section 7.3. In this section we will see
how the flow around an expanding bubble depends only on $v_{\text{w}}$ and
$\alpha_{\text{n}}$.

The fluid flow is also of importance for baryogenesis scenarios, which can be
explained through a first-order electroweak phase transition [60]. A
relativistic wall speed increases the chances of detecting gravitational
imprints, but decreases the efficiency with which baryon number is generated
[61].

###  7.1 Wall speed

For a bubble to expand, the interior pressure must be larger than the exterior
one. The pressure is given by

$p=\frac{\pi^{2}}{90}g_{\text{eff}}T^{4}-V_{T}(\phi),$ (7.1)

where $g_{\text{eff}}$ are the effective relativistic degrees of freedom. The
bubble expands if $V_{T}(\phi_{\rm b})<V_{T}(0)$, as in this case the internal
pressure exceeds the external one. This can only happen below the critical
temperature.

The bubble wall (which marks the boundary between the phases) must therefore
move outward through the plasma. The dynamics of the field $\phi$ driving the
phase transition interacting with the plasma of particles is given by

$\Box\phi-V^{\prime}_{T}(\phi)=-\eta_{T}(\phi)\,u^{\mu}\partial_{\mu}\phi\,,$
(7.2)

where $\eta_{T}=\tilde{\eta}{\text{d}m^{2}}/{\text{d}\phi}$, cf. Eq. (5.26).
This is just the standard Klein Gordon equation equipped with a friction term
due to the interaction of the bubble with the plasma.

Let us now assume that the wall-speed is constant in the rest frame of the
universe. Furthermore, we assume that the bubbles are macroscopic objects such
that the bubble wall can be approximated to be planar. Let the wall move into
the $z$-direction, and let us assume that it sets up a flow which settles down
to a steady state. In the frame of the wall, the fluid moves with speed
$u^{\mu}(z)$, and the field profile is $\phi=\phi(z)$. Let us also assume that
the fluid speed changes little as it moves across the wall, so that

$u^{\mu}\simeq\gamma_{\rm w}\,(1,0,0,v_{\text{w}})\,.$ (7.3)

We can then rewrite Eq. (7.2) as

$\partial^{2}_{z}\phi-V^{\prime}_{T}(\phi)=-\eta_{T}(\phi)\gamma_{\rm
w}v_{\text{w}}\partial_{z}\phi\,.$ (7.4)

Multiplying both sides of the above equation with $\partial_{z}\phi$ and
integrating over $z$, we obtain

$\int\text{d}z\,\left[\frac{1}{2}\partial_{z}(\partial_{z}\phi)^{2}-\frac{\text{d}V_{T}}{\text{d}z}\right]=-\int\text{d}z\,\eta_{T}(\phi)\gamma_{\rm
w}v_{\text{w}}\left(\partial_{z}\phi\right)^{2}\,.$ (7.5)

After performing the integral on the left hand side of the equation we find

$\Delta V_{T}\simeq\gamma_{\rm
w}v_{\text{w}}\int\text{d}z\,\eta_{T}(\phi,\gamma_{\text{w}})\left(\partial_{z}\phi\right)^{2}\,,$
(7.6)

where $\Delta V_{T}=V_{T}(z=+\infty)-V_{T}(z=-\infty)$ is the pressure
difference across the wall. To obtain the wall speed, one needs to solve this
equation for $v_{\text{w}}$. In practice, it is difficult to obtain a
solution. In some cases a constant-$v_{\text{w}}$ solution may not even exist,
in which case the bubble wall must continue to accelerate. The case where
$v_{\text{w}}\rightarrow 1$ is called runaway solution; the wall speed becomes
ultra-relativistic [62, 63].

However, provided the product $\gamma_{\text{w}}\eta_{T}$ does not decrease
with $\gamma_{\text{w}}$, a constant-$v_{\text{w}}$ solution always exists.
This seems to be the generic expectation for a phase transition in a gauge
theory like the Standard Model [63, 64].

Figure 10: Due to thermal fluctuations, bubbles where the Higgs is in the
stable ground state occur that expand into regions where the Higgs is still in
the metastable ground state. At early times, the bubbles do not overlap, but
later they combine until the entire fluid is in the stable ground state.

In the following discussion, we assume that the transition completes in much
less than a Hubble time, so that we can neglect the expansion of the universe.

In terms of $v_{\text{w}}$ the radius of a bubble at time $t$ that nucleated
at time $t^{\prime}$ is given by $R=v_{\text{w}}(t-t^{\prime})$. Consequently,
the volume $\mathcal{V}$ of this bubble is

$\mathcal{V}(t,t^{\prime})=\frac{4\pi}{3}v_{\text{w}}^{3}(t-t^{\prime})^{3}\,.$
(7.7)

In particular, this allows us to determine the fractional volume of the
universe in the broken phase occupied by all bubbles. Ignoring overlaps and
collisions between bubbles, the fractional volume is given by the volume of a
bubble nucleated at time $t^{\prime}$ multiplied by the number density of
nucleated bubbles that nucleated during $(t^{\prime},t^{\prime}+dt^{\prime})$,
see Fig. 10. The latter is given by $\differential
n(t^{\prime})=\Gamma(t^{\prime})\differential t^{\prime}/{\mathcal{V}}$, such
that the fractional volume in the metastable phase is

$h(t)=1-\int_{t_{\text{c}}}^{t}\frac{4\pi}{3}v_{\text{w}}^{3}(t-t^{\prime})^{3}\frac{\Gamma(t^{\prime})}{\mathcal{V}}\differential
t^{\prime}$ (7.8)

when we do not take into account overlaps. When we include overlaps, the
fractional volume in the metastable phase is given by [56, 65, 66]

$h(t)=\exp[-\int_{t_{\text{c}}}^{t}\text{d}t^{\prime}\,\frac{4\pi}{3}v_{\text{w}}^{3}(t-t^{\prime})^{3}\frac{\Gamma(t^{\prime})}{\mathcal{V}}\Bigg{]}\,.$
(7.9)

At early times, when the exponent is small, we recover the linear law as in
Eq. (7.8). At late times, $h(t)$ tends rapidly to zero as expected.

Next, we consider the behaviour of the bubble nucleation rate per unit volume
$\Gamma(t^{\prime})/\mathcal{V}$, which governs the progress of the phase
transition. It is most sensitive to the dimensionless combination $S=E_{\rm
c}/T$, which decreases rapidly from infinity for $T<T_{\text{c}}$. Hence
$\Gamma(t^{\prime})/\mathcal{V}$ grows very rapidly from zero.

As a first estimate, let $t_{H}$ be the time where one bubble is nucleated per
Hubble volume per Hubble time. This means that at $t_{H}$ the nucleation rate
per volume $\mathcal{V}$ is given by

$\frac{\Gamma(t_{H})}{\mathcal{V}}=H^{4},$ (7.10)

where $H$ is the Hubble rate. The time $t_{H}$ can be thought of as the moment
when the phase transition starts.

As a rough estimate, suppose there is a phase transition with
$T_{\text{c}}=100\,\text{GeV}$. We estimated previously that
$\Gamma(t^{\prime})/\mathcal{V}\sim
T_{\text{c}}^{4}\exp(-E_{c}/T_{\text{c}})$, where we have dropped all
dimensionless constants. The Friedmann equation tells us that the Hubble rate
is $H(T_{\text{c}})^{2}\sim GT^{4}$. Hence the first bubble nucleates when $S$
has dropped to

$S_{H}\sim\log(G^{2}T_{\text{c}}^{4})\sim 100\,.$ (7.11)

The logarithm means that the value of $S$ needed for the first bubble to
nucleate in a Hubble volume is insensitive to all the constants we have
dropped.

As the universe continues to cool, the first bubbles grow, and new bubbles
appear, which convert the old phase into the new one. Let us consider a later
reference time, which we denote by $t_{f}$, such that $h(t_{f})=1/e$ is
satisfied. That means that when $t=t_{f}$ roughly $64\%$ of the universe is
converted to the Higgs phase.

We now define the transition rate parameter $\beta$ as

$\beta\equiv\frac{\differential}{\differential
t}\log\left(\frac{\Gamma(t)}{\mathcal{V}}\right)\Bigg{|}_{t=t_{f}}\,.$ (7.12)

where the derivative is to be evaluated at $t_{f}$. Around this time, we can
now write Eq. (6.41) after a Taylor expansion of $\log\Gamma$ as

$\frac{\Gamma}{\mathcal{V}}=\frac{\Gamma_{f}}{\mathcal{V}}e^{\beta(t-t_{f})}$
(7.13)

with $\Gamma_{f}\equiv\Gamma(t_{f})$. We are now able to perform a saddle-
point approximation to the integral Eq. (7.9), as

$\displaystyle-\log h(t)$ $\displaystyle\simeq$
$\displaystyle\int^{t}\differential t^{\prime}\frac{4\pi}{3}v_{\rm
w}^{3}(t-t^{\prime})^{3}\frac{\Gamma_{f}}{\mathcal{V}}e^{\beta(t^{\prime}-t_{f})}$
(7.14) $\displaystyle=$ $\displaystyle\frac{4\pi}{3}v_{\rm
w}^{3}\frac{\Gamma_{f}}{\mathcal{V}}\frac{3!}{\beta^{4}}e^{\beta(t-t_{f})}\,.$
(7.15)

Using the definition of $t_{f}$, we find that

$8\pi\frac{v_{\text{w}}^{3}}{\beta^{4}}\frac{\Gamma_{f}}{\mathcal{V}}=1\,,$
(7.16)

which can be used to compute $t_{f}$, once the form of $\Gamma(t)$ is known.
Thus we find that $h(t)$ takes a very simple form, namely

$h(t)=\exp[-e^{\beta(t-t_{f})}\Big{]}\,.$ (7.17)

In Fig. (11) we show the fractional volume that is in the metastable ground
state, i.e. $h(t)$, as a function of time $t$. We also plot the derivative of
that quantity, which quantifies the rate at which the physical volume of the
metastable ground state is changing. It turns out that the reference time
$t_{f}$ admits an alternative interpretation. Namely, the rate $\dot{h}$ at
which the universe is converted from the old to the new ground state has a
maximum at time $t_{f}$. This can be seen from the second time derivative of
the fractional volume,

$\displaystyle\frac{\differential}{\differential t}\left(\frac{\differential
h(t)}{\differential t}\right)$
$\displaystyle=\beta^{2}\left(1-e^{\beta(t-t_{f})}\right)\exp\left(-e^{\beta(t-t_{f})}+\beta(t-t_{f})\right)\,.$
(7.18)

Hence, for $t=t_{f}$ the above expression vanishes due to the first term in
the brackets. This is also demonstrated in Fig. 11. The time555In order to
find a relation between $t_{H}$ and $t_{f}$, we note that Eq. (7.13) reads
$\Gamma_{H}=\Gamma_{f}\exp[\beta(t_{H}-t_{f})]$. Using Eqs. (7.10) and (7.16),
and solving for $t_{H}$ yields  $t_{H}=t_{f}+\beta^{-1}\ln\left(8\pi v_{\rm
w}^{3}\frac{H^{4}}{\beta^{4}}\right)\,.$ The logarithm evaluates to negative
values (much) larger than unity. Therefore, the time between $t_{H}$ and
$t_{f}$ is some orders of magnitude larger than $\beta^{-1}$. $t_{H}$ that we
introduced first, is always smaller than $t_{f}$, i.e. $t_{H}<t_{f}$. At the
time $t_{H}$ basically the entire universe is still in the metastable phase as
can be seen in the plot. Instead, the phase transition completes very rapidly
around $t_{f}$ within a duration set by $\beta^{-1}$. To draw the plot in Fig.
11 we use the realistic values of $T_{H}=105.2\,{\rm GeV}$ and
$T_{f}=104.9\,{\rm GeV}$ for the temperatures and $\beta=8950H$, and $v_{\rm
w}=0.9$. This implies that $t_{H}=2.66\cdot 10^{-10}\,{\rm s}$ and
$t_{f}=2.72\cdot 10^{-10}\,{\rm s}$ such that $t_{f}-t_{H}\sim 10^{-13}\,{\rm
s}$. Shortly after the nucleation rate $\dot{h}$ is at its maximum at $t_{f}$,
essentially the entire universe is in the stable ground state.

Figure 11: The time-evolution of the fractional volume in the metastable
ground state, $h$, (solid blue) and its derivative $\dot{h}$ normalized to the
maximum value (dashed red) is shown. At the time $t_{f}$ the change of the
fractional volume is at its maximum, i.e. the phase transition is very fast in
this moment. Shortly after, essentially the entire universe is in the new
ground state. The time $t_{H}$ is very close before $t_{f}$. Their relation
depends on the values of $v_{\rm w}$, $\beta$ and $H$. To draw the plot, we
use the realistic values of $T_{H}=105.2\,{\rm GeV}$, $\beta=8950H$ and
$v_{\rm w}=0.9$. Then the phase transition has completed at a temperature
$T_{f}=104.9\,{\rm GeV}$. This choice implies that $t_{f}-t_{H}\sim
10^{-13}\,{\rm s}$.

The bubbles can nucleate only in the metastable phase. Hence, using Eq. (7.13)
we can deduce that the bubble number density $n_{\text{bubble}}(t)$ in the
universe is given by

$\displaystyle n_{\text{bubble}}(t)$ $\displaystyle=$
$\displaystyle\int^{t}\text{d}t^{\prime}\,\frac{\Gamma(t^{\prime})}{\mathcal{V}}h(t^{\prime})=\frac{\Gamma_{f}}{\mathcal{V}}\int^{t}\text{d}t^{\prime}\,e^{\beta(t^{\prime}-t_{f})}h(t^{\prime})$
(7.19) $\displaystyle=$
$\displaystyle-\frac{\Gamma_{f}}{\mathcal{V}}\frac{1}{\beta}\int^{t}\text{d}t^{\prime}\frac{\differential
h(t^{\prime})}{\differential
t^{\prime}}=\frac{\Gamma_{f}}{\mathcal{V}}\frac{1}{\beta}\left(1-h(t)\right)\,,$
(7.20)

where we used Eq. (7.17) in the second line. At late times when
$t\rightarrow\infty$ and upon using Eq. (7.16), this expression simplifies to

$n_{\text{bubble}}(\infty)=\frac{\Gamma_{f}}{\mathcal{V}}\frac{1}{\beta}=\frac{\beta^{3}}{8\pi
v_{\text{w}}^{3}}\equiv\frac{1}{R_{*}^{3}},$ (7.21)

which is the final bubble density, and defines the mean bubble centre
separation $R_{*}=\left(n_{\text{bubble}}(\infty)\right)^{-1/3}$. As we will
see later, the mean bubble centre separation is linked to the position of the
peak of the gravitational wave power spectrum produced by the first-order
phase transition.

Note that for deflagrations, one should really take into account the heating
of the fluid in front of the bubble wall, which will reduced the net
nucleation rate. In the extreme case of a slow wall and large latent heat,
nucleation can stop altogether, and the interior of the bubbles can reheat to
the critical temperature. In this case, the universe is in a mixed phase, and
the transition takes about a Hubble time to complete [67].

###  7.2 Relativistic combustion

Figure 12: Fluid velocities right ahead of and behind the wall in the case of
a supersonic deflagration. The velocities just ahead the wall are denoted by
$+$ while the velocities just behind the wall are denoted by a $-$ subscript.
The left panel shows the velocities in the universe rest frame in which the
wall moves with the speed $v_{\text{w}}$. The right panel shows the velocities
in the wall rest frame where $v_{\text{w}}=0$ by definition and we indicate
the fluid velocities in this frame by a tilde, $\tilde{v}_{\pm}$. The enthalpy
is denoted by $w$. Figures taken from Ref. [66].

Now we will assume that the wall velocity $v_{\rm w}$ is known, and discuss
solutions to the hydrodynamic equations following Refs. [66, 68]. Let
$u^{\mu}$ be again the 4-velocity of the fluid and we introduce the enthalpy
$w=e+p$ of the fluid. First, we work in the wall rest frame where
$v_{\text{w}}=0$, see Fig. 12, and introduce a subscript $\pm$ to distinguish
between quantities evaluated just ahead of the wall ($+$) and just behind the
wall ($-$). From the conservation of energy-momentum at the wall it follows
that

$w_{-}\tilde{\gamma}^{2}_{-}\tilde{v}_{-}=w_{+}\tilde{\gamma}^{2}_{+}\tilde{v}_{+}\,,\qquad
w_{-}\tilde{\gamma}^{2}_{-}\tilde{v}_{-}^{2}+p_{-}=w_{+}\tilde{\gamma}^{2}_{+}\tilde{v}_{+}^{2}+p_{+}\,,$
(7.22)

where $\tilde{v}_{\pm}$ are the fluid velocities in the rest frame of the wall
and $\tilde{\gamma}_{\pm}=(1-\tilde{v}_{\pm}^{2})^{-1/2}$. The enthalpy
density $w$ and pressure $p$ are scalars, and so they are the same in the wall
frame and the universe rest frame. These are the bubble wall junction
conditions. They can be rearranged to give

$\tilde{v}_{+}\tilde{v}_{-}=\frac{p_{+}-p_{-}}{e_{+}-e_{-}}\,,\qquad\frac{\tilde{v}_{+}}{\tilde{v}_{-}}=\frac{e_{-}+p_{+}}{e_{+}+p_{-}}$
(7.23)

For convenience, let us introduce the trace anomaly

$\theta=\frac{1}{4}(e-3p)\,,$ (7.24)

which is proportional to the trace of the energy-momentum tensor. The trace
anomaly should vanish for an ultra-relativistic equation of state. We denote
the difference of the trace anomaly just ahead and behind the wall by
$\Delta\theta=\theta_{+}-\theta_{-}$. From here, one defines the transition
strength $\alpha_{+}$ and the enthalpy ratio $r$ as

$\alpha_{+}=\frac{4\Delta\theta}{3w_{+}}\,,\qquad r=\frac{w_{+}}{w_{-}}$
(7.25)

In terms of these quantities, Eq. (7.23) can be rearranged as

$\tilde{v}_{+}\tilde{v}_{-}=\frac{1-(1-3\alpha_{+})r}{3-3(1+\alpha_{+})r}\,,\qquad\frac{\tilde{v}_{+}}{\tilde{v}_{-}}=\frac{3+(1-3\alpha_{+})r}{1+3(1+\alpha_{+})r}$
(7.26)

These equations can be solved for $\tilde{v}_{+}$ in terms of $\tilde{v}_{-}$
(or vice versa). The result is

$\tilde{v}_{+}=\frac{1}{1+\alpha_{+}}\left(\frac{\tilde{v}_{-}}{2}+\frac{1}{6\tilde{v}_{-}}\pm\sqrt{\left(\frac{\tilde{v}_{-}}{2}-\frac{1}{6\tilde{v}_{-}}\right)^{2}+\frac{2}{3}\alpha_{+}+\alpha_{+}^{2}}\right)$
(7.27)

which depends only on $\tilde{v}_{-}$ and $\alpha_{+}$ but not on $r$. As we
explain in more detail later, the upper sign is to be taken for
$\tilde{v}_{-}>1/\sqrt{3}$, while the lower sign is to be taken if
$\tilde{v}_{-}<1/\sqrt{3}$ in order to ensure the solution to be physical.
This implies that both velocities are either subsonic (deflagrations) or
supersonic (detonations). Further, $\tilde{v}_{+}$ must be positive because
the fluid must flow through the wall from the outside to the inside of the
bubble. This requires $\alpha_{+}<1/3$. In Fig. 13, $\tilde{v}_{+}$ as a
function of $\tilde{v}_{-}$ is shown for different values of $\alpha_{+}$.

Figure 13: The fluid velocity just ahead the wall $\tilde{v}_{+}$ as a
function of the fluid velocity just behind the wall $\tilde{v}_{-}$ in the
rest frame of the wall for different transitions strength parameters
$\alpha_{+}$. For physical solutions to exist, both velocities have to be
subsonic (deflagrations) or supersonic (detonations) as indicated by the black
lines. Non-physical solutions are indicated by the gray lines. The fluid
velocity has an extremum at $\tilde{v}_{+}=c_{\rm s}$ with the speed of sound
set to $c_{\rm s}=1/\sqrt{3}$ in this plot. The velocity profiles at the blue
dots are presented in Fig. 15. Figure taken from Ref. [66]

So far, from the conservation of energy-momentum, we obtained a relation
between the fluid velocities ahead and behind the wall which also allows one
to compute the enthalpies $w$. We can obtain two more independent equations by
projecting the conservation equation onto the fluid 4-velocity
$u^{\mu}=\gamma(1,\vec{v})^{\mu}$ and the space-time orthonormal vector
$\bar{u}^{\mu}=\gamma(v,\vec{v}/v)^{\mu}$. These vectors satisfy
$\bar{u}_{\mu}\bar{u}^{\mu}=1$ and $u_{\mu}\bar{u}^{\mu}=0$. Projecting the
conservation equation yields

$\displaystyle 0=u_{\mu}\partial_{\nu}T^{\mu\nu}$ $\displaystyle=$
$\displaystyle-\partial_{\mu}(wu^{\mu})+u^{\mu}\partial_{\mu}p\,,$ (7.28)
$\displaystyle 0=\bar{u}_{\mu}\partial_{\nu}T^{\mu\nu}$ $\displaystyle=$
$\displaystyle
w\bar{u}^{\nu}u^{\mu}\partial_{\mu}u_{\nu}+\bar{u}^{\mu}\partial_{\mu}p\,,$
(7.29)

which are referred to as continuity equations.

In order to simplify these equations, we assume that the bubbles are
spherically symmetric. The radius of the bubble is given by $R=v_{\text{w}}t$
where we set the nucleation time to $t^{\prime}=0$. Since there is no length
scale involved in the problem, the differential equations should exhibit
similarity solutions that depend on the dimensionless coordinate $\xi=r/t$
only. We can write the fluid velocity as
$\vec{v}=v(r,t)\vec{r}=v(\xi)\vec{r}$, where $\vec{r}$ is a unit radial
vector. The continuity equations can be rearranged to

$\displaystyle\frac{\differential v}{\differential\xi}$ $\displaystyle=$
$\displaystyle\frac{2v(1-v^{2})}{\xi(1-\xi v)}\left(\frac{\mu^{2}}{c_{\rm
s}^{2}}-1\right)^{-1}\,,$ (7.30) $\displaystyle\frac{\differential
w}{\differential\xi}$ $\displaystyle=$ $\displaystyle w\left(1+\frac{1}{c_{\rm
s}^{2}}\right)\gamma^{2}\mu\frac{\differential v}{\differential\xi}\,.$ (7.31)

Here, $c_{\rm s}^{2}=\differential p/\differential e$ is the speed of sound
and

$\mu=\frac{\xi-v}{1-\xi v}$ (7.32)

is the fluid velocity at $\xi$ in a frame that is moving outward at speed
$\xi$. Solving these equations requires numerical techniques.

In general, the speed of sound $c_{\rm s}$ depends on the temperature which
changes across the wall. The situation simplifies if one assumes an
ultrarelativistic equation-of-state in both phases such that $c_{\rm
s}^{2}=1/3$ everywhere. An example of such an equation of state is the “bag”
model, where

$p_{\text{s}}=a_{\text{s}}T^{4}-V_{\text{s}},\quad
p_{\text{b}}=a_{\text{b}}T^{4},$ (7.33)

where $a_{\text{s}}$, $a_{\text{b}}$ and $V_{\rm s}$ are positive constants
with $a_{\text{s}}>a_{\text{b}}$. The subscripts s and b refer to the
symmetric (i.e., where $\phi=0$) and broken (i.e., where $\phi=\phi_{\rm b}$)
phases respectively. In this case, the fluid speed equation can be integrated
separately. More realistic equations of state have also recently been
considered in Ref. [69].

Solutions are obtained by integrating the continuity equations starting at the
position of the wall, $\xi_{\rm w}=v_{\text{w}}$, while the wall is assumed to
be infinitesimally thin. The boundary conditions at the wall read
$v\rightarrow v_{\pm}=\mu(\xi_{\rm w},\tilde{v}_{\pm})$ as
$\xi\rightarrow\xi_{\rm w}^{\pm}$ where $\xi_{\rm
w}^{\pm}=v_{\text{w}}\pm\delta$ and $\delta$ infinitesimally small. On the
other hand, the fluid velocity must vanish, $v=0$, at the center of the bubble
$\xi=0$ due to spherical symmetry and at $\xi=1$ due to causality because we
assume the fluid to be undisturbed until a signal from the wall arrives.

The solutions can be classified by how the boundary conditions are satisfied.
Specifically, there are only two possibilities to smoothly approach $v=0$.
Either one starts with $v=0$ or one starts in the region $\xi>c_{\rm s}$ and
$\mu(\xi,v)>c_{\rm s}$ implying that $\differential v/\differential\xi>0$ and
then integrating backwards in $\xi$. The only other way to meet the boundary
conditions is by a discontinuity, i.e., by a shock. Hence, we can identify the
following three classes of solutions:

  * •

Subsonic deflagrations: In a subsonic deflagration the wall moves at a
subsonic speed, $v_{\text{w}}<1/\sqrt{3}$, and the fluid is at rest everywhere
inside the bubble, $v_{-}=0$. In the wall rest frame we thus find
$\tilde{v}_{-}=v_{\text{w}}$. The fluid velocity just ahead the wall in the
universe frame is $v_{+}=\mu(\tilde{v}_{+}(\alpha_{+},v_{\text{w}}),\xi_{\rm
w})$. In order to ensure $v_{+}>0$ one has to choose the negative sign in Eq.
(7.27) as anticipated. Further one finds that ahead of the wall the fluid
velocity decreases until a shock occurs outside which the fluid is at rest.
This situation is depicted schematically in the left panel of Fig. 14. The
velocity profile is shown in the upper left panel of Fig. 15.

  * •

Detonations: A detonation is characterized by a fluid exit speed in the wall
frame $\tilde{v}_{-}>1/\sqrt{3}$, with the fluid being at rest everywhere
outside the bubble, $v_{+}=0$. In the wall rest frame we thus have
$\tilde{v}_{+}=v_{\text{w}}$. The condition on $\tilde{v}_{-}$ means there is
a minimum for $\tilde{v}_{+}$ called the Chapman-Jouguet speed $v_{\rm CJ}$.
It is given by

$v_{\rm
CJ}(\alpha_{+})=\frac{1}{\sqrt{3}}\left(\frac{1+\sqrt{\alpha_{+}+3\alpha_{+}^{2}}}{1+\alpha_{+}}\right)\,.$
(7.34)

In order to ensure that $\tilde{v}_{-}(\alpha_{+},v_{\text{w}})>c_{\rm s}$ and
hence $\differential v/\differential\xi>0$ the positive sign has to be taken
in Eq. (7.27) as mentioned before. As result, the fluid velocity smoothly
decreases as one moves further inside the bubble until it is at rest at
$\xi=c_{\rm s}$. The schematic visualization can be found in the right panel
of Fig. 14 while the velocity profile is shown in the upper right panel of
Fig. 15.

  * •

Supersonic deflagrations (hybrids): There is a hybrid between the two
aforementioned classes of solutions with the wall moving at supersonic speed,
$v_{\text{w}}>1/\sqrt{3}$, which occurs for $\tilde{v}_{-}=1/\sqrt{3}$. A
physical solution for $\tilde{v}_{+}(\alpha_{+},1/\sqrt{3})$ exists provided
that it is larger than the wall speed to ensure a positive $v_{+}$. In front
of the wall, the fluid behaves in exactly the same manner as for a subsonic
deflagration. The hybrid is schematically shown in the middle panel of Fig. 14
while the velocity profile is shown in the upper middle panel of Fig. 15.

In the lower panels of Fig. 15 we present the enthalpy profiles for all three
classes. Note that the solutions are found once $\alpha_{+}$ is given, which
requires knowing the energy density and pressure just in front of the wall.
However, one normally specifies these quantities far in front of the wall,
where the undisturbed plasma is at the nucleation temperature. For
detonations, this presents no problem, but for deflagrations shooting method
is needed to find the correct value of $\alpha_{+}$ which matches to the
asymptotic energy and pressure.

For more details on these solutions we refer to Refs. [66, 68, 69, 70].

Figure 14: Sketch of the three different cases of relativistic combustion. A
subsonic deflagration (left) occurs detonation when the fluid is at rest
inside the bubble and the wall moves at subsonic speed. The opposite case is a
detonation (right) where the fluid outside the bubble is at rest and the wall
moves at supersonic speed. In the hybrid case (mid) the wall speed is
supersonic and the fluid is moving both ahead and behind the wall. Credit: D.
Cutting.

Figure 15: The fluid velocity $v$ and the enthalphy $w$ as a functions of the
dimensionless coordinate $\xi=r/t$ for the three classes of solutions is
depicted. The wall is located at $\xi_{\rm w}$. The left panel shows the
velocity profile for a subsonic deflagration, where the fluid velocity
smoothly decreases from the wall until it reaches a shock ahead of which the
fluid is at rest. The middle panel shows supersonic deflagration (hybrid). The
right panel shows a detonation, where the velocity smoothly increases until
the location of the wall, while the fluid is at rest outside the bubble. The
three cases are represented by the blue dots in Fig. 13. Figures taken from
Ref. [66].

###  7.3 Energy redistribution

Roughly speaking, the expanding bubble converts potential energy of the scalar
field into kinetic energy and heat. The kinetic energy fraction of a single
bubble is a reasonable estimate of the kinetic energy fraction of the entire
fluid flow [71, 66]. The kinetic energy fraction of a single bubble can then
be used to estimate the power of the gravitational wave signal. Hence, making
the above statement quantitatively precise is crucial.

To remind ourselves, the spatial components of fluid energy-momentum tensor
are given by

$T^{i}_{\ j}=w\gamma^{2}v^{i}v_{j}+p\delta^{i}_{j}\,$ (7.35)

with $\vec{v}$ the fluid velocity, $w=e+p$ the enthalpy, $e$ the total energy
and $p$ the pressure of the fluid. The kinetic energy of the fluid can be
obtained as the trace of the energy-momentum tensor minus the trace in the
rest frame of the fluid, and integrating over space. The kinetic energy
fraction is thus given by

$K=\frac{1}{\mathcal{V}\bar{e}}\int\differential^{3}x\,w\gamma^{2}v^{2}\,.$
(7.36)

where $\mathcal{V}$ is the averaging volume, and $\bar{e}$ is the mean energy
density.

The kinetic energy fraction of a single bubble describes how much of the
initially available energy contained in the bubble is converted into kinetic
energy, which can

$K_{1}=\frac{3}{\xi_{\rm w}^{3}e_{\rm
s}}\int\differential\xi\,\xi^{2}w\gamma^{2}v^{2}\,,$ (7.37)

where $e_{\rm s}$ is the mean energy density in the symmetric phase.
Conservation of energy, and the rapidity of the transition compared with the
Hubble rate, implies that the mean energy density in the broken phase $e_{\rm
b}=e_{\rm s}$. Our estimate of the kinetic energy fraction in the broken phase
is therefore $K=K_{1}$.

It is also useful to define the enthalpy-weighted root-mean-square 4-velocity
of the fluid $\bar{U}_{\text{f}}$, through

$\bar{U}_{\text{f}}^{2}=\frac{3}{\xi_{\rm
w}^{3}\bar{w}}\int\differential\xi\xi^{2}w\gamma^{2}v^{2}\,,$ (7.38)

which can be related to the kinetic energy fraction via

$K=\Gamma\bar{U}_{\text{f}}^{2}\,,$ (7.39)

where $\Gamma=\bar{w}/\bar{e}$ is the adiabatic index of the fluid in the
broken phase.

In the previous section we introduced the transition strength parameter
$\alpha_{+}$, which is defined in terms the enthalpy and the trace anomalies
just ahead and behind the wall. These quantities are perturbative by
definition and thus $\alpha_{+}$ is difficult to compute in practice. Instead,
it would be nice to have a quantity that measures the transition strength in
terms of background variables only. As before, we denote variables that
correspond to the symmetric phase (i.e., where $\phi=0$) with a subscript
${\rm s}$ while variables that correspond the broken phase (i.e., where
$\phi=\phi_{\rm b}$) with a subscript ${\rm b}$. In terms of these background
variables, one defines the transition strength parameter as

$\displaystyle\alpha_{\rm n}=\frac{4}{3}\frac{\Delta\theta}{w_{\rm
s}}\Bigg{|}_{T=T_{\rm n}}\,,$ (7.40)

where now $\Delta\theta=\theta_{\rm s}-\theta_{\rm b}$. The two definitions of
the transition strength coincide only in the case of detonations within the
bag model. The transition strength parameter is defined at the nucleation
temperature $T_{\rm n}$. Note that in the literature there exist various
definitions of the transition strength parameter $\alpha$ that are related in
nontrivial and model-dependent ways.

It is convenient to define the efficiency factor $\kappa$ that quantifies how
much of the available energy is converted into kinetic energy. The available
energy is determined by the trace anomaly $\theta$ introduced in Eq. (7.24).
It is related to the potential energy of the scalar field $V_{T}(\phi)$ as

$\theta=V_{T}(\phi)-\frac{1}{4}T\frac{\partial V_{T}}{\partial T}\,.$ (7.41)

This implies that not all of the potential energy can be converted into
kinetic energy and heat. Therefore, one defines the efficiency factor as666In
the bag model, it is assumed that $\theta_{\rm b}=0$ such that $\Delta\theta$
coincides with the so-called bag constant $\epsilon$.

$\kappa=\frac{3}{\xi_{\rm
w}^{3}\Delta\theta}\int\differential\xi\,\xi^{2}w\gamma^{2}v^{2}\,,$ (7.42)

where $\Delta\theta$ is again the difference of the trace anomaly between the
symmetric and broken phase. Consequently, efficiency factor, transition
strength and kinetic energy fraction are related as

$K=\frac{\kappa\alpha_{\rm n}}{1+\alpha_{\rm n}+\delta_{\rm n}}\,,$ (7.43)

where $\delta_{\rm n}=4\theta_{\rm b}/(3w_{\rm s})$. The efficiency factor can
be computed numerically as a function of $\alpha_{\rm n}$ and $\xi_{\rm w}$ as
in Ref. [68]. To give an explicit example, in the case of detonations and
using the bag model, where $\theta_{\rm b}=0$, they find the approximate
relation

$\kappa\simeq\frac{\alpha_{\rm n}}{0.73+0.083\sqrt{\alpha_{\rm n}}+\alpha_{\rm
n}}$ (7.44)

for $\xi_{\rm w}\to 1$. For $\alpha_{\rm n}<1$, one can generally take
$\kappa\sim\alpha_{\rm n}$, except for low wall speeds, and wall speeds near
the Chapman-Jouguet speed, where the parametric dependence is closer to
$\kappa\sim\sqrt{\alpha_{\rm n}}$ [68]. In addition, relations between the
transition strength and the inverse duration of phase transition were recently
established for a variety of BSMs in Ref. [72, 73].

###  7.4 Sound waves

Radial perturbations of the radial velocity field are longitudinal in
character, and so we can think of the bubble expansion as an explosion,
generating a compression wave which, once the bubbles have disappeared,
propagate through the fluid as sound waves. We conclude this section by
studying sound waves in a relativistic fluid.

Let us work in the plane wall approximation and let the $z$-direction be
orthogonal to the wall. Then, the wall moves only in the $z$-direction and
also the fluid is perturbed in the $z$-direction only. Let us write the
perturbed energy density as $e=\bar{e}+\delta e$ and the perturbed pressure as
$p=\bar{p}+\delta p$ with $\\{\delta e,\delta p,v^{z}\\}\ll 1$ with $v^{z}$
the fluid velocity perturbation. The components of the fluid energy-momentum
tensor simplify to

$T^{tt}=w\gamma^{2}-p\,,\quad T^{tz}=w\gamma^{2}v^{z}\,,\quad
T^{zz}=w\gamma^{2}(v^{z})^{2}+p\,,$ (7.45)

while the other components vanish identically. From the $t$-component of the
energy-momentum conservation equation we find

$\partial_{t}\delta e+\bar{w}\partial_{z}v^{z}=0$ (7.46)

and from the $z$-component

$\bar{w}\partial_{t}v^{z}+\partial_{z}\delta p=0\,.$ (7.47)

Note that both $\delta e$ and $\delta p$ depend on temperature as $\delta
p=(\frac{\partial p}{\partial T}/\frac{\partial e}{\partial T})\delta e=c_{\rm
s}^{2}\delta e$. Therefore, the $t$\- and $z$-components of the conservation
equation can be combined to

$(\partial_{t}^{2}-c_{\rm
s}^{2}\partial_{z}^{2})v^{z}=0\,,\quad(\partial_{t}^{2}-c_{\rm
s}^{2}\partial_{z}^{2})\delta p=0\,.$ (7.48)

These equations describe sound waves that travel at the speed $c_{\rm s}$
through the fluid. Sound waves are a collective mode of the fluid velocity $v$
and temperature $T$. They are longitudinal as the fluid velocity varies along
the direction of travel.

##  8 Gravitational Waves

Figure 16: The spectrum of gravitational waves with possible sources and
detectors. The gravitational wave spectrum from first-order phase transitions
in the early universe are expected to peak in the sensitivity range of LISA.
Credit: NASA.

Gravitational waves were predicted by Einstein in 1916 [74, 75], although it
took about forty years for them to be understood as physical, rather than
coordinate artefacts (see e.g. [76]). According to the general theory of
relativity, they are perturbations of space-time that travel at the speed of
light. They are generated by an accelerating asymmetric mass distribution;
more precisely, a distribution of energy-momentum with a time-dependent
quadrupole moment. The strongest astrophysical sources of gravitational waves
are compact binary systems: combinations of neutron stars, black holes or
white dwarfs. Besides these astrophysical sources there are also cosmological
sources, in particular the early Universe (see e.g. [22]). In this section we
give a rough overview about the production of gravitational waves by first-
order phase transitions. For a detailed discussion we refer to the review
paper by the LISA working group [12]. Good textbooks are [77, 78]. Fig. 16
gives an overview over the spectrum of gravitational waves and possible
sources.

###  8.1 Introduction

Gravitational waves yield expansion of spacetime in one direction and
contraction in the other, in the plane perpendicular to the propagation
direction. Modern detection methods are mostly based on measuring the
variation of the distance between two test masses by interferometry, as
sketched in Fig. 17. The first direct detection was performed in 2015 by the
LIGO/VIRGO science collaborations [10].

The results of the first measurement [10] showed that two black holes with
masses $M_{1}=36\pm 6M_{\odot}$ and $M_{2}=29\pm 4M_{\odot}$ produced the
signal. According to numerical simulations of such an event, an energy
equivalent of $M_{\rm gw}\simeq 3M_{\odot}$ must have been radiated away as
gravitational waves.

Name | Location | Arm length
---|---|---
GEO | Germany |  $0.6$km
(a)LIGO | USA (2) |  $4$km
(a)VIRGO | Italy |  $3$km
KAGRA | Japan |  $4$km
LIGO-India∗ | India |  $4$km

Figure 17: Sketch of the gravitational wave detector LIGO. Table 3: Current
and planned (*) ground-based gravitational wave detectors using laser
interferometry.

After this first direct detection, in 2017 gravitational waves of the merger
of a binary neutron star system were observed [79, 80]. This signal was
accompanied by the detection of a gamma ray burst (GRB) which allows to put
constraints on the travel speed of gravitational waves $c_{\rm gw}$. It was
found that gravitational waves travel at speed of light to a large precision,

$|c_{\rm gw}^{2}-1|\lesssim 10^{-5}\,.$ (8.1)

This puts tight constraints on modifications to general relativity [81, 82]. A
total of 11 secure detections are reported from the first two observing runs
(O1 and O2), with 56 further candidates for further analysis from O3, which
ended in March 2019. Further exciting results are expected from the expanding
network of detectors. In Tab. 3 we summarize existing and planned ground-based
interferometric gravitational wave detectors.

We are interested in gravitational waves that are produced in the early
Universe. As a rough estimate of the frequency, take an event at time $t$ that
produces gravitational waves. The minimal frequency is given by $f=t^{-1}\sim
H$ with $H=\dot{a}/a$ the Hubble rate of that time and $a$ the scale factor.
Due to cosmic expansion, the minimal frequency that we can observe today is
redshifted as $f_{0}=a(t)/a(t_{0})f$, where $t_{0}$ is the time now. In Tab. 4
we summarize the minimal frequencies from different events in the early
Universe.

Event | Time [s] | Temperature [GeV] | Frequency [Hz]
---|---|---|---
QCD phase transition | $10^{-3}$ | $0.1$ | $10^{-8}$
EW phase transition | $10^{-11}$ | $100$ | $10^{-5}$
? | $10^{-25}$ | $10^{9}$ | $100$
End of inflation | $\geq 10^{-36}$ | $\leq 10^{16}$ | $\geq 10^{8}$
Table 4: Overview about the minimal frequencies of gravitational waves from
different potential sources in the early universe.

Before going on, let us be a bit more quantitative in our outline of
gravitational waves. We start by considering perturbations about Minkowski
background $\eta_{\mu\nu}$ as

$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}\,,$ (8.2)

where $h_{\mu\nu}\ll 1$ is the metric perturbation. A gravitational wave
$h_{ij}(t,\vec{x})$ is a propagating mode of the transverse
($\partial_{i}h_{ij}=0$) and traceless ($h_{ii}=0$) part of the metric
perturbation, satisfying the equation

$\ddot{h}_{ij}-\nabla^{2}h_{ij}=16\pi GT_{ij}^{\rm TT}\,,$ (8.3)

where $T^{\rm TT}_{ij}$ is the transverse and traceless part of the energy-
momentum tensor, and $G$ is Newton’s constant. The constraints on $h_{ij}$
reduce the number of physical propagating modes to two, called ‘plus’ ($+$)
and ‘cross’ ($\times$) polarisations.

Provided their amplitude is small, gravitational waves can themselves be a
source of energy-momentum, described by the tensor [78]

$T_{\mu\nu}^{\rm gw}=\frac{1}{32\pi
G}\langle\partial_{\mu}h_{ij}\partial_{\nu}h_{ij}\rangle\,.$ (8.4)

We aim at computing the power spectrum of the energy density of gravitational
waves. From the $tt$-component of the energy-momentum tensor, we find the
energy density,

$\rho_{\rm gw}=\frac{1}{32\pi G}\langle\dot{h}_{ij}^{2}\rangle\,.$ (8.5)

After a spatial Fourier transform of the metric, one can consider the
contribution of a frequency interval $\differential f$ to the total energy
density, where the frequency of the gravitational wave $f$ is related to the
wavenumber of the Fourier transform $k$ by $f=k/2\pi$. More often, one studies
the energy density per logarithmic frequency interval $\differential\ln f$,
which can be written

$\frac{\mathrm{d}\rho_{\rm gw}}{\mathrm{d}\ln f}=\frac{\pi}{4G}f^{3}S_{h}(f),$
(8.6)

which can be treated as a definition of the one-sided strain power spectral
density $S_{h}(f)$, defined for positive frequencies, $f>0$. From here, one
can indicate the amplitude of a gravitational wave spectrum at frequency $f$
by the root power spectral density $h(f)=\sqrt{S_{h}(f)}$, which has dimension
$\sqrt{{\rm Hz}^{-1}}$, or the (dimensionless) characteristic strain
$h_{c}(f)=\sqrt{fS_{h}(f)}$.

A convenient measure in cosmology is the fractional density in gravitational
waves

$\Omega_{\rm gw}=\frac{\rho_{\rm gw}}{\rho_{\rm tot}}$ (8.7)

where $\rho_{\rm tot}$ is the total energy density of the universe, related to
the Hubble rate $H$ by the Friedmann equation,

$H^{2}=\frac{8\pi G\rho_{\rm tot}}{3}.$ (8.8)

With these quantities, one can define the fractional gravitational wave energy
density per logarithmic frequency,

$\frac{\mathrm{d}\Omega_{\rm gw}}{\mathrm{d}\ln f}=\frac{1}{\rho_{\rm
tot}}\frac{\mathrm{d}\rho_{\rm gw}}{\mathrm{d}\ln
f}=\frac{2\pi^{2}}{3H^{2}}f^{3}S_{h}(f)\,.$ (8.9)

often referred to as the gravitational wave (power) spectrum. The
characteristic strain and the fractional density in gravitational waves are
related through

$h_{c}(f)=\frac{H_{0}}{f}\sqrt{\frac{3}{2\pi^{2}}\frac{\mathrm{d}\Omega_{\rm
gw}}{\mathrm{d}\ln f}}\,.$ (8.10)

Going from Minkowski space to the Friedmann-Lemaître-Robertson-Walker metric
describing an expanding universe, one finds that the gravitational waves
behave just like electromagnetic radiation: the frequency is redshifted, and
the energy density decreases as the fourth power of the scale factor [77].

The metric perturbations $h_{ij}$ cause changes in lengths in directions
perpendicular to the propagation direction of the wave, which can be measured
by interferometry. Let $l_{i}$ and $m_{i}$ be components of unit vectors along
the arms of an interferometer located at $\vec{x}$. The strain, or relative
change in length, in the detector, in the limit that the arm length is much
less than the wavelength of the gravitational wave, is

$h(t)=\frac{1}{2}h_{ij}(t,\vec{x})(l_{i}l_{j}-m_{i}m_{j})\,,$ (8.11)

whose Fourier transform is

$\tilde{h}(f)=\frac{1}{2}\int_{-\infty}^{\infty}{\rm d}t\,e^{-2\pi
ift}h_{ij}(l_{i}l_{j}-m_{i}m_{j})\,.$ (8.12)

We define the one-sided strain power spectral density of the signal in the
interferometer $S^{I}_{h}(f)$ as

$\langle\tilde{h}(f)\tilde{h}^{*}(f^{\prime})\rangle=\frac{1}{2}S^{I}_{h}(f)\delta(f-f^{\prime})\,.$
(8.13)

This is not the same quantity as the strain power spectral density in a
gravitational wave. In general they are related by a response function
$\mathcal{R}^{I}(f)$,

$S^{I}_{h}(f)=\mathcal{R}^{I}(f)S_{h}(f)\,.$ (8.14)

For an interferometer with perpendicular arms, whose arm length is much less
than the wavelength so that $fL\to 0$, $\mathcal{R}^{I}(f)\to 1/5$ [83].

Figure 18: Visualization of the orbit of LISA. LISA is a space-based laser
interferometer designed to measure gravitational waves with frequencies of