Field theoretic derivation of bubble wall force

We derive a general quantum field theoretic formula for the force acting on the 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, the force grows as the Lorentz factor squared, such that the bubbles cannot run away. We apply our formalism to a massive real scalar field, the standard model and its simple portal extension. Next, we compare with the existing literature. In particular we discuss the differences and similarities of our work with that of B\"odeker and Moore [arXiv:0903.4099]. We find that a bubble can run away if scatterings are negligible across the wall (ballistic limit), but it always attains a finite Lorentz factor if scatterings are efficient across the bubble. 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-of-equilibrium) 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 at 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 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]. 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 energy-momentum tensor of the bubble-plasma system. When applied to the limit when an approximate Lorentz symmetry 1 holds, we find that the friction force scales as (γ 2 − 1), excluding the possibility of a runaway. 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 the bubbles can 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 one-loop 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 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 (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 (1) simplifies to, Several remarks are in order: • The energy-momentum tensor in Eqs. (2-3) is a composite operator which diverges, so to make it finite it ought to be regularized and renormalized; • The expectation values in (3) ought to be calculated in a state that includes the bubble, which approximates it well enough; • Eqs. (2)(3) simplify further in the bubble frame, in which the terms containing time derivatives drop out. 2 In the bubble frame, one can then integrate equations (2)(3) across the bubble to obtain, where ∆ T p,b µν denotes the change of the µν components of the energymomentum tensor of the plasma or bubble across the bubble. Since Eq. (4) contains one time and one spatial derivative acting on a field, it must vanish. Namely, its vacuum part is Lorentz covariant and hence proportional to ∂ i ∂ 0 ∆x 2 (x; x )| x →x = 0, where ∆x 2 (x; x ) = −(t − t ) 2 + x − x 2 , while the thermal contribution is proportional to the derivative ∂ i ∂ 0 acting on an even function of ∆t = t−t (cf. Eq. (21)), which vanishes when evaluated at coincidence, t → t. Therefore, Eq. (5) contains all of the relevant information about the phase transition dynamics. The bubble wall speed is determined by the balance of the two terms in (5).
When calculating the plasma contribution in (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 , T b zz = P b , 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 force F per volume V on an expanding bubble, where we used, P = P p +P b and ∆P b = ∆ T b zz . The relation (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. By making use of (7) one can determine the bubble speed. This simple observation constitutes the principal result of this work.
Before we proceed to applications, in what follows we discuss applicability of formula (7).

1.
To arrive at (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 on-shell across the bubble interface. To estimate when the quasiparticle approximation is reliable, recall that according to the uncertainly 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) , where we assumed that E ∼ T . 4 . 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 (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. Figure 1: The propagation of a particle (horizontal solid black line) in presence of field dependent mass insertions, m(x) = λ/2φ 0 (x) (vertical dashed lines). 4 The condition (8) is not the most general one. Namely, in view of Einstein's relation, E = p 2 + m 2 , the on-shell condition is most restrictive for highly infrared particles, for which p 2 < m 2 , such that E ≈ m and the on-shell condition (8) becomes γv < γ < Lm. Since typically the most massive particles dominate the bubble friction, it is often the case that m ∼ T and (8) applies. If, however, the heaviest particle's mass is much smaller than T , then (8) ought to be replaced by the more stringent condition, γv < γ < Lm When bubbles are very fast (or very thin) and the criterion (8) is violated, quantum off-shell effects can become significant. To capture that, when constructing the energy-momentum 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. [15]. 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 [16,17,18]. While progress has been made in applying such a formalism in baryogenesis/leptogenesis scenarios [19,20,21,22] (for reviews see [4,23]) and in cosmic inflation [24,25,26,27,28], little or no progress has been made in studying cosmological phase transitions.
2. While our conclusions based on the consideration of the energymomentum tensor operator are general, the form of Eq. (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 if 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, where Here η and ζ denote the shear and bulk viscosity, respectively, and (αβ) means that the indices are symmetrized. Notice that τ µν is orthogonal to u µ , u µ τ µν = 0 = τ µν u ν , and that its trace is proportional to the bulk viscosity, Tr[τ µν ] = g µν τ µν = 3ζ∇ · u, where we used, u 2 = −1 and u µ (u · ∇)u µ = 0. If the plasma is in thermal equilibrium, then (10) vanishes. 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 (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. (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 (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 (11) by dV to obtain its local form, 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 5 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. (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 oneloop approximation and in local thermal equilibrium, ∆P = −∆V eff , where V eff is the one-loop effective potential. 6 Therefore, we have ∆s = − ∆V eff . It would be worth investigating whether the two approaches are equivalent in more general situations.
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 (7) for the phase transition dynamics, in what follows we first consider a simple real scalar field model and then more general models.

Real scalar field
In this section we calculate the one-loop energy momentum tensor of a selfinteracting, massless scalar field in the presence of an expanding spherical 6 To see that, note that from Eq. (30) we can read off the thermal contribution to the pressure, P = −(6π 2 β 4 ) −1 [z −1 ∂ z J B (6, z)] z=βm , where J B (6, z) is given in Eq. (23). By a suitable partial integration one can show, z −1 ∂ z J B (6, z) = −3J B (4, z), from which we conclude, P = (2π 2 β 4 ) −1 J B (4, βm) = −V eff , which completes the proof.
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 (13) is more stringent than the on-shell condition (8). If (13) 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 and then discuss how our results are affected if ballistic approximation represents a more appropriate description. The free scalar field action is given by, 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 (14) will generate in the presence of a condensate φ = φ 0 a field dependent mass, 7 If the background field φ 0 = φ 0 (x) varies in spacetime, the mass (16) will follow the suit. As long as the variation is slow, it can be treated as an adiabatically varying quantity.
Varying (14) with respect to g µν results in the energy-momentum tensor, Here we are interested in evaluating the expectation value of the one-loop energy-momentum tensor, The contributing one-loop diagram is shown in figure 2, where the dashed line denotes the scalar propagator and the cross (×) indicates the energymomentum tensor insertion.

Local thermal equilibrium approximation
Since at the one-loop order the energy-momentum tensor (18) 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 (13) 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 (19) obeys the equation of motion, where m = m(φ 0 ) is the field dependent scalar field mass. The thermal propagator in (20) can be thought of as the inverse of the operator, √ −g ( − m 2 )× δ D (x−y), where = g µν ∇ µ ∇ ν denotes the d'Alembertian operator with thermal boundary conditions imposed. In practice this can be done by inserting the thermal density operatorρ th in (18). 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, 8 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 (21) 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 the propagator in (21), 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 (21), such that the vacuum part is divergent in D = 4 and ought to be regularized. The vacuum contribution in (24) was evaluated by noting that the Bessel function K ν (z) in Eq. (21) 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 (25) 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 (18) we need, where we introduced the usual T * time ordering which is defined to commute with the two external derivatives in (26). The vacuum part of (26) is −2η µν times the linear coefficient in ∆x 2 of the propagator in (21), which is easily extracted from the integer series in (25), When taken together with the vacuum part in (24) this then implies, T [L] vac = 0, such that, which has the form of a cosmological constant. The thermal contribution to (26) is given by the integral over p µ p ν of the coincident thermal integral in (21). 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 (30) are the thermal contributions to the one-loop stress-energy tensor.
To complete the calculation, we still ought to renormalize the vacuum contribution (28), 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 (16), 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 (28) around D = 4, where µ is an arbitrary scale and m 2 = λφ 2 0 /2. Comparing (33) with (32) we see that, removes the divergence from (33), resulting in the following renormalized energy-momentum tensor, From the point of view of the bubble force calculation, the change in the entropy density (6) across the bubble is what determines the force on the bubble interface, which is determined by the second line in (35) (Lorentz covariant contributions have a vanishing entropy density), Inserting this into Eq. (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. of a real scalar field thermal plasma across the nucleated bubble as a function of the scalar mass m/T . Since the scalar mass m = λ/2φ 0 (x) increases across the bubble as the scalar condensate increases, the entropy in the broken phase decreases. The maximum amount by which the entropy density can change is 2π 2 T 3 /45, which is formally reached when m → ∞, when the entropy density inside the bubble tends to zero (horizontal dashed).  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. Therefore, we conclude that, if local thermal equilibrium is maintained across the bubble, it cannot run away. 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.

Ballistic approximation
If the thermalization rate is not large enough to satisfy Eq. (13) 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 (37) can be written as a general function of m 2 + p 2 z and p ⊥ (or, equivalently, of E), 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 (40) will 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 (40) will be absent. When the boundary condition (39-40) is imposed on (38) one obtains the general ballistic solution, which is conveniently broken into three parts as follows.
Case A. Transmission (t + ) from the symmetric phase: Case C. Transmission (t − ) from inside the bubble: where we temporarily introduced an index 0 on m to denote the mass deep inside the bubble. Note that in the reflected solution (42) 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 (41)(42)(43) 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 γ in the plasma frame quantities in Eq. (7), making thus the γ dependence of the bubble force more complex than what Eq. (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 that is the pressure in the direction in which the bubble expands. For simplicity in Eq. (44) we took account of the plasma contribution only (since the vacuum contribution is unchanged). The form of the pressure in (44) can be obtained from the zz part of Eq. (26), provided the thermal distribution function 1/(e βEp − 1) in (21) gets replaced by the ballistic distribution function (41)(42)(43). This ballistic approximation is valid as long as the semiclassical on-shell condition (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 (44) 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. (7). To evaluate the integrals in (44), it is convenient to first integrate over E, in which case the ranges of integration are, E ∈ [ p 2 z + m 2 , ∞), p z ∈ (−∞, ∞). By making use of EdE = p ⊥ dp ⊥ and integrating e.g. the transmitted contribution (41) over E one gets, (45) 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 (45) becomes, where v ≡ v(γ) = 1 − (1/γ 2 ). An analogous procedure gives for the re-flected contribution (42) and for the transmitted contribution (43), The total pressure is then simply the sum of the three contributions, and the pressure difference across the bubble, There is a subtlety involved in evaluating (49) in the reflected contribution (47), 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 (46), (47) and (48) one obtains for the pressure difference, 9 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 (51) 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 (51) 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. (51) is to be compared with the pressure difference across the interface obtained assuming local thermal equilibrium (lte) close to the interface, for which f = 1/[e βγ(E−vpz) − 1] and hence, ∆P z,lte (βm, γ) = π 2 90β 4 4(γ 2 −1)+1 From Eqs. (41)(42)(43) we see that, in the ballistic approximation, all particles ascending onto the interface (Cases A and B) slow down, such that f > f lte , 9 To arrive at Eq. (51) we made use of the integral, 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. 10 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 (51) 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 γ, implying that in the ballistic limit the bubbles can run away. The maximum force reached when γ → ∞ can be quite straightforwardly evaluated from Eqs. (51-52) as follows. In the limit when γ 1, and the pressure difference ∆P z in Eq. (51) can be estimated as, In this limit our estimate (55) agrees with the result of Ref. [1], Eq. (2.4).
To get the result (55) we have used the small argument expansion of the J B integral in the first line of Eq. (51). The terms in the second line can be neglected, since in that limit the polylogarithm functions are exponentially suppressed as ∼ e −2γβm , which is easily seen from their small argument expansion. Finally, the third line in (51) is exponentially suppressed at least as ∼ e −γβm , 11 and hence does not contribute at the leading order to (55).
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, 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. , and m = 0.8T (solid blue). As expected, the force is smaller for smaller masses.
upon which the integral naturally splits into two parts, I t− = I 1 + I 2 , where 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 I t− is suppressed at least as ∼ e −γβm .
To illustrate how the bubble force depends on the particle mass, in figure 7 we plot ∆P z in (51) 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 by how much the bubble force calculated in the ballistic approximation is smaller from that in the local thermal equilibrium (lte), in figure 8 we plot ∆P z (51) 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. Finally, in figure 9 we plot the natural logarithm of ∆P z (51) as a function of the Lorentz factor γ for a fixed mass m. Left panel shows, from 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 producing 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 can run away. In the opposite limit, when particle scatter efficiently and can be consequently approximated by a local thermal equilibrium, the entropy production and the bubble force maximize. As long as a local thermal equilibrium applies, entropy production grows as the Lorentz factor squared and the bubbles cannot run away. Instead they reach a terminal velocity which can be accurately estimated from the one-loop the energy-momentum tensor. We emphasize that formula (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, both the ballistic and 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 longstanding 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 energymomentum conservation law (3). When written in the bubble wall frame (5) it can be recast as (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. (58) indicates that the force on the bubble is present only when the bubble is moving.
The precise form of the energy-momentum tensor in (58) 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 one can use the thermal part of the massive scalar propagator (21) to calculate the force in (58). When integrated across the bubble, Eq. (58) can then be recast as, 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. (59) we see two unimportant differences. Firstly, the force in Eq. (59) 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. 12 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 on both sides of the wall is imposed as was done above, Eq. (3.3) from Ref. [1] and our result in Eq. (59) do agree up to the minor differences already discussed.
The above analysis shows that, starting with our main result (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] when studying rapidly expanding bubbles?
First thing to observe is that in Ref. [1] Eq. (3.3) was not the one used to obtain the final result, it was instead, Eq. (3.5), which follows from Ref. [9]. Motivated by the large γ limit considered, it was assumed that the distribution does not change while crossing the wall. This enables treating the momenta belonging to different sides of the wall under a single integral, enabling the use of the energy-momentum conservation, which in turn breaks the Lorentz structure of the expressions. An approximate Lorentz invariance was crucial for deriving our Eq. (7), and thus this can be viewed as the source of the difference.
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. (37)). 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 (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 goe beyond the semiclassical treatment.

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 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.
Standard model. Even though it is known that the transition in the standard model is a crossover [30], 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. (7) it is the change in the entropy density across the bubble, T ∆s = ρ p + P p (see (36)) that determines the force per unit area on the bubble. With this in mind and from Eq. (87) 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. (191) and (35-36)), (62) Based on the same Abelian Higgs model analysis (see Eq. (191)) we can conclude that W ± and Z bosons contribute as two and one massive gauge bosons, respectively, Note that in Eqs. (62-64) we have assumed that a local thermal equilibrium is attained. The total change in the entropy density in a model with the field content of the standard model is then obtained by simply summing the four contributions in Eqs. (61-64), T ∆s SM = T ∆s t + ∆s H + ∆s W + ∆s Z .
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 bubble -is, 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].  include portal models [31] (an important class of which are conformal portal models, see e.g. Refs. [32,33,34,35] 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 [35]) -studying how the additional field content affects the transition dynamics can still be useful. In this case, we ought to add to (61-64) the contribution of the singlet, 13 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 -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. [35]. Even though the dynamics of such multistage transitions can be analyzed with Eq. (7), such transitions are beyond the scope of this paper.

Summary and discussion
We derive a general quantum field theoretic formula (7) for the force acting on expanding bubbles of a first order phase transition. The formula has simple physical interpretation. If local thermal equilibrium is attained across the bubble interface, the bubble force is proportional to the entropy increase across the bubble, with the proportionality constant being the Lorentz factor squared, γ 2 − 1 = v 2 /(1 − v 2 ). 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 (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 (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 (7) to two simple cases. First we assume thermal equilibrium on both sides of the bubble wall. Our results show that the force grows as the Lorentz factor-squared, ∝ (γ 2 − 1), which qualitatively differs from the well-known result of Ref. [1], which found that, for sufficiently strong transitions and in a large γ limit, the force saturates to a finite value, such that bubbles can exhibit a runaway behavior. That claim was relaxed in a subsequent paper [10], where a higher-loop mechanism was found which prevents the runaway, but still permits a large Lorentz factor. Our analysis of the ballistic limit, in which the plasma particles are in thermal equilibrium outside the bubble and interact semiclassically but so fast that they do not thermalize, reveals a similar behavior -the force saturates for large γ factors.
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. These 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 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 (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.

Appendix: The energy-momentum tensor
In this appendix we present some details of the calculations of the one-loop energy-momentum 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 [36,37,38,39], 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 [40,41], 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 (67) one easily gets the equation of motion for the fermionic Feynman propagator, 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 (69) 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 (73) 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. (67) 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. 14 We wish to calculate (74) 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 energy-momentum insertion, which is obtained by varying (74) 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γ and 1 is the unity matrix in the spinor space, one finds that the free thermal fermionic propagator (70) can be expressed in terms of a 'scalar' propagator i∆ F (x; x ), where with E p = p 2 + m 2 . The principal difference between i∆ F (x; x ) and the scalar propagator i∆(x; x ) in (21) is in the thermal part, where the Bose-Einsten distribution function, n BE = 1/(e βEp − 1) in Eq. (21) is replaced by the Fermi-Dirac distribution, n FD = 1/(e βEp + 1) in (76) 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 (77) in D = 4. The result (77) is to be compared with (24) and (23), where the difference is in the form of the fermionic integral (78), but also in the sign of the coincident thermal propagator in (77). 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 energymomentum tensor (74). 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 Eq. (75) in mind, it is clear that only the term containing two γ µ 's contributes. By recalling that, Tr , we can rewrite Eq. (80) 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 (74) contributes equally, such that we have where Tr[1] = 2 D/2 . When Eqs. (82) and (79) are compared with the analogous results for the real scalar field of section 3, we see that -up to the factor Tr[1] = 2 D/2 which counts the number of degrees of freedom of a Dirac fermion -the fermionic one loop energy-momentum tensor can be obtained from the real scalar one by replacing i∆(x; x ) by i∆ F (x; x ). This then immediately implies that, the expectation value of the Lagrangian vanishes, i.e. T [L Ψ ] = 0 and that (cf. Eqs. (28) and (33)) Just as in the case of the real scalar, the vacuum part of (83), which equals, can be renormalized by the counterterm action (31). 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 (84) reads (see Eqs. (32)(33)(34)), 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.
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 (88) 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 (88) 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 (90), 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, where The solution of (93) can be written as, where i∆ 0 (x; x ) is the massless limit of the thermal scalar propagator (21), which equals (see the n = 0 term of the second series in Eq. (25)), x ) is given in Eq. (22) and, for simplicity, we have evaluated the thermal part (97) 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 (97) is not there to cancel the principal part (P) of the vacuum propagator (96) 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 energymomentum tensor, which for the total photon action, consisting of Eq. (88) plus the gauge fixing part (90), reads, where we also took an expectation value of the operator-valued energymomentum tensor. Next, the derivatives in (98) can be taken out of the expectation values, provided one replaces the T with the T * time ordering, such that the energy-momentum tensor (98) can be recast as, This expression can be simplified by making use of the tensor structure of the photon propagator (95), where we took account of the antisymmetry property of the derivatives, ∂ α i∆ 0 (x; x ) = −∂ α i∆ 0 (x; x ). The first line (101) originates from the gauge field action (99), while the second line (102) comes from the gauge fixing term contribution (100). Notice that both contributions in (101-102) do not dependent on the gauge parameter ξ, and combine into, It is now clear that in dimensional regularization the vacuum part of i∆ 0 (x; x ) in (96) does not contribute to (101). However, the thermal part does contribute and the result is, where 15 The result (104) 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 (90) which made all four polarizations of the photon dynamical (cf. Eq. (91)), thus explaining (104). Even though Eq. (104) 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, where c = c(x) is a complex Grasmannian (anticommuting) scalar field and c = c * (x). The corresponding canonical momenta are, (the minus in Π c = δS ghost /δ∂ 0 c comes from the anticommuting of δ/δ∂ 0 c withc) which imply the following nonvanishing canonical quantization rules, 16 from where one obtains the ghost propagator equation, with Had we defined the ghost propagator with the ghost field ordering as indicated after the second equality in (110), we would have obtained −iδ D (x−x ) on the right hand side of (109). Ghosts are complex, bosonic, anticommuting fields, so their propagator equation is simply related to that of a massless The argument presented in footnote 17 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, whose determinant, 17 One can consider the positive frequency thermal Wightman function for the ghost field, Tr e −βĤ gh , 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, from where we infer, 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.
can be represented as a path integral over the anticommuting Faddeev-Popov complex ghost fields, On the other hand, Eq. (118) implies, such that its contribution to the effective action is of the form, which equals minus twice that of a real scalar field, which contribute as, . 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 17, 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. 18 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 (122) which governs the dynamics of the ghosts and which, upon a partial integration, can be recast as, (123) 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 [43,44,45]. 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, 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 (126) 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. 19 The action (126) is not problem free however. Indeed, upon varying (126) 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 (127), which are equivalent to the equations obtained from the original action (122-123). 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 (123), which is feasible at the one-loop level and that is what we do next. Varying the action (123) with respect toφ gh andφ gh gives, where and for simplicity we assumed in (129) a Minkowski background. Eq. (131) 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 23 below), where i∆ th M (x; x ) denotes the thermal part of the massive scalar propagator (21). In Eq. (132) 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 (132) in practical calculations.
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, where we set, → ∂ 2 and → ∂ 2 . We can now use the ghost propagator (132) to show that the vacuum part of the energy momentum tensor (133) 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. (133) contributes (on-shell), which equals minus twice the contribution of a massless scalar field (cf. Eqs. (105) and (116)), resulting in, Upon adding this ghost contribution to (104) one finally obtains the gauge independent and physically correct one-loop energy-momentum tensor for Figure 14: The one-loop Feynman diagrams contributing to the one-loop energy-momentum tensor of a massless gauge field. The photon contribution is the diagram with wiggly lines and the ghost diagram is dashed and oriented.
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 (137), 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 the the correct results 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, where F µν = ∂ µ A ν − ∂ ν A µ is the field strength and M is its mass. The Proca field has D − 1 dynamical degrees of freedom and no gauge symmetry, such that the calculation of the energy-momentum tensor is rather simple since it is not complicated by gauge symmetry. From (138) one easily obtains the equation of motion (in Minkowski space) for the Proca field, Upon acting ∂ µ on (139) 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, 20 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 (141) dictates the form of the Proca propagator, where i∆ M (x; x ) is the real scalar propagator (21) whose mass is m = M .
We are now ready to calculate the energy-momentum tensor of the Proca field. Varying the action (138) with respect to the inverse metric gives, 20 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. (115-116).
Upon pulling the derivatives out of (143) and making use of the transverse form of the Proca propagator (142) one obtains (cf. Eq. (101)), The vacuum contribution can be extracted from Eq. (21) and the integer series of (25), The thermal contribution is a homogeneous solution that satisfies the on-shell condition, (∂ 2 − M 2 )i∆ th M (x; x ) = 0, such that can be evaluated where we made use of, (1/∂ 2 )i∆ th M (x; x ) = (1/M 2 )i∆ th M (x; x ). Comparing this with (26) we conclude that a massive Proca field has D − 1 degrees of freedom (three degrees in D = 4), which is in agreement with the expectation. Namely, one out of the D degrees of freedom of the massive vector field A µ gets removed by the transversality condition (141).
The vacuum contribution (145) diverges and thus must be renormalized. If M is a tree-level mass, it is easy to see that the counterterm action is that a cosmological constant, 21 whose energy-momentum tensor is, By comparing (148) with (145) one can see that the minimal subtraction demands, When (148) is added to (145) 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 (35)), 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 (151) is invariant under the following local field transformations, For example, the covariant derivative (155) transforms multiplicatively under (154), When the mass parameter µ 2 > 0 in (152), then µ is the tree level scalar mass. In this symmetric phase, the model (151-152) 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. (35) with mass m → µ and the energy-momentum tensor of a massless gauge field is given in Eq. (137).
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 (154)). The gauge field A µ is assumed not to develop a condensate. Then ϕ 1 is the Higgs field and ϕ 2 is the Goldstone boson of the symmetry 'broken' by the condensation. 22 Upon inserting (156) into (152) 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 (156) via the term M ∂ µ A µ ϕ 2 . The fields decouple in the 't Hooft gauge [46], 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 (158), δ θ A µ = ∂ µ θ and δ θ ϕ 2 = −gθ(φ 0 + ϕ 1 ). From (160) 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 (161). 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. 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.
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 (159) we can which, in the massless limit M → 0, reduces to that of a massless gauge field in Fermi gauge (93). A formal solution of (163) can be written as (cf. Eq. (95)), where i∆ M (x; x ) is the massive scalar propagator (21) of mass M . The solution (95) can be recast in a more convenient form (for ξ = 1 and ξ ≥ 0) as, 23 This form of the solution is convenient since the propagator (168) splits naturally into a gauge independent transverse part and a gauge dependent longitudinal part.
We are now ready to calculate the corresponding energy-momentum tensor. By varying Eq. (159) with respect to g µν one gets (cf. Eqs. (99-100) and (143)), Acting on the transverse and longitudinal tensor structures of the propagator (168) yields, 23 When writing the solution (168) we made use of the fact that the particular solution of a sourced differential equation of the form, can be written as, In the singular (Feynman) gauge when ξ = 1 the solution of (165) is a parametric derivative, We leave the analysis of this case as an exercise to the reader.
By making use of (21) and (25)