Field-theoretic derivation of bubble-wall force

We derive a general quantum field theoretic formula for the force acting on expanding bubbles of a first order phase transition in the early Universe setting. In the thermodynamic limit the force is proportional to the entropy increase across the bubble of active species that exert a force on the bubble interface. When local thermal equilibrium is attained, we find a strong friction force which grows as the Lorentz factor squared, such that the bubbles quickly reach stationary state and cannot run away. We also study an opposite case when scatterings are negligible across the wall (ballistic limit), finding that the force saturates for moderate Lorentz factors thus allowing for a runaway behavior. We apply our formalism to a massive real scalar field, the standard model and its simple portal extension. For completeness, we also present a derivation of the renormalized, one-loop, thermal energy-momentum tensor for the standard model and demonstrate its gauge independence.


Introduction
Modeling the bubble-wall velocity of expanding bubbles at a first order electroweak transition is of an essential importance for an accurate modeling of the gravitational wave production [2] and for baryogenesis at the electroweak scale [3,4]. In spite of a lot of progress, we still lack a reliable first principle calculation which -based on an (out-ofequilibrium) quantum field theoretic framework -provides a formula which can be used to get a quantitatively reliable information about the phase transition dynamics. In this paper we present a first principle (quantum field theoretic) derivation of the force acting on expanding bubbles that reached stationary state during a strong first order transition, providing thus an important step towards that noble goal.
To determine the dynamics of expanding bubbles one ought to know the friction force exerted on them by the plasma. On the other hand, the bubbles back-react on the plasma, thereby changing its properties. One way to find a general solution of this complex problem is to solve the Boltzmann equations for the relevant species in the presence of expanding bubbles [5][6][7][8], which is a formidable task. For that reason, in many papers the bubble velocity is treated as a free parameter whose value is assumed or roughly estimated.
An important question from the point of view of gravitational wave production is whether the bubble can run away, i.e. permanently accelerate, asymptotically reaching the speed of light. If such a situation is possible the latent heat of the transition is pumped into the scalar field, resulting in a characteristic gravitational wave spectrum with a very strong peak amplitude. An important step towards solving this puzzle was made by Bödeker and Moore in ref. [1], where they considered the bubble force in the relativistic limit and at the JHEP01(2021)070 leading order in the relevant coupling constants, following an earlier work by Arnold [9]. They found that, in the limit of a large Lorentz factor γ, the force does not depend on γ and thus concluded that the bubbles can run away if there is enough latent heat released. In their follow-up paper [10], the authors considered the next-to-leading order effects, finding a γ dependence in the friction force, which in principle precludes the runaway scenario, but still allows for highly relativistic walls. Some of the relevant works that discuss how fast the bubbles can expand include [11][12][13][14]. 1 In particular, ref. [12] refines some of the arguments presented in [1,10].
The approach of the present paper is based on covariant conservation of the energymomentum tensor of the bubble-plasma system. When applied to the limit when an approximate Lorentz symmetry 2 holds, we find that the friction force scales as (γ 2 − 1). Therefore, for the most of physically interesting cases, the wall will reach stationary state for moderate velocities and thus will not run away. This is the case when scatterings are efficient and local thermal equilibrium is enforced. If, on the other hand, scatterings are inefficient, a ballistic approximation [8] better describes the actual situation. In that case the bubble force saturates and there is no obstacle for the bubbles to run away. In this work our formalism is applied to a toy model with one real scalar field in a heat bath, as well as to a model with the standard-model field content featuring a first order phase transition and to its simple extension.
The paper is organized as follows. In section 2 we derive our main formula for the bubble force and discuss its applicability and possible generalizations. In sections 3 and 5 we show how to apply our formalism to a real scalar field and to the standard model and its simple portal model extension. In section 4 we compare our formalism with the literature, and in section 6 we conclude. In an extensive appendix we calculate the oneloop energy-momentum tensor of the standard model in local thermal equilibrium. A particular attention is devoted to renormalization and to showing gauge independence of the renormalized energy-momentum tensor.

The bubble force from the renormalized energy-momentum tensor
When taken together, the energy-momentum tensor of the plasma T p µν and expanding bubbles T b µν must be covariantly conserved, For simplicity we shall assume that the bubbles are large and nearly spherical, such that their front can be approximated by nearly planar walls propagating through the plasma. Because the surface tension of the bubbles is typically quite large, this approximation is justified. Moreover, we shall assume that the following hierarchy of scales holds, L, D, τ int = 1/Γ R H , where L is the bubble wall thickness, D is the relevant diffusion time (recall that JHEP01(2021)070 c = 1), Γ is the scattering rate of the relevant plasma species, and R H = 1/H is the Hubble radius characterizing the expansion of the Universe. Then the covariant derivative in (2.1) can be approximated by an ordinary derivative and the expansion of the Universe can be, to the leading order in adiabatic expansion, encoded by the temperature dependence on the scale factor of the Universe a(t). Taking account of these, equation (2.1) simplifies to, Several remarks are in order: • The energy-momentum tensor in eqs. (2.2)-(2.3) is a composite operator which diverges, so to make it finite it ought to be regularized and renormalized; • The expectation values in (2.3) ought to be calculated in a state that includes the bubble, or which approximates it well enough; • Eqs. (2.2)-(2.3) simplify further in the bubble frame, in which the terms containing time derivatives drop out. This applies when the bubble reaches stationary state and that is where our analysis holds. The regime of nonstationary bubbles is worth investigating in a separate work.
In the bubble frame, one can then integrate equations (2.2)-(2.3) across the bubble to obtain,

4)
∆ T p zz + ∆ T b zz = 0 , (2.5) where ∆ T p,b µν denotes the change of the µν components of the energy-momentum tensor of the plasma or bubble across the bubble. In what follows we focus on eq. (2.5) since the bubble-wall speed is determined by the balance of the two terms in (2.5) which encapsulate the driving force of the vacuum energy and the friction force from interactions with the plasma.
When calculating the plasma contribution in (2.5), it is convenient to calculate it in the plasma frame, in which (T 00 ) plasma ≡ ρ p and (T zz ) plasma ≡ P p . As T µν is a tensor, Lorentz boosting it to the bubble frame results in, T p zz = (γ 2 − 1)(ρ p + P p ) + P p , 3 where in the last relation we used ρ b + P b = 0 (since the bubble possesses no entropy). Taking account of the thermodynamic relation for the entropy density s, we arrive at the following expression for the friction force F per volume V on an expanding bubble, F V ≡ −∆P = (γ 2 − 1)T ∆s , (2.7)

JHEP01(2021)070
where we used, P = P p + P b and ∆P b = ∆ T b zz . The relation (2.7) tells us that the change in the pressure is balanced by the change in the entropy density of the plasma across the bubble and that the effect grows quadratically with the Lorentz factor γ = 1/ √ 1 − v 2 , unless some compensating γ dependent terms arise from the plasma frame ∆s, which we discuss later. Since in the derivation of eq. (2.7) we assumed stationary bubbles, making use of (2.7) one can determine the terminal bubble speed. This simple observation constitutes one of the principal results of this work.
Equation (2.4) in the plasma frame reads ∆(γ 2 vT s) = 0 and thus implies that (with nonzero variation of entropy across the bubble) the temperature or velocity of the plasma must change across the bubble. Here we focus on the bubble-wall force, eq. (2.7), therefore studying the effect of eq. (2.4) is beyond the scope of the present work. 4 Before we proceed to applications, in what follows we discuss applicability of formula (2.7).

1.
To arrive at (2.7) we took an expectation value of the energy-momentum tensor. This can be exacted by making use of the full quantum formalism, which gives accurate answers, but it is hard to implement because it involves quantum field theory in an out-of-equilibrium setting. A much simpler procedure is to take a semiclassical limit, according to which the plasma can be described as a collection of quasiparticles that remain approximately onshell across the bubble interface. To estimate when the quasiparticle approximation is reliable, recall that according to the uncertainity principle particles can be off-shell for short periods of time satisfying, ∆t 1/∆E ∼ 1/E, where E = p 2 + m 2 is the energy of the particle. Conversely, when the wall passage time ∼ L/(γv) is longer than the bubble wall thickness, L/γ (which gets Lorentz-contracted in the plasma frame), then particles will be approximately on-shell, i.e.

γv < γ < LT
(on − shell condition) , (2.8) where we assumed that E ∼ T . 5 A typical bubble wall thickness L at the electroweak transition is of the order L ∼ 10/T [8], implying that when γ 10 the plasma quasiparticles will be on-shell and the semiclassical kinetic formalism applies. When the on-shell condition (2.8) is satisfied, that does not yet mean that local thermal equilibrium is reached. In fact, an additional condition -efficient scatterings on the bubble interface -must be met for a local thermal equilibrium to be reached. We elaborate more on that below. The adiabatic approximation we employ here can be regarded as a semiclassical approximation, in which the effects of varying backgrounds are modeled by mass insertions along a particle's trajectory, as illustrated in figure 1. This is a kinematic effect enforced by the energy conservation in the bubble frame, and the approximation holds as long as the quantum off-shell effects are small. On top of this, particles can interact and scatter off each other. Even though scattering effects can be important for a complete understanding of the phase transition dynamics, we postpone their study for future work since they require the inclusion of higher-loop effects. In this work we investigate only two limits, namely very rapid and very slow scatterings. When bubbles are very fast (or very thin) and the criterion (2.8) is violated, quantum off-shell effects can become significant. To capture that, when constructing the energymomentum tensor, one ought to use exact mode functions in the moving bubble background, which is a much harder endeavor. For an example of how that can be done for fermionic fields and without loop effects, see e.g. ref. [17]. To consistently include the quantum loop effects on top of such a tree level treatment, one would have to solve the corresponding out-of-equilibrium problem using a perturbative quantum field theoretic framework such as the Schwinger-Keldysh formalism [18][19][20]. While progress has been made in applying such a formalism in baryogenesis/leptogenesis scenarios [21][22][23][24] (for reviews see [4,25]) and in cosmic inflation [26][27][28][29][30], little or no progress has been made in studying cosmological phase transitions.

2.
While our conclusions based on the consideration of the energy-momentum tensor operator are general, the form of eq. (2.7) is based on the assumption that the (expectation value of the) energy-momentum tensor is well approximated by a perfect fluid form, T µν = (ρ + P)u µ u ν + g µν P. To get a better understanding of the limitations of this approximation, recall that the perfect fluid form can be viewed as the leading order approximation in a gradient expansion. Including the first order gradient corrections yields the well known expression, (2.10) Here η and ζ denote the shear and the bulk viscosity, respectively, and (αβ) means that the indices are symmetrized. Notice that τ µν is orthogonal to u µ , u µ τ µν = 0 = τ µν u ν , and that
To see that recall that in the plasma (bubble) frame, u µ = (1, 0, 0, 0) (u µ = (γv, 0, 0, v)), such that (if one neglects the expansion of the Universe) the covariant derivatives acting on u µ give zero. Therefore, one gets a nonvanishing contribution from the viscosity part of the energy-momentum tensor (2.10) only if both (a) plasma velocity is perturbed from its thermal equilibrium form and (b) the viscosity coefficients do not vanish. In a perturbative treatment the leading order contribution to the energy-momentum tensor comes from the one-loop approximation. Since the one-loop contribution is non-dissipative, the viscosities (whose nature is dissipative) acquire nonvanishing contributions only at two and higher loops. Therefore in weakly coupled theories the dominant contribution to the bubble force is captured by the one-loop calculation and the bubble dynamics can be obtained from eq. (2.7). That does not mean that an accurate answer for the bubble dynamics can be obtained just from a one-loop analysis, as higher loops may be essential for determining the accurate form of the state with respect to which one takes the expectation value of the energy-momentum tensor.

3.
Even though formula (2.7) was obtained based on non-equilibrium considerations, it has a deceptively similar form to the fundamental thermodynamic law, where E denotes the energy, S is the entropy, P is the pressure, V is the volume, µ the chemical potential and N the particle number of the system. To see that let us divide (2.11) by dV to obtain its local form, ρ = T s − P + µn , (2.12) where ρ = dE/dV is the energy density, s = dS/dV is the entropy density and n = dN/dV is the number density. Next, close enough to thermal equilibrium µ 0 6 and the contribution of µn can be neglected and one obtains a standard thermodynamic relation, that the change in energy density plus pressure across the bubble equals to the change in the entropy density times the temperature, ∆(ρ + P) = T ∆s, which was used in the derivation of eq. (2.7).
Let us now look more closely at how our approach compares with the usual description according to which the transition dynamics is governed by the latent heat release = ∆ρ and by the change in the effective potential across the bubble, see e.g. ref. [8]. Our approach here is instead based on the change in the entropy density across the bubble, ∆s = ∆ρ + ∆P = + ∆P, which naturally arises from the energy-momentum conservation. In the one-loop approximation and in local thermal equilibrium, ∆P = −∆V eff , where V eff is the one-loop effective potential. 7 Therefore, we have ∆s = − ∆V eff . It would be worth investigating whether the two approaches are equivalent in more general situations. 6 The condition µ 0 does not mean that eq. (2.7) does not apply away from thermal equilibrium in which chemical potentials are appreciable. It just means that we have subsumed all the relevant effects that contribute to the bubble friction into the entropy increase. 7 To see that, note that from eq. (3.18) we can read off the thermal contribution to the pressure, P = −(6π 2 β 4 ) −1 [z −1 ∂zJB(6, z)] z=βm , where JB(6, z) is given in eq. (3.11). By a suitable partial integration one can show, z −1 ∂zJB(6, z) = −3JB(4, z), from which we conclude, P = (2π 2 β 4 ) −1 JB(4, βm) = −V eff , which completes the proof.

JHEP01(2021)070
The usual intuition from statistical physics tells us that, in the limit of a thick bubble interface, entropy does not increase across the interface if local thermal equilibrium is maintained, i.e. plasma particles crossing a thick interface constitutes an isoentropic process. This is true only if a single particle species forms the plasma. However, in realistic applications many particle species are present. Then the system naturally splits into active species, which exert a significant force on the interface, and passive species, whose mass remains approximately constant or zero across the interface, such that they do not exert a significant force on the interface. The passive species form a heat reservoir. The standard thermodynamic picture then applies. The entropy density in the active part of the system reduces significantly, thus exerting the force on the bubble. The heat reservoir absorbs heat, thus heating up. If the heat capacity of the reservoir is large enough, the total temperature of the system plus reservoir does not change much, and the local thermal equilibrium description, which assumes equal temperature on both sides of the interface, can be considered as the leading order approximation. In realistic situations, in which the standard model contains most of the degrees of freedom, the active particles are the top quark, the weak gauge bosons and the Higgs particle and they comprise about 20% of the relativistic degrees of freedom; the remaining 80% constitute a large heat reservoir.
In order to illustrate how to use (2.7) for the phase transition dynamics, in what follows we first consider a simple real scalar field model and then more general models. For the scalar theory we will analyze two opposite cases: the local thermal equilibrium, where we assume that the interactions in the wall are very efficient; and the ballistic approximation in which the interactions in the wall are inefficient. The force that we find in both cases will display a very distinct behavior.

Real scalar field
In this section we calculate the one-loop energy momentum tensor of a self-interacting, massless scalar field in the presence of an expanding spherical bubble at a first order phase transition. For simplicity we shall first consider the case in which the scalar is in a local thermal equilibrium (lte), which approximates well the scalar field state if the thermalization time scale, τ th = 1/Γ th is smaller than the time ∆t b = L/(γv) it takes the bubble to pass by an observer in the plasma frame, i.e.
Since Γ th is controlled by a coupling constant, which is typically smaller than one, τ th > 1/T , the condition (3.1) is more stringent than the on-shell condition (2.8). 8 If (3.1) is not met, then particles that move across the bubble partially thermalize, and in the extreme case do not thermalize at all, in which case a ballistic approximation applies. In what follows we firstly calculate the energy-momentum tensor by assuming local thermal equilibrium 8 We do not verify whether this condition is met for the toy model that we consider here, since in this model no phase transition is present and we treat it only as an illustration of the method. More realistic models are considered in section 5. and then discuss how our results are affected if ballistic approximation represents a more appropriate description. The free scalar field action is given by,

JHEP01(2021)070
where g µν is the inverse of the metric tensor g µν , the signature of the metric is mostly plus, g = det [g µν ] and m is the scalar mass. The origin of the mass can be either a tree-level mass, m 0 , or it can be generated by a field condensate. For example, adding a scalar self-interaction, to the Lagrangian in (3.2) will generate in the presence of a condensate φ = φ 0 a field dependent mass, 9 If the background field φ 0 = φ 0 (x) varies in spacetime, the mass (3.4) will follow the suit. As long as the variation is slow, it can be treated as an adiabatically varying quantity. Varying (3.2) with respect to g µν results in the energy-momentum tensor, Here we are interested in evaluating the expectation value of the one-loop energymomentum tensor, The contributing one-loop diagram is shown in figure 2, where the dashed line denotes the scalar propagator and the cross (×) indicates the energy-momentum tensor insertion.

Local thermal equilibrium approximation
Since at one-loop order the energy-momentum tensor (3.6) can be evaluated by making use of the field propagator, the next natural step is to construct the scalar field propagator, which is in general defined by, Hereρ denotes the density operator, which contains the information about the scalar field state and T denotes the time ordering operator. In local thermal equilibrium (3.1) holds and the density operator can be approximated by its local equilibrium form,ρ →ρ th = e −βĤ φ /Tr[e −βĤ φ ], whereĤ φ is the Hamiltonian operator.
In general the free scalar propagator (3.7) obeys the equation of motion, where m = m(φ 0 ) is the field-dependent scalar field's mass. The thermal propagator in (3.8) can be thought of as the inverse of the operator, Alembertian operator with thermal boundary conditions imposed. In practice this can be done by inserting the thermal density operatorρ th in (3.6). Furthermore, since time scales that govern the phase transition dynamics are typically much shorter than the Hubble time, the expansion of the Universe can be considered as adiabatic and in the leading order adiabatic approximation the d'Alembertian reduces to the usual wave operator, → ∂ 2 = η µν ∂ µ ∂ ν (the expansion of the Universe is then captured by taking T = T (t)). That means that for our purposes it suffices to determine the thermal propagator for a massive scalar in Minkowski space in slowly varying background fields (the metric tensor and the scalar field). That propagator is well known and in the plasma frame it reads, 10 where K ν (z) is the Macdonald function (the modified Bessel function of the second kind), E p = p 2 + m 2 and, is the invariant distance function on Minkowski space where the appropriate ı prescription for the Feynman propagator is also indicated. Note that in the free propagator (3.9) the vacuum (∝ K (D−2)/2 ) and thermal contributions neatly split and that the thermal contribution is finite, and thus can be evaluated in D = 4. While it is not generally possible to write in a closed form the thermal part of JHEP01(2021)070 the propagator in (3.9), its coincident and near coincident limits are possible to express in terms of the bosonic thermal integral, based on which one can obtain the thermal contribution to the one-loop effective potential of a bosonic degree of freedom in n space-time dimensions. For example, in n = 4 and for a particle of mass m, the one-loop thermal effective potential is V . With these remarks in mind we can write a closed form expression for the coincident propagator (3.9), such that the vacuum part is divergent in D = 4 and ought to be regularized. The vacuum contribution in (3.12) was evaluated by noting that the Bessel function K ν (z) in eq. (3.9) can be expanded around the lightcone as a sum of two series (ν = (D − 2)/2), and then used the fact that in dimensional regularization, by a clever use of the suitable analytic extension, one finds that the series with D dependent powers in (3.13) does not contribute at coincidence. Since all power-law divergences get subtracted in this way, this feature became known as the automatic subtraction. In order to evaluate (3.6) we need, where we introduced the usual T * time ordering which is defined to commute with the two external derivatives in (3.14). The vacuum part of (3.14) is −2η µν times the linear coefficient in ∆x 2 of the propagator in (3.9), which is easily extracted from the integer series in (3.13), When taken together with the vacuum part in (3.12) this then implies, T [L] vac = 0, such that, which has the form of a cosmological constant.

JHEP01(2021)070
The thermal contribution to (3.14) is given by the integral over p µ p ν of the coincident thermal integral in (3.9). Recalling that p 0 = E p , the 00 and ij contributions ought to be evaluated separately (the 0i contribution vanishes), such that T [L] th = 0, and the contributions in (3.18) are the thermal contributions to the one-loop stress-energy tensor.
To complete the calculation, we still ought to renormalize the vacuum contribution (3.16), which for a constant bare mass can be done by adding a cosmological constant counterterm. However, here we are primarily interested in a mass generated by a field condensate (3.4), and the suitable counterterm action is of the form, which contributes to the energy-momentum tensor as, We shall use the minimal subtraction scheme, and to that purpose expand (3.16) around D = 4, where µ is an arbitrary scale and m 2 = λφ 2 0 /2. Comparing (3.21) with (3.20) we see that, removes the divergence from (3.21), resulting in the following renormalized energymomentum tensor, From the point of view of the bubble force calculation, the change in the entropy density (2.6) across the bubble is what determines the friction force on the bubble interface, which is determined by the second line in (3.23) (Lorentz covariant contributions have a vanishing entropy density),   Inserting this into eq. (2.7) gives an expression for the bubble speed as a function of the change in the pressure across the bubble −∆P due to the bubble nucleation and the change of the plasma entropy density. Since −∆P depends on the amount of supercooling before bubbles start nucleating and on the detailed form of the effective potential, the toy model Lagrangian considered in this section cannot be used for a meaningful estimation of ∆P.
For that reason we shall not attempt to estimate it here, but instead we shall treat it as a free parameter of the transition.
In figure 3 we show how the entropy density changes across the bubble as a function of the scalar mass. The relativistic plasma limit, s = 2π 2 T 3 /45 (dashed) is reached in the limit when the mass in the broken phase m → ∞, because then the entropy density inside the bubble tends to zero. The main application of eq. (2.7) is to determine the terminal velocity (or Lorentz factor) of a bubble in stationary state. Figure 4 shows the bubble's Lorentz factor γ as a function of the scalar mass and the strength of the transition, expressed as ∆P across the bubble for a moderately strong transition, ∆P = −0.1/β 4 (solid black), for a strong transition ∆P = −0.5/β 4 (solid orange) and for a very weak transition ∆P = −0.01/β 4 (dashed green). For each choice of ∆P there is a minimum bubble Lorentz factor, The maximum γ is reached when m → 0, in which case the bubble runs away (since in that limit there is no force). However, for every m 2 > 0 and ∆P < 0, no matter how small they may be, a finite γ is reached. Figure 5 shows the bubble speed for the same choice of parameters as in figure 4. Notice that, independently of ∆P, all curves begin at v = c for m = 0. The horizontal dashed lines indicate the lowest attainable bubble speed for a given phase transition strength ∆P.
We see that in the case of local thermal equilibrium, the bubble velocity stabilizes at relatively low values. Therefore, even based on our analysis which is limited to moderate velocities with γ < 10, one can conclude that the bubbles cannot run away.

Ballistic approximation
If the thermalization rate is not large enough to satisfy eq. (3.1) then the ballistic approximation [8] is the more appropriate one to model the state of the field. In this case, one assumes that particles are in thermal equilibrium in the symmetric phase far in front of the expanding bubble, and they move across the wall so fast that they interact semiclassicaly, but the time is so short that they do not reach a local thermal equilibrium. To solve for the force acting on a bubble in this case, one can solve the Liouville equation in the bubble frame, By observing that, ∂ z = (∂ z m 2 )∂ m 2 , one sees that the general solution of (3.25) can be written as a general function of m 2 + p 2 z JHEP01(2021)070 which is equivalent to saying that E and p ⊥ are conserved in the bubble frame. On the other hand, in front of the bubble, where m = 0 (inside the bubble where m = 0), the solution is given by, respectively. Notice that the negative p z branch of the distribution function (3.28) exists only if thermalisation inside the bubble takes place. Since the bubbles at the electroweak transition grow large before they start colliding, this condition will be satisfied for typical bubbles. If not, the negative branch (3.28) will be absent. When the boundary condition (3.27)-(3.28) is imposed on (3.26) one obtains the general ballistic solution, which is conveniently broken into three parts as follows.
Case A. Transmission (t + ) from the symmetric phase: where we temporarily introduced an index 0 on m to denote the mass deep inside the bubble. Note that in the reflected solution (3.30) we included both the particle incoming on the wall (−m 0 < p z < 0) and the reflected particles (0 < p z < m 0 ), because they both contribute to the force on the wall. The ballistic solution (3.29)-(3.31) tells us that quasiparticles are still on-shell, but they are not thermally distributed. As a consequence of this departure from thermal equilibrium Lorentz symmetry is violated. This then introduces a dependence on the Lorentz factor γ JHEP01(2021)070 in the plasma frame quantities in eq. (2.7), making thus the γ dependence of the bubble force more complex than what eq. (2.7) suggests.
To obtain the bubble force one must calculate the change in the plasma pressure across the bubble interface. The pressure in the bubble frame reads, where the index z indicates that P z is the pressure in the direction in which the bubble expands. For simplicity in eq. (3.32) we took account of the plasma contribution only (since the vacuum contribution is unchanged). The form of the pressure in (3.32) can be obtained from the zz part of eq. (3.14), provided the thermal distribution function 1/(e βEp − 1) in (3.9) gets replaced by the ballistic distribution function (3.29)-(3.31). This ballistic approximation is valid as long as the semiclassical on-shell condition (2.8) holds true. Since the pressures in the other two spatial directions do not change by the expanding bubble, we do not need to consider them here. The quantity P z in (3.32) represents the total pressure in the z direction in the wall frame exerted by the plasma, and it is therefore equal to T p zz = (γ 2 − 1)T s + P p , where s and P p are the plasma frame quantities, and therefore -up to an unimportant term P p -it equals the differential bubble force in eq. (2.7). To evaluate the integrals in (3.32), it is convenient to first integrate over E, in which case the ranges of integration are, . By making use of EdE = p ⊥ dp ⊥ and integrating e.g. the transmitted contribution (3.29) over E one gets, The integral over p z cannot be done analytically. To simplify it, it is convenient to introduce dimensionless variables, M ± = βγ(1 ± v)m, x = βγ(1 ± v)p z , and the integral (3.33) becomes, . An analogous procedure gives for the reflected contribution (3.30) and for the transmitted contribution (3.31),

JHEP01(2021)070
The total pressure is then simply the sum of the three contributions, (3.37) and the pressure difference across the bubble reads There is a subtlety involved in evaluating (3.37) in the reflected contribution (3.35), which is simply equal to, ∆P r z = P r z . This is because the contribution of the reflected particles to ∆P r z is given by P r z in front of the wall, minus the pressure at the turning point. But this contribution is zero because the phase space at the turning point is zero (all particles at the turning point have p z = 0, such that the integral over p z vanishes). Upon combining all of the contributions (3.34), (3.35) and (3.36) one obtains for the pressure difference, 11 where Li s (z) = ∞ n=1 (z n /n s ) denotes the polylogarithm function and we dropped the index 0 on the mass. The first term in (3.39) is the pressure at the vanishing mass, which is precisely of the form, (γ 2 − 1)T s 0 + P 0 , with s 0 = (ρ 0 + P 0 )/T , ρ 0 and P 0 being the entropy density, energy density and pressure of an ultrarelativistic plasma with one degree of freedom. The third line in (3.39) comes from particles penetrating the interface from inside the bubble. As noted above, this contribution will be present only if thermalization inside the bubble takes place. Eq. (3.39) is to be compared with the pressure difference across the interface obtained assuming local thermal equilibrium (lte) close to the interface, To arrive at eq. (3.39) we made use of the integral,

JHEP01(2021)070
From eqs. (3.29)-(3.31) we see that, in the ballistic approximation, all particles ascending onto the interface (Cases A and B) slow down, such that f > f lte , implying that P z > P z,lte , from which we infer ∆P z < ∆P z,lte , meaning that the ballistic approximation yields a smaller bubble force and thus also faster bubbles. 12 This should not come as a surprise, since the ballistic approximation completely neglects particle interactions on the interface, it minimizes entropy production across the bubble, thus underestimating the bubble force.
To get an idea by how much, in figure 6 we plot ∆P z defined in (3.39) as a function of the Lorentz factor γ for m = 0.1T (orange dashed), m = 0.2T (black dashed), m = 0.4T (solid red), and m = 0.8T (solid blue). As expected, the force is larger for larger masses. Unlike in the case of the local thermal equilibrium force, the ballistic force saturates with γ. Interestingly, from the point of view of applicability range of our approach, the force saturates already for moderate values of γ < 10, suggesting the possibility of runaway in the ballistic limit. Motivated by that, we evaluate the maximum force reached when γ → ∞ from eqs. (3.39)-(3.40) as follows. In the limit when γ 1, 42) and the pressure difference ∆P z in eq. (3.39) can be estimated as, hence ∆Pz > ∆P z,lte , contribute sub-dominantly to the force. This must be the case because in the wall frame the distribution function of the descending particles is always suppressed when compared with that of the ascending particles, implying that they also sub-dominantly contribute to the bubble force. 13 To see this, let us rewrite the integral in the third line in (3.39) , where we used the variables x = γx, M = βm and we included the prefactor −1/γ 4 in the definition of the integral. Next, it is convenient to introduce Lorentz boost variables, upon which the integral naturally splits into two parts, where the last estimate holds in the limit γ → ∞. The second integral cannot be evaluated exactly, but it can be bounded from above by the simpler integral obtained by replacing √ x 2 + M 2 → x in the exponent, which can be evaluated, from where we conclude that It − is suppressed at least as ∼ e −γβm .  Notice that, if there are heavy particles in the plasma, m T , the bubble force is ∝ m 2 and thus gets saturated at much larger values and for large Lorentz factors, γ m/T . Namely, even though the number of heavy particles is exponentially suppressed in the plasma frame, their energy in the wall frame gets boosted by the Lorentz factor γ, thereby reducing their suppression. An important lesson to take from this analysis is that, if a phase transition is strong and bubbles are relativistic, then the existence of very heavy particles (with a mass m T ) can be of a crucial importance for the correct determination of the terminal bubble-wall velocity. In particular, heavy particles can determine whether the bubbles run away or not.
To illustrate how the bubble force depends on the particle mass, in figure 7 we plot ∆P z in (3.39) as a function of m/T for a fixed Lorentz factor γ. The values of γ (starting from bottom up) are γ = 1.1 (orange dashed), γ = 1.5 (black dashed), γ = 2 (solid red), and γ = 5 (solid blue), respectively. For small masses the force in the ballistic approximation increases rather slowly (left panel), but then it continues increasing and eventually saturates for very large masses. This is in contrast with the local thermal approximation result in figure 3, which shows that the lte force already saturates for quite modest masses. In order to get a better insight into how much the bubble force calculated in the ballistic approximation is smaller than in the local thermal equilibrium (lte), in figure 8 we plot ∆P z (3.39) as a function of m/T for a fixed Lorentz factor γ. In the left panel we show, from bottom up, γ = 1.1 (lte is solid blue curve, ballistic is dashed blue) and γ = 2 (lte is solid red curve, ballistic is dashed red), while in the right panel the force for γ = 5 (lte is solid blue curve, ballistic is dashed blue) and γ = 2 (lte is solid red curve, ballistic is dashed red) are shown. Notice that the force in the lte approximation is always larger, but for sufficiently large masses the two forces asymptote to the same value. This means that very massive particles (m T ) contribute more in the ballistic approximation than in the lte approximation. This observation can be important for the correct determination of the phase transition dynamics, particularly in systems with very heavy degrees of freedom. . Notice that the force calculated in the local thermal equilibrium approximation can be orders of magnitude larger than that calculated in the ballistic approximation. Notice further that, as γ increases the ratio of the two forces increases, which can be explained by recalling that, as γ increases, the ballistic force saturates.
To conclude, the analysis in this section shows that the entropy production generated by a passing bubble (and thus also the bubble force) minimizes in the limit when particle interactions are negligible and can be treated in a ballistic approximation, in which case the bubbles' runaway is expected. In the opposite limit, when particles scatter efficiently and can be consequently approximated by a local thermal equilibrium, the entropy production and the bubble force maximize. As long as the local thermal equilibrium approximation applies, entropy production grows as the Lorentz factor squared and quickly reach terminal velocity which can be estimated from the one-loop energy-momentum tensor. We emphasize that formula (2.7) derived in section 2 applies in both of those limits as well as in more general situations, provided one takes proper care of the state of the field. Furthermore, JHEP01(2021)070 both the ballistic and the local thermal equilibrium approximation can be understood as the leading order approximation in a more general treatment, in which they represent the leading (zeroth) order approximation to an expansion in powers of the perturbative parameter ΓL/v for the ballistic approximation and v/(ΓL) for the local thermal equilibrium approximation, where L is the thickness of the bubble interface, v its speed and Γ the relevant thermalization rate.

Comparison with existing results
Since the problem of bubble dynamics at first order transitions is long-standing and there exists rich literature of the topic, it is useful to clarify where we differ when compared with the existing work. Our principal statement is that the force on an expanding bubble can be extracted from the energy-momentum conservation law (2.3). When written in the bubblewall frame (2.5) it can be recast as (2.7), from which we get a differential force on the bubble, where ρ p and P p are plasma frame quantities and T p zz is calculated in the bubble frame. Note that P p is just T p zz in the plasma frame, such that eq. (4.1) indicates that the force on the bubble is present only when the bubble is moving.
The precise form of the energy-momentum tensor in (4.1) is complicated and it depends both on the state chosen and on the accuracy at whichT p zz (x) is evaluated. For simplicity of the argument here we work in the one-loop approximation and assume a state in which the bubble profile is an adiabatic function of z. Then for a real scalar field in thermal equilibrium (outside the bubble) one can use the thermal part of the massive scalar propagator (3.9) to calculate the force in (4.1). When integrated across the bubble, eq. (4.1) can then be recast as,

JHEP01(2021)070
where we made use of the bubble frame relations, and of the differential chain rule, d(p 2 z ) = dz∂ z (p 2 z ) = −dz(dφ 0 /dz)(dm 2 /dφ 0 ). If we compare eq. (3.3) from ref. [1] to the result in eq. (4.2) we see two unimportant differences. Firstly, the force in eq. (4.2) is two times larger, and secondly, eq. (3.3) in ref. [1] does not subtract the term needed to make the force vanish in the static limit. 14 Moreover, eq. (3.3) takes into account that the distribution function depends on the z coordinate and is different at each point. Once a simplifying assumption of local thermal equilibrium in front of the bubble is imposed as was done above, eq. (3.3) from ref. [1] and our result in eq. (4.2) do agree up to the minor differences already discussed.
The above analysis shows that, starting with our main result (2.7), and up to the well understood differences, our general quantum expression for the bubble force reduces to the standard semi-classical expression in refs. [1,7,8]. The question that then naturally arises is: why did we reach a different conclusion from [1] and from the result above in eq. (4.2) when studying rapidly expanding bubbles in thermal equilibrium?
The answer lies in the use of the energy conservation across the wall as in eq. (4.3). By applying it, we impose a non-equilibrium distribution inside the bubble and thus we violate Lorentz symmetry which was crucial for obtaining the (γ 2 − 1) scaling of the friction force in the lte case.
In tracing individual particles and applying energy-momentum conservation the approach of ref. [1] is closer to our classical treatment of the ballistic regime. We derived the distribution inside the bubble, assuming thermal equilibrium outside and applying Liouville equation (see eq. (3.25)). If highly relativistic wall velocities are considered typical momenta in the wall frame are much larger than the masses involved. Then, the distribution inside and outside the bubble are approximately the same. This brings our analysis of the ballistic limit close to the analysis of ref. [1] and, moreover, explains why in that limit the result -that the force saturates for large γ -is similar. Nevertheless, it is important to emphasize that our treatment of the ballistic limit is not limited to the highly relativistic regime. Moreover, our general formula (2.7) allows one to calculate the bubble force not just in local thermal equilibrium or in the ballistic approximation, but also in more general settings that go beyond the semiclassical treatment. The limited, one-loop computation of this work does not include the transition radiation which was treated in refs. [10], however, higher-order computations should account also for that effect.

Standard model and its extensions
In this section we utilize the results from sections 2, 3 and appendix A to study the phase transition dynamics in the standard model and some of its popular extensions. To keep our analysis as simple as possible we do not analyze here the bubble nucleation (from which one 14 In ref. [1] it is mentioned that the force coming from the difference in vacuum potential on two sides of the wall should be taken into account. Then, a static limit is discussed and a conclusion is reached that a wall can be static only at the critical temperature.

JHEP01(2021)070
can infer the latent heat of the transition, bubble surface tension of the transition, etc.), neither we address the non-equilibrium aspects of the transition. Instead, we compute the force on the expanding bubbles, from which one can then infer the dynamics of the transition if one knows the pressure difference across the bubble interface. Below we apply the local thermal equilibrium approximation. To verify whether it holds one should evaluate the condition given by eq. (3.1) using some estimates for thermalization rates, see e.g. refs. [7,8].
Standard model. Even though it is known that the transition in the standard model is a crossover [32], it is useful to analyze it since its particle content is verified by experiments and the first microscopic analysis of the dynamics of the electroweak phase transition was presented in [7,8], where the authors assumed that the transition is first order. It is well known that not all fields of the standard model are relevant for the dynamics of the electroweak transition. Namely, only those fields which significantly contribute to the thermal effective potential and which exert a large force on the bubble are important, and these are the fields whose mass is of the order of the temperature. For the standard model that selects: the top quark (12 relativistic degrees of freedom, m t = 173 GeV), the Higgs boson (1 relativistic degree of freedom, m H = 125 GeV), W ± bosons (6 relativistic degrees of freedom, m Z = 91.2 GeV) and the Z boson (3 relativistic degrees of freedom, m Z = 91.2 GeV). All other particles are much lighter and do not significantly contribute. For example, the next heaviest particle is the bottom quark, whose mass is about 5 GeV. Unless there is a large supercooling such that nucleation occurs at a very low temperature, comparable with or lower than 5 GeV, the bottom quark and other particles of the standard model are not important for the phase transition dynamics.
According to eq. (2.7) it is the change in the entropy density across the bubble, T ∆s = ρ p + P p (see (3.24)) that determines the force per unit area on the bubble. With this in mind and from eq. (A.21) we can extract the top contribution, From the analysis of the Abelian Higgs model in its condensate phase in appendix A we infer that the Higgs boson after symmetry breaking contributes as one massive scalar field (see eqs. (A.125) and (3.23)-(3.24)), Based on the same Abelian Higgs model analysis (see eq. (A.125)) we can conclude that W ± and Z bosons contribute as two and one massive gauge bosons, respectively, To get an impression of how strong the bubble force is in a theory with the matter content of the standard model, in figure 10 we plot the increase in the entropy density across the bubble wall as a function of βm t (when plotting figure 10 we made use of βm i = (m i /m t )βm t with i = H, W, Z). When compared with the single scalar field in figure 3, we see that the entropy increase -and thus also the force on the bubbleis, as expected, much larger in the standard model than in the real scalar case. This is so because there are many more heavy degrees of freedom in the standard model. To be precise, twenty two in the standard model vs one in the real scalar field.
In figures 11 and 12 we show the bubble Lorentz factor and the corresponding expansion speed for the standard model as a function of the top condensate in units of the temperature. When compared with the real scalar field in figures 4 and 5, we see that the bubbles become non-relativistic already for reasonably strong transitions for which βm t 1. Given that the top Yukawa y t 1, this is equivalent to φ 0 /T 1. Recalling the Shaposhnikov's baryon washout criterion for the strength of the transition, φ 0 /T ≥ 1, we infer that the bubbles at a strongly first order phase transition with the standard model matter content are typically subsonic, which broadly speaking agrees with the results of the more detailed microscopic analyses of Moore and Prokopec presented in refs. [7,8].

Standard model with a singlet.
There are several types of extensions of the standard model in which the standard model is extended by a scalar field which is a singlet under the standard model gauge group. Examples include portal models [33] (an important class of which are conformal portal models, see e.g. refs. [34][35][36][37] and references therein) and supersymmetric models, which include scalar singlet fields, a notable example being the NMSSM.
Quite generically both classes of models lead to strongly first order transitions. Even though in portal extensions the transition often proceeds in two stages -the nucleation  along the scalar singlet direction is followed by a rolling along the Higgs direction [37] studying how the additional field content affects the transition dynamics can still be useful. In this case, we ought to add to (5.1)-(5.4) the contribution of the singlet, 15 where m s and N s denote the singlet mass and its number of relativistic degrees of freedom, respectively (for example, N s = 2 if the singlet is a complex scalar). As a general remark, adding more massive degrees of freedom generally increases the entropy production across the expanding bubbles, which in turn increases the force on the bubbles, thus slowing them down. On the other hand, if the character of the phase transition is changed -as it is for example in the aforementioned conformal models - 15 If the phase transition occurs along the scalar direction and then followed by rolling in the Higgs direction, then only the particles acquiring mass in the first stage are relevant for the pressure calculation.

JHEP01(2021)070
strong super-cooling can be present, resulting in a large latent heat release, which in turn accelerates the walls. It is also worth noting that, if there is more than one scalar field in the theory, then multistage transitions are possible, see e.g. ref. [37]. Even though the dynamics of such multistage transitions can be analyzed with eq. (2.7), such transitions are beyond the scope of this paper.

Summary and discussion
We derive a general quantum field theoretic formula (2.7) for the terminal velocity of expanding bubbles of a first order phase transition. The formula has a simple physical interpretation. If local thermal equilibrium is attained across the bubble interface, the friction force acting on a bubble is proportional to the entropy increase across the bubble, with the proportionality constant being the Lorentz factor squared, Our formula is applicable to quantum field theories both in and out-of equilibrium and can be applied at the one or higher-loop level. In this paper we show how to apply (2.7) at the one-loop level in the toy model with one scalar field which exhibits spontaneous symmetry breaking (in section 3), as well as to the standard model and its simple extensions (in section 5). Our analysis applies to bubbles with a moderate Lorentz factor, γ 10, and when scatterings are efficient, such that the local thermal equilibrium approximation applies and the relevant two-point functions are approximately thermal. A particular attention is devoted to how to obtain gauge independent results, the details of which are given in appendix A. Our formula (2.7) generalizes the previous known semi-classical formula of refs. [7,8], and reduces to it in the adiabatic one-loop approximation. However, our analysis reveals an important dependence of the force on the relativistic γ factor of the wall which was not known before.
We apply our formula (2.7) to two simple cases. First we assume local thermal equilibrium. Our results show that the force grows as the Lorentz factor-squared, ∝ (γ 2 − 1), and thus the bubbles quickly reach their terminal velocity. This allows to expect, even though our computation formally does not apply to γ > 10, that runaway is not possible in this scenario. The other limit that we analyzed -the ballistic approximation in which the particles do not thermalize efficiently -displays behavior closer to the runaway scenario featuring friction force that saturates already for moderate γ values. These two limiting cases -the local thermal equilibrium which overestimates the force and the ballistic case which underestimates it -can be used to estimate the true force acting on a bubble, which needs to be between the two extrema. Moreover, the force in more realistic scenarios can be estimated as perturbations around one of our two solutions.
While the analysis presented in this work is at places simplistic, it can capture the leading order contribution to the bubble force in a broad range of situations. Nevertheless, there are situations in which our analysis is not accurate enough. For example, when a transition is not very strong, the force can be significantly altered by higher-loop effects which can induce a moderate, or even large, departure from equilibrium, whose effects can be transported away from the bubble interfaces. This type of corrections, as well as their gauge dependence, should be carefully investigated and one should reassess how the analysis presented here is affected when these effects are included. Furthermore, the latent heat JHEP01(2021)070 released by the transition can induce significant effects on the plasma, which can propagate and dissipate in the form of sound waves and turbulence, both of which can heat up and induce a large scale motion of the plasma, thus affecting the force on the bubbles, which should also be taken into account. These are just some of the unaccounted-for effects, which can influence the dynamics of the transition and which are captured by our formula (2.7). In order words, this paper provides an important step towards an accurate modeling of the dynamics of the electroweak phase transition, which can be of paramount importance for a quantitative understanding of the gravitational wave production and baryogenesis at the electroweak scale.

Acknowledgments
We would like to thank Guy Moore for insightful comments which allowed to significantly improve this article, to Dietrich Bödeker for critically reading the manuscript and to Yann Gouttenoire for a fruitful correspondence. This work was in part supported by the D-

A The energy-momentum tensor
In this appendix we present some details of the calculations of the one-loop energymomentum tensor for a massive fermionic field, massless and massive gauge fields and the Abelian Higgs model all in thermal state in Minkowski space. The real scalar field is considered in the main text in section 3. Most of the material covered in this appendix can be found in the literature [38][39][40][41], however not in a single source. Dirac Fermion. The action for a Dirac fermionic field Ψ(x) that suffices for the one-loop calculation is of the form, is the field dependent mass, y is a Yukawa coupling and ∇ µ is the covariant derivative acting on a spinor field [42,43], Here Γ µ is the spin(or) connection, e µ a (x) is the tetrad field, which lifts tensors into the tangent space. For example, g µν (x)e µ a (x)e ν b (x) = η ab and |g| = |det[e νb ]| ≡ |e|, where η ab is the (flat) tangent space metric at a spacetime point x µ . From (A.1) one easily gets the equation of motion for the fermionic Feynman propagator,

JHEP01(2021)070
where α, β, γ are spinor indices, and andρ in is the fermionic density operator (in Heisenberg picture). The Dirac matrices γ µ (x) build a Clifford algebra, and obey the standard relations on spacetime and on the tangent space, respectively, such that γ µ (x) = e µ a (x)γ a , where the Latin letters a, b, c denote the flat tangent space indices. The right hand side of (A.3) follows from the canonical quantization relation, which follows from the form of the canonical momentum of Ψ, Π Ψ = −ia D−1 Ψ * . The factor a D−1 in (A.7) originates from specifying to a homogeneous cosmological spacetime, in which the metric is diagonal with √ −g = a D and 1/a comes from the tetrad e µ b = a −1 δ µ b that projects γ µ (x) onto the tangent space according to, γ µ = e µ a γ a . The stress-energy tensor is obtained by varying the action in eq. (A.1) with respect to the tetrad field according to T Ψ µν = −|e| −1 e a(µ δS δe ν) a and its expectation value reads, which is covariantly conserved, ∂ µ T Ψ µν = 0. 16 We wish to calculate (A.8) at the one-loop order, which corresponds to the Feynman diagram in figure 13, where the solid oriented line is the free thermal fermionic propagator and the cross indicates the fermionic energymomentum insertion, which is obtained by varying (A.8) with respect to Ψ(y) andΨ(y ).
In adiabatic regime (see section 2), when the effects of the expansion can be neglected, and by making use of (iγ µ ∂ µ − m)(iγ ν ∂ ν + m) = (∂ 2 − m 2 )1, where ∂ 2 = η µν ∂ µ ∂ ν , ∂ µ ∂ ν = ∂ ν ∂ µ , and 1 is the unity matrix in the spinor space, one finds that the free thermal fermionic propagator (A.4) can be expressed in terms of a 'scalar' propagator i∆ F (x; x ), with E p = p 2 + m 2 . The principal difference between i∆ F (x; x ) and the scalar propagator i∆(x; x ) in (3.9) is in the thermal part, where the Bose-Einsten distribution function, n BE = 1/(e βEp − 1) in eq. (3.9) is replaced by the Fermi-Dirac distribution, n FD = 1/(e βEp + 1) in (A.10) and an overall minus sign in the fermionic thermal part (which can be traced back to the anticommuting nature of the fermionic fields). This affects, for example, the coincident fermionic propagator through, where we introduced a fermionic thermal integral, and we should keep in mind that ∂ z J F (n, z) < 0. Notice that we evaluated the thermal contribution to the coincident fermionic propagator (A.11) in D = 4. The result (A.11) is to be compared with (3.12) and (3.11), where the difference is in the form of the fermionic integral (A.12), but also in the sign of the coincident thermal propagator in (A.11). This sign difference can be traced back to the fact that the fermionic propagator is formally a one-loop quantity, and each fermionic loop contributes with a minus sign when compared with a bosonic loop.
We have now all the elements needed to evaluate the one-loop energy-momentum tensor (A.8). For that we firstly need the contribution of the mass term, where Tr[1] = 2 D/2 is the number of fermionic degrees of freedom in D spacetime dimensions, which reduces to four in D = 4 (two chiralities and two helicities). Next, where, as before when we considered the real scalar field in section 3, the T * product indicates that the derivative ∂ x ν commutes with the time ordering operator T * . With

JHEP01(2021)070
eq. (A.9) in mind, it is clear that only the term containing two γ µ 's contributes. By recalling that, , we can rewrite eq. (A.14) as, where the last equality follows from the symmetry of the propagator under the exchange, x ↔ x . Analogous considerations show that the second term in (A.8) contributes equally, such that we have Just as in the case of the real scalar, the vacuum part of (A.17), which equals, can be renormalized by the counterterm action (3.19). If the fermion mass m is generated by a real scalar field condensate, where y is a Yukawa coupling and φ 0 denotes a scalar condensate that may be adiabatically varying in space and time, then in the minimal subtraction scheme the counterterm coupling required to renormalize (A.18) reads (see eqs. (3.20)-(3.22)), Upon adding the counterterm contribution to the energy-momentum tensor one obtains the sought-for renormalized one-loop energy-momentum tensor for the Dirac fermion, This result implies that the renormalized vacuum contribution to the energy density of a Dirac fermion is negative and four times as large as that of a real scalar field. Its thermal contribution harbors four fermionic degrees of freedom which obey Fermi-Dirac statistics.

JHEP01(2021)070
Abelian gauge field model. The action for an Abelian gauge field reads, where F µν = ∂ µ A ν − ∂ ν A µ is the field strength and A µ is an Abelian gauge field, a prime example of which is the photon of the quantum electrodynamics (QED). The action (A.22) possesses an Abelian gauge symmetry, under which the field transforms as, where Λ(x) is an arbitrary function of space and time. This means that A µ contains redundant (unphysical) degrees of freedom which make the kinetic operator for A µ derived from (A.22) non-invertible. Without going into the details of the gauge fixing procedure, this is resolved by adding a gauge fixing term to the action, which makes the kinetic operator invertible, but does not change any physical quantity. A legitimate gauge fixing is the Fermi gauge, which is convenient since it has one gauge parameter ξ which can be used to control gauge dependence of the results. The corresponding gauge fixing action and Lagrangian are, The canonical momentum of the theory is obtained by a variation of the total action, S tot = S EM + S Fermi with respect to ∂ 0 A µ (x), and in Minkowski background it reads, such that e.g. Π 0 A = ξ −1 η ρσ ∂ ρ A σ does not vanish. Due to the added gauge fixing term (A.24), all of the canonical momenta become dynamical, such that the canonical commutation relation in the Fermi gauge is simple, The Feynman propagator equation in the Fermi gauge is therefore, The solution of (A.27) can be written as, where i∆ 0 (x; x ) is the massless limit of the thermal scalar propagator (3.9), which equals (see the n = 0 term of the second series in eq. (3.13)), where ∆t = t − t , r = x − x , ∆x 2 (x; x ) is given in eq. (3.10) and, for simplicity, we have evaluated the thermal part (A.31) in D = 4 (because it is finite in D = 4 and therefore does not need to be regularized). The first term in the second line (A.31) is not there to cancel the principal part (P) of the vacuum propagator (A.30) in D = 4, but instead it is needed to get the thermal contribution to vanish in the zero temperature limit, β → ∞.
We are now ready to proceed with the calculation of the one-loop energy-momentum tensor, which for the total photon action, consisting of eq. (A.22) plus the gauge fixing part (A.24), reads, where we also took an expectation value of the operator-valued energy-momentum tensor. Next, the derivatives in (A.32) can be taken out of the expectation values, provided one replaces the T with the T * time ordering, such that the energy-momentum tensor (A.32) can be recast as, This expression can be simplified by making use of the tensor structure of the photon propagator (A.29), It is now clear that in dimensional regularization the vacuum part of i∆ 0 (x; x ) in (A.30) does not contribute to (A.35). However, the thermal part does contribute and the result is, The result (A.38) cannot be correct, since it suggests that the photon has four degrees of freedom, instead of two of the physical photon (the two transverse polarizations with helicities, h = ± ). The reason is the gauge fixing term (A.24) which made all four polarizations of the photon dynamical (cf. eq. (A.25)), thus explaining (A.38). Even though eq. (A.38) does not depend on the gauge parameter ξ, it is not correct because we did not take account of the contribution from the Faddeev-Popov ghosts. It is well known that in the vacuum the contribution from ghosts in Abelian gauge theories vanishes, for more general states such as thermal states however, ghosts do contribute and thus have to be taken into account. To show that, consider the Fadeed-Popov ghost action and Lagrangian associated with an Abelian gauge field, from where one obtains the ghost propagator equation, Had we defined the ghost propagator with the ghost field ordering as indicated after the second equality in (A.44), we would have obtained −iδ D (x − x ) on the right hand side of (A.43). Ghosts are complex, bosonic, anticommuting fields, so their propagator 17 The small argument expansion of (A.31) is, i∆ th 0 (x; x ) ∼ 1 12β 2 − π 2 180β 4 (r 2 + 3∆t 2 ).
The argument presented in footnote 19 shows that, as a consequence of treating the ghosts as Grassmannian fields in thermal equilibrium, implies that they should obey a Fermi-Dirac statistic.
The reason why ghosts are treated as anticommuting fields is that, upon integrating them out, one reproduces the correct Faddeev-Popov O FP (x; x ) determinant. This determinant is obtained by a functional variation of the gauge fixing condition with respect to the gauge parameter, which ensures the correct field-dependent integration measure along gauge orbits, and which for the problem at hand reads, can be represented as a path integral over the anticommuting Faddeev-Popov complex ghost fields, On the other hand, eq. (A.52) implies, such that its contribution to the effective action is of the form, 19 One can consider the positive frequency thermal Wightman function for the ghost field, Tr e −βĤ gh , (A. 46) whereρ th is the thermal density matrix andĤ gh is the ghost Hamiltonian. By taking the cyclic properly of the trace, c(t, x) = e −iĤ gh c(0, x)e iĤ gh and the non-commuting nature of the ghost fields, one obtains the following Kubo-Martin-Schwinger condition for the ghost, which reads in momentum space, i∆ + (p α ) = −e −βp 0 i∆ − (p α ). This suggests that the thermal part of the ghost propagator should obey a Fermi-Dirac statistic.

JHEP01(2021)070
which equals minus twice that of a real scalar field, which contribute as, Γ φ = (i/2)Tr [ln (O φ (x; x ))]. This observation has led to the development of the ghost field formalism, according to which ghosts are -just like fermions -anticommuting fields, but they obey a Bose-Einstein statistic. However, as we argue in footnote 19, that is inconsistent with the notion that a thermal state is defined in terms of the thermal density operator. Since ghosts fields are commonly viewed as 'unphysical', the field practitioners have gotten used to this inconvenience, and declared by fiat that ghosts fields obey a Bose-Einstein statistic. 20 Our take on this is that, once degrees of freedom are added to the Hamiltonian, they are degrees of freedom and have to be dealt with as such.
Here we take a different take on ghosts. Note that the Faddeev-Popov determinant can be also written as, where φ gh is now defined as a complex (commuting) scalar ghost field andφ gh = φ * gh . The price to pay is the non-local action (A.56) which governs the dynamics of the ghosts and which, upon a partial integration, can be recast as, The principal disadvantage of this action is that it is non-local. One should not be scared by non-local actions however, as they have successfully been dealt with in the literature on dark energy [45][46][47]. One may wonder whether the nonlocality can be dealt with by introducing auxiliary fields, χ = −1 φ gh , in terms of which the ghost action can be rewritten into an on-shell (weakly) equivalent (S gh ≈ S gh ) form as, where λ = λ(x) andλ =λ(x) are Lagrange multiplier fields. Upon solving for the constraint fields λ andλ, one gets the original action. However, varying with respect to χ andχ yields,

JHEP01(2021)070
Now λ andλ can be eliminated in favor of χ andχ, to obtain yet another on-shell equivalent ghost action, This action is local, and that is desirable. However the price we paid to get (A.60) is in that χ andχ are ghost scalars, i.e. they have a negative kinetic term. As long as we do not include dynamical gravity, this need not be fatal for the theory, and we can work safely with it. 21 The action (A.60) is not problem free however. Indeed, upon varying (A.60) with respect toχ, χ,φ gh and φ gh we get, Obviously, these equations cannot be the correct equations for the ghost sector. The correct equations are obtained by acting 1/ on the first and third equation in (A.61), which are equivalent to the equations obtained from the original action (A.56)-(A.57). The above exercise is instructive, as it teaches us that one ought to be extra careful when making use of the procedure which is known to be valid for constrained systems and whose dynamics is described by local Hamiltonians. In order to avoid such pitfalls, we shall work here with the non-local version of the theory (A.57), which is feasible at the one-loop level and that is what we do next. Varying the action (A.57) with respect toφ gh andφ gh gives, and for simplicity we assumed in (A.63) a Minkowski background. Eq. (A.65) implies that the ghosts φ gh andφ gh are bosonic fields which obey a Bose-Einstein statistic, and the solution with thermal boundary conditions can be written as (see footnote 25 below), where i∆ th M (x; x ) denotes the thermal part of the massive scalar propagator (3.9). In eq. (A.66) we shifted the poles of the massless propagator, (k 0 ) 2 = k 2 by δM 2 = M 2 , which regulates the would-be singular behavior of its thermal part. In what follows, it will become clear how to use (A.66) in practical calculations.

JHEP01(2021)070
The next step is the one-loop energy-momentum tensor from the ghost fields, T gh µν = −(2/ √ −g)(δS ghost /δg µν ), whose expectation value is, which, upon extracting the derivatives, can be recast as, x →x (A.68) where we set, → ∂ 2 and → ∂ 2 . We can now use the ghost propagator (A.66) to show that the vacuum part of the energy momentum tensor (A.67) vanishes in dimensional regularization. The thermal part, on the other hand, gives a non-trivial contribution. Indeed, by taking account of, one sees that only the first term in eq. (A.67) contributes (on-shell), which equals minus twice the contribution of a massless scalar field (cf. eqs. (A.39) and (A.50)), resulting in, Upon adding this ghost contribution to (A.38) one finally obtains the gauge independent and physically correct one-loop energy-momentum tensor for the photon in thermal equilibrium, The Feynman diagrams contributing to the one-loop energy-momentum tensor of a massless gauge field (the photon) are shown in figure 14, where both the photon loop (wiggly) and the ghost loop (dashed) are shown. From the result (A.71), one can read off the energy density and pressure of an ideal gas of photons, ρ = π 2 T 4 /15, P = π 2 T 4 /45, such that P = (1/3)ρ, which is the correct result for an ultrarelativistic plasma with two degrees of freedom. The corresponding entropy density is, s = k B (P + ρ)/T = k B 4π 2 T 3 /45. The important lesson to learn from this calculation is that, even though we obtained a result in which any dependence on the gauge parameter ξ dropped out even without including the ghosts, the result was in a subtle way incorrect, and only when we included the ghosts' contribution we obtained the physically correct result. Massive gauge field. We include in this appendix a massive Proca gauge field A µ , since it appears in the calculation of the energy-momentum tensor of the Abelian Higgs model and consequently of the standard model. The Proca action is given by,

JHEP01(2021)070
Upon acting ∂ µ on (A.73) we see that the Lorentz condition, η µν ∂ µÂν = 0 is automatically imposed by the equation of motion, which kills one degree of freedom. Because of this, it is natural to require that the Feynman propagator obeys the equation of motion, 22 such that the Proca propagator is transverse on both legs, In the language of gauge theories this can be formally viewed as the Landau gauge of the massive gauge field propagator of the Abelian Higgs model, whose detailed analysis can be found below in this appendix. The transversality property (A.75) dictates the form of the Proca propagator, where i∆ M (x; x ) is the real scalar propagator (3.9) whose mass is m = M . 22 The precise meaning of the operator, G(x; x ) ≡ (1/∂ 2 )δ D (x − x ), is revealed upon acting with ∂ 2 on both sides of the equation, which gives ∂ 2 G(x; x ) = δ D (x − x ), which means that iG(x; x ) = i∆0(x; x ) is the massless scalar propagator, whose precise form is subject to boundary conditions. Since here we are interested in the thermal propagator on Minkowski space, iG(x; x ) is the thermal propagator for a massless scalar and its spacetime dependence is given in eqs. (A.49)-(A.50).

JHEP01(2021)070
When (A.82) is added to (A.79) one obtains the renormalized vacuum contribution to the energy-momentum tensor. Taking account of the thermal contribution as well (which is three times that of the real scalar field given in (3.23)), one obtains the renormalized one-loop energy-momentum tensor of a massive vector field in the thermal state, Abelian Higgs model. This model, also known as scalar quantum electrodynamics (SQED), consists of one complex scalar Φ and one Abelian gauge field A µ and its action is given by, where the gauge-covariant derivative is defined as, where g is the gauge coupling constant, µ 2 is the scalar mass parameter and λ Φ its quartic self-coupling. The model possesses an Abelian gauge symmetry. This means that the Lagrangian (A.85) is invariant under the following local field transformations, When the mass parameter µ 2 > 0 in (A.86), then µ is the tree level scalar mass. In this symmetric phase, the model (A.85)-(A.86) contains one massive complex scalar with mass m = µ and one massless gauge field. The corresponding contributions to the one-loop energy-momentum tensor have already been calculated: the energy-momentum tensor of a complex scalar field is two times that of a real scalar given in eq. (3.23) with mass m → µ and the energy-momentum tensor of a massless gauge field is given in eq. (A.71).
When, on the other hand, µ 2 < 0, the scalar field Φ exhibits a spontaneous symmetry breaking and acquires a condensate by the famous Brout-Englert-Higgs (BEH) mechanism. To study the mechanism, it is convenient to decompose Φ as, where we take the field condensate φ 0 = Φ to be real (this can be always achieved by a suitable global gauge transformation (A.88)). The gauge field A µ is assumed not to develop JHEP01(2021)070 a condensate. Then ϕ 1 is the Higgs field and ϕ 2 is the Goldstone boson of the symmetry 'broken' by the condensation. 24 Upon inserting (A.90) into (A.86) we get (up to boundary terms) for the quadratic part of the Lagrangian, where we have defined M 2 = (gφ 0 ) 2 and m 2 H = 2λφ 2 0 and for simplicity calculated L (2) on Minkowski background, on which g µν → η µν . Note that the scalar and gauge perturbations couple in the gauge (A.90) via the term M ∂ µ A µ ϕ 2 . The fields decouple in the 't Hooft gauge [48], in which the total quadratic gauge fixed Lagrangian reads (with g µν → η µν ), where we also included the Faddeev-Popov ghost Lagrangian, which can be obtained from the Faddeev-Popov determinant, which is obtained by considering how the gauge fixing condition varies with respect to infinitesimal gauge transformations, Λ(x) → θ(x), where M 2 c = ξM 2 and we made use of (A.92), δ θ A µ = ∂ µ θ and δ θ ϕ 2 = −gθ(φ 0 + ϕ 1 ). From (A.94) it follows that the corresponding Lagrangian can be written in terms of Grassmannian ghost fields as, where L (2) ghosts is the quadratic part of (A.95). That means that the model consists of two massive scalars with masses m 2 1 = m 2 H = 2λφ 2 0 and m 2 2 = ξM 2 = ξ(gφ 0 ) 2 , one massive gauge field with M 2 = g 2 φ 2 0 and a massive ghost with M 2 c = ξM 2 , where we assume ξ ≥ 0 for stability. All of these particles contribute to the energy-momentum tensor, and the corresponding Feynman diagrams are shown in figure 15. 24 In fact, since the Ward identities generated by the symmetry transformation (A.88) are respected both in the symmetric and in the condensate phase, the gauge symmetry is never really broken. Nevertheless, because the scalar condensate generates a mass for the gauge field, and a massive Proca theory does not possess a gauge symmetry, the condensate phase is often -but not correctly -referred to as the broken phase. In section 3 we analyzed a massive scalar field, and in this appendix a massless gauge field, a massless ghost and a massive gauge field, but the massive gauge field was a Proca field which does not possess any gauge symmetry.

JHEP01(2021)070
This means that we still ought to analyze the massive gauge field and the massive ghost. Let us start with the massive gauge field. From (A.93) we can read off the propagator equation of motion, We leave the analysis of this case as an exercise to the reader.