Fluid equations for fast-moving electroweak bubble walls

The cosmological electroweak phase transition can be strongly first order in extended particle physics models. To accurately predict the speed and shape of the bubble walls during such a transition, Boltzmann equations for the CP-even fluid perturbations must be solved. We point out that the equations usually adopted lead to unphysical behavior of the perturbations, for walls traveling close to or above the speed of sound in the plasma. This is an artifact that can be overcome by more carefully truncating the full Boltzmann equation. We present an improved set of fluid equations, suitable for studying the dynamics of both subsonic and supersonic walls, of interest for gravitational wave production and electroweak baryogenesis.


I. INTRODUCTION
The electroweak phase transition in the early universe is known to be a smooth-crossover within the standard model (SM), given the measured value of the Higgs boson mass [1,2]. The addition of new particles coupling to the Higgs can turn it into a strongly first order phase transition, proceeding by the nucleation of bubbles of the true, electroweak symmetry breaking vacuum, in the initially symmetric plasma. This possibility has been widely studied because of its potential for providing electroweak baryogenesis (EWBG) [3][4][5], and gravity waves that might be observable in the upcoming LISA experiment [6,7].
An important parameter for the efficiency of baryon or gravitational wave production is the terminal speed v of the bubble walls, with baryogenesis generally favoring slower walls, while faster walls tend to produce stronger gravity wave signals. To determine v and other relevant properties of the bubble wall, within a given particle physics model, one must self-consistently solve for the perturbations to the fluid induced by the wall; these are needed to determine the frictional force acting on the wall, that brings it to a state of steady expansion.
In previous literature on this subject, quantitative study of fast-moving walls has been hampered by an apparent singularity of the fluid equations occurring at the sound speed c s = 1/ √ 3, that we will explicitly demonstrate below. This makes a microscopic calculation of the friction in such cases problematic, motivating phenomenological estimates for the friction [8][9][10][11], or else leaving aside supersonic walls altogether [12]. Complementary approaches have been used to study the ultrarelativistic limit [13][14][15]; in this work we are primarily interested in velocities v c s rather than v ∼ = 1. We argue that the apparent sound barrier is an artifact of a particular truncation of the Boltzmann equations for the fluid perturbations, and that sensible solutions exist for wall speeds up to v = 1 by making a better choice.
A similar observation was recently made in ref. [16] in the context of the CP-odd fluid perturbations that are needed to compute the source terms for EWBG, but the analogous study for the CP-even perturbations, relevant to determining the bubble wall properties, has not been done. It requires more work because the perturbation in the local temperature δτ = δT /T (not needed for the EWBG source terms) must now be included in the network. The optimal way of doing this turns out to be somewhat subtle, as we will discuss.
We start by reviewing the standard approach in section II and the pathology of the perturbations it predicts for supersonic walls. We derive improved fluid equations in section III, and in section IV the solutions of the old and new formalisms are compared for a typical background wall profile, as a function of the wall velocity v. These results are used in section V to compute the predictions for the friction term in the Higgs field equation of motion, that determines the bubble wall shape and speed. There we highlight the problems with the old approach and their absence in the new one. Conclusions are given in section VI. Formulas for the coefficients of the new fluid equations are presented in Appendix A, and the results of refined estimates for the collision terms are explained in Appendix B.

II. OLD FORMALISM (OF)
We begin by recapitulating the method that has been used in previous literature for computing the plasma perturbations [9,12,[17][18][19][20][21]. These are the deviations of the distribution function f for a given particle away from its equilibrium form, that have been parametrized as [17,18,22] where β = 1/T , γ = 1/ √ 1 − v 2 and the equilibrium part, with δX = 0, is expressed in the rest frame of the bubble wall, taken to be planar and moving to the left. µ is the dimensionless chemical potential (in units of temperature) and u is the velocity perturbation. The wall frame is convenient for expressing the Boltzmann equation since the solutions in this frame are stationary, where δf = −(df v /dX) δX ≡ −f v δX is the perturbation, (m 2 ) = dm 2 /dz for a particle whose mass depends on the background Higgs field h(z) (and possibly other fields like a singlet scalar) in the wall, and f v is the equilibrium distribution in the wall frame.
To approximately solve eq. (2), three moments are taken, by integrating over momenta d 3 p with the respective weight factors [23] 1, γ(p z − vE) and γ(E − vp z ), giving three coupled ordinary differential equations for the perturbations q ≡ (µ, u, δτ ) , that can be written in the 3 × 3 matrix form with a rate matrix Γ from the moments of the collision term C and a source S ∼ vβ 2 (m 2 ) from the Liouville operator L in (2) acting on f v . The A v matrix depends on v in such a way that A −1 v Γ becomes singular at v = c s , and has only positive eigenvalues for v > c s . By constructing a Green's function to solve eq. (3) [24], one can see that this implies that the perturbations q must strictly vanish in front of the wall for v > c s . Ref. [16] has argued that this kind of behavior is unphysical, since the fluid equations (3) describe particle diffusion, which is a physically distinct process from the propagation of sound waves. There is no reason why diffusion should be suddenly quenched in the vicinity of a supersonic wall, since some fraction of particles in front of the wall can still travel fast enough to get ahead of it.

III. IMPROVED FLUID EQUATIONS (NF)
In this section we propose a new formalism (NF) for the fluid equations, motivated by the recent paper [16]. In that work, the problem of artificial suppression of diffusion for supersonic walls was overcome, following a longestablished method of dealing with the velocity perturbation u [25][26][27]. The adoption of a specific form for u is known to lead to unphysical results, that can be avoided by instead writing the perturbations in the form where now δX omits the velocity perturbation u, which is instead encoded through δf u in such a way that To deal with other integrals involving δf , one makes a factorization ansatz for any quantity Q. This procedure was shown in ref. [25] to lead to nonsingular diffusion in front of supersonic walls, so long as one carefully evaluates the full v-dependence of A v , rather than linearizing it in v, and weighting the Boltzmann equation by the moments 1, p z /E. However ref. [25] only considered the case of CP-odd perturbations, where δτ plays no significant role and hence was omitted. Our purpose in this work is to extend those results to include δτ , whose value is needed for the full solutions to the field equations determining the shape and speed of the bubble walls. To determine this additional perturbation, a third moment is needed. We find that by choosing the weighting factor E, in addition to 1 and p z /E (all defined in the wall frame), the resulting A v matrix becomes where the dimensionless functions C m,n v and D m,n v are defined as With this choice, det A v has no singularity for wall speeds between v = 0 and 1, and it gives the desired behavior in which diffusion ahead of the wall only gets suppressed in the limit v → 1. The source term becomes In previous literature, the coefficients corresponding to C m,n v and D m,n v were usually calculated in the limit of vanishing mass (as well as only leading order in v), but we find that the variation of m 2 (z) for the relevant particles within the wall can have a significant impact on the shape of the solutions. We thus retain the full mass-and v-dependence of those functions. Moreover, it is possible to analytically determine the v-dependence by boosting to the plasma frame (see Appendix A).
We have also updated the components of the collision matrix Γ to account for the new choice of moments. The calculation of ref. [17] is improved by correcting some errors pointed out in ref. [28] and by using a Monte Carlo algorithm to compute more accurately the collision integrals. The new values of the collision terms are given in Appendix B.

IV. SOLUTIONS FOR A STANDARD MODEL-LIKE PLASMA
Next we apply the improved fluid equations to a SMlike plasma in the context of a first order electroweak phase transition. The species that couple most strongly to the Higgs boson are the top quark t and the electroweak gauge bosons. The W and Z bosons are approximated as having the same distribution functions, and we will refer to them collectively as W bosons. The remaining particles form a background fluid which is assumed to be in thermal equilibrium (µ bg = 0) at a z-dependent temperature T + δτ bg (z) [17]. Even if they are not driven out of chemical equilibrium by the phase transition, these lighter fields still play an important role in the dynamics of the bubble wall. One might also expect the Higgs boson distribution to be perturbed, but its small number of degrees of freedom makes its contribution negligible compared to that of t or W . It is therefore included in the background fluid (and similar reasoning could also be applied to additional fields not present in the SM, e.g., a singlet scalar). The complete set of matrix equations for the t, W and background components is where A t , A W , S t and S W are given in eqs. (7) and (9), using the appropriate equilibrium distribution functions. The A matrix for the background fluid is with N f and N b respectively the fermionic and bosonic number of degrees of freedom included in the background fluid (N f = 78 and N b = 19 in the SM). We evaluate A t and A W at m = 0 because all the particles in the background fluid are approximately massless. Energy and momentum conservation fixes Γ bg,t = −12Γ t and Γ bg,W = −9Γ W [9], and Γ t and Γ W are evaluated in Appendix B. To solve the system (10), one can eliminate q bg using the third equation; however the fact that µ bg = 0 makes one of the three "bg" equations redundant. We have chosen to keep the first and third "bg" component equations (corresponding to the weighting factors 1 and p z /E), since this leaves A bg nonsingular for v ∈ (0, 1). The result is whereÃ −1 bg is the inverse of the 2×2 A bg matrix, projected onto the 1,3 columns and 2,3 rows of a 3 × 3 matrix. It can be written in terms of the 3 matrices  (11): The six remaining equations take the form with To compare the new and old formalisms (denoted by NF and OF in the following) for a generic first order phase transition, we model the bubble wall using a tanh ansatz for the background Higgs field, where h 0 is the VEV of the Higgs in the broken phase and L is the wall thickness. As an example we solve eqs. (14) within the OF and NF for T = 100 GeV, h 0 = 150 GeV, L/γ = 5/T [29] and several wall velocities, using the collision rates given in [17] for the OF and the ones evaluated in Appendix B for the NF. We include a factor γ in L in order for the wall to have a constant thickness in the plasma frame. The solutions are shown in Figure 1, for a series of increasing wall velocities. One can notice that within the NF, the perturbations in front of the wall (z < 0) vanish only in the limit v → 1, as required by causality. This is not the case in the OF, whose solutions always vanish in front of the wall for v > 1/ √ 3. As argued in ref. [16], this behavior is unphysical, since there is no reason for particles not to be able to diffuse in front as long as their v z velocity component is higher than v.
As a consistency check, we observe that the linearization of the Boltzmann equation in δX and u is justified, since all the perturbations are generally well below unity in magnitude. We have tested that this condition holds for most wall parameters; the linearization starts to break down only in the extreme cases of very fast (v 0.95) and thin walls (L 1/T ).

V. CONSEQUENCES FOR WALL FRICTION
An important application is the calculation of the friction term F in the Higgs equation of motion multiplied by h = dh/dz [9],  where V eff | T is the finite-temperature potential evaluated at the unperturbed background temperature, and Here the sum is over the species t and W , and N i is the corresponding number of degrees of freedom. An exact solution to eq. (18) exists only for a specific wall velocity and shape, and so the accurate estimation of F is important for determining the wall properties. An ansatz such as (17) can give a rough approximate solution, where v and L are determined by demanding that two moments of eq. (18) vanish [9,12,17], for example We plot F (z) constructed from the OF and NF solutions in the bottom row of Figure 1. At small v, the friction predicted by NF is ∼ 50% larger, leading us to expect the NF to predict a smaller wall velocity than the OF. This difference is mainly due to our improved calculation of the collision integrals and the fact that we keep the full mass dependence of the C m,n v and D m,n v functions. In this very coarse grid on velocity space, v = 0.2, 0.5, 0.7, 0.95, the friction appears to be qualitatively different only for v = 0.95. The negative values of F on the leading side of the wall in the NF indicate a force that will tend to make the wall thicker, as one can infer from eq. (18). But the larger positive contribution means the net friction is positive, and the wall will not necessarily run away.
To better understand the differences at large v, we study F (z) region on a finer velocity grid v = 0.8, 0.85, 0.9, 0.95 in Figure 2, while also varying the wall thickness as LT /γ = 5 (see note [29]), 10, 20 between successive rows. The negative contribution to F in the NF first appears at v = 0.95 for LT /γ = 20, while for narrower walls it becomes evident at lower speeds. These are also the points where the NF friction generally starts to become lower than that predicted by the OF, which could make runaway solutions [13] more likely.
So far, the differences between the OF and NF predictions seem to be quantitative, but not dramatic. However in the vicinity of v = c s , the OF has a pathological feature that foils numerical algorithms for determining the wall parameters, and casts doubt on its physical validity. This is illustrated in Figure 3(a), where we show that the profile for the friction evolves in a discontinuous fashion as v passes through the sound barrier, and it becomes negative everywhere at v = 1/ √ 3 ∼ = 0.58. This is even more evident in Figure 3 On the other hand, the corresponding solutions in the NF are smooth, and hardly evolve over the same narrow range of v. We have argued that the smooth behavior is the one expected on physical grounds, since diffusion should not care about the sound speed. Even if the solutions of the OF fluid equations happen to give a reasonable result away from v = c s , any numerical method for determining wall parameters by searching for simultaneous solutions of the system (20) is likely fail whenever the trial value of v crosses v = c s . We find that the second moment M 2 also has a similar discontinuity.

VI. CONCLUSION
In this work we have pointed out a shortcoming at high wall speeds (v c s ) with the fluid equations that have been used, since their introduction in ref. [17], to calculate the friction F on electroweak bubble walls. We have also proposed a modification to these equations that solves the problem. It is reassuring that the two approaches give results that are not too different from each other at low wall speeds-and there the difference arises mainly because we have improved estimates of the collision rates, rather than the changes in formalism that become important at high v. Near the sound barrier and above, the differences are more significant, with our new results evolving continuously as a function of v, whereas the old ones exhibit a discontinuity in F at v = c s . The new system predicts lower friction at high v > c s compared to the old one, which is likely to lead to faster walls. At low v the opposite is true. Application of these methods to a realistic model is underway [30].
The new elements in our treatment are a different choice of weighting factors for taking moments of the Boltzmann equation, and a different treatment of the velocity perturbation. The latter has long been recognized and recently highlighted in the high-v context in ref. [16]. While there are strong theoretical motivations for the velocity perturbation, the choice of weighting factors is more arbitrary, and cannot be justified a priori.
Instead we have made a phenomenological determination, by finding a set of moments that give the expected behavior for the fluid perturbations as a function of v. One could characterize it as an educated guess, that should be validated by finding a more exact solution of the full Boltzmann equations. There are several ways one could imagine doing this. Instead of three moments and three perturbations, one could increase this number to N and look for convergence of a physical quantity like the friction with increasing N . Alternatively, one could approximate the distribution function f by taking N bins in momentum space and seeking convergence with growing N . This is an investigation we hope to undertake in future work.
Acknowledgment. We thank A. Friedlander, K. Kainulainen and D. Tucker-Smith for useful discussions. This work was supported by NSERC (Natural Sciences and Engineering Research Council, Canada) and FRQNT (Fonds de recherche Nature et technologies, Québec). The coefficients appearing in the A matrix generally depend on the local particle masses m(z)/T and the wall velocity v. They can be evaluated numerically directly from their definition (8), but it is also possible to analytically calculate their v-dependence, by making the substitution E → γ(E + vp z ) and p z → γ(p z + vE) to boost the integration variables to the plasma frame. This transforms f v to f 0 , the equilibrium distribution function evaluated at v = 0, and leaves the combination d 3 p/E invariant.
In this way, the C m,n v and D m,n v functions can be expressed as a sum (finite or infinite) of C m,n 0 and D m,n 0 , the corresponding functions evaluated at v = 0. One can show that (henceforth omitting the subscript 0) With these, it is sufficient to compute the required C m,n and D m,n at only a few values of m/T and use interpolation to quickly compute them for any m/T . The infinite series are all well-behaved: they are exact at v = 0 and v = 1 using only the first term of the series, and an accuracy of less than 1% for all v ∈ [0, 1] is achieved using a small number of terms.

Appendix B: Evaluation of the collision rates
We discuss here the calculation of the collision integrals by a corrected and improved version of the method used in ref. [17]. The collision term for a given particle species is where the sum is over all the relevant processes listed in Table I, p is the momentum of the incoming particle whose distribution is being computed, N p is its number of degrees of freedom, k is the momentum of the other incoming particle, and p , k are the momenta of the outgoing particles. |M i | 2 is the squared scattering amplitude, summed over the helicities and colors of all the external particles. The distribution functions appearing in P are Fermi-Dirac or Bose-Einstein depending the respective external particles, and the ± is + for bosons and − for fermions. P can be simplified by expanding it to linear order in the perturbations. Using the definition (1) of the distribution function with δX(p) = µ+βγδτ (E p −vp z )−δf /f v , one can show that P becomes where the sum is over the external particles not in equilibrium and the ± in front of δX is + for incoming particles and − for outgoing particles. The quantities needed for the fluid equations are the moments of C[f ]. These have the general form where we boosted to the plasma frame to get the second line. Using the substitution (6), the perturbations become in that frame Following the treatment of ref. [17], the calculation of the collision rates has been done to leading log accuracy, where it is justified to neglect the masses of all the external particles, which implies E p = p. One can also neglect s-channel contributions and the interference between diagrams because they are not logarithmic. To account for thermal effects, we use propagators of the form 1/(t − m 2 ) or 1/(u − m 2 ), where m is the exchanged particle's thermal mass. It is given by m 2 g = 2g 2 s T 2 for gluons, m 2 q = g 2 s T 2 /6 for quarks, m 2 W = 5g 2 w T 2 /3 for W bosons and m 2 l = 3g 2 w T 2 /32 for leptons [31]. The top quark collisions are dominated by their strong interactions; we include only contributions to |M| 2 of order g 4 s for t interactions. For the W bosons, we include terms of order g 2 s g 2 w and g 4 w . The relevant processes are shown with their corresponding |M| 2 in Table I [32].
To evaluate the integrals in (B4), one can first use the delta function and the symmetry of the integrand to analytically perform five of the twelve integrals. This can be done efficiently using the parametrization detailed in refs. [28,33,34]. The remaining seven integrals can be evaluated analytically using several approximations, justified to leading log accuracy. However, we have found that it is more precise to numerically compute these integrals, which can be done with a Monte Carlo algorithm. One can use a stratified sampling algorithm or VEGAS to reduce the variance, but this is generally not necessary since it only takes a few seconds to get an accuracy of ∼ 1% in most cases.
With the linearization of P[f ] made in (B3), the moments of the collision term can be written as linear combinations of the three perturbations: T Γ    τ,t = (2.26v + 4.82v 2 − 9.32v 3 + 6.54v 4 ) × 10 −3 Γ (3) τ,W = (2.48v + 6.27v 2 − 11.9v 3 + 8.12v 4 ) × 10 −3 Our results and differ from those of [17] by factors of O(1). Even taking account of the errors previously mentioned, our results are still roughly 2 times smaller. As discussed in ref. [12], this discrepancy is due to the various leading log approximations made in [17] in order to analytically evaluate the collision integrals. Either procedure is valid to leading accuracy, which gives an estimate of the theoretical uncertainty associated with this approximation. It may be worthwhile (though quite laborious) to include subleading contributions for future studies relying upon these fluid equations.