Phase transitions in the early universe

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 firstorder 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. Copyright M. Hindmarsh et al. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 25-09-2020 Accepted 08-02-2021 Published 15-02-2021 Check for updates doi:10.21468/SciPostPhysLectNotes.24


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. [17,21,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.

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.

Basic thermodynamics -the bosonic harmonic oscillator
The basic object of thermodynamics is the partition function Z(T ) = Tr e −βĤ , (2.1) whereĤ is the Hamiltonian operator and β = 1/T the inverse temperature with T the temperature. The free energy, entropy and energy of the system are given by 3) (2.4) respectively. First, we will study a single bosonic harmonic oscillator, the simplest case. To compute its partition function we consider 〈n|e −βĤ |n〉 , (2.5) where the |n〉 are the eigenstates of the HamiltonianĤ of the harmonic oscillator with angular frequency ω, together satisfyingĤ For the partition function and the free energy this yields

7)
F bho (T, ω) = 1 2 ω + T ln 1 − e −βω , (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φ( x, t) and decompose it into its Fourier modeŝ where we postulate that the operatorsâ,â † satisfy the commutation relation The equation of motion for the free field is the Klein-Gordon equation, (2.12) which in terms of the Fourier modes reads 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 | k|. The partition function of a bosonic field (indicated by the subscript B) thus factorizes into Z B = k Z bho (T, ω 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 Again, the sum over all momenta k is defined over a finite volume V. In the infinite volume limit V → ∞, the sum is replaced by an integration as The free energy density, i.e., the free energy normalized to the volume V, then becomes 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 counterterm. The standard renormalisation convention takes the zero-temperature ground state free energy to be zero. +O m T 6 .
Here γ E ≈ 0.57721 is the Euler-Mascheroni constant. While the first two terms follow in a relatively simple way using the ζ-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 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 after using Eq. (2.6). Consequently, the free energy for fermions is given by In order to generalize the above expression to fermionic fields we introduce a Dirac spinor fieldΨ α ( x, t), which creates and destroys massive fermions. Every such spinor has four components, denoted by the index α. 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 Ψ α . Analogously to the bosonic case, the free fermionic field is a collection of independent harmonic oscillators, four for each momentum mode | k|. The partition function for a fermionic field F hence factorizes into where Z fho is the partition function of a single fermionic harmonic oscillator, cf. Eq. (2.23). This leads to the fermionic free energy 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 The vacuum energy density V 0,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. In analogy to the bosonic case the free energy density can be written as In Fig. 2  and bosons. In the high temperature limit the free energy can be approximated as Compared to Eq. (2.21), 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.28) together give the free energy of a single Dirac fermion field.

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].

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 φ, 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, φ = v EW 246 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, 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 O(1).
For the other fields of the standard model, c i 1. The Higgs particle is a quantised fluctuation around the ground state, with mass The Higgs field is unique in that its mass does not in general depend linearly on φ. 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 Φ, but for completeness we mention that 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  Figure 3: This figure shows the effective number of relativistic degrees of freedom g eff of a Standard Model plasma as a function of temperature, taking into account interactions between particles, with both perturbative and lattice methods [24].
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 Here, g eff is the effective number of relativistic degrees of freedom, given by (3.6) where N F is the number of Dirac fermions, N V is the number of massive vectors, N V0 is the number of massless vectors and N 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.21) or Eq. (2.28), resp. For the Standard Model at high energies this value is g eff = 106.75. 1 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 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 eff . The second, field-dependent term in Eq. (3.5), V T (φ), is called the thermal effective Higgs potential. According to Eqs. (2.21) and (2.28) it is given by [23,[25][26][27] 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. [27,31].
For high temperatures, the thermal effective potential can be approximated by an expansion in φ/T yielding (3.8) where A, D are constants and λ T depends only logarithmically on the temperature. In the Standard Model we have where we have dropped the logarithmic dependence of λ 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 T c , the potential has a minimum at φ = 0, which is the only ground state or equilibrium state of the system. We will refer to the ground state where φ = 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 which is well-defined only if 2A 2 < 9λD. This case will be of particular interest to us, as in this case the two minima at φ = 0 and φ(T c ) = 2AT c /3λ 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 φ = 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 φ = v EW 246 GeV.

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. For large temperatures T T c (red) the potential has a minimum at φ = 0 and the ground state is symmetric. Below the temperature T 1 > T 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.
In an interacting theory, one splits the HamiltonianĤ =Ĥ 0 +Ĥ 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 (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 1. We can then try to compute the partition function Z = Tr e −β(Ĥ 0 +Ĥ I ) , (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 with f the phase space density. For a boson, which approaches T /ω k for frequencies low compared with the temperature, ω k T . In this limit, the expansion parameter reads = g 2 T ω k , (3.18) which is greater than unity for k g 2 T . The expansion parameter diverges as | k| → 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 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. 2 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 where c is a theory-dependent constant, and m 0 is the mass of the field at zero temperature. Therefore, the expansion parameter ε for massless gauge bosons (such as the photon) with m 0 = 0 is of the order of the coupling constant, ε ∼ g 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 ε 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.

Beyond weak coupling and the Standard Model
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 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 φ 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 firstorder 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.

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.

Distribution function
In the presence of several interacting harmonic oscillators, we need to also include a number operatorN in the definition of the partition function where µ is the chemical potential. The number operator is obtained aŝ 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 p: From now on we switch to the notation ω p → E p such that the dispersion relation reads E 2 p = p 2 + m 2 . The particle number densities are Hence, let us introduce the 1-particle distribution function as where η = +1 for bosonic fields and η = −1 for fermionic fields. The µ = 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, describes the average number of particles of momentum p in a 3-phase space volume element ( x, x + d x) × ( p, p + d 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 We can combine both quantities and define the particle current Furthermore, we define the energy density e(x), the 3-momentum density Π i (x) and the 3- We can combine the last three quantities to form the energy-momentum tensor, where T 00 = e, T 0i = Π i = T i0 , and T i j = Π i j . The integral measure transforms as a scalar under Lorentz transformations because we can rewrite it as where θ (p 0 ) ensures positivity of the energy. Hence, for the energy-momentum tensor T µν to transform as a 2-tensor under Lorentz-transformation (and the particle current j µ as a 1tensor), the distribution function f (p, x) has to be a Lorentz-scalar.

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 µ (τ), p µ (τ)), which is parametrised by proper time τ. The position and momentum after a small proper time interval dτ hence change as where F µ describes an external 4-force. This is sketched in Fig. 6, which shows the phase space at two time steps τ and τ + dτ. The particle distribution functions at time τ and τ + dτ must be the same because the particle number is conserved in an infinitesimal phase space volume. This leads to A Taylor expansion around dτ = 0 yields 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 µ p µ = 0 or (b) We introduce the on-shell condition by hand to arrive at the collisionless relativistic Boltzmann equation x p dx dx dp dp Figure 6: This figure depicts how the phase space volume occupied by particles within (x, x +dx) and (p, p+dp) changes in a time step dτ to (x , x +dx ) and (p , p +dp ). 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 µ = 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 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 dτ because scattering can remove and add particles to the phase space volume element (dp, dx) around (p, x). Let R (R) describe a scattering in time dt which removes (adds) an initial (final) particle with momentum p at position x. The Boltzmann equation including collisions reads Let us quantify R andR 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 where W is called the scattering function. Here, we introduced the short-hand notation d 3 p = d 3 p/(2π) 3 . We can write W in terms of the cross section σ as in terms of the Mandelstam variables s = (p 1 + p 2 ) 2 , t = (p 1 − p 1 ) 2 , and the scattering angle θ is obtained from cos θ = 1 − 2t/(s + 4m 2 ). Let us introduce the collision function C[ f ] =R − R, where the notation reminds us that R andR depend on the distribution function. We now have derived the relativistic Boltzmann equation 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.

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 func- as a consequence of particle number and momentum conservation. The function ψ is a collision invariant. Now we replace the collision function using the Boltzmann equation (4.25). For the case b µ = 0, the above integral then implies 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 µ = 0, we can set a = 0 to arrive at where we took the partial derivative out of the integral and used the definition in Eq. (4.13). This establishes conservation of the tensor, which is the energy-momentum tensor of the system of particles with distribution function f .

Local equilibrium and perfect fluid
The fluid is in local equilibrium in the presence of collisions when R( (4.22) and (4.23) we find that the fluid is in local equilibrium when, ∀ p 1 , p 2 , p 1 , p 2 , . 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 (4.32) in local equilibrium. We rewrite a(x) and b µ (x) in a suggestive notation, , with u 2 = −1, and we recover the classical equilibrium distribution function We then see that β(x) = 1/T is the inverse temperature, µ(x) the chemical potential, and u µ (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, 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 The additional factors of 1± f ( ) i 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 (4.36) Separating primed and unprimed variables reveals that ln As for the classical case, this implies that the distribution function in equilibrium can be written as 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 µ, β and u yields the following expression for the distribution function, In the early universe the chemical potential is negligible, µ 0. Then, the energy-momentum tensor in local equilibrium reads where all quantities depend on x. This is the energy-momentum tensor for a perfect fluid with energy density e and pressure 3 p. Despite an ambiguity in the notation, the pressure is not to be confused with 4-momentum. 3 The pressure is defined as

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 µ . 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 φ) is the space-time coordinate of a particle, parameterised by the particle's proper time τ, and g µν is the space-time metric.
Varying this action we obtain the equation of motion From this equation we identify the 4-momentum p µ = m dx µ / dτ, and the force 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 energymomentum assuming there are no external sources. Now, consider instead the Boltzmann equation with collisions and external forces, where we have introduced the on-shell condition again. We multiply both sides with p ν , integrate over momenta, and use Eq. (4.26) to find 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 ∂ p ν /∂ p µ = δ ν µ . Using Eq. (5.3), this yields where the 4-momentum is constrained such that p 0 = E 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 where λ 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.

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 µν f and one for the Higgs T µν φ , as where we label the fluid pressure as p 1 for now. As an example we use the symmetry breaking potential for the scalar field. Note that there are some field-dependent terms in the fluid pressure, which we define as 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) The full thermally corrected potential then reads In the high-temperature expansion, V 1 (φ, T ) contains all T -dependent terms of the potential in Eq. (3.7).
We demand that the full energy-momentum tensor is conserved, and hence non-conservation of T µν f has to appear in T µν φ . Using Eq. (5.9) this implies 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 where prime denotes the derivative with respect to the argument φ. 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, 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, Consequently the equation for small departures from fluid equilibrium (5.19) becomes to first order in δ f . Let us repackage the contribution from the effective potential into the fluid by defining the fluid pressure as p(φ, which yields for its conservation equation. The right hand side is generated by departures from the equilibrium phase-space distribution f 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, where we use the fluid 4-velocity u µ to contract the open index in order to respect isotropy. The proportionality factorη might in general depend on the field value, the temperature, and the fluid γ-factor, orη =η(T, φ, γ). The overall mass dimension ofη is zero, so it must depend on the ratio φ/T . This leaves us with the equation of motion 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), where w = e + p = Ts is the enthalpy density, and s = dp/dT is the entropy density. From this we can derive a conservation equation for the entropy current S µ = su µ , 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·∂ φ. In this model entropy is generated by the interaction of the field and the fluid.

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. 4 The order parameter is a scalar field φ, and we will take the values of the field in the two phases to be φ = 0 in the high temperature phase and φ = φ b in the low temperature phase.
The transition proceeds via bubble nucleation: individual bubbles of the new ground state φ = φ b appear due to thermal fluctuations. If the temperature of the universe is below the critical value T c , these bubbles grow if they exceed a critical radius R 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 φ at temperature β = 1/T , which is the functional integral with the Hamiltonian where π is the conjugate momentum to the field φ. We will use the partition function to calculate the rate per unit volume of fluctuations over a thermal barrier in the potential V T .

Transition rate in 0 dimensions
To get started, we discuss a simple toy model without spatial dimensions [51]. Therefore, spatial gradients ∇φ do not exist, and the integrals over functions π and φ become ordinary integrals. The theory describes a particle with coordinate φ in a potential V T as depicted in Fig. 8. The partition function for this simple system is We can immediately perform the integral over π yielding To carry out the remaining integration over φ we use the saddle point approximation. In principle the potential V T (φ) has three stationary points, two minima, a local one at φ = 0 and a global one at φ = φ b as well as a maximum at φ = φ 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 φ = 0. To infer the transition rate of particles fluctuating to the stable ground state at φ = φ 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 φ = 0. We hence expand as φ = 0 + δφ to obtain The second stationary point is at φ = φ m with V T (φ m ) < 0. As the second derivative at the stationary point is negative here we need to deform the integration contour into the complex Figure 8: The effective thermal potential V T (φ) for a temperature below the critical temperature, T < T c . The potential has two minima φ = 0 an φ = φ b , and a maximum at φ = φ 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.
plane, which results in the integral acquiring an imaginary part. At this second stationary point we obtain, for the imaginary part, This expression is imaginary because of the negative sign for V T . 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 β = Z 0 β + Z 1 β . Computing the free energy we obtain where is the barrier height in the potential, and we expanded the logarithm for β∆V T 1. The two signs arise due to the choice of the integration contour. Next, we compute the the probability flux Γ across the potential barrier given by The δ-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 Comparing this equation to the expression for the free energy (6.7) establishes the relation The probability flux density is thus proportional to the imaginary part of the free energy. To be explicit, it is given by 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 τ m = 1/ V T (φ m ) . If the mean free time between interactions is much less than τ m , then it is more appropriate to study the diffusion of particles across the barrier, which can be modelled by the Langevin equatioṅ Here, γ is the diffusion constant, and ξ(t) represents the forces exerted by the heat bath. It is a Gaussian random variable, obeying 〈ξ(t)ξ(t )〉 = 2γT δ(t − t ), meaning that it is uncorrelated from one instant to the next, but has a mean amplitude 2γT . In the case γτ m 1, the rate of particles crossing over the barrier is This is a factor 1/γτ 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].

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.

The critical bubble
The starting point is the path integral which for a thermal field theory is given by where the Hamiltonian H is given by Eq. (6.2) now including spatial gradients of the scalar field. We can perform the integration over π, as it is Gaussian, leading to 15) where 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 δE T [φ]/δφ = 0, leading to the following differential equation 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 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 φ = 0 at a temperature just below T c . The asymptotics of the solution are therefore: as r → ∞, φ → 0, and as r → 0, φ (r) → 0, where the prime denotes derivative with respect to r. A trivial solution φ(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λ = 9λD 2A 2 , (6.18) provided one scales the potential with the vacuum energy difference, the field with its stable phase value φ b , and radial distance with the effective mass M = D T 2 − T 2 0 . The results for selected values ofλ 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φ(r).
At the critical temperature,λ → 1, and the potential is approximately This leads to the following approximate solution to the equation of motion, where R c is the radius of the critical bubble, w = 1/M (T c ) the thickness of the bubble wall, and R c w . This solution is plotted in the right panel of Fig. 9 for different values ofλ. We refer to this approximate solution as a thin-wall bubble.
The energy of a critical bubble can be computed from the energy functional 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φ and the energy of the metastable state φ = 0, or 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.
is the free-energy difference between the two phases, there exists an approximate solution for E c as a function of radius R, given by Here, σ denotes the surface tension given by with A and λ the constants appearing in the thermal Higgs potential (3.8). It is possible to introduce a critical radius R 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.
such that the critical radius is given by Plugging this expression into Eq. (6.23), we find the energy of the critical bubble as In the thin-wall approximation, the energy of the critical bubble can be computed as a function of temperature T as [56] E c (T ) = 16π 3 where L is the latent heat. It is given by L = T c (p s − p b ) with p s and p b the fluid pressure in the symmetric and in the broken phase, resp. Hence, the energy E c diverges at the critical temperature.
Note the different uses of the word "critical" in critical bubble, with radius R c and energy E c , and the critical temperature T c .

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 [φ]. In the case at hand the two previously mentioned classical solutions to the equations of motion, the trivial solution φ(r) = 0 and the critical bubble solution φ(r) =φ(r), extremise E T [φ]. 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 φ = 0 and φ = φ m , we now evaluate the partition function at two classical solutions out of the field configuration space of φ. So the partition function will be a sum of two terms Z 0 and Z 1 evaluated around the extrema φ = 0 and the critical bubble φ =φ, respectively.
Let us first expand around the classical solutionφ as φ =φ + δφ. Performing the actual saddle-point approximation we find To evaluate the functional integral it is useful to diagonalise the operator − ∇ 2 + V T (φ). In order to do so we consider the following eigenvalue equation with δφ = α c α ψ α , where the ψ α are normalized so that β d 3 x|ψ α | 2 = 1. The measure then takes the following form For the non-zero eigenvalues, we have to perform Gaussian integrals, resulting in 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φ, given by Taking the spatial gradient yields These are the translational zero modes, as they correspond to the change in the field when performing a small translation For these modes, the normalised integration variables are from which it follows that where V is the volume of the system. Furthermore, also the trivial solution φ = 0 is an extremum. Therefore, if we solve Eq. (6.15) at φ = 0 we obtain The solution φ = 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, Z Z 0 + 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 λ − , linked to an unstable mode (i.e. expansion or contraction of the bubble).
We can already see that the ratio Z 1 /Z 0 contains an exponentially small factor exp(−β E c ). Therefore, we can expand the logarithm in the free energy and approximate its imaginary part as 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 V is then given by where we have introduced a standard notation for the eigenvalue products. The symbol det 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 V T (0) 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 V T (0) ∼ g 2 T at T 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 λ 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 γ −1 ∼ g 4 ln(1/g) [59], bringing an extra factor of order g 2 ln(1/g).

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 w . Having good knowledge of the wall speed is important, as the flow setup 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 α n , introduced in Section 7.3. In this section we will see how the flow around an expanding bubble depends only on v w and α 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].

Wall speed
For a bubble to expand, the interior pressure must be larger than the exterior one. The pressure is given by where g eff are the effective relativistic degrees of freedom. The bubble expands if , 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 φ driving the phase transition interacting with the plasma of particles is given by where η T =ηdm 2 /dφ, 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 µ (z), and the field profile is φ = φ(z). Let us also assume that the fluid speed changes little as it moves across the wall, so that We can then rewrite Eq. (7.2) as Multiplying both sides of the above equation with ∂ z φ and integrating over z, we obtain After performing the integral on the left hand side of the equation we find  where ∆V T = V T (z = +∞)−V T (z = −∞) is the pressure difference across the wall. To obtain the wall speed, one needs to solve this equation for v w . In practice, it is difficult to obtain a solution. In some cases a constant-v w solution may not even exist, in which case the bubble wall must continue to accelerate. The case where v w → 1 is called runaway solution; the wall speed becomes ultra-relativistic [62,63]. However, provided the product γ w η T does not decrease with γ w , a constant-v 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].
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 w the radius of a bubble at time t that nucleated at time t is given by R = v w (t − t ). Consequently, the volume V of this bubble is 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 multiplied by the number density of nucleated bubbles that nucleated during (t , t + d t ), see Fig. 10. The latter is given by dn(t ) = Γ (t )dt /V, such that the fractional volume in the metastable phase is 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] 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 Γ (t )/V, which governs the progress of the phase transition. It is most sensitive to the dimensionless combination S = E c /T , which decreases rapidly from infinity for T < T c . Hence Γ (t )/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 V is given by 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 c = 100 GeV. We estimated previously that Γ (t )/V ∼ T 4 c exp(−E c /T c ), where we have dropped all dimensionless constants. The Friedmann equation tells us that the Hubble rate is H(T c ) 2 ∼ GT 4 . Hence the first bubble nucleates when S has dropped to 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 β as 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 Γ as with Γ f ≡ Γ (t f ). We are now able to perform a saddle-point approximation to the integral Eq. (7.9), as Using the definition of t f , we find that  the reference time t f admits an alternative interpretation. Namely, the rateḣ 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, d dt 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 time 5 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 β −1 . To draw the plot in Fig. 11 we use the realistic values of T H = 105.2 GeV and T f = 104.9 GeV for the temperatures and β = 8950H, and v w = 0.9. This implies that t H = 2.66 · 10 −10 s and t f = 2.72 · 10 −10 s such that t f − t H ∼ 10 −13 s. Shortly after the nucleation rateḣ is at its maximum at t f , essentially the entire universe is in the stable ground state.
The bubbles can nucleate only in the metastable phase. Hence, using Eq. (7.13) we can 5 In order to find a relation between t H and t f , we note that Eq. (7.13) reads Γ H = Γ f exp β(t H − t f ) . Using Eqs. (7.10) and (7.16), and solving for t H yields 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 β −1 . deduce that the bubble number density n bubble (t) in the universe is given by where we used Eq. (7.17) in the second line. At late times when t → ∞ and upon using Eq. (7.16), this expression simplifies to which is the final bubble density, and defines the mean bubble centre separation R * = (n bubble (∞)) −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].

Relativistic combustion
Now we will assume that the wall velocity v w is known, and discuss solutions to the hydrodynamic equations following Refs. [66,68]. Let u µ 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 w = 0, see Fig. 12, and introduce a subscript ± to distinguish between quantities evaluated just ahead of the wall (+) and just behind the wall (−). From the conservation of energymomentum at the wall it follows that (7.22) whereṽ ± are the fluid velocities in the rest frame of the wall andγ ± = (1 −ṽ 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ṽ For convenience, let us introduce the trace anomaly θ = 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 ∆θ = θ + −θ − . From here, one defines the transition strength α + and the enthalpy ratio r as In terms of these quantities, Eq. (7.23) can be rearranged as These equations can be solved forṽ + in terms ofṽ − (or vice versa). The result is which depends only onṽ − and α + but not on r. As we explain in more detail later, the upper sign is to be taken forṽ − > 1/ 3, while the lower sign is to be taken ifṽ − < 1/ 3 in order to ensure the solution to be physical. This implies that both velocities are either subsonic (deflagrations) or supersonic (detonations). Further,ṽ + must be positive because the fluid must flow through the wall from the outside to the inside of the bubble. This requires α + < 1/3. In Fig. 13,ṽ + as a function ofṽ − is shown for different values of α + . 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 µ = γ(−1, v) µ and the space-time orthonormal vectorū µ = γ(−v, v/v) µ . These vectors satisfyū µū µ = −1 and u µū µ = 0. Projecting the conservation equation yields (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 w t where we set the nucleation time to t = 0. Since there is no length scale involved in the problem, the differential equations should exhibit similarity solutions that depend on the dimensionless coordinate ξ = r/t only. We can write the fluid velocity as v = v(r, t) r = v(ξ) r, where r is a unit radial vector. The continuity equations can be rearranged to Here, c 2 s = dp/de is the speed of sound and is the fluid velocity at ξ in a frame that is moving outward at speed ξ. Solving these equations requires numerical techniques. In general, the speed of sound c 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 2 s = 1/3 everywhere. An example of such an equation of state is the "bag" model, where where a s , a b and V s are positive constants with a s > a b . The subscripts s and b refer to the symmetric (i.e., where φ = 0) and broken (i.e., where φ = φ 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, ξ w = v w , while the wall is assumed to be infinitesimally thin. The boundary conditions at the wall read v → v ± = µ(ξ w ,ṽ ± ) as ξ → ξ ± w where ξ ± w = v w ± δ and δ infinitesimally small. On the other hand, the fluid velocity must vanish, v = 0, at the center of the bubble ξ = 0 due to spherical symmetry and at ξ = 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 ξ > c s and µ(ξ, v) > c s implying that dv/dξ > 0 and then integrating backwards in ξ. 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 w < 1/ 3, and the fluid is at rest everywhere inside the bubble, v − = 0. In the wall rest frame we thus findṽ − = v w . The fluid velocity just ahead the wall in the universe frame is v + = µ(ṽ + (α + , v w ), ξ 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ṽ − > 1/ 3, with the fluid being at rest everywhere outside the bubble, v + = 0. In the wall rest frame we thus haveṽ + = v w . The condition onṽ − means there is a minimum forṽ + called the Chapman-Jouguet speed v CJ . It is given by v CJ (α + ) = 1 3 1 + α + + 3α 2 In order to ensure thatṽ − (α + , v w ) > c s and hence dv/dξ > 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 ξ = c 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 w > 1/ 3, which occurs forṽ − = 1/ 3. A physical solution forṽ + (α + , 1/ 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 α + 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 α + which matches to the asymptotic energy and pressure. For more details on these solutions we refer to Refs. [66,[68][69][70].

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 [66,71]. 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 with 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 where V is the averaging volume, andē 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 where e 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 b = e 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 (7.38) which can be related to the kinetic energy fraction via where Γ =w/ē is the adiabatic index of the fluid in the broken phase.
In the previous section we introduced the transition strength parameter α + , 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 α + 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 φ = 0) with a subscript s while variables that correspond the broken phase (i.e., where φ = φ b ) with a subscript b. In terms of these background variables, one defines the transition strength parameter as It is convenient to define the efficiency factor κ that quantifies how much of the available energy is converted into kinetic energy. The available energy is determined by the trace anomaly θ introduced in Eq. (7.24). It is related to the potential energy of the scalar field V T (φ) as This implies that not all of the potential energy can be converted into kinetic energy and heat. Therefore, one defines the efficiency factor as 6 where ∆θ 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 = κα n 1 + α n + δ n , (7.43) where δ n = 4θ b /(3w s ). The efficiency factor can be computed numerically as a function of α n and ξ w as in Ref. [68]. To give an explicit example, in the case of detonations and using the bag model, where θ b = 0, they find the approximate relation κ α n 0.73 + 0.083 α n + α n (7.44) for ξ w → 1. For α n < 1, one can generally take κ ∼ α n , except for low wall speeds, and wall speeds near the Chapman-Jouguet speed, where the parametric dependence is closer to κ ∼ α 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].

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 zdirection only. Let us write the perturbed energy density as e =ē + δe and the perturbed pressure as p =p + δp with {δe, δp, v z } 1 with v z the fluid velocity perturbation. The components of the fluid energy-momentum tensor simplify to while the other components vanish identically. From the t-component of the energy-momentum conservation equation we find ∂ t δe +w∂ z v z = 0 (7.46) and from the z-componentw Note that both δe and δp depend on temperature as δp = ( Therefore, the tand z-components of the conservation equation can be combined to These equations describe sound waves that travel at the speed c 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.

Gravitational Waves
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.

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 ± 6M and M 2 = 29 ± 4M produced the signal. According to numerical simulations of such an event, an energy equivalent of M gw 3M must have been radiated away as gravitational waves. 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 gw . It was found that gravitational waves travel at speed of light to a large precision, 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 ∼ H with H =ȧ/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.
Before going on, let us be a bit more quantitative in our outline of gravitational waves. We start by considering perturbations about Minkowski background η µν as g µν = η µν + h µν , (8.2) where h µν 1 is the metric perturbation. A gravitational wave h i j (t, x) is a propagating mode of the transverse (∂ i h i j = 0) and traceless (h ii = 0) part of the metric perturbation, satisfying the equationḧ where T TT i j is the transverse and traceless part of the energy-momentum tensor, and G is Newton's constant. The constraints on h i j reduce the number of physical propagating modes to two, called 'plus' (+) and 'cross' (×) polarisations.
Provided their amplitude is small, gravitational waves can themselves be a source of energymomentum, described by the tensor [78] We aim at computing the power spectrum of the energy density of gravitational waves. From the t t-component of the energy-momentum tensor, we find the energy density, After a spatial Fourier transform of the metric, one can consider the contribution of a frequency interval d 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π. More often, one studies the energy density per logarithmic frequency interval d ln f , which can be written A convenient measure in cosmology is the fractional density in gravitational waves where ρ tot is the total energy density of the universe, related to the Hubble rate H by the Friedmann equation, With these quantities, one can define the fractional gravitational wave energy density per logarithmic frequency, dΩ gw often referred to as the gravitational wave (power) spectrum. The characteristic strain and the fractional density in gravitational waves are related through 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 i j 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 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 whose Fourier transform is We define the one-sided strain power spectral density of the signal in the interferometer S I h ( f ) as 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 R I ( f ), (8.14) Figure 18: Visualization of the orbit of LISA. LISA is a space-based laser interferometer designed to measure gravitational waves with frequencies of about 10 −2 Hz. The interferometer consists of three satellites that are orbiting the Earth in a triangle. Figure taken from Ref. [9].
For an interferometer with perpendicular arms, whose arm length is much less than the wavelength so that f L → 0, R I ( f ) → 1/5 [83]. For gravitational waves sourced by phase transitions, the relevant mission is Laser Interferometer Space Antenna (LISA). LISA is a space-based interferometric gravitational wave detector that works with three satellites orbiting the Earth, see Fig. 18. These are equipped with lasers and photodetectors such as to detect small changes in the separation of the satellites by measuring the time delays in signals sent between them, through the technique of time delay interferometry [84][85][86]. The interferometer arms have a length of 2.5 Gm such that LISA is most sensitive at frequencies in the range 10 −3 -10 −2 Hz. The target sources are binary white dwarfs, merging supermassive black holes at the centres of distant galaxies, and early-universe physics at the TeV-scale such as first-order phase transitions. The planned launch year is 2034. In Fig. 19 we show the projected sensitivity of LISA and the amplitudes of expected astrophysical sources. In addition we have added the signal from a first-order phase transition, which we discuss further below. As can be seen, LISA is expected to measure gravitational waves coming from phase transitions, but other sources produce a signal in the same frequency interval as well. For details on LISA, the expected sources of gravitational waves, and its relevance for probing fundamental physics we recommend Refs. [9,12,87]. Other planned space-based gravitational missions include DECIGO [88], Taiji [89] and TianQin [90].

GWs from first-order phase transitions
In this section we discuss the expected gravitational waves signal produced by a first-order phase transition. The details and derivation are quite involved and require numerical techniques which is beyond the scope of this course. Instead, we explain the key quantities that can be inferred from the gravitational wave power spectrum and their relation to the results of the previous sections. For a review we refer to Ref. [12] that discuss the predicted power spectra for various BSMs.
The general picture is the following. Initially, when temperatures are high, T T c , the thermal Higgs potential (free energy) has a minimum only at φ = 0. We recall that minima Figure 19: The expected sensitivity of LISA and possible astrophysical sources of gravitational waves. The phase transition signal, indicated by the blue line, is predicted for a transition with nucleation temperature T n = 500 GeV, strength parameter α = 0.3, transition rate to Hubble rate ratio β/H n = 100, and wall speed v w = 0.4. Background figure taken from [9].
of the free energy correspond to equilibrium states. Therefore, the Higgs is everywhere in the symmetric phase with thermal fluctuations around the minimum. As the temperature drops, the thermal Higgs potential develops a second minimum and the minima are separated by a potential barrier. At a critical temperature T = T c , both minima are degenerate. Below the critical temperature, it is thermodynamically preferred for the Higgs to occupy the new minimum. However, it has to cross the potential barrier. This it can do either by quantum tunneling [49,91] or by thermal fluctuations [50]. Here, we will study the thermal case.
Thermal fluctuations can drive the Higgs into the symmetry-breaking phase in limited regions of space. The most probable shape of the regions is a spherical bubble, inside which the Higgs is already in the new phase, while outside the bubbles the Higgs is still in the metastable phase. Large enough bubbles expand into the fluid such that the rest of the universe enters the symmetry-broken phase.
We have already discussed the rate at which bubbles appear, and how they expand, in sections 6 and 7. Once the bubbles start colliding, we must rely on a mixture of numerical simulations (see e.g. Ref. [71]) and modelling to describe the motion of the fluid and the production of gravitational waves.
We recall that near equilibrium, the system can be described as an ultrarelativistic fluid coupled to a Higgs field, according to equations (5.26) and (5.24). For concreteness, we state the equations of motion as used in Ref. [71] ready for numerical simulations. As before, the fluid four-velocity is given by u µ = γ(−1, v) µ . The Higgs field satisfies the relativistic equation which follow from Eq. (5.24) together with Eq. (5.25). Here, E = γe is the fluid energy density and Z i = γ(e + p)u i the fluid momentum density with u i = γv i . The uncertainty in the precise form of η turns out not to matter much: different choices affect only the precise profile of the domain wall as it moves through the fluid, and the hydrodynamics ensures that the fluid flows around expanding bubbles are the same for a given wall speed v w and phase transition strength α n . Finally, the metric perturbation is described by Eq. (8.3). In practice, one evolves an unconstrained symmetric tensor u i j with the tensor source (8.16) and projects out the transverse traceless part h i j (t, k) from the Fourier transform u kl (t, k) with the appropriate projector. In the bubbles of a thermal transition, the scalar field part of the source is confined to a microscopically thin wall, while the fluid source is distributed over a significant volume. The ratio of the energies can be estimated as w /R, where w is the wall thickness, and R is the bubble radius. The wall thickness w is very roughly the inverse temperature, and given that the total energy density is e ∼ T 4 , the ratio can be estimated as w In general, Ω sw is thought to yield the dominant contribution, unless the transition is very strong, or the bubbles are as large as the cosmological horizon.
The key quantities which determine these power spectra are the following: • T n , the cosmological phase transition nucleation temperature; • α n , the phase transition strength parameter at the nucleation temperature, which is related to the scalar potential energy released during the phase transition and given in Eq. (7.40); • v w , the bubble wall speed.
• β, the transition rate parameter (7.12), which can be thought of as the inverse phase transition duration, and in combination with v w determines the mean bubble centre separation R * through Eq. (7.21).
The belief is that only these four parameters determine the power spectrum, and therefore that the gravitational wave power spectrum is potentially a precise probe of these quantities, which are in principle computable from the Lagrangian of a specific theory. 7 Hence, there is great excitement about the potential of LISA to probe physics beyond the Standard Model 8 .
One can estimate that the gravitational wave density parameter after a phase transition is [71,96] Ω gw ∼ (H n τ v )(H n τ ac )K 2 , (8.18) where τ v is the lifetime of the shear stresses, τ ac is the autocorrelation time of the fluid flow, and K is the kinetic energy fraction, defined in Eq. (7.36). The timescales depend on the properties of the fluid flow through its characteristic length scale L f and a characteristic flow speed V f K, which both ultimately depend on the parameters R * , α n and v w . The Hubble time H −1 n is a maximum time scale, as the expansion of the universe reduces the factor H 2 n in (8.18) and therefore the effective power of the source. Initially, one can estimate the characteristic length scale as L f ∼ R * . The kinetic energy fraction K(v w , α n ) of an individual bubble is a good estimate of the global kinetic energy fraction [71,96], unless the transition is strong and proceeds by deflagrations [97]. In that case, there is significant suppression of the fluid kinetic energy, which is associated with the slowing of the bubble walls when they encounter a hot region around another bubble.
For sound waves the characteristic speed is the sound speed c s 1/ 3, and the autocorrelation time is therefore τ ac ∼ R * /c s . The decay of a fluid flow is a non-linear process. From the non-linear terms in the fluid equations (8.16) one can form the non-linearity timescale The approximate lifetime of the gravitational wave source is then 7 Recent work has demonstrated that a fifth parameter, the sound speed, should be added to the list [69,94]. Holographic calculations indicate that at phase transitions in strongly coupled theories the sound speed can be quite different from 1/ 3 [95] 8 While the four parameters listed above allow one to fit the gravitational wave spectra obtained from numerical simulations, it has been found that two of these parameters, the strength of the phase transition α, and the inverse phase transition duration β, are not entirely independent from each other [72,73]: Stronger phase transitions, characterized by larger values of α, take more time to complete and hence also predict a smaller value of β. Conversely, weaker phase transitions, predicting smaller values of α, proceed rather quickly. In simple words, the underlying reason for this relationship is rooted in the fact that the strength of the phase transition is related to the potential difference between the metastable and the stable ground state, ∆V . As the temperature drops, the difference between the metastable and the stable ground state becomes larger. Strong phase transitions are a result of a sufficient amount of supercooling, i.e., the phase transition proceeds well below the critical temperature, where the two minima are degenerate, and hence take more time to complete. A more careful calculation of gravitational wave production in an expanding universe [98] shows that a better approximation for the effective source lifetime is (8.19) which assumes that the RMS fluid velocity stays constant until τ nl , when it immediately vanishes. Hence we can estimate which goes as (H n R * ) 2 K 3/2 for fluid flow lifetimes much less than the Hubble time. We recall from Eq. (7.43) that the kinetic energy fraction is where κ(α n , v w ) is the efficiency with which Higgs potential energy is turned into kinetic energy, discussed around Eq. (7.44).
A more detailed survey of the expected lifetimes and correlation times of the various sources can be found in Ref. [12]. It is important to note that numerical simulations and modelling show that the dimensionless constant of proportionality in Eq. (8.20) is O(10 −2 ) and not of order unity [71].
The shape of the gravitational wave spectrum is best understood for acoustic production, which has been explored by numerical simulations [71] and are now described by a precise theoretical framework, the sound shell model [66]. This is only the most recent development in a long history of modelling of gravitational wave production at a thermal phase transition [68,[99][100][101][102].
As sketched in Fig. 20 (left) the power rises as k 9 to a domed peak at kR * ∼ 10, and then decreases as k −3 . If the wall speed is close to the speed of sound, v w ∼ c s , then the domed peak becomes a slow k 1 increase towards a peak wavenumber determined by the sound shell thickness ∆R * . At very low wavenumbers, where the signal is generically small, different power laws may emerge [103,104]. For wall speeds away from the sound speed, the shape of the domed peak can be approximated by a simpler broken power law rising as k 3 and falling as k −4 [71].
The turbulent stage, however, is far less well understood, as simulations of turbulent flows are very challenging [97,107], and consequently the modelling [102,105,106] relies on untested assumptions. The main assumptions that go into the modelling are how the flow autocorrelation time depends on wavenumber, and also how the global quantities such as kinetic energy and correlation length evolve with time. The differences between these assumptions explain the different predictions for the power laws at low and high wavenumber, sketched in Fig. 20 (right). Only dedicated numerical simulations can resolve these issue. The simulations of turbulent vortical flow [107] seem to disagree with all the predictions, having a slow rise of around k 1 to the peak, and a k −8/3 decay at high wavenumber.
Other unresolved questions are how the vortical component of the flow, which leads to turbulence, is generated in the first place, how important it is relative to the the compressional part, and how efficient they are at generating gravitational waves. On the first question, recent simulations [97] show that vorticity is generated by bubble collisions in strong phase transitions. For example, in the simulation shown in Fig. 21), the vorticity is associated with regions where bubbles collide, and about 20% of the fluid kinetic energy is in vortical flow. The figure also shows the above-mentioned heating of the symmetric phase, here confined to shrinking asymmetrical droplets. There is clearly much work to be done before we are able to accurately compute a gravitational wave power spectrum for a phase transition of any strength. At the moment, this means regions where simulations have been carried out [71] and models [66] have reasonable agreement in the range 0.4 v w 0.9 and α n 0.1.
For ultra-strong transitions (α n 1), the energy density of the universe becomes dominated by the potential energy of the scalar field, and the expansion of the universe starts to accelerate. This is the original idea behind cosmic inflation [4], and it was quickly realised that the nucleation rate parameter cannot be much smaller than the Hubble rate, as otherwise the majority of the Universe would still be inflating [65,108]. Close to this boundary, the bubble separation R * would be of order the Hubble length, and so the gravitational signal would be maximal, motivating recent study [92,109,110]. In particular, bubble dynamics in the near-vacuum state is rather different from the thermal energy dominated case we consider in these lectures [92]. Recent numerical simulations involving just the scalar field are helping to develop a picture of what happens when the fluid plays no role at all [55,93,111]. One old idea which can be tested at the same time is that vacuum bubble collisions could have generated primordial black holes [108,[111][112][113].
The gravitational wave power spectrum at much larger scales than R * is also not fully understood. There are predicted to be some characteristic power laws in the gravitational wave power spectrum at long wavelengths connected with expanding shells of shear stress [103,104,114,115], which in the strong supercooling case would not be dominated by the sound shell model signal. The fluid flow length scale is also expected to grow due to nonlinearities in the fluid equations, which may also contribute to the large-scale power [116].

Comparison with GW observations
Finally for this section, we see how calculations of the gravitational wave signal match up to the potential for observations. First, we study how the gravitational wave frequency and intensity change between their generation and their observation.
The frequency scale of the spectrum when it is produced can be taken to be f * = R −1 * , which is also an estimate of the peak frequency. The frequency today f * ,0 is determined by the redshifting of the gravitational waves since the time of the phase transition. Using the Fried- Hence LISA will be most sensitive to electroweak-scale (100 -1000 GeV) transitions with bubble separations between 10 −2 and 10 −3 of the Hubble length. The intensity of the gravitational wave signal is also affected by the expansion of the Universe. While the energy density of the Universe is dominated by radiation, the density parameter of gravitational waves Ω gw remains constant, as the energy densities of both components redshift the same way. However, once radiation gives way to non-relativistic matter as the dominant component, the density parameter of gravitational waves Ω gw decreases, until today it is reduced by a factor F gw,0 = (3.57 ± 0.05) × 10 −5 100 g eff 1 3 , (8.23) where the uncertainty comes mainly from the measurement of the Hubble rate today, taken here to be the Planck 2015 best-fit value H 0 = 67.8 ± 0.9 km s −1 Mpc −1 [117].
We can now present a simple model for the gravitational wave power spectrum from firstorder phase transitions, which captures the parametric understanding of the power outlined in this section, gives a simple functional form based on numerical simulations [12,71], and includes an explicit attenuation factor due to the decay of the flow [98]. The contribution to the fractional density in gravitational waves in a logarithmic frequency interval d ln f from a phase transition with mean bubble centre separation to Hubble length ratio H n R * generating a kinetic energy fraction K is dΩ gw,0 d ln f = 2.061F gw,0Ωgw (H * R * ) 2 K + H * R * K 2 C( f / f p,0 ) , (8.24) whereΩ gw = 0.012 is a numerically determined constant, representing the approximate efficiency with which kinetic energy is turned into gravitational waves, C(s) is the function C(s) = s 3 7 4 + 3s 2 7/2 , (8.25) describing the gravitational wave power around the peak, and f p,0 = 10 f * ,0 is an estimate of the peak frequency, again from the numerical simulations. The number 2.061 is an approximation to 3/ ∞ 0 d ln(s)C(s). In this very simple model the power at small s is over-estimated, and at large s under-estimated, compared to the sound shell model, and has significant deviations around the peak for wall speeds near the speed of sound. Its domain of reliability is α n 0.1, 0.4 v w 0.5 and v CJ (α n ) 0.9, where v CJ is the minimum speed of a detonation, given in Eq. (7.34). This represents the parameter range where gravitational wave power spectra have been extracted from numerical simulations. The model does not get the shape right for deflagrations with v w near or above the speed of sound (see Fig. 20).
We are now in a position to estimate the expected density fraction in gravitational waves from a phase transition, or equivalently the expected characteristic strain. As the power spectrum has a definite peak, the estimate represents both the total power and the peak power. Taking T n = 500 GeV, α n = 0.3, β/H n = 10 2 , and v w = 0.4, we find from Eqs. (8.22), (8.20), (8.21), and (7.44) that Ω gw,0 ∼ 10 −12 , (8.26) and that the peak frequency estimate is f * ,0 10 −3 Hz. At this frequency, the characteristic strain of a signal with a gravitational wave density parameter of order 10 −12 is h c ∼ 10 −21 , (8.27) where we have used Eq. (8.10). Estimates of the power spectrum, in the simplified broken power law approximation in Eq. (8.24), can be generated by the PTPlot tool [12,118]. In Fig. 19 we have plotted the signal from a phase transition predicted by the more refined sound shell model [66,119], with the parameters given above. While the strength of the transition is in the range where the amplitude and shape of the signal are merely indicative, the amplitude agrees with our order-of-magnitude estimate above, and such a transition would be detected very clearly by LISA, despite the competing astrophysical signals. The most significant of these will probably be the foreground from white dwarf binaries in our galaxy, of which there are thought to be tens of millions. This foreground will always be present, unlike the transient signals from merging supermassive black holes in distant galaxies. However, the anisotropy of the galactic white dwarf foreground should make it vary in strength throughout the year as LISA orbits the sun and the relative orientation of the plane of the satellite constellation and the galaxy changes. This will make it distinguishable from the isotropic background expected from phase transitions (see [83] for a review of gravitational wave signals and detection methods).

Summary
Gravitational waves are an extraordinary new tool for astronomy and cosmology, and can be used to directly observe processes which happened in the early universe. In particular, firstorder phase transitions in the early universe give rise to gravitational waves with characteristic power spectra. The source of the gravitational waves is the shear stresses in the fluid set up by the expansion and collision of bubbles of the low-temperature symmetry-broken phase during the phase transition. The shear stresses are initially overlapping pressure waves, or sound, which may be strong enough to generate turbulence when they collide.
The dynamics of the phase transition can be described in terms of only four properties of the phase transition: the bubble nucleation temperature T n , the strength parameter α n , the transition rate parameter β, and the bubble wall speed v w . These are all in principle computable from an underlying particle physics model, and in these lectures we have gained some insight into the methods. We have also seen how relativistic combustion theory explains how kinetic energy is generated in the fluid, which eventually shows up as the sound waves and turbulent motion.
The current excitement is that information about these phase transition parameters is contained in the power spectrum of the gravitational waves, and that the space mission LISA could detect these gravitational waves if the phase transition took place at a time around t = 10 −11 s [12], when the temperature was in the range 100 GeV to 1 TeV. This is the scale of electroweak symmetry-breaking. Hence LISA is a particle physics experiment, as well as an astrophysical observatory. A first-order electroweak phase transition is possible only if there is physics beyond the Standard Model, perhaps in the form of extra Higgs fields. LISA has the potential to discover new fundamental physics. This is a dynamic and developing field. In the last section we already mentioned the current uncertainty about the hydrodynamics of strong phase transitions (α n 0.1). It is also not yet known how accurately LISA can measure the phase transition parameters, but even a 20% measurement of (say) the wall speed is not yet matched by the accuracy of the calculations from underlying particle physics models. Furthermore, the framework of homogeneous nucleation theory is taken from condensed matter, where nucleation by thermal fluctuation has never been observed to take place. There is a long-standing mystery about the lifetime of the metastable superfluid 3 He-A phase towards decay into 3 He-B. The theory described in these lectures predicts lifetimes of experimental samples of superfluid 3 He-A of order 10 1000000 seconds, while in the laboratory, even when taking great care over impurities, 3 He-A lasts only a few minutes (see e.g. [120]). There is much to be done to realise the goal of making LISA into a particle physics experiment to complement the Large Hadron Collider.