From Boltzmann equations to steady wall velocities

By means of a relativistic microscopic approach we calculate the expansion velocity of bubbles generated during a first-order electroweak phase transition. In particular, we use the gradient expansion of the Kadanoff-Baym equations to set up the fluid system. This turns out to be equivalent to the one found in the semi-classical approach in the non-relativistic limit. Finally, by including hydrodynamic deflagration effects and solving the Higgs equations of motion in the fluid, we determine velocity and thickness of the bubble walls. Our findings are compared with phenomenological models of wall velocities. As illustrative examples, we apply these results to three theories providing first-order phase transitions with a particle content in the thermal plasma that resembles the Standard Model.


Introduction
Strong first-order cosmological phase transitions predict a variety of interesting phenomena: gravitational waves [1,2,3,4,5,6], baryogenesis [7], magnetic fields [8] and many more. The thermodynamic features of a phase transition, as e.g. its critical temperature, latent heat and order parameter, can be easily determined using standard techniques [9,10,11]. Its (out-ofequilibrium) dynamic properties are instead more difficult to predict. Among these properties, for a first-order phase transition the speed of the expanding bubbles and their wall thickness are probably the most relevant.
In particular, little is known about the wall velocity in most models. Determining it hinges on quantifying the friction that is exerted by the fluid on the bubble wall. This requires a framework that captures out-of-equilibrium features of the plasma. On a technical level this can be achieved by solving Boltzmann equations. This route has been followed for the Standard Model (SM) [12,13,14,15] and for the Minimal Supersymmetric SM (MSSM) [16].
A second way to quantify the friction is to use a phenomenological approach [17,18,19,20,21,22,23,24]. In this case, the friction is modeled by an additional dissipative term in the Higgs equation of motion. This involves a free friction coefficient that is inferred by matching to the full Boltzmann treatment. Even though this approach gives reasonable results in the small wall velocity limit, it has its limitations. For example, it is not clear whether these results can be extrapolated to supersonic wall speeds. Extensive numerical simulations of the phase transitions have followed this approach [17,25,26].
In the present work we attempt to fill the gap between the full Boltzmann treatment and the phenomenological approach. Rather than adding an ad hoc term to the Higgs equation of motion, we solve numerically the Boltzmann equations. Subsequently, we fit the obtained friction in terms of the parameters characterizing a first-order phase transition. The interpolations we provide can be easily applied to models with first-order phase transitions.
From a technical point of view, the parametrization we provide assumes a given set of particle species contributing to the friction. We indeed consider a SM-like framework, in which the friction is dominated by the electroweak gauge bosons and top quarks [15]. In extensions of the SM, however, any particle that is not too heavy and is strongly coupled to the Higgs contributes. For instance, in the parameter region of the MSSM suitable for electroweak baryogenesis [27,28] (but in tension with LHC date [29,30,31] and possible magnetogenesis [32]), also stops participate in the friction [16]. Our parametrization then underestimates the friction in the MSSM. On the other hand, it well applies to models beyond the SM with not too many new degrees of freedom coupled to the Higgs. The gauge-singlet extension belongs to this class of theories. It is weakly constrained by present collider measurements [33,34,35,36] and can provide very strong two-stage phase transitions if the singlet acquires a vacuum expectation value (VEV) before the electroweak symmetry breaking [37].
The paper is organized as follows. From section 2 to section 4 we rederive the fundamental Boltzmann and Higgs equations in the Schwinger-Keldysh formalism. Compared to the semi-classical derivation in Moore and Prokopec's paper [15] (called M&P in the following), the resulting Kadanoff-Baym equations have the advantage to allow for a systematic inclusion of quantum corrections. This should be relevant to develop a unifying framework to determine the wall velocity and baryogensis [38,39,40,41]. It also permits to derive the transport equations for relativistic wall velocities. In particular, we assume that the system is close enough to equilibrium such that the flow ansatz we consider can be linearized in terms of deviations from equilibrium. This approximation does not automatically imply a small wall velocity at all. For example, if all particles are weakly coupled to the Higgs, the wall speed is large although the fluid remains close to equilibrium. In section 5 we compare our results with phenomenological approaches to the wall velocity. In section 6 we apply our Boltzmann approach to several models. We start with the SM with a (experimentally excluded) small Higgs mass in order to facilitate the comparison between our results and M&P. Then, we study the SM with a low cutoff (including additional φ 6 operators) and a singlet extension of the SM. Conclusions are given in section 7.

Equations of motion
The dynamics of the particles in the plasma is described by the Kadanoff-Baym equations whose gradient expansion reduces to the usual Boltzmann equations. The advantage of Kadanoff-Baym equations over Boltzmann equations is two-fold. First, the Kadanoff-Baym equations naturally provide the forces that act on the particles due to the Higgs background. Second, the Kadanoff-Baym equations allow for accurate treatment of spin in the case of fermionic particles.
In the context of the Kadanoff-Baym formalism, certain two-point functions encode the dynamics of the system. In particular, the Wightman function G < encodes the particle distribution functions. At leading order, the Kadanoff-Baym equations for a scalar degree of freedom in the gradient expansion read where the term "coll" summarizes the collision contribution. The first equation is the so-called constraint equation and encodes the fact that the Wightman function G < can be expressed in terms of the particle distribution function f ( p, x): Using this ansatz and the identity eq. (2) leads to the relation The first term corresponds to free floating of the particles. In the nonrelativistic limit it reduces to the usual kinetic term of the Boltzmann equation: The second term describes the force acting on the particles. To understand its effect, we can imagine a Higgs background that is constant in space and only depends on time. In this case, we expect the three-momentum of the particles to be conserved. So the energy p 0 has to change in order to ensure the t-dependent on-shell condition p 2 = m 2 (t). This behavior is well reflected by the force term as it admits solutions of the form f ( p, In the analysis of the wall velocity, it is the force contribution that drives the plasma out of equilibrium. The complexity of the problem however lies in the collision terms. They depend on the interactions between the particle species and will be discussed in more detail in section 4.

Symmetries and conservation laws
Before studying the effect of the collision terms, it is useful to discuss some symmetries of the problem. To this aim, we consider the energy-momentum tensor and the particle current. By the conservation of the former we deduce the equation of motion of the background.

Four momentum
The spatial variation of the classical background can be seen as a bubble wall separating the inner (electroweak broken) and external (electroweak unbroken/symmetric) phases. Ultimately, we are interested in the velocity of this wall once the growing bubble reaches a steady expansion regime. In order to quantify this speed, one needs the equation of motion of the Higgs in the plasma 1 . A simple way of achieving it is to use the energy-momentum conservation of all particles in the plasma and the Higgs background (the expansion of the Universe can be neglected during the phase transition): with n running over each species in the plasma. The energy momentum tensor of the plasma of each species n can be expressed as 2 1 Hereafter we assume only the Higgs field to acquire a VEV and to act as a classical background. See section 6.3 for a more general discussion. 2 In the following we focus on bosonic degrees of freedom in the plasma. However, the conclusion does not change for fermions.
On the other hand, the energy momentum tensor of the classical field background is Hence, the divergence of this energy-momentum tensor reads Finally, by plugging eqs. (9) and (11) into (7), one obtains the Higgs equation of motion Note the collision terms are absent. They indeed cancel out when the sum over all species is performed. This is a consequence of the energy-momentum conservation in the decay and scattering amplitudes and can be checked explicitly once the model is specified. In thermal equilibrium, the last term in eq. (12) can be easily related to the thermal contribution to the Higgs finite-temperature effective potential V (φ, T ). This can be verified by using the relations where "pressure" stands for the pressure of a bosonic gas in the plasma and f (E) does for the Boltzmann distribution (i.e. f ( p, x) → f (E) in the equilibrium limit). Therefore, splitting the distribution functions into an equilibrium part plus some deviations δf n yields The last contribution is the so-called friction term. We remark that this result, which is based on the Kadanoff-Baym equations, reproduces the finding in M&P where the Higgs equations was obtained in the WKB approximation.

Charges
The system conserves electric charge and also iso-spin in the symmetric phase. This should be reflected in the equations. In particular, the particle current should be conserved in the sum over species (weighted by the corresponding charges). Notice that in the conservation equation of the charge, the force term does not enter. Indeed, after partial integration, it turns out that This reflects the physical picture that the force modifies the trajectory of the quasi-particles but does not change their charge. The collision terms on the other hand contain decay and annihilation processes that change the individual particle numbers but also conserve the total charge.
4 Transport equations of the plasma components

The fluid approximation
In order to solve the equation of motion (15) one needs to determine the friction contribution. This requires to identify the correct particle distribution function. As in M&P, we consider the flow ansatz where the four-velocity u µ (x), the chemical potential µ(x) and the inverse temperature β(x) are space dependent. In the limit of negligible space dependence, it reproduces the usual Boltzmann distribution in the frame boosted by the four-velocity u µ . Contrarily to ref. [15], in which this ansatz was first used, here we do not require small fluid velocity u µ , although we still assume small spatial dependence (the consistency of this assumption will be checked a posteriori). We can hence use for the individual particle species and linearize in the following in the fluctuations δτ , δu and δµ when necessary. Temperature changes and the chemical potential are encoded in the dimensionless quantities δτ and δµ in units of the temperature. Changes in the fluid velocity are encoded in δu that fulfills u µ δu µ ≃ 0 in order to achieve the correct normalization for the four-velocities u µ and u µ + δu µ . The space-independent part of each quantity is fixed at its value far outside the bubble wall. In particular, the constant part of the chemical potential can be neglected [15]. In the following we need two different types of averages and we define Using these definitions, the out-of-equilibrium densities can be expressed in terms of fluctuations and equilibrium densities. For example for the fourcurrent one finds in leading order of the fluctuations where the zero subscript denotes the equilibrium quantities with a distribution function f for fixed background values (δu = δτ = δµ = 0). These functions are still space-time dependent due to their mass dependence. The equations of motion of the system can in principle be obtained from (25). However, using the properties under Lorentz transformations of the different functions and dimensional analysis, they can be brought to a form that allows for a more intuitive interpretation (details are given in appendix A). Subsequently, the divergences of the four-current and the energy-momentum tensor turn into and We solve these equations in the planar wall approximation. We also work in the wall frame (with the z-axis orthogonal to the wall and oriented towards the broken phase) in which the plasma velocity and its fluctuation are u µ = γ(1, v w ) and δu µ = δvū µ = δvγ(v w , 1) 3 . As we focus on the steady velocity regime, the substitutions m 2 (x) → m 2 (z), u µ ∂ µ → γv w ∂ z andū µ ∂ µ → γ∂ z apply. The linearized eqs. (26) and (27) can hence be expressed as 4 A · q ′ + coll = S , where q = (δµ, δτ, δv) and the prime denotes the dimensionless derivative q ′ = γβ∂ z q. The matrix A and the source S have the form where the coefficients c i depend on the spin statistics of the species we are dealing with. Working at lowest order in m/T as in M&P, for bosons [having p = π 2 T 4 /90 and n = ζ(3)T 3 /π 2 ] one finds 5 whereas for fermions [with p = 7π 2 T 4 /720 and n = 3ζ The linearized fluid equation (28) agree with the system obtained in M&P in the limit of non-relativistic wall velocities.
Notice that in the limit v w → 0 the matrix A has one vanishing eigenvalue. On the other hand, the collision terms do not go to zero. It is then ensured that in this limit, q = 0 is the unique solution. 4 Linearizing the equations is justified (among other conditions) when the change in mass is small compared to the temperature, m 2 T 2 .
5 One might wonder if it is feasible to neglect the mass dependence in those coefficients. After all, the leading coefficient only corresponds to the mean-field approximation while the phase transition relies on the interplay between the mean-field and higher contributions. However, these terms are only comparable because the zero-temperature contribution to the φ 2 operator almost cancels the mean-field contribution close to the critical temperature. In the matrix A only the finite temperature contributions are relevant. No cancellation occurs and higher orders can be neglected.

Standard Model-like plasma content
In principle, in order to determine the friction, one has to solve a system of differential equations (28) (as many as the number of species) coupled to each other via the collision terms. Nevertheless, we are interested in theories in which the particle content of the thermal bath resembles the SM one. The species that are relevant during the phase transition are therefore the electroweak gauge bosons (simply called W bosons hereafter) and top quarks. The remaining particles are not driven out of equilibrium and act as a background 6 . Quarks and gluons are the large portion of the background particles. They are strongly coupled and hence share the same plasma fluctuations. Furthermore, their chemical potential vanishes since the gluons quickly equilibrate. The equations of motion of the background can be deduced from the two relations arising from energy-momentum conservation.
After linearizing the collision terms [15], the fluid equation (28) applied to the W and top fields reads whereas for the background particles it leads to The quantities A W and A t are given by and A f are defined as the matrix A with the coefficients (30) and (31), respectively. Considering a SM-like background (with decoupled right handed neutrinos), one has A bg = 19A b + 78A f . The values of Γ W and Γ t are summarized in appendix B. Since the total energy momentum is conserved, it follows that N W Γ W + Γ bg,W ∝ (1, 0, 0) and N t Γ t + Γ bg,t ∝ (1, 0, 0), with N W = 9 and N t = 12. This system of differential equations can be more easily solved by removing q bg from (32) and (33) and then determining the background equation by integration. This shows that the fluctuations q W and q t vanish by construction far inside and far outside the bubble, where both sources S W and S t vanish. On the other hand, the background fluctuations q bg cannot vanish on both sides. As previously mentioned, we chose to match the solution in the symmetric phase in front of the wall. Notice that due to energy-momentum conservation, the absolute change of q bg along the wall cannot depend on the wall shape. In our approximation it can be determined by knowing the change of the W-boson and top masses in units of the temperature.

Higgs equation of motion
The equation of motion of the Higgs (15) can be linearized as well. From the expansion of the fluid ansatz (18) one obtains The system of differential equations (32), (34) and (35) constitutes the basis of our numerical analysis. We solve it by means of the two-parameters wallshape ansatz where φ 0 and L are respectively the VEV of the Higgs in the broken phase and the wall thickness, both during the bubble expansion. This ansatz seems particularly appropriate for weak phase transitions. In this case, the profile of the tunneling bounce (i.e. the instanton solution connecting the two phases [51]) is very similar to eq. (36) and such a shape is expected to be kept during the bubble evolution. Instead for very strong phase transition the bounce profile may qualitatively differ from eq. (36). Nevertheless, it seems reasonable that the bubble wall acquires the above configuration once it approaches the steady velocity regime. This is also seen in recent simulations [25] and we assume this shape in our analysis. In order to implement the constraint (35), we take the moments These relations have a physical interpretation (cf. section 3.1). Eq. (37) declares that in the steady velocity regime, the total pressure on the wall vanishes. Its equilibrium part is the potential difference

The shock front in the deflagration mode
We have seen in the last section that by minor modifications, the equations obtained in M&P are also valid in the relativistic regime as long as the phase transition is weak enough. However, there is a further reason why only slow walls are considered in M&P. When v w is equal to the sound speed c s , the sign of one eigenvalue in the fluid system is flipped, and therefore the whole dynamics changes. Additionally, the linearization for the background fields ceases to be valid at v w ≈ c s . In the deflagration mode, a shock wave builds up in front of the expanding bubble [43]. Accordingly, the fluid velocity and temperature in front of the wall differ from those of the symmetric phase. Indeed, the shock wave sets the fluid in motion and heats the plasma. This decreases the pressure difference experienced by the Higgs and reduces the latent heat released in the plasma. This effect can be strong enough to dominate the dynamics of the bubble expansion, whose description can be inferred only from hydrodynamic considerations [44]. To deal with this issue, we follow the procedure of ref. [17,45] which matches the plasma velocity and enthalpy in front and behind the wall. Unlike the analysis in M&P, we solve the non-linear equations and do not rely on small fluid/wall velocities in this step.

Phenomenological approaches
In this section we discuss phenomenological approaches to the bubble wall friction. In this kind of study, the Higgs equation (15) is assumed to be effectively described as The effective friction η may involve an explicit dependence on φ and/or v w . Typically, it is deduced either by a matching to the existing results in Boltzmann treatments [23] or by the relaxation time approximation [22]. Often, it is supplemented by a further equation that sets the temperature variation in the wall and that may be derived from the energy-momentum conservation of the plasma (assumed to be in local equilibrium [17,23]). Depending on the parameter region and the level of sophistication, eq. (39) can reproduce almost all features of the full Boltzmann treatment. It has however its limitation. The most striking is that the friction force scales with the wall thickness as 1/L. Moreover, the dependence on the wall velocity can be quite involved. For instance, in the highly relativistic regime, v w → 1, the term u µ ∂ µ φ in eq. (39) is enhanced by a Lorentz factor, and leads to finite wall velocities even for extremely strong phase transitions. It is however known that bubble walls can enter a runaway regime [46]. Phenomenologically, this can be cured by introducing a 1/γ w factor in η(φ, v w ) [23] (or an even more complicated dependence [24]). This seems quite ad hoc and one can instead wonder whether the discrepancy has deeper origins. As we will see, the dependence on the wall velocity already starts to become quite non-trivial nearby the speed of sound.
In the following we demonstrate in which cases the Boltzmann treatment agrees with the phenomenological approach. We consider a system of equations obtained by integrating the moments (37) and (38). The nonequilibrium part contains a first contribution coming from the fluctuations δµ and δτ , and a second contribution coming from the background fields δτ bg . These scale differently in terms of v w , φ 0 /T and L, and we hence treat them differently. Our parametrization then reads where the quantities f f l , g f l and f bg , g bg are the fluid and background functions (depending on v w , φ 0 /T and L T ) and W is given by The functions f f l , f bg , g f l , g bg heavily simplify for thick bubble walls. More specifically, under the condition the kinetic term in the fluid equations (32) and (33) can be neglected and the background equation (34) yields (in the last relation we neglected the z-dependence in c 1 ). Hence, f f l and g f l scale as 1/L T , whereas f bg and g bg are independent of it. Their dependence on L T , v w and φ 0 /T are shown in figures 1-3 for several phase transitions parameters. By inspecting the eigenvalues of the scattering terms Γ, one can observe that the criterion (42) amounts to L T ≫ 20. This is well reproduced by the numerical results in figure 1 that show how the friction terms approach above scaling. Generally, the friction terms of the background fields display the proportionalities f bg ∝ m 4 ∝ φ 4 and g bg ∝ φm 4 ∝ φ 5 . The friction terms of the fluid have an additional dependence on φ via the mass-dependent coefficient c 1 in (29), which is absent in the background due to its vanishing chemical potential δµ bg (see the discussion on the background in section 4.1). A fit to the numerical data yields f f l ∝ φ 7/2 and g f l ∝ φ 9/2 . This behavior is shown in figure 2.
Using the above proportionalities, a reasonable fit for the parameterizations (40) turns out to be (for L T ≫ 20) The rather complicated dependence on the wall velocity can be disentan- gled and traced back to different origins. The quadratic corrections of the friction components in terms of wall velocity come from the dependence of the eigenvalues of the system. Also the factor 1/(c 2 s − v 2 w ) in the background field arises from the eigenvalues of the matrix A in (43). On the other hand, the √ γ enhancement in the fluid functions is due to the suppression of the collision terms. This enhancement suggests a divergent friction, but for very fast walls, γ ≫ 10, the friction approaches a constant value, which is just given by the "free fluid solution" (i.e. Γ → 0). In contrast, the background functions are suppressed by an additional factor 1/γ compared to the fluid functions. This can be deduced from the fluid equations, which imply that in the ultra-relativistic limit the background fields have to be space-independent due to vanishing source and collision terms. In this case equilibration to the true temperature and fluid velocity only happens far behind the bubble wall. Interestingly, the background contribution to the friction is negative for subsonic wall velocities. In fact, this term encodes the impact of the temperature variation on the Higgs field, namely The leading contribution to d 2 V T /dτ dφ comes from the mean-field term V T ∝ m 2 T 2 and is positive. For the deflagration mode, the temperature drops across the wall and makes this term (49) negative. For supersonic wall velocities, the sign changes and this term acts as an additional friction, hindering the wall expansion.  For small velocities, the usual friction dominates, but the contribution from the background grows with an additional factor 1/(c 2 s −v 2 w ). One curious consequence of this behavior is that there is a wall velocity with maximal friction and hence a maximal velocity in the deflagration mode. This is a dynamically effect and not related to the considerations about entropy increase in [17]. For example, we find the numerical values v w < 0.37 or v w > 0.74 This effect seems not be present in the phenomenological approach (39). In this case, even if the change in temperature is accounted for, the terminal velocity in the deflagration mode seems to be the speed of sound and not significantly below it [23]. Likewise, there is a minimal velocity for the detonation mode. Hence, there results a gap (in terms of pressure difference ∆V ) for which no solution exists to the linearized Boltzmann equations. So, even if the friction is well represented in the phenomenological approach eq. (40), the contribution from the background is quite different. Some examples for the fluctuations are shown in the figures 4 to 6. The first two plots show a deflagration and a detonation in the thick wall regime. The last example is a detonation with a relatively thin wall. The first important point is that the fluid and background fluctuations are small in all cases, what justifies the linearization of the equations in (32). This is even true for wall velocities close to the speed of sound or supersonic wall velocities. Next, we see that depending on the parameters, the profiles in the wall can be quite different than in the relaxation time approximation (43). In particular, the background fields do not need to be monotonic.
In some baryogenesis analyses, the wall thickness is not derived from the Higgs equation (39) but taken from the tunneling bounce profile. This amounts to neglecting the friction term in eq. (40). It is clear that this approximation breaks down for very thick walls, since in this regime the friction dominates over the Higgs kinetic terms. Depending on the velocity of the bubble wall, this procedure can lead to thicker or thinner walls than in dynamical treatments including friction.
In conclusion, the phenomenological model (39) is only a good description in the regime where the bubble walls are thick [cf. eq. (42)] and the wall velocity is much below the speed of sound. One can use eqs. (40) and (45)- (48) as an improved phenomenological model. To this aim, one can proceed by: i) computing the nucleation temperature of the phase transition, e.g. via the bounce analysis; ii) guessing the value of v w ; iii) determining the shock front and the temperature in front of the wall; iv) calculating the wall thickness from the second constraint in eq. (40) (the result is rather insensitive to v w ); v) checking whether the first equality in eq. (40) is satisfied and, if not, repeating the procedure from ii .

Runaway regime
In [23] it has been argued that the analysis of runaway walls can be used to deduce the friction coefficient η and that this procedure leads to very similar results as the matching to the solutions of the Boltzmann equations. In order to find finite friction in the limit v → 1, this paper assumed an additional    explicit factor 1/γ in η that cancels the factor γ present in the four-velocity u µ . The runaway regime results when the pressure difference from the fluid is too low to compensate for the pressure difference from the Higgs field in the wall. In the highly-relativistic regime, the pressure difference from the fluid can be readily evaluated [46]. It is equal to the free energy difference in the mean-field approximation (evaluated using the temperature in front of the wall). Hence, the friction approach (39) can lead to the same runaway criterion only if the contributions to the finite temperature potential beyond mean-field equals the friction term.
Interestingly, the Boltzmann approach leads to a finite friction in the v w → 1 limit. This can be seen by inspecting eqs. (28) and (29). For γ w → ∞, their kinetic and source terms are linear in γ w , whereas the collision terms are not. Thus, one can neglect the collision terms in the Boltzmann equations. Naively, one may wonder whether this depends on our notation since we absorbed a factor γ w in the definition ofū µ . It has instead to be noticed that the velocity fluctuations do not enter the Higgs equations, and this is why the friction is finite in this limit. Nevertheless, our fluid ansatz is not justified for such a regime. It is not guaranteed that the friction calculation actually leads to the results in [46] that uses the proper particle distribution functions of the highly-relativistic limit. In fact, the two results scale quite differently. For example, in the SM the leading terms beyond mean-field are the thermal cubic contributions. Only bosonic degrees of freedom contribute to them which is quite opposite to the friction terms, where fermions yield numerically even larger contributions. In this light, it seems plausible that the agreement found in reference [23] is specific to the considered model, namely the SM with low cutoff (see section 6.2).

Applications to models
In this section we apply the method we have previously discussed, to calculate the wall dynamics in simple models providing first-order phase transitions. The first model we consider is the SM with a small (and experimentally excluded) Higgs mass. The analysis of this scenario allows to compare our approach with the original calculation in M&P. It also permits to point out some peculiarities of those first-order phase transitions that rely solely on a temperature-induced cubic term in the free energy.
The second model we analyze is the SM with a low cutoff. In this framework we will check the consistency of our results with those of ref. [23] where a phenomenological approach is employed.
Finally, we discuss the phase transition in the gauge-singlet scalar extension of the SM. This framework has enough free parameters to disentangle the effect of the pressure due to the finite-temperature potential from that of the phase transition strength. Qualitative behaviors not emerging in the previous two models will be highlighted.
The numerical results we present are based on the following procedure. We use the bounce method to determine the bubble action S(T ) and we define the nucleation temperature T n such that S 3 (T n )/T n = 140 [51]. This also provides L nucl , the thickness of the wall at the nucleation time. (L nucl is defined such that the integrations of the bounce profile and of the function (36) with L = L nucl , are equal). The wall thickness satisfying the fluid constraints is dubbed L dyna in the following. Depending on the wall velocity, the temperature in front of the wall varies. We calculate this temperature using the methods discussed in section 4.4 and denote it by T w .

The Standard Model with light Higgs
For a first model, we consider the SM with a Higgs mass m h ≤ 70 GeV and compare with the findings of M&P. The comparison has of course only illustrative purposes: for a Higgs mass in agreement with the LHC measurements [47,48] the electroweak phase transition in the SM is a crossover [49,50].
We implement the high-temperature expansion of the one-loop effective potential [51] where the Coleman-Weinberg corrections are included as follows: with v 0 = 246 GeV, a f ≃ 14 and a b ≃ 223. For such a potential, we obtain the numerical results presented in figure  7. As expected, the thicknesses L nucl and L dyna (upper right panel) are closer for weak phase transitions (cf. central left panel). Moreover, the wall velocity (upper left panel) is rather constant, but the thickness shrinks with stronger phase transitions. This behavior of v w is due to the fact that the dependence of the friction on φ 0 /T and the wall thickness LT is almost identical to the one of the pressure difference along the wall. Concerning the phenomenological approach in the SM, the wall thickness is sufficiently large and the wall velocity is sufficiently small to make the phenomenological approaches (39) or (40) feasible.
Finally, our numerical findings are close to the results found by M&P but are not identical. The reasons for this discrepancy can be traced back to the determination of the shock front (which we treat non-linearly unlike M&P), the mass dependence in c 1 (that is neglected in M&P) and the slightly different potential.

Standard Model with a low cutoff
As a second example we chose a simple extension of the SM. It contains the SM supplemented by new physics coming into play at a scale M and producing an effective φ 6 operator at low energy. This framework allows for strong first-order phase transitions with a Higgs mass compatible with present LHC data [55]. The additional content is chosen such that it affects the Higgs potential but does not contribute to the friction.
In order to compare our results to those obtained via a purely hydrodynamic approach to wall velocities, we consider a similar framework as analyzed in ref. [23]. In this case the high temperature expansion of the Higgs effective potential is given by where g 1 , g 2 are the electroweak gauge couplings, h t is the top-Yukawa coupling, and Q is the renormalization scale fixed at Q = m t . The zerotemperature part together with the renormalization conditions is used to determine the µ and λ parameters.
Our numerical results are displayed in figure 8. They are very similar to those obtained in ref. [23] by means of the phenomenological approach based on eq. (39). The discrepancy is indeed smaller than a few percent, as table 1 shows. In comparison, the wall velocity we obtain is slightly smaller (larger) for light (heavier) new physics, namely, M ∼ 900 GeV (M ∼ 700 GeV). In particular, the wall thickness is small and the wall velocity is rather close to the speed of sound, which reduces the friction. At the same time, we find slightly weaker phase transitions in our potential, what reduces the wall velocity.

Singlet model
The third model we consider is the SM with an additional scalar. The extra scalar is an electroweak singlet and thus is only coupled to the Higgs. In order to reduce the free parameters we take a model with a manifest Z 2symmetry. A peculiar feature of this model is that it allows for very strong phase transition already in the mean field approximations [37]. The Higgs potential of the model can be parametrized as [37] V sing (h, s, Its temperature-dependent contribution is taken in the mean field limit, and it is absorbed into the quadratic couplings: In this parametrization, the strength of the phase transition at the critical temperature T c turns out to be Two of the five free parameters in the potential are fixed by imposing φ 0 (T = 0) = v 0 and m h = 125 GeV. The remaining three parameters can be expressed as a function of the singlet mass at zero temperature m s , the ratio φ 0 (T c )/T c and the coupling λ m .
We restrict ourselves to the parameter space where the singlet acquires a VEV before the electroweak phase transition. Because of baryogenesis, we also require the transition to be rather strong and the wall velocity to be subsonic. In spite of these constraints, the allowed parameter space is still too large to allow for an extensive analysis. We then focus on some benchmark points where the free energies in the broken and unbroken phases at T = T n are not much different [T n is determined by the two-dimensional bounce in the (φ, s) plane]. This automatically avoids those configurations leading to runaway bubble expansions [46]. At the same time, it makes the bubble walls relatively thick.
Our benchmark points are listed in table 2. The corresponding numerical findings are also reported. As expected, since we are implicitly working in the thick bubble regime, the bubble wall thickness L dyna deviates from its initial value L nucl (see section 4.3). We have also quoted the effective friction η calculated in the phenomenological approach (40). In this model, the latent heat and the wall thickness are more or less independent parameters. Hence, we can use the model to test a regime where the wall velocity is not too large but the wall thickness is relatively small. As expected from the discussion in section 5, we find that in this regime the non-trivial dependence on the wall thickness can lead to an effective friction coefficient that differs by up to 50% from the one obtained in the SM (see figure 7).

Conclusion
Considerable progress has been made recently to determine the asymptotic wall velocity of bubbles generated during a first-order phase transition. Two different procedures have been developed to address this issue: the full Boltzmann treatment [15,16] and the so-called phenomenological approach [17,18,19,20,21,22,24,23]. In the former case, the dynamics of the fluid components is determined close to the interface between the two plasma phases and used in the equation of motion of the Higgs. In the latter case the equation of motion of the Higgs is supplemented by a phenomenological friction term without considering its microscopic origin. In the present paper we have described how to extend the regime of applicability of the Boltzmann approach and how to make contact to the phenomenological approach from first principles. Concerning the Boltzmann treatment, we have adopted the Schwinger-Keldysh formalism to rederive the fluid and background equations of motion including the correct relativistic behavior. In contrast to former work in the literature [15], we have not assumed small wall velocities but only small deviations from thermal equilibrium. Small deviations do not necessitate small wall velocities, as we have explained in the main text and showed in some examples. The Boltzmann treatment with these new equations can hence be applied also to models leading to supersonic detonation fronts. Subsequently, the new fluid and background equations have been solved numerically assuming a thermal bath populated by a Standard Model-like particle content at electroweak scales. In the Boltzmann approach, the Higgs equation contains two qualitatively different contributions from the plasma.
The first contribution comes from the particle species that are driven out-of-equilibrium due to interactions with the wall. It depends parametrially only on the Higgs wall thickness, the strength of the phase transition (in terms of φ/T ) and the wall velocity. This term is parametrized in the phenomenological approach. It turns out that usually the phenomenological approach is well justified as long as the (Lorentz-contracted) wall thickness is much larger than the mean free path of the particles in the plasma. In this case the relaxation-time approximation can be used to solve the Boltzmann equations, what justifies the phenomenological approach. Above criterion amounts for a SM-like particle content to a constraint on the wall thickness in terms of the temperature, L T ≫ 20. If the phase transition produces thinner walls, the friction can be reduced considerably. For example we find cases with 30% less friction at L T ≃ 10. However, overall this first term in the Boltzmann equation is reproduced quite well.
The second term comes from the majority of particle species that is not driven out of equilibrium directly but nevertheless their temperature and velocity changes due to the latent heat that is released into the plasma and finally distributed under all degrees of freedom. In the phenomenological approach, this contribution is determined using local energy-momentum conservation but ignoring out-of-equilibrium effects. For small wall velocities, these background fields are not important but their impact increases if the wall velocity approaches the speed of sound. In this regime, the system changes from deflagrations to detonations. We have found that the out-ofequilibrium effects can have a large impact on the background fields. In extreme cases and depending on the latent heat, the Boltzmann approach can lead to no static bubble wall solutions at all. One main difference to the phenomenologial approach is that this gap is quite substantial.
In order to benchmark the described improvements, we have analyzed some specific models. The case of the Standard Model with a Higgs mass below 70 GeV shows that, in the non-relativistic regime with weak phase transitions, our Boltzmann treatment agrees with former results using an expansion in small wall velocities [15]. Some minor differences arise due to our non-linear treatment of the shock front of the wall. Next, the Standard Model with a low cutoff allows for a comparison between our outcomes and previous results obtained by means of the phenomenological approach [23]. Also here, no substantial discrepancies emerge. This is due to the fact that the wall velocity in these two models is relatively low. Besides, the latent heat and the strength of the phase transition are connected to a common scale (namely the Higgs mass and the cutoff, respectively). When these two quantities are decoupled from each other, the effective friction coefficient can vary more strongly even for relatively small wall velocities. This happens, if the model parameters allow to make the Higgs bubble wall thinner (without approaching the speed of sound). This can occur for instance in the singlet extension of the Standard Model, that we discussed as a last model.
In summary, the Boltzmann approach to bubble wall velocities was so far limited to small wall velocities. In this regime phenomenological models to bubble wall friction lead to quite good results (once the friction coefficient is known and the bubble wall is not too thin, L T ≫ 20). In the present work, we generalized the Boltzmann approach to larger wall velocity. In this regime, the phenomenological approach does not perform so well. The main reason is that the latent heat of the phase transition is first released into the particle species that couple strongly to the Higgs and then distributed under all remaining degrees of freedom by scatterings. The latter process is not represented well in the phenomenological approaches. Therefore the full Boltzmann treatment has to be used when the wall velocity approaches the speed of sound.

A Lorentz properties and dimensional analysis of the kinetic equations
The relation (25) displays certain relations between the equilibrium densities that are not explicit in their definitions. Assume that u µ is a four vector that is not constraint to the normalization u 2 = 1. Then one has At the same time, Lorentz invariance implies for J µ the form J µ = u µ n. The equilibrium distributions are a function of the combination u µ β only, such that derivatives with respect to u µ can be written in terms of derivatives with respect to T . For example Furthermore This is consistent with (25), since δτ encodes just the fluctuations in the temperature. Likewise one finds and βu µJ µ 0 = −T ∂ T N .
The density M is by construction symmetric in the three indices what implies the obvious relation T ∂ T p = ω. Also in this case, contraction with the velocity reproduces the derivative with respect to temperature βu λM µνλ 0 = −u µ u ν T ∂ T ω + g µν T ∂ T p = −T ∂ T T µν 0 .
In fact, there are more consistency relations hidden related to the mass dependence that can be made explicit by analyzing the dimensionality of the various functions. For example, the four current is of dimension three, so it fulfills the relation Analogously, one finds for the energy-momentum tensor the relation Actually, these relations hold also for the different components of the densities, namely the enthalpy ω and the pressure p (dimension 4) as well as the densities n (dimension 3) and N (dimension 2). These relations are important to establish a local equilibrium in case of a static wall. The equation for the current then reads ∂ µ J µ 0 = ∂ µ m 2 u µ ∂ m 2 n = coll = 0 , which is automatically fulfilled for static walls due to ∂ µ m 2 u µ = 0. On the other hand, the corresponding relation for the energy-momentum tensor reads with ∂ µ T µν 0 = ∂ µ m 2 (u µ u ν ∂ m 2 ω − g µν ∂ m 2 p) .
The first term vanishes again automatically for a static wall. The second cancels against the term involving N only thanks to (73) and the relation g µν T µν = m 2 N. These relations are not very illuminating and also to remove derivatives with respect to m using these equations are not really simplifying the system. Nevertheless, these relations are important for two reasons. First, if one wants to achieve explicit energy momentum conservation, these relations help to guide which terms are important. Second, these relations also lead to cancellations in the non-equilibrium case. The easiest way to see this is by looking at the original equations for the Wightman functions (1). If the fluid ansatz is used, all the derivatives acting on the mass (which only show up in the on-shell delta function) are canceled by corresponding momentum-derivatives acting on the on-shell delta-function. The only terms that remain involve derivatives acting on X.

B Linearized collision terms
The matrices Γ W , Γ t containing the collision terms appearing in the fluid equations in the wall frame are given by where the index q stands for the particle species. The following linearized collision terms were take from [15]. The numerical values for bosons and fermions are Γ µ f 1 = 0.00899 T, Γ µ b1 = 0.00521 T , Γ µ f 2 = 0.01752 T, Γ µ b2 = 0.01012 T , Γ δT f 2 = 0.06906 T, Γ δT b2 = 0.03686 T , The corresponding diagrams have been calculated in the leading-log approximation. In the case of the W-bosons this approximation works very well since the gauge coupling is small. The uncertainties of the top quark collision terms are much larger due to their color charge. The errors were estimated to be up to 50% [15]. In order to check the impact on the friction we multiplied the matrix Γ t by a factor of χ = 0.5 − 1.5. The friction from the fluctuations in the thick wall regime scales with approximately 1+χ 2 . The reason for this is that the contribution to the friction from top quarks and W-bosons is of the same order. The impact on the background contribution is not as easy to assess, but much smaller than the effect on the fluctuations. From this we conclude that even if the collision terms of the top quarks are off by 50% the resulting correction to the wall velocity is just of order 25%.