Bubble Expansion and the Viability of Singlet-Driven Electroweak Baryogenesis

The standard picture of electroweak baryogenesis requires slowly expanding bubbles. This can be difficult to achieve if the vacuum expectation value of a gauge singlet scalar field changes appreciably during the electroweak phase transition. It is important to determine the bubble wall velocity in this case, since the predicted baryon asymmetry can depend sensitively on its value. Here, this calculation is discussed and illustrated in the real singlet extension of the Standard Model. The friction on the bubble wall is computed using a kinetic theory approach and including hydrodynamic effects. Wall velocities are found to be rather large ($v_w \gtrsim 0.2$) but compatible with electroweak baryogenesis in some portions of the parameter space. If the phase transition is strong enough, however, a subsonic solution may not exist, precluding non-local electroweak baryogenesis altogether. The results presented here can be used in calculating the baryon asymmetry in various singlet-driven scenarios, as well as other features related to cosmological phase transitions in the early Universe, such as the resulting spectrum of gravitational radiation.


Introduction
Cosmological phase transitions in the early Universe are interesting for a variety of reasons. They can produce observable gravitational radiation [1][2][3][4], seed primordial magnetic fields [5,6], affect the abundance of thermal relics [7,8], and otherwise play an important role in the cosmological history of the Universe [9]. Perhaps most notably, such phase transitions can give rise to a viable mechanism of baryogenesis, provided the transition is first-order.
A first-order phase transition can occur in the early Universe when two vacua of the theory coexist for some range of temperatures. If this is the case, an energy barrier exists between the two and the system can transition to the state with lower free energy via quantum tunneling or thermal fluctuations [10][11][12][13]. Physically, this corresponds to the formation of spherical bubbles in the ambient metastable vacuum. These bubbles grow and can reach a steady state expansion velocity. Inside the bubble, some subset of the scalar fields have non-zero vacuum expectation values (VEVs), while outside they do not. If the non-vanishing condensate breaks the SU (2) L gauge symmetry of the Standard Model (SM), as at the electroweak phase transition (EWPT), non-perturbative sphaleron transitions, which violate B + L, will be quenched inside the bubble and active outside. These processes can act on chiral charge currents diffusing in front of the wall to source a net baryon asymmetry. If the sphaleron rate is significantly suppressed in the broken electroweak phase, the asymmetry can be frozen in once captured by the expanding bubble. Roughly, this requires [14,15,116] h T 1, where T is a temperature associated with the phase transition and h is some combination of fields charged under SU (2) L usually identified with the Standard Model-like Higgs. This picture is known as "non-local" or "transport-driven" electroweak baryogenesis (EWB) and is an elegant explanation for the origin of the observed baryon asymmetry of the Universe [17][18][19][20][21][22]. Clearly this mechanism relies on several requirements beyond sphaleron suppression in the broken phase. One must ensure a significant amount of CP -violation to source the chiral charge currents in front of the wall. This is well-known, has been studied extensively, and has motivated several experiments in search of CP -violating signatures, such as permanent electric dipole moments [23]. There is another, somewhat less-appreciated, requirement for successful baryogenesis through this mechanism: bubbles must expand slowly enough for sphalerons to convert a significant fraction of the CP -asymmetry to a net baryon density [19][20][21]24]. In other words, the diffusion of the chiral plasma excitations must be efficient in front of the bubble wall 1 .
Precisely how slowly the wall must move in this scenario depends on several factors, including the amount of CP -violation and the details of diffusion in front of the bubble.
However, it is generally the case that bubble expansion must at least be slower than the speed of sound in the plasma 2 , c s ∼ 1/ √ 3 ≈ 0.58. Otherwise, diffusion in front of the wall will be very inefficient. Even if the wall moves subsonically, the predicted value of the steady-state wall velocity, v w , is an important input into any non-local electroweak baryogenesis calculation. Previous studies have found that the predicted baryon asymmetry typically peaks around v w ∼ 0.01, falling off as ∼ 1/v w or faster for larger values [26][27][28][29] (for velocities much smaller than 0.01, the sphalerons can begin to equilibrate the asymmetry). In some cases the dependence of the predicted asymmetry on the wall velocity can be less severe [30][31][32], as the scaling hinges on the form of the dominant CP -violating source (see e.g. Refs. [33][34][35]). Nevertheless, a determination of the wall velocity is often important even for rough estimates of the baryon asymmetry when the wall is expected to move quickly.
The electroweak phase transition in the Standard Model is not first-order [36,37]. Thus, electroweak baryogenesis necessarily requires some physics beyond the SM. There have been many such scenarios proposed in the literature. One of the most popular and straightforward involves augmenting the Standard Model Higgs sector by a gauge singlet scalar field [38][39][40][41][42][43][44][45][46][47][48][49][50][51], providing a tree-level cubic term in the effective potential. Such a cubic term can easily give rise to the barrier required for a first-order transition, and the singlet nature of the new state(s) in many cases can ensure compatibility with current phenomenological constraints [47][48][49][50][51], including the observation of the 125 GeV Standard Model-like Higgs at the LHC [52,53]. This class of models, though popular and simple, is often expected to produce fast-moving bubble walls when the singlet field VEV changes appreciably during the electroweak phase transition [44]. This is simply because the additional field direction contributes to the pressure difference between the phases, which drives the expansion of the bubble, but does not experience a substantial drag force from the plasma 3 . Although this fact was recognized several years ago [44], studies of EWB typically focus on the strength of the phase transition and assume CP -violation can be added in separately without considering the effects of the bubble wall dynamics on generating the baryon asymmetry. In fact, while many studies have since considered baryogenesis in these scenarios with a changing singlet VEV, the wall velocity has never been directly calculated 4 , nor has it been shown that the resulting bubble walls can propagate subsonically as required for successful EWB. Our aim here is to fill this gap.
The bubble wall velocity is an important quantity to compute apart from baryogenesis considerations. For example, models with additional gauge singlets that predict a strong first-order transition can source gravitational waves through bubble collisions and turbulence (see e.g. Refs. [1,4,[55][56][57][58][59][60][61][62][63][64][65][66][67][68]). These scenarios may be effectively probed by upcoming gravitational wave experiments, such as eLISA [69], or Big Bang Observer [70]. However, in order to connect these observations to an underlying theory, one must be able to reliably calculate the velocity of the expanding bubble, as well as the other properties of the phase transition.
A detailed calculation of the wall velocity is rather involved. Building on previous work [26,[71][72][73][74][75], Moore and Prokopec were the first to calculate the velocity for the Standard Model case microphysically in Refs. [76,77] (Ref. [54] also recently revisited this calculation). Five years later, John and Schmidt performed an analogous study in the minimal supersymmetric Standard Model (MSSM) with a light scalar top quark (stop) [78]. Around the same time, the effects of infrared gauge bosons on the wall velocity were calculated in Ref. [79]. Most recently, Ref. [54] extended the results of Moore and Prokopec to other SM-like scenarios, including that in which the VEV of a singlet scalar is approximately stabilized during the EWPT. To date, these remain the only full microphysical calculations of the wall velocity existing in the literature 5 . Recent years have seen progress in matching models onto these existing results [54,[80][81][82] and in the hydrodynamic considerations associated with bubble wall expansion (see e.g. [83][84][85][86][87][88][89]).
Scenarios in which an the VEV of an additional singlet scalar field changes appreciably during the transition merit separate consideration 6 . This is because one must account for the friction on the singlet field. Neglecting these contributions can lead to a drastic overestimate of the wall velocity. A proper treatment requires computing several new classes of interaction rates in the plasma, which can be rather involved. Also, the additional field direction complicates the equations of motion for the condensates. Nevertheless, the calculation can be done and is especially important given the current status of electroweak baryogenesis in light of collider searches. Until recently, the MSSM light stop scenario [90,91] was considered by many to be the most plausible setting for electroweak baryogenesis beyond the SM. Now light stops are in severe tension with both direct LHC searches [92,93] and measurements of the Standard Model-like Higgs couplings [94][95][96][97]. Similar conclusions hold true for many different models relying on large thermally-induced cubic terms to strengthen the phase transition [94,98,99]. This situation has led to a renewed interest in singlet-driven scenarios, since they can be much more elusive at colliders [50]. An analysis of the wall dynamics would mark an important step forward in understanding electroweak baryogenesis in these models. 5 By 'microphysical', we mean calculations explicitly computing the friction exerted by the plasma on the wall, as opposed to those using a phenomenological viscosity parameter. The former involves determining the various interaction rates in the electroweak plasma and solving for the deviations from thermal equilibrium around the bubble wall, as we describe in detail below. 6 In the remainder of this study we will take 'singlet-driven scenarios' to refer to those in which the singlet VEV changes non-negligibly during the transition. This can occur even in models with a discrete Z2 symmetry at T = 0. Examples in which the singlet VEV is approximately stabilized during the EWPT can be treated by the techniques developed for the Standard Model (or MSSM, if the singlet contributions to the finite-T cubic term are large) since only the Higgs field is involved in the transition (see e.g. Refs. [54,82] for an application of this approach).
The goal of this study will be to demonstrate how the electroweak bubble wall velocity can be calculated and to extract some general features of the wall dynamics in singlet-driven scenarios . We will generally follow the strategy and techniques developed in Refs. [54,[76][77][78][79] but modified to account for the singlet field direction. The methods described will be applicable to many different models, although for the sake of simplicity we will frame our discussion in the real singlet extension of the SM, sometimes known as the 'xSM' [43].This scenario should encapsulate the most relevant features of models with singlet-driven firstorder phase transitions. Notably, the xSM does not feature any new sources of CP -violation, which could in principle significantly alter the wall dynamics when included. We comment further on this below and the reader should keep this in mind as we proceed.
We will focus on two different schemes for calculating the wall velocity. The first is explicitly gauge-independent and neglects the contributions to the effective potential and friction from the SU (2) L gauge bosons, while the second includes them. For slower bubble walls (such that the friction on the singlet field is large), both calculations yield similar results, while for faster walls the gauge boson contributions become increasingly important in slowing down the expansion.
The remainder of this study is structured as follows. In Section 2, we introduce the singlet-driven scenario and the finite temperature effective potential which will be used throughout our study. We then move on to computing the wall velocity. Following Ref. [77], the calculation can be broken down into several parts. First, the phase transition properties must be computed and the temperature near the bubble inferred from hydrodynamic considerations, as discussed in Sec. 3. There we also discuss the equations of motion (EOMs) for the bubble wall and consider the simple case of wall velocities in the ultra-relativistic limit. Possible values for the steady state wall velocity are those such that the equations of motion for the wall are satisfied. The EOMs depend on the deviations from thermal equilibrium of all plasma excitations in front of the wall. These are discussed in Sec. 4. We next move on to solving the system of equations for the deviations from equilibrium and the equations of motion in Sec. 5. The calculation is then applied to the parameter space of the xSM consistent with all phenomenological constraints in Sec. 6. We find that sufficiently strong phase transitions may possess no subsonic solutions, and that v w 0.2 for points with h /T n ≥ 1 in the parameter regions considered. We infer that bubbles may expand slowly enough for singlet-driven electroweak baryogenesis, but only in certain portions of the parameter space. Results for the bubble wall profiles are also presented. Finally, our main findings and conclusions are summarized in Sec. 7. We also provide a brief appendix which compares the interaction rates we have calculated with those appearing elsewhere in the literature.
Before proceeding it is important to note that there are other electroweak baryogenesis scenarios that do not rely on diffusion in front of the wall, and hence that do not require slow bubble walls. Local EWB [100][101][102][103], in which the baryon number and CP -violation occur in the same region at the bubble wall boundary, is one such example 7 . However, this typically leads to a highly suppressed total baryon asymmetry relative to the non-local case, 7 Cold Electroweak Baryogenesis [104][105][106] also falls into this category. since the sphaleron transitions turn off near the outer edge of the bubble wall [19]. More recently, an interesting scenario was presented in Ref. [107], in which bubbles can expand quickly enough to significantly reheat the plasma inside the bubble. Secondary bubbles can then nucleate near which transport-driven baryogenesis can occur. This is an intriguing possibility, however it requires a substantial amount of reheating which only occurs for very strong phase transitions. The reader should bear these alternative scenarios in mind as we proceed.

A Singlet-Extended Higgs Sector
To study the dynamics of singlet-driven electroweak phase transitions, we will work in the real singlet extension of the Standard Model. This simple scenario has been studied in depth in the literature, from the standpoint of electroweak baryogenesis, dark matter, LHC signatures, and more (see e.g. Refs. [42,43,45,49,[108][109][110][111][112][113] and references therein). The arguments and methods we discuss here can be straightforwardly applied to other singlet extensions of the SM, such as the next-to-minimal supersymmetric Standard Model (nMSSM [31,32,41] or NMSSM [32, 38-40, 47, 48]) and other scenarios with real or complex gauge singlets [51].
The tree-level potential is taken to be 8 where S is a real scalar singlet under the Standard Model gauge groups and H is a complex SU (2) L doublet. Both the singlet and CP -even neutral component of H are assumed to obtain vacuum expectation values during electroweak symmetry breaking. Throughout our discussion, we will also assume that both VEVs vanish in the high-temperature phase for simplicity (this is discussed further below). Non-vanishing VEVs correspond to minima of the effective potential for non-zero background field values φ h , φ s . These classical fields are those relevant for computing the properties of the phase transition. At a given temperature, we can expand H and S about the background fields, GeV. The zero temperature singlet VEV, φ s (T = 0) ≡ v s can vary. Throughout our study we will identify h with the Standard Model-like Higgs discovered at the LHC [52,53] and take s to be a pure singlet with no mixing at tree-level. The phenomenology of this setup and our choices for the various parameters are detailed in Section 6.1 below.

The Effective Potential
For a homogeneous background field configuration φ(x) ≡ φ, the ground state of the theory corresponds to a minimum of the effective potential V eff (φ). At one loop, V eff is given by the tree-level potential (expanded around the background fields), modified by additional Coleman-Weinberg terms.
At finite temperature and density the physical ground state of the theory is altered by the interactions of the scalar field φ with the ambient plasma. The vacua of the theory can then be determined from the finite-temperature effective potential, V eff (φ, T ). In the simple case involving one background field, it is given by [14] where the plus and minus signs correspond to the bosonic and fermionic contributions, respectively, and the N i are the associated number of degrees of freedom for the species i. This expression generalizes straightforwardly to the case of more than one background field. The functions J ± are given by In the high-temperature limit they admit a useful expansion, given by with a b = 16π 2 e 3/2−2γ E , a f = π 2 e 3/2−2γ E , and γ E the Euler-Mascheroni constant. Note that the thermal contributions above correspond to momentum integrals of equilibrium distribution functions for all species in the plasma coupled to φ [14,44].

Gauge-Invariance
The finite temperature effective potential is only gauge-invariant at its extrema [114,115]. Thus, tunneling calculations depending on the potential away from the local minima are in general gauge-dependent. This will result in a gauge-dependent determination of the nucleation temperature for the phase transition, T n , and ultimately the wall velocity. To avoid this as much as possible, our primary analysis will only consider terms in the effective potential which are explicitly gauge-invariant. Thus, we will not include the T = 0 Coleman-Weinberg corrections, or the finite temperature cubic and tadpole terms in the high-temperature effective potential (gauge-dependence in the tadpole may enter at higher perturbative order [49]). This is precisely the strategy followed by Ref. [49] in analyzing the phase transition properties of the xSM. The finite-temperature effective potential in this case becomes (2.6) Although morally satisfying, the gauge-invariant approach has the disadvantage of sometimes neglecting numerically important contributions to the effective potential. While it should capture the physics we are interested in, namely the singlet contributions to the potential and the friction on the expanding bubble, it is important to consider the effects of the neglected terms, especially the gauge boson cubic term 9 . To this end, we will also show results for the wall velocities including the gauge boson cubic term and friction in Landau gauge. In this case the effective potential in Eq. 2.6 is modified by the additional term We will find that these contributions (and the friction from gauge bosons) are numerically significant in some cases, especially for faster moving bubble walls.

Preliminaries: Phase Transitions, Hydrodynamics, and the Wall Equations of Motion
For a given point in the model parameter space with a first-order electroweak phase transition, we are interested in determining the steady-state velocity of the bubble wall separating the electroweak-symmetric and broken phases. To do so we must first determine the properties of the phase transition (most importantly its characteristic temperature), as well as the equations of motion for the scalar field condensates. Our discussion will roughly follow that of Ref. [77], which we draw from frequently throughout the remainder of this study.

Bubble Nucleation
Using the effective potential in Eq. 2.6, first-order transitions can occur when two local minima coexist for some range of temperatures. The background fields can then "tunnel" from the origin to the new vacuum, in which φ h , φ s = 0. This can begin to occur below the critical temperature, T c , at which the two relevant vacua are degenerate. Bubbles begin to nucleate efficiently at the nuceation temperature, T n , determined by requiring the expectation value for one bubble to nucleate per Hubble volume to be ∼ O(1). At finite temperature, the nucleation probability is determined by the O(3)-symmetric instanton interpolating between the metastable and and true vacua with the lowest ratio of threedimensional Euclidean action 10 to temperature, S 3 /T [12][13][14]. The tunneling probability 9 The inclusion of the tadpole term acts primarily to shift the high-temperature minimum away from φs = 0. We have computed the wall velocities for several scenarios with this term included and obtain values similar to those found neglecting the tadpole. For simplicity we will not consider this term further in this study, although our methods can be straightforwardly modified to include it. 10 For all the cases we consider, the nucleation temperatures are much larger than the inverse radii of the instantons, and so the O(3) bounce is indeed the relevant quantity to consider.
per unit volume is then given by where the pre-factor A(T ) is only weakly temperature-dependent. Using dimensional analysis to estimate A(T ) and assuming typical electroweak temperatures, one finds that the nucleation temperature is approximately determined by S 3 (T n )/T n ≈ 140 [14]. We adopt this definition in the rest of our study. Inside the nucleated bubble, both the VEV of the Higgs and singlet fields will be nonzero by assumption. From the standpoint of electroweak baryogenesis, only the value of φ h is important for sphaleron suppression. We will therefore define a strongly first-order phase transition, occurring at the nucleation temperature T n , by [116] φ h (T n ) T n 1.
This condition is free from explicit gauge-dependence in our primary setup, since we have neglected all gauge-dependent terms in the effective potential. It should be noted that the baryon number preservation condition above still contains several implicit assumptions, discussed in detail in Ref. [15].

Temperature Variations
The nucleation temperature defined above is that of the ambient plasma when bubbles begin to form efficiently. Once formed, however, the temperature is no longer homogeneous. The phase transition releases latent heat into the plasma and the expansion of a subsonic bubble heats up the medium in front of it. The temperature in the broken phase will thus differ from that immediately outside the bubble, which in turn is not the same as the typical nucleation temperature of the bubble. To relate these various quantities requires a treatment of the plasma hydrodynamics. These changes in temperature can have large effects on the expansion of the bubble [77], and so we must take them into account. Far away from the bubble, the relevant temperature is that at which bubble nucleation occurs, T n . We wish to obtain the temperature in the vicinity of bubble wall. To do so, let us consider the wall-plasma system, with the plasma modeled as a perfect relativistic fluid. Hydrodynamic equations can be obtained by requiring conservation of the wall-fluid stress-energy tensor [71], We define the 'fluid -' or 'plasma frame' such that the fluid is at rest far from the bubble and in its center. This is the frame which we use to define the wall velocity v w and the wall profile parameters. Solutions to the fluid equations in the plasma frame can typically be classified as either 'detonations', in which the bubble velocity exceeds the sound speed c s in the plasma, or 'deflagrations', in which v w < c s 11 . Successful subsonic electroweak baryogenesis typically requires a deflagration solution, since otherwise diffusion in front of the bubble is inefficient. We will restrict ourselves to this case. Consider an expanding bubble with free energy V eff (φ − , T − ) inside and V eff (φ + , T + ) immediately outside ('±' subscripts will correspond to quantities outside/inside the bubble). The equations of state (EoS) for the two phases can be written as where p ± and ρ ± are the pressure and energy density of the fluid in either phase, and The above form for the equations, taken from Ref. [83], are inspired by the so-called 'Bag EoS', but involves the temperature-dependent quantities a ± (T ), ± (T ). Fortunately, we can safely neglect the temperature dependence in a ± , ± , using their values at T = T n . This is because the free energy (and hence a ± , ± ) are dominated by light degrees of freedom, which contribute a constant term to the free energy in each phase that does not vary significantly between T n and T c for the cases we consider. We find that using a ± (T n ), ± (T n ) in Eq. 3.5 reproduces the full result for the pressure and energy density within a few percent. This is fortunate as it allows us to avoid several issues arising for more complicated temperature dependence in the EoS, such as the variation of the sound speed in the plasma [87]. Ultimately we will take T eq (x) ≡ T + +δT bg (x) to be the (space-time-dependent) temperature entering the equilibrium distribution functions for the various particles in the plasma. We thus need to determine T + , the temperature just outside the bubble. The pressure and energy density can be related to the fluid velocities on either side of the phase boundary by integrating Eq. 3.3 across the wall. This yields expressions for the fluid velocities v ± in the wall frame which depend on T ± : These velocities can be simply transformed to their analogs in the fluid frame Note that for a subsonic deflagration, v − = 0 and so v w = −v − , as discussed in e.g. Refs. [71,83].
One then needs to relate the temperature T + to T n . In the subsonic deflagration case, the bubble wall is preceded by a shock front moving with velocity v sh in the fluid frame. The temperatures T 1,2 and fluid velocities v 1,2 on either side of the shock front (in its rest frame) will be different; the equation of state is however the same (again neglecting small temperature variations in a + (T ), + (T )). We will use the subscripts 1,2 to denote quantities inside and outside the shock front, respectively. One can again integrate across the interface and use the fact that the fluid is at rest beyond the shock front (i.e. v 2 = 0 → v 2 = −v sh ) with temperature T 2 = T n . This yields an expression for v 1 in terms of T 1 and T n : Again, the corresponding fluid velocity in the fluid frame v 1 is simply given by velocity addition, v 1 = 3v 2 1 −1 2v 1 . Throughout our calculation we neglect the curvature of the bubble wall. In this approximation, the temperatures and fluid velocities (in the fluid frame) between the bubble wall and the shock wave are simply constant [84], and so one can set and solve for T + in terms of T n , v w . Previous studies suggest that using the planar approximation instead of the full solutions to the spherical hydrodynamic equations can reproduce the full result for the wall velocity to within a few percent [82].
With the temperature T + and the static properties of the phase transition determined in this way, we can now consider the asymptotic behavior of the bubble after its formation.

Wall Equations of Motion
The main object for our analysis will be the bubble wall equations of motion corresponding to the set of scalar fields φ i = φ h , φ s . These can be derived by requiring conservation of the energy-momentum tensor for the scalar field condensates computed in a WKB approximation [77], or directly from the Kadanoff-Baym equations [54]. We are interested in the stationary limit of the equations of motion in the plasma frame; that is, we want to investigate the bubble wall once it has reached its terminal velocity (if it exists), with the pressure driving the expansion precisely counterbalanced by the drag force exerted on the bubble by the plasma. This is illustrated in Fig. 1.
Neglecting the curvature of the bubble, in the rest frame of a stationary (non-accelerating) bubble wall all functions will be depend only on z, the distance from the phase boundary. Consequently, in the plasma frame, all functions depend only on the coordinate x ≡ z +v w t, where v w is the wall velocity in the plasma frame and we have assumed that the wall is moving to the left. In the stationary wall limit, the equations of motion then simplify to where primes indicate differentiation with respect to x. Here the sum is over all fields coupling to the scalar field φ i , E j is the (space-time-dependent) energy of the particle j, , and δf j is the deviation from the equilibrium distribution function for the species j.
Solutions to the above equations of motion typically only exist for one subsonic value of the constant v w . This is the quantity we wish to determine. To do so, one must find profiles φ i (x) such that Eq. 3.9 is satisfied, which in turn requires solving for the deviations from equilibrium of the various species in the plasma. These deviations, along with the equilibrium contributions, are responsible for the drag force on the bubble wall. Unfortunately, the δf j depend non-trivially on v w and the bubble profile, so Eq. 3.9 represents a set of integro-differential equations.

Aside: Runaway Bubbles and Tree-level Cubic Terms
Before moving on to the case of non-relativistic bubbles (relevant for electroweak baryogenesis), we can begin by considering the wall dynamics in a simple limit: that of ultrarelativistic, "runaway" bubbles [44], with Lorentz factor γ 1. In this case, the friction on the bubble from the plasma in the large-γ limit is too small to counterbalance the pressure difference between the vacua, which drives the expansion. Ref. [44] showed that this situation is common in singlet-driven transitions, so it is important to review this case before moving on to the non-relativistic regime.
Following Ref. [44], a runaway solution to the equations of motion exists provided (3.10) at the nucleation temperature. Here, f 0 is the equilibrium distribution function of the species i, and φ ± are the field values at the minima of the potential. In the high-T limit, there is a simple interpretation of this criterion in terms of the high-temperature expansion of the thermal effective potential: a runaway solution will exist if it is energetically favorable to tunnel to the broken phase in the 'mean-field' potential, obtained by retaining only the T 2 terms in Eq. 2.6. In other words, (3.11) The above expression indicates that all points found with a first-order phase transition in our gauge-invariant approach (retaining only the quadratic finite-T terms) would feature an ultra-relativistic wall solution if there were no other contributions to the effective potential. This may appear incompatible with our goal of determining subsonic solutions to the equations of motion but it is not. First of all, including the finite temperature cubic term inevitably changes the the transition temperature and the effective potential at that temperature. This can cause the same parameter space point to instead feature V no cubic eff (φ + , T n ) < V no cubic eff (φ − , T n ) , and hence no runaway solution. We indeed find this to be the case for most points considered when including the gauge boson cubic term in our parameter scans. Even if a runaway solution exists for the EOMs including the full finite-T effective potential, there is another important caveat. The criterion in Eq. 3.11 assumes that the bubble is in the ultra-relativistic regime to begin with. However it is instead possible for the friction to prevent the bubble from ever reaching such large velocities required for Eq. 3.10 to be valid. In fact, hydrodynamic effects alone can obstruct the wall from expanding ultra-relativistically [117]. Thus, even if a particular parameter space point admits a runaway solution, it may not be realized if a subsonic stationary solution exists. On the other hand, even if no runaway solution exists, one with v w > c s might. The reader should thus bear in mind that our approach will find subsonic solutions to the equations of motion, not guarantee that they are realized. This is also true of previous studies [54,[76][77][78].
Equation 3.10 shows that the friction force acting on the wall takes a very simple form in the γ 1 limit. This is not the case for the subsonic walls we are interested in. Determining the wall velocity in the γ ∼ 1 regime requires a careful calculation of the various deviations from equilibrium in the plasma. This is what we discuss in the following section.

Setup
With the temperature T ≡ T + inferred from hydronamic considerations, the first step towards solving the bubble wall equations of motion in the non-relativistic (γ ≈ 1) case is determining the distribution functions f i for the various excitations appearing in Eq. 3.9. To do so, we will primarily utilize a perturbative effective kinetic theory approach [118,119], as in previous studies [76][77][78] (we will take a somewhat different approach for the corresponding gauge boson friction, which should be modeled classically as discussed below). This treatment applies to weakly coupled excitations with local interactions and short wavelengths compared to the length scale of the bubble wall in the plasma frame, i.e.
where L w is the wall width. Typical momenta are of order p ∼ T , but softer excitations will be present in the plasma as well. We will assume that the kinetic theory description is viable in the range p gT , which is reasonable for the particles we will be interested in given the values we find for the wall widths. Here and throughout this section g represents a generic dimensionless coupling of the theory 12 that is assumed to be small. Infrared (IR) excitations with momenta p T will not be captured by this treatment, since their interactions cannot be properly described by a local collision term. These contributions can be important for the bosonic species [79], but the perturbative effective kinetic theory should provide an adequate estimate of the damping force on the bubble wall, provided that very infrared excitations are equilibrated quickly [77], as we will assume for most of the species we are interested in 13 .
In the effective kinetic theory we consider, the quasiparticle distribution function for the species i satisfies the Boltzmann equation in the fluid frame, where C[f ] i is a local collision integral. The collision term involves all interactions of the species i with all other excitations in the plasma. It can be written as where the sum is over all 4-body processes ij → mn, with the momenta labeled as p, p , k , and k moving clockwise around the diagram starting with particle i. The matrix elements include finite-temperature effects (discussed below) and are summed over helicities and colors of all four external quasiparticles, then divided by the number of degrees of freedom corresponding to species i, N i (N h = 1, N t = Nt = 6) 14 . The population factor is with the upper (lower) signs corresponding to bosons (fermions) and f a the appropriate Bose-Einstein or Fermi-Dirac distribution function for particle a, which we assume to take the form In Eq. 4.3, the prefactor of 1/2 takes care of both the symmetry factor when identical particles are present in the final state, and the double counting that occurs from the unrestricted sum over m and n. The Boltzmann equations above apply to all quasiparticles in the plasma satisfying Eq. 4.1 with sufficiently high momentum. However, examining Eq. 3.9, we see that only the distribution functions of field excitations with significant couplings to the relevant scalar fields involved in the phase transition are required. Since these particles have significant couplings to the Higgs and singlet scalar fields, we will refer to them as 'heavy'. Also, δf i = δf i (p, x) has some space-time-dependence, arising in part from the spatial variation of the background fluid temperature and velocity across the bubble wall, as discussed in Sec. 3. The background fluid is in local thermal equilibrium and comprises all 'light' effective degrees of freedom. Note that quasiparticles with large field-independent masses will be irrelevant for our purposes, since their distribution functions feature significant Boltzmann suppression. Also, precisely which fields should be considered 'heavy', 'light', or irrelevant depends on the given model. For the singlet-driven scenarios we are concerned with here, the heavy fields will be the top quarks, gauge, Higgs, and singlet bosons.
To find approximate solutions to the Boltzmann equations for the heavy species and background, we will utilize the 'fluid ansatz' [77], in which case the perturbations are assumed to take the form Here µ j , δT j , δv j are the chemical potential, temperature perturbation, and velocity perturbation of the species j, respectively, in the plasma frame. We have assumed that the fields with small couplings to the scalar condensates φ h,s are in thermal equilibrium at a common space-time-dependent temperature T + + δT bg (x) and velocity v bg (x) with vanishing chemical potential, as in Ref. [77]. The assumption that µ bg ≈ 0 is valid whenever the total background particle destruction rate is larger than that for the heavy particles, as will be the case here (all pure gluon rates are enhanced by the large color factors and Bose statistics). The space-time-dependence in δT bg , v bg arises from the change in masses of the corresponding particles moving from the φ i = 0 phase inside the bubble to the φ i = 0 vacuum outside. Throughout this study, we will work to linear order in the perturbations, which are assumed to be small (µ j /T , δT j /T , δT bg /T , δv j , v bg 1). This should be the case for moderately strong phase transitions, and we verify the validity of this assumption a posteriori. It should be noted that this treatment can be extended to accommodate large fluid velocities in front of the wall [54], although this will not be necessary for any of the transitions we consider. As a result, we set all Lorentz γ factors to 1 throughout our calculation.
With the above definitions, the population factor P is given to linear order in the perturbations by where the '0' subscript indicates the corresponding equilibrium distribution function. Note that the background temperature and velocity perturbations do not enter the collision integrals to linear order.
To determine µ i , δT i , and δv i we follow Refs. [76][77][78] and take three moments of each equation, multiplying by 3 and solve the resulting expressions for the perturbations. For a given heavy species, the relevant three equations are given in the plasma frame by where an ingoing particle of the relevant species has momentum p and where Further details can be found in Ref. [77]. The resulting collision terms for each heavy field i can be written as The background excitations also satisfy a set of Boltzmann equations, which arise from Eq. 4.8 with µ bg ≈ 0. The sum above is over all background species, with c 4 ≡ c 4 the heat capacity of the plasma. As for the heavy quasiparticles, the collision terms can be written as Although δT bg and v bg do not enter the collision integrals, the perturbations corresponding to the heavy excitations do. The convention for evaluating the matrix elements is the same as for the heavy particles, with all background excitations treated as one species. Thus, every heavy particle process involving the background excitations will contribute to Eqs. 4.12. We will calculate all of the contributions relevant for singlet-driven transitions in the next subsection.

Relevant Excitations and Interaction Rates
In the SM and its singlet extensions, the relevant heavy species to consider in Eq. 3.9 above are typically the top quarks, SU (2) L gauge, Higgs, and singlet bosons, with the Higgs and singlet excitations being the dominant source of friction on the singlet field condensate. We will consider two different sets of contributions to the total friction. For the gaugeinvariant calculation, we include only the top quark, Higgs, and singlet contributions. When incorporating the gauge boson cubic term, we will also account for the friction arising from the SU (2) L gauge bosons. We will not include the Goldstone friction contribution, since we drop the corresponding finite-T cubic term from the effective potential. This should be a reasonable approximation as the Higgs excitations will only make up about 20% of the total friction on the wall 15 .
The friction from each species enters the bubble wall EOM (Eq. 3.9) through the derivative of the corresponding mass squared. For the top quarks, the effective mass squared is with the corresponding thermal self-energy correction neglecting the subdominant thermal SU (2) L and U (1) Y contributions. The Higgs and singlet require slightly more care. Throughout the remainder of this study we will neglect all mixing effects between the SM-like Higgs and singlet excitations. As discussed in Sec. 6.1, we will choose the parameters of the T = 0 Lagrangian such that the mixing vanishes in the broken electroweak phase. At finite temperature and across the bubble wall this will no longer be the case. However, in the high temperature limit, the effective neutral scalar mass matrix is diagonal, since off-diagonal thermal corrections are proportional to dimensionful parameters and vanish as T → ∞. For temperatures around the electroweak phase transition, the thermal masses still dominate the mixing matrix, and so this should be a decent approximation. The relevant field-dependent masses, including the leading thermal corrections, are then with the thermal masses 15 We do include the Goldstones, φ 0 , φ ± in the various interaction rates. They only appear as external legs of the diagrams we consider, and the corresponding matrix elements are gauge-independent. They are treated as a background species.
where we have left out the light fermion Yukawa contributions. Since we are neglecting the finite-temperature tadpole contribution to the effective potential we also drop the a 1 terms in the masses above. This is required for consistency, since these terms are precisely those that give rise to the finite-temperature tadpole. Finally, as in Ref. [77], we treat the transverse SU (2) L gauge bosons as a single species W with field-dependent mass squared Transverse excitations do not acquire a thermal mass at leading order in the couplings.
Longitudinal modes obtain an effective thermal (Debye) mass at leading order, corresponding to the inverse screening length of electric potentials in the plasma [120]. This is given by in the Standard Model. Since the gauge boson friction is dominated by very infrared excitations, only the transverse contributions will be relevant.
Our strategies for dealing with each of these types of excitations will differ. As we will see below, the top quark and Higgs interaction rates are typically sizable, and so the collision term plays an important role in the corresponding Boltzmann equations. This is not expected to be the case for singlet quasiparticles at high temperature. Contrary to the tops and Higgs, we will assume that the singlet interactions are slow. In this case, the collision term can be neglected. The corresponding Boltzmann equation decouples from the rest of the system and can be solved exactly. We discuss this further in Secs. 4.2.2 and 5.1. Finally, the gauge boson contributions are dominated by infrared dynamics and require a classical treatment, which has been worked out in Ref. [77] and discussed in Sec. 4.3 below.
Let us first consider the interactions involving the top quark, Higgs, and background excitations.

Top, Higgs, and Background Excitations
Solving the Boltzmann equations for the perturbations µ t,h , δT t,h , δT bg , δv t,h , and v bg requires computing the collision integrals corresponding to all the four-body interactions involving t, h, and the background fields. This task is rather daunting due to the sheer number of allowed processes. However, the dominant interactions will be of O(α 2 s ) for the top quarks, and O(α s α t ), O(α 2 t ) for the Higgs bosons, where α s = g 2 3 /4π, α t = y 2 t /4π. We will therefore focus on these interactions, neglecting, for example, contributions involving a factor of α w , which are numerically small compared to the Yukawa-type contributions for the Higgs bosons 16 .
To estimate the relevant interaction rates, we will work at leading order in all couplings in the high-T , weak coupling limit, neglecting all terms of O(m 2 /T 2 ) (here m should be understood as either a zero-temperature or thermal mass). This is the approximation used in all previous microphysical studies of the wall velocity [76][77][78], as well as in the context of plasma properties in arbitrary high-temperature gauge theories [121,122]. This approximation can begin to break down inside the bubble wall for the top quarks and scalars and could in principle be improved upon in the future. Nevertheless, it should reproduce at least the correct parametric dependence of the full leading order result. In this limit we can neglect the effect of the space-time-varying masses on the interaction rates. Hard excitations dominate the phase space for the relevant top quark and Higgs collision integrals, and can be characterized by massless dispersion relations to leading order in the couplings and in the high-temperature limit.
Although the external quasiparticles can be treated as massless in this approximation, infrared excitations appearing as mediators in tand u-channel diagrams naively result in logarithmic IR divergences in the Higgs and top quark scattering amplitudes. These divergences are cut off by the interactions of the mediator with the plasma. For longwavelength excitations, the corrections comprise so-called 'hard thermal loops' (HTLs), and result in a breakdown of the perturbative expansion. The corrections can be resummed into a thermal self-energy correction to the propagator, valid in the low momentum limit 17 . The self-energy is typically of order gT and so these processes can produce sizable logarithmic enhancements of the corresponding matrix elements, scaling as ∼ g(T ) 4 log 1/g(T ) at high temperatures. This provides a useful way of categorizing the most important diagrams contributing to M i in C[f ] in the high-T limit.
A full leading-order determination of the effective scattering rates in the plasma is possible [122,123], though computationally more involved and beyond the scope of this work. It would be interesting to revisit in the future. We will instead work in a 'leading logarithm' expansion, keeping only contributions of order ∼ g(T ) 4 log 1/g(T ), which are typically the largest. In this approximation, only 4-body rates with tand u-channel diagrams contribute. For further details on this approximation, see Refs. [77,121].
Another subtlety arises in computing scattering rates in the high-T limit involving soft tor u-channel exchange. The thermal self-energies involved in the propagators are generally momentum-dependent [122]. Previous studies of the wall velocity neglected these contributions, simply replacing them with the corresponding Debye (thermal) masses for the corresponding gauge bosons (fermions). However, including the momentum-dependent self-energies, which enter at leading order in the couplings, can have a significant effect and has been shown to often provide better agreement between the leading log and full leading order results for plasma transport coefficients in high-temperature gauge theories [122]. Consequently, we will use the full momentum-dependent HTL-resummed propagators [122,124,125] in computing the various collision integrals for the top quark and Higgs excitations.
The relevant processes and their associated vacuum matrix elements are listed in Table 1. All terms of O(m 2 /T 2 ) have been dropped. The leading log matrix elements are summed over the helicities and colors of the external particle (but not particle-antiparticle). : tt ↔ hg, φ 0 g:  Table 1. Relevant 4-body processes and their corresponding matrix elements in the leading log approximation. The matrix elements are summed over the helicities and colors of all four external states (as well as flavors and quark -anti-quark for tq → tq). The excitation appearing on the internal propagators in each case is listed in the right-hand column. Note that other tand u-channel processes exist, but do not contribute logarithmically to the collision integrals or are suppressed by powers of couplings small compared to g 3 , y t .
These contributions will also be divided by the number of degrees of freedom of the species under consideration when entering the various Γ i k,j . The vacuum matrix elements must be modified to include the medium-dependent effects discussed above. At leading order, this amounts to inserting the momentum-dependent HTL self-energies on the internal lines. To translate the vacuum matrix elements above to their finite-temperature analogs, we can use the results of Refs. [119,122]. For diagrams with an exchanged fermion in the leading log approximation, this amounts to the replacement u t − s t → 4 Re(p · q k · q * + s q · q * ) | q · q| 2 (4.20) with q µ = p µ − p µ − Σ µ (p − p ) and Σ µ (q) the fermionic HTL self-energy function [122,124,125] with m f (T ) the leading order fermion thermal mass, given approximately by m f (T ) ≈ 1/6g 3 T for quarks (see Eq. 4.16).
For the gluon exchange diagrams, we must replace with the retarded thermal equilibrium gluon propagator D µν (q) given by with m D (T ) the Debye mass of the gauge boson (= √ 2g 3 T for the gluon).
With the above replacements, we can now perform the collision integrals numerically. We do so using the phase space parametrization discussed and detailed in Refs. [119,122,126] and the Vegas Monte Carlo routine included in the Cuba package [127]. The integrals converge reasonably quickly for most cases on a standard desktop computer to reasonable numerical precision.
The results of this numerical evaluation are given in Equations 4.25-4.29 below. These values can then be plugged into Eqs. 4.8 and 4.11 for the perturbations. For the Higgs bosons, the rates (in the leading log approximation) are  while the corresponding contributions to the top quark distributions are

(4.26)
For the top quarks, while their contributions to the Higgs distributions are  Note that the background rates for the top quarks tend to appear larger than their counterparts in Eq. 4.27 above. This is simply because in the background rates we have summed over all contributions, while the rates in Eqs. 4.25-4.28 are the average values per degree of freedom (e.g. divided by N t = 6 in the top quark case). The latter rates will be multiplied by the appropriate N i factors when they enter the bubble wall equation of motion. Also, note that the contributions in Eqs. 4.26 and 4.28 are negative because they arise from diagrams with the relevant species on the outgoing legs of the Feynman diagrams. In previous work, the above integrals were performed analytically using several approximations and without incorporating the (momentum-dependent) self-energies. The different computational methods used here change the rates by O(1) factors relative to the results in Ref. [77]. There were also some algebraic errors in the results of Ref. [77], as pointed out in Ref. [121], that contribute to the discrepancy. Although our treatment is still formally at the same order as that of Ref. [77], the HTL-improved calculation in many cases is expected to more closely reproduce the leading order result (see e.g. Fig. 1 of Ref. [122]). Nevertheless, modulo the algebraic mistakes in Ref. [77], our interaction rates are no more accurate than those of Moore and Prokopec in approximating the full leading order results; they should simply be thought of as arising from a different set of approximations. We comment further on the differences between the rates found in Ref. [77] and those reported above in Appendix A. The reader should bear in mind that the predicted wall velocity will tend to be higher if the rates computed in Ref. [77] are used instead of ours. This is because the former are larger and thus result in faster equilibration for the various perturbations. Note also that the collision integrals computed and listed above depend only on the Standard Model degrees of freedom, and as such are quite general. They can be used in various applications beyond those considered in this work.
Before moving on, some comments regarding the higher order contributions neglected in the fluid approximation are in order. The assumed form for the perturbations is that of a perfect fluid and can be thought of as a truncated expansion in powers of momentum. That is, the fluid ansatz assumes that the effects of higher angular moments p Y m ( p) in the distribution functions are negligible [77] (the Y m are spherical harmonics). For the top quarks this is a good approximation, since we find that the velocity perturbations typically satisfy δvT /δµ 0.1, while the contributions from higher moments should scale roughly like (δvT /δµ) [77]. On the other hand, the Higgs bosons have smaller interaction rates than the tops, and so the fluid approximation begins to break down for strong phase transitions. For this reason we will restrict ourselves to moderately strong phase transitions with φ h (T n )/T n 1.1 in our consideration of the xSM in section 6. Already in this regime some points will be found to possess no subsonic solutions. Further details on the limitations of the fluid approximation can be found in Appendix B of Ref. [77].

Singlet Contributions
Excitations of the singlet field will also contribute to the friction on the bubble wall. The corresponding collision integral for singlet quasiparticles is dominated by scattering processes involving four external scalars. At high temperatures, the resulting effective interaction rates are typically suppressed relative to those for the processes involving external fermions. To see this, note that processes with t-channel diagrams involving two external scalars and two external fermions schematically contribute to the first moment of the Boltzmann equation in the small-q limit (here q ≡ |p − p |). The logarithmic divergence is cut off by the thermal self-energy of the exchanged quasiparticle. The upper limit q ∼ T corresponds to the breakdown of the small q approximation. For processes involving four external scalars with a scalar exchanged in the t-channel, the integrals of the Bose-Einstein distribution functions are also IR sensitive. Cutting off the distribution functions with a parameter with mass dimension 1 such that f 0 (p/T ) → f 0 (p/T + /T ), we find the corresponding contribution to be in the small q and /T regime (a is a cubic coupling with mass dimension 1). The divergence is now nominally quadratic (again cut off by the self-energy of the exchanged scalar), and the integral of over the distribution functions is cut off by the thermal masses of the external scalars in the infrared. This suggests that a rough leading order estimate of the scalar quasiparticle scattering rates should include thermal masses in the Bose-Einstein distribution functions. Computing the dominant contributions involving the cubic and quartic couplings and performing the resulting integrals, we find interaction rates that are significantly smaller than those for the tops and Higgs across the range of couplings and temperatures we consider, despite the nominally more severe divergence structure. This is because for large temperatures, ∼ gT and so the schematic rate in Eq. 4.31 is suppressed by 1/T 3 . This is expected, since at high temperatures all dimensionful parameters of the zero temperature theory should be irrelevant [123]. Meanwhile the quartic interactions do not contribute a small q divergence at leading order. We therefore expect the singlet qusiparticle collision term to be small. Assuming this is the case, the fluid approximation is likely to provide a rather poor estimate of the corresponding friction. Instead, we will make a 'free particle' approximation [72,73] for these excitations, dropping the collision term for the singlet. In this case, the Boltzmann equation can be solved exactly, without taking moments. The solution is given by Eq. 5.3 of Ref. [77], to lowest order in v w , and is reproduced below in Eq. 5.1. We will include the corresponding contribution to the equations of motion when computing the friction. Note that the presence of non-negligible interactions would decrease the friction and increase the predicted wall velocity.
This treatment assumes that the dominant singlet excitations are well described by our perturbative kinetic theory, with a local collision term. This should be true for the hard excitations. For much softer excitations, a classical treatment with the short wavelength fluctuations integrated out is likely more appropriate [79,128]. Unlike classical Yang-Mills fields (discussed in the next subsection), the classical scalar field is not overdamped [129], and so infrared excitations are likely to equilibrate quickly. We thus neglect the effect of infrared singlet modes on the friction. This approximation is likely rather crude and should be revisited in the future. Including the IR contributions would increase the friction and potentially slow the wall down. The reader should bear this in mind as we proceed.

Gauge Boson Contributions
Finally, for our calculations incorporating the finite-temperature gauge boson cubic term in Eq. 2.7, we will include the friction from the SU (2) L gauge bosons. In contrast with the top quarks, Higgs, and singlet excitations, the friction in this case is dominated by infrared degrees of freedom, which can be treated approximately as over-damped classical fields [129] as opposed to the perturbative approach utilized for the other species. The distribution functions in the classical limit can be shown to satisfy [129] where N is a noise term. This equation serves as the analog of Eq. 4.2. Further discussion of its derivation and applicability can be found in Ref. [129]. Note that hard gauge boson excitations also exert a drag force on the bubble wall, as computed in Ref. [77]. However, we have verified that these contributions are substantially suppressed relative to that from the IR gauge bosons, as found in Ref. [129]. We do not include them.

Solving for the Wall Velocity
With the collision terms evaluated, we can now solve the Boltzmann equations to determine the perturbations for a given field profile and wall velocity. The goal is then to find the value of v w and the profile (and hence the perturbations) such that the equations of motion are satisfied. We will describe how this can be done below. First, let us consider solutions to the Boltzmann equations given a particular profile and value of v w .

Exact Solution for the Singlet Excitations
As mentioned above, the singlet equation can be decoupled from the rest of the system and solved exactly in the free-particle limit. The result been discussed in detail previously [72,73,77], and so we simply quote it here. The integral appearing in the equations of motion 3.9, to lowest order in v w , is where the upper (lower) sign is for fermions (bosons). The function Q is defined as 2) with m 0 s (T ) the singlet mass (including thermal contributions) in the broken phase (i.e. at z → ∞). This integral is multiplied by ∂m 2 s (φ h , φ s , T )/∂φ h in the Higgs EOM, and by ∂m 2 s (φ h , φ s , T )/∂φ s in the singlet equation.

Exact Solution for the IR Gauge Contributions
We can also solve for the classical gauge boson contribution in Eq. 4.32. The result is [129] where the quantity x * solves m W [φ h (x * )] = 1/L h , with L h the SM-like Higgs wall width. For smaller x, the WKB description used to derive Eq. 4.32 breaks down. For more discussion on this point, see Ref. [129]. Note that this value cuts off the IR divergence of Eq. 5.3.

Solving the top-Higgs System
It remains to solve the equations for the top quark, Higgs, and background excitations.
Here we follow the methods found in Refs. [77,78] with some slight modifications.
Since we are interested in static solutions to the equations of motion in the wall frame, all quantities depend only on x and so the derivatives in Eqs. 4.8 can be re-written as ∂ t → v w d/dx, ∂ z → d/dx. The Boltzmann equations in the static limit are therefore a set of linear ordinary inhomogeneous differential equations.
To solve them, the equations for the background temperature and velocity, Eqs. 4.11, can be used to eliminate T bg and v bg from the top quark and Higgs equations. Defining a vector of perturbations Eqs. 4.8 can then be written as and the source vector (5.8) The field-dependent masses are given by Eq. 4.15.
The system of equations can be solved by simple Green's function techniques. Following Ref. [77] we define the matrix χ such that where λ k are the eigenvalues of A −1 Γ. It is then straightforward to define the vector Green's function in terms of which the perturbation δ i is given by These solutions can then be inserted into the equation for the variation in the background temperature, defined with respect to T + , its value far ahead of the bubble in the shock front.

Approximate Solutions to the Equations of Motion
With the perturbations determined, we can now try to identify solutions to the wall equations of motion. In terms of the perturbations δ j , Eq. 3.9 reads, for the gauge-invariant case, where the last term is given in Eq. 5.1. If the gauge boson contributions are included, the RHS of Eq. 5.3 should be added to the LHS of the above expression. The boundary conditions are φ h,s (x → ∓∞) = φ h,s;± (T + ) and φ h,s (x → ±∞) = 0. This system of equations will typically admit a solution for certain values of v w and profile φ h,s (x). Our strategy will be to vary the profile and scan over values of v w consistent with a deflagration bubble, looking for parameters such that the equations of motion (and Boltzmann equations) are simultaneously satisfied. All parameter space points we consider have at most one deflagration solution. Eq. 5.13 represents a set of integro-differential equations for Φ ≡ (φ h , φ s ) T (in the following discussion it will be useful to use explicit vector notation). To find its approximate solutions, we will follow the strategy of Refs. [77,78] and use an ansatz for the field profiles which depend on only a few parameters. Of course in using an ansatz it is unlikely that the full equations of motion will be satisfied exactly. However, we can reasonably approximate a solution by scanning over the ansatz parameters and imposing physical constraints. For a given choice of parameters, the Boltzmann equations can be solved exactly and the results inserted into the EOM. A set of parameter values such that all constraints are simultaneously satisfied corresponds to an approximate solution to the original equations. This strategy has been employed in previous calculations of the wall velocity [54,76,77] and we expect the results obtained in this way to be a decent approximation to the full numerical solution.
Before analyzing Eq. 5.13 further, some useful insight can be gained from solving the corresponding field equations with the assumption of constant friction of the form subject to the boundary conditions Φ i (x → ∓∞) = Φ i,± , Φ i (x → ±∞) = 0, where ∇ φ ≡ (∂/∂φ h , ∂/∂φ s ) T and F is the same for both field directions. Clearly this is not a realistic case, but we will improve on it below. These equations are much simpler than the full integro-differential equations of Eq. 5.13. As detailed in Ref. [48], the solution to the equations of motion can be found numerically via path deformations, and corresponds to a limit in which all the friction on the wall is parallel to the trajectory in the (φ h , φ s ) field space within the wall. This can be seen by noting that dΦ/dx is a "velocity" vector in field space, so FdΦ/dx always acts parallel to the trajectory. In all cases we consider, the solution to these equations of motion is well-fit by a tanh ansatz with parameters in the fluid frame. Here L i = L h , L s are the wall widths and δ i are offsets allowing for a good fit to the numerical solution. We can take δ h = 0 without loss of generality. Note that if we had allowed φ s = 0 in the electroweak-symmetric phase, we could have instead used φ s (x) = φ 0 s + ∆φ s /2 1 + tanh x−δs Ls for the ansatz, with φ 0 s the singlet VEV in the symmetric phase and ∆φ s the change in VEV across the wall. The remaining analysis would proceed in the same way.
How does the situation change when including a realistic friction term? The friction is no longer proportional to dΦ/dx and so will have some component perpendicular to the field space trajectory, acting as an effective "normal force" along the path. However, there is a fortunate simplification we can make if the friction perpendicular to the field space trajectory found by solving Eq. 5.14 is negligible. In this case, the field will not be significantly deformed from its field space path found using the constant friction equations of motion, although the field profile in physical space will change. In other words, if we write the solution to Eq. 5.14 as Φ(s) where s = s(x) is some parameter such that |dΦ/ds| = 1, the effect of a change in the friction parallel to Φ(s) will only be to alter s(x). Meanwhile, a change in the friction normal to the profile would result in a change of Φ(s) itself. Applying this reasoning to the tanh ansatz (which we find to be a good fit to the constant friction solution), the effect of altering only the friction parallel to the trajectory and neglecting that normal to the path will simply be an overall simultaneous re-scaling of all the widths and offsets: L h,s → aL h,s , δ s → aδ s . (5.16) This is the only change in the tanh profile that will not deform the path in field space. Then, starting from the constant friction solution, the problem can be reduced to finding the values of v w and a such that the pressure and pressure gradient in the wall vanish: (5.17) The above constraints will only be satisfied for values of v w and a such that the wall is not accelerating or expanding/contracting, as required for the steady-state solution we are seeking. This is a simple generalization of the strategy used in the SM case in Refs. [54,76,77], where v w and L w are varied.
In what follows, we will assume that the friction force normal to the field space path determined from Eq. 5.14 is negligible, such that the discussion of the above paragraph applies. The validity of this assumption can be checked a posteriori (which we do), but there is an intuitive reason why it should often be reasonable. At very high temperatures the potential is stabilized at the origin (in our approximation) by the effective thermal masses ∼ gT of all the scalars. Around T ∼ 0, the potential is necessarily stabilized at a minimum away from the origin. The temperature of the phase transition is such that two minima exist simultaneously and are (nearly) degenerate. Provided that the tree-level contribution to the barrier is not too large, this can only occur if there is a significant cancellation between the zero-temperature and finite-temperature corrections in some direction of the Φ field space. This approximate cancellation will be largest along the field space trajectory Φ(s) found by solving Eq. 5.14 such that, schematically, The resulting ridge in the finite-temperature effective potential is precisely that along which the cubic term becomes relevant. We can insert the solution to Eq. 5.14 into Eq. 5.13 and see how we expect the solution to change when going to the full EOM. With the approximate cancellation of Eq. 5.18 in effect, the full equation of motion parallel to Φ(s) is then schematically along Φ(s). The (approximate) cancellation has made the contribution from the friction the leading effect. The change in the friction term from Eq. 5.14 to 5.13 will alter s(x) but leave the field space trajectory Φ(s) unchanged. In contrast, the cancellation in Eq. 5.18 is not expected to hold in the perpendicular direction along the path (∝ d 2 Φ/ds 2 ), unless the minimum of the potential away from the origin lies in the bottom of a shallow valley. In the absence of such an approximate continuous symmetry, the full EOM perpendicular to Φ(s) is and so the effect of the friction in this direction is perturbatively small, resulting in a negligible correction to the perpendicular component of Eq. 5.14 and hence to Φ(s).
The above discussion suggests the following strategy for finding approximate solutions to the equations of motion for a given parameter space point: 1. Compute the phase transition properties, namely the order parameter and T n . We do this using the CosmoTransitions package [130].
2. Solve for the constant friction profile from Eq. 5.14. This can be done using path deformations [48] or otherwise. Fit the solution to the tanh ansatz Eq. 5.15.
3. Solve the hydrodynamic relations to obtain T + for various values of v w .

4.
Vary the values of v w and a. For each pair, solve the Boltzmann equations as discussed above, using T = T + and L i = aL 0 i , δ s = aδ 0 s with L 0 i and δ 0 s obtained from the numerical solution of Eq. 5.14.
5. Insert the solutions for the perturbations (and background temperature profile) into Eq. 5.13, then compute the constraints in Eq. 5.17. The values of v w and a satisfying Eq. 5.17 can be found by interpolating between the results of the scan.
This method generalizes that of Refs. [54,76,77] to accommodate the additional singlet field direction.
The results for v w and a obtained in this way will still produce a residual 'normal force' perpendicular to the trajectory in field space when inserted back into Eq. 5.13 due to the neglect of the friction in the direction ∝ d 2 Φ/ds 2 . Defining s(x) = |Φ(x)|, the tangent and normal unit vectors to the field space path Φ(s) arê The 'normal force' along Φ(s) is given by A full solution to the equations of motion would guarantee that N(x) = 0 for all x. This will not be true for our approximate solutions.
To check that this residual normal force is indeed negligible, we can deform the profile to eliminate it. The deformation can be performed along the lines suggested by Ref. [130] for computing the critical bubble profile. It typically results in small changes to the original field space profile, which in turn have very little effect on the perturbations and constraints in Eq. 5.17, since the curvature perpendicular to the path is typically significant. This suggests that the wall velocity and profile found in the way outlined above should indeed provide a reasonable approximation to those obtained from the full solution of the equations of motion, at least in the cases we consider. There may be exceptions elsewhere in the parameter space.
Note that the procedure outlined in Steps 1-5 above is quite general, and can be adapted beyond singlet models to other scenarios with multiple field directions, provided Eq. 5.18 is approximately satisfied.

Parameter Space and Phenomenology
We now turn to the parameter space of the real singlet extension of the Standard Model as an application. Our goal is not a comprehensive analysis of this model. Instead, we focus on a sample of the parameter space consistent with current observations and a strongly first-order phase transition.
The authors of Ref. [49] recently performed a detailed analysis of the electroweak phase transition in this setup and so we use their study as a guide. Recall that we have identified the excitation h as the Standard Model-like Higgs with m h 125 GeV. The couplings of the discovered Higgs are very close to those expected in the Standard Model [131,132]. Therefore we will assume that the doublet H couples precisely as in the Standard Model and that there is no mixing between the singlet and Higgs at tree-level 18 . This corresponds to the choice a 1 + 2a 2 v s = 0 (6.1) and immediately fixes λ in terms of the observed Higgs mass, m 2 h = 2λv 2 = (125 GeV) 2 . Note that Ref. [49] showed that deviations from this no-mixing limit are allowed by current LHC measurements and so this requirement can be relaxed. In our scans, we will vary both the cross-quartic coupling, a 2 , and the zero-temperature singlet VEV, v s . Then Eq. 6.1 can be used to determine a 1 , given our choice for v s and a 2 .
We will also assume b 3 = 0. The barrier required for a first-order electroweak phase transition will then arise primarily from the tree-level mixed cubic coupling a 1 . Again, this is not required by the phenomenology, and this choice can be altered without significantly affecting any of our arguments. Note that much of the NMSSM parameter space compatible with a strongly first-order electroweak phase transition lies close to the corresponding limit |κA κ | |λA λ | [47,48]. Equating the minima of the tree-level potential with the VEVs v = 246 GeV and v s yields the conditions which allows us to solve for µ 2 and b 2 . We also take the tree-level singlet mass squared (in the zero mixing limit), as an input, and use the above expression to solve for b 4 . The free parameters are thus a 2 , v s , and m 2 s . There are some additional requirements on the theory that allow us to hone in on phenomenologically viable values for these quantities. First of all, stability of the T = 0 potential requires b 4 > 0, which in turn limits the values of m 2 s we can consider via Eq. 6.3. Also, if m 2 s is too small, the decay h → ss would cause large deviations in the width of h which are not observed experimentally. On the other hand, if m s > 2m h , di-Higgs production at the LHC can place significant constraints on the model [133,134]. For simplicity, we avoid this regime and choose m h /2 < m s < 2m h .
Further insight can be gained from considering the expected behavior of the electroweak phase transition strength and wall velocity as a function of these free parameters. In particular, larger values of a 2 are favorable from the standpoint of small wall velocities. This is because increasing a 2 corresponds to lowering T c , as can be seen by noting that T c corresponds approximately to the temperature such that m 2 h (φ h , φ s , T ) ∼ 0. Lowering T c results in lower values of φ c satisfying φ c /T c 1. Smaller field values in turn lower the total pressure difference between the vacua. This pressure difference is what drives the expansion of the bubble and which must be compensated for by the friction. Furthermore, larger values for a 2 correspond to larger m 2 s (φ h , φ s ), which in turn increases the friction from the singlet on the wall.
With the above reasoning in mind, we choose two sets of parameters likely to be promising for electroweak baryogenesis and across which we can compute the bubble wall velocity. These are For both sets of points we vary v s , which corresponds to varying the strength of the electroweak phase transition. This is clear, since higher values of v s correspond to larger |a 1 | via Eq. 6.1 and hence a larger contribution to the barrier separating the electroweak minimum from the origin at finite temperature. We vary v s up to values such that either the fluid approximation breaks down, or subsonic solutions no longer exist. This corresponds roughly to values of v s between 30-100 GeV for Sets 1 and 2. Note that with our choices of parameters, every point in Sets 1 and 2 will satisfy all current phenomenological constraints. Specifically, the electroweak vacuum is absolutely stable for all points considered while the absence of s − h mixing ensures that both h and the new singlet-like state are compatible with current observations and limits.

Results
Finally we arrive at our results for the xSM. The wall velocities computed in a parameter scan for the two sets of points (Set 1 and 2) described above are shown in Fig. 2. The critical bubble profile and nucleation temperature are computed using CosmoTransitions [130]. The solid lines depict the outcome of the gauge-invariant method. Stronger transitions correspond to faster moving bubble walls. The perturbative fluid approximation becomes worse and breaks down for stronger phase transitions, and so we cut off our scans above φ h (T n )/T n ∼ 1.1. Wall velocities for stronger phase transitions will only be larger than those shown. The dashed lines depict the resulting wall velocities when including the gauge boson cubic terms and friction. The values of v w are smaller in this case, although for Set 2 the gauge boson contribution makes a less significant difference. This suggests that our gauge-invariant treatment should provide a reasonable, though rough, estimate of the wall velocity when v w is not too large. smaller pressure difference between the phases. The friction on the bubble wall also tends to be enhanced for larger thermal masses.
Interestingly, for strong first-order phase transitions, we find that subsonic solutions to the equations of motion may not exist. This is because as v w → c s , the background temperature contribution begins to dominate in the Higgs and singlet field EOMs (it is proportional to 1/(c 2 s − v 2 w )). As pointed out in Ref. [54], the background terms typically enter with a relative sign to those from the heavy species, thus reducing the total friction for subsonic deflagrations. This behavior is seen for Set 1 in Fig. 2: no subsonic solution exists for the gauge-invariant case with φ h (T n )/T n 1. Including the gauge-dependent terms, subsonic solutions can extend up to φ h (T n )/T n ∼ 1.1, but not higher. We conclude that viable non-local electroweak baryogenesis in singlet-driven models is incompatible with very strong first-order phase transitions, at least in some cases. This can be at odds with sphaleron suppression inside the bubble, as seen for Set 1.
Even if a subsonic solution exists, the bubbles tend to expand rather quickly from the standpoint of successful EWB. For example, previous studies of CP -violating sources in the MSSM [27][28][29] suggest that electroweak baryogenesis tends to be most efficient for v w ∼ 0.01, while Fig. 2 indicates that v w > 0.2 for most points featuring a strongly firstorder phase transition. Viable bayogenesis in singlet-driven scenarios may thus require substantially more CP -violation than in models with slow walls (such as the MSSM with light stops) to overcome the suppression arising from large v w .
Our methods also allow us to determine the wall widths and offset for the subsonic configurations. These quantities are important inputs for microphysical calculations of the baryon asymmetry. The resulting bubble wall profiles for Sets 1 and 2 are shown in Fig. 3. The offset can change sign, with the singlet field lagging behind that of the SM-like Higgs for stronger phase transitions. For φ h (T n )/T n 1, the wall widths are typically ∼ O(5/T ). This is substantially smaller than typical values in Standard Model-like cases and consistent with the findings of Ref. [48] in the NMSSM. Thin walls follow from the large pressure difference due to the changing singlet VEV during the transition. This is in fact promising for electroweak baryogenesis, since in many cases the CP -violating sources scale as ∼ 1/L w [28,135].
One may ask to what extent we should expect similar results beyond the minimal real singlet extension of the Standard Model. After all, the xSM is known to be incomplete from the standpoint of electroweak baryogenesis, since it contains no new source of CPviolation. However, the model can be modified slightly to incorporate CP -violation by e.g.
complexifying the singlet and adding CP -violating Higgs-singlet couplings, or by including additional higher dimension CP -violating operators, as in Ref. [112], in cases where the singlet VEV vanishes at T = 0. Neither possibility should significantly alter the friction on the bubble wall. We expect similar conclusions in other CP -violating extensions of the xSM. Our findings followed primarily from the form of the friction, which is dominated by the top quarks and gauge bosons for the SM-like Higgs field and the singlet and Higgs excitations for the singlet field. As long as this is the case, the results beyond the minimal model should be qualitatively similar to those we have found here. In fact, this is not unreasonable: it would be difficult for new states to couple as strongly to the Higgs field as the top quark without violating existing phenomenological constraints, for example. Regardless, the methods and ingredients presented in sections 4-5 can be used to determine the wall velocity beyond the minimal xSM, although this may require computing additional interaction rates involving the excitations of the new species in the plasma.

Summary and Conclusions
In this study, we have seen how to compute the electroweak bubble wall velocity at singletdriven first-order phase transitions. This extends previous work which applied to the Standard Model-and MSSM-like cases. For concreteness, we framed our discussion in the real singlet extension of the Standard Model, or xSM, although our methods can be used in other models involving singlets at the electroweak phase transition.
Some of the key findings of this study are: • As anticipated, bubbles tend to expand rather quickly at first-order phase transitions driven by tree-level cubic terms in which the singlet vacuum expectation value changes appreciably. We have found v w 0.2 for all points with φ h (T n )/T n ≥ 1 in the xSM. These wall velocities may be compatible with electroweak baryogenesis in some cases, provided a sufficiently strong source of CP -violation. One should bear in mind, however, that the free-particle approximation made for the singlet excitations may overestimate the corresponding friction, and thus lead to an underestimate of v w .
• The most promising parameter space for slower bubble walls, and hence for electroweak baryogenesis, features larger thermal masses for the singlet and SM-like Higgs field. This translates into larger values of a 2 and b 4 in the xSM.
• Strong phase transitions may exhibit no subsonic solution and hence not allow for viable transport-driven electroweak baryogenesis. For example, considering a particular set of parameters in the xSM, we found no points with v w < c s for φ h (T n )/T n 1.1.
• A gauge-invariant estimation of the bubble wall velocity is possible and should approximate the full solution rather well for the slowest walls. The (gauge-dependent) SU (2) L gauge boson contributions become important for faster moving bubble walls.
• Wall widths are typically of order ∼ 5/T for strong first-order phase transitions as required for electroweak baryogenesis. These values are considerably smaller than their Standard Model analogs which are often used in the literature.
Along with the above points, our treatment of the friction on the bubble wall can be useful in various applications related to the electroweak plasma. For example, the interaction rates computed for the top quarks and Higgs bosons can be used to extract diffusion constants for these species, valid at leading log order and including the effects of hard thermal loops, which are important in transport calculations for electroweak baryogenesis. A rough estimate along the lines of Ref. [21,77] suggests D h ∼ 13/T , D t ∼ 2/T .
Our results are promising from the standpoint of observable gravitational radiation. The peak amplitude for the stochastic gravity wave background produced at a phase transition is enhanced for faster moving bubbles and larger pressure differences. The ingredients presented in this study can be used to more precisely compute the resulting gravity wave spectrum in concrete models involving singlets. Recent work [67,68] suggests that the peak amplitude of the signal from a strong electroweak-scale phase transition may in fact be significantly larger than previously realized. It would be interesting to analyze singlet driven phase transitions given these new hydrodynamic insights along with our predictions for the wall velocity in concrete models. An observable gravity wave signal could provide exciting (indirect) evidence for a first-order phase transition, and possibly electroweak baryogenesis, in the early Universe.
There are several ways to improve over the methods presented this study. One might hope to move beyond the simple fluid approximation to be able to study stronger phase transitions. Also beneficial would be a full leading-order determination of the quasiparticle interaction rates entering the Boltzmann equations for the various perturbations. Other improvements include accounting for the effects of the spherical bubble geometry on the hydrodynamics, formulating a gauge-independent treatment incorporating the gauge and Goldstone bosons (which is a difficult problem), and considering full numerical solutions to the bubble wall equations of motion rather than utilizing an ansatz. These improvements, required for more precise determinations of the wall velocity, would become much more important if the LHC or a future collider were to unearth direct evidence for a singlet-extended Higgs sector. In the interim, we expect our methods to provide a decent approximation of the bubble wall velocity in singlet-driven scenarios, which remain a particularly compelling setting for electroweak baryogenesis in light of current experimental constraints.

A Comparing Effective Interaction Rates
The effective interaction rates we have computed and listed in Equations 4.25-4.29 differ by O(1) factors from those appearing in the classic references [20,77]. In this appendix, we show that this discrepancy can be explained by the slightly different set of approximations made in evaluating the integrals analytically in previous works. Technically, both the methods used in Refs. [20,77] and in this work are valid 'leading log approximations' to the full first-order result, in that only processes contributing logarithmically to the collision integrals are considered. They differ primarily in their treatments of the non-logarithmic pieces of the various momentum integrals. The discrepancies can therefore be understood to demonstrate the uncertainties associated with the leading log approximation and the importance of performing a full leading-order calculation for more precise results in future studies.
To understand the different approaches to evaluating the collision integrals, let us consider the process tt → gg. First of all, Ref. [77] does not include the symmetry factor for the corresponding matrix elements, resulting in a factor of two discrepancy before evaluating any integrals [121]. Including the symmetry factor, the leading-log matrix element is ≈ 64/9g 4 3 u/t, where we have averaged over the top quark degrees of freedom. The resulting contribution to Γ t µ1,t , evaluating the integral numerically and including the top quark momentum-dependent self-energy, is including the correct symmetry factor. Starting from the same matrix element, the different methods for evaluating the integrals results in almost a factor of 4 difference in the result. Simply evaluating the integral in Ref. [77] numerically, but including the thermal mass in the propagator instead of the HTL momentum-dependent self-energy, we find ∆Γ t µ1,t ≈ 1.5 × 10 −3 T, (A. 3) suggesting that the simple propagator replacement over-estimates the integral by about 40%, but does not account for the whole discrepancy. In fact, most of the difference comes from the various approximations made to arrive at the analytic result in Eq. A.2. In particular, all non-logarithmic contributions are dropped, while numerically evaluating the integrals includes all of the various contributions. For example, the final result for ∆Γ t µ1,t in Ref. [77] contains an integral over the plasma frame angle θ between p and k, Moore and Prokopec drop the constant piece, keeping only the logarithm. However, the contribution of the constant piece is numerically comparable, and of opposite sign, to the logarithmic term. Performing the remaining integrals over |p|, |k|, the contribution without the constant term is ∼ 3.1 × 10 −3 T , while including it yields ∼ 1.8 × 10 −3 T , which is significantly closer to the results we have obtained. Similar approximations are made in performing the other integrals, and for the other rates. It is worth reiterating that neither approach includes all processes contributing at leading order in the gauge couplings. Our interaction rates should simply be viewed as a slightly different approximation to the full leading order result. Future studies of the wall velocity in different models may find it beneficial to compare the results obtained from our interaction rates and those of Ref. [77] to assess the uncertainty expected from the leading log approximation.