Lower bound on the electroweak wall velocity from hydrodynamic instability

The subsonic expansion of bubbles in a strongly first-order electroweak phase transition is a convenient scenario for electroweak baryogenesis. For most extensions of the Standard Model, stationary subsonic solutions (i.e., deflagrations) exist for the propagation of phase transition fronts. However, deflagrations are known to be hydrodynamically unstable for wall velocities below a certain critical value. We calculate this critical velocity for several extensions of the Standard Model and compare with an estimation of the wall velocity. In general, we find a region in parameter space which gives stable deflagrations as well as favorable conditions for electroweak baryogenesis.


Introduction
A first-order electroweak phase transition may explain the observed baryon asymmetry of the universe (BAU). Indeed, such a phase transition would provide all the Sakharov conditions, namely, baryon number violation, C and CP violation, and a departure from thermal equilibrium. For a quantitatively successful electroweak baryogenesis (EWB) an extension of the Standard Model (SM) is needed, such that there is enough CP violation as well as a sufficiently strong first-order phase transition. In the standard mechanism for EWB (see [1] for a recent review), the departure from equilibrium acts in two different ways. On the one hand, the expansion of bubbles of the broken-symmetry phase builds up non-equilibrium particle densities in front of the bubble walls. These densities are asymmetric for left handed particles and their antiparticles due to CP violating interactions with the wall. This asymmetry is transported to the unbroken-symmetry phase, where it biases the weak sphaleron processes which violate baryon number. The generated baryon asymmetry reenters the broken-symmetry phase. As a result, the bubble walls leave behind a net baryon number density. On the other hand, before this baryon number density recovers the equilibrium, the baryon number violating processes should be turned off. Otherwise the generated BAU would be washed out. Such a suppression of the sphaleron processes indeed occurs inside the bubble, as long as the Higgs background field φ b in the broken-symmetry phase satisfies the well known condition where T is the temperature. The ratio φ b /T plays the role of an order parameter, and the condition (1) expresses the baryogenesis requirement of a strongly first-order phase transition.
Although for a Higgs mass as large as 125GeV the electroweak phase transition is a smooth crossover, many extensions of the SM give strongly first-order phase transitions. Most investigations of EWB concentrate in the value of the order parameter φ b /T and the sources of CP violation for specific models. Since the computation of the velocity v w of bubble walls is too involved, a specific value is often assumed (typically v w = 0.1) to obtain a result for the BAU 1 . However, the generated BAU has also an important dependence on v w . Indeed, for very small velocities thermal equilibrium is restored and a small baryon asymmetry is generated. On the other hand, if the wall velocity is too large the diffusion of left-handed density perturbations is not efficient, and the resulting baryon number density is again small. In other words, a departure from equilibrium is needed, but such a departure should not be too strong. As a consequence, the generated baryon asymmetry has a maximum for a certain wall velocity v w = v peak . The value of v peak depends on the time scales associated to particle diffusion and baryon number violation, and is in general in the range 10 −2 < v peak < 10 −1 (see, e.g., [3,4,5]). A sizeable BAU is more easily obtained if v w is close to v peak . Moreover, any model which gives supersonic velocities is in conflict with the standard EWB mechanism [6].
Subsonic wall velocities are possible due to the friction with the plasma, which generally causes the walls to reach a terminal velocity. This velocity is given by the balance between the driving force, which depends on the pressure in the two phases, and the friction force, which depends on the microscopic interactions of the bubble wall with plasma particles. The driving force is very sensitive to hydrodynamics. As a consequence of nonlinear hydrodynamics, there are different kinds of stationary solutions for the propagation of the wall [7]. The solutions which can be realized in a cosmological phase transition (see, e.g., [8,9,10,11,12,13,14]) are weak deflagrations, which are subsonic, Jouguet deflagrations, which are supersonic, and Jouguet or weak detonations, which are supersonic too. Hence, the case of interest for baryogenesis is that of weak deflagrations.
It is well known that the stationary propagation of a weak deflagration front may be unstable [7]. For the case of a relativistic equation of state, the stability of deflagrations was first studied in Ref. [15]. The result was that deflagrations are always unstable under perturbations above a certain wavelength. This analysis was improved in Ref. [16]. The main improvement was to take into account the dependence of the stationary velocity on the temperature. The main result of Ref. [16] was that the deflagration is stable for wall velocities above a certain critical value v crit . Numerical simulations [17] agree with such a stabilization. In Ref. [18], the results of [16] were improved by taking into account temperature fluctuations on both sides of the wall, as well as the fact that the reheating due to the release of latent heat depends on the wall velocity. For small amounts of supercooling, the results of Ref. [18] agree with those of Ref. [16].
In Ref. [16], the stability analysis was applied to the electroweak phase transition for the minimal Standard Model with unrealistic values of the Higgs mass, which gives a strong enough phase transition for EWB. For Higgs masses higher than m H = 40GeV , the critical velocity below which deflagrations are unstable was found to be v crit 0.07. This result was compared with the wall velocity calculations [19,20], which gave v w 0.1. Therefore, the result of Ref. [16] indicated that the electroweak deflagration is stable. However, both v w and v crit depend on the model and should be recalculated for each extension of the SM.
The calculation of v w is more involved than that of v crit and depends on more details of the model. Indeed, the value of the critical velocity depends only on thermodynamical parameters which can be derived from the free energy density. In contrast, the actual value of the stationary velocity depends (besides thermodynamics and hydrodynamics) on the friction of the wall with the plasma. The computation of the friction force involves considering Boltzmann equations for the out-of-equilibrium particle densities in front of the wall. For the SM, a thorough calculation (including reheating effects) [21] gave wall velocities in the range 0.36 < v w < 0.44 (for 0 < m H < 90GeV). A similar calculation for the Minimal Supersymmetric Standard Model (MSSM) [22] gave smaller velocities, v w = (5 − 10) × 10 −2 , due to the larger particle content of this model (essentially, due to the contribution of top squarks). To our knowledge, these two results constitute the only detailed microphysics calculations for specific models. The reason for this is the difficulty of computing the collision terms for the Boltzmann equations. In spite of this, many investigations of the friction were performed. In particular, a study of the overdamped evolution of gauge fields [23] showed that infrared boson excitations generally increase the friction and, consequently, cause smaller wall velocities than previous studies. In particular, for the SM the estimated wall velocity was v w 0.01 for m H ≃ 80GeV and v w ≃ 0.1 for m H ≃ 45GeV.
In this paper we shall investigate the possible instability of the electroweak wall velocity for several extensions of the SM. The main motivation for this is the fact that the deflagration instability may affect the baryogenesis scenario. Indeed, notice that the value of v crit obtained in Ref. [16] lies within the optimal range for EWB. Moreover, given the general uncertainties and large errors in the estimations of v w , the value of v crit , which is much easier to calculate, provides a lower bound for v w which may be important to constrain baryogenesis. It is worth mentioning also that an instability of the stationary wall propagation may have several cosmological consequences, such as the generation of magnetic fields [24] or gravitational waves [18].
The plan of the paper is the following. In Sec. 2 we review the hydrodynamics of a wall which propagates as a deflagration and we discuss the stability of such a stationary solution as a function of thermodynamic parameters. In Sec. 3 we calculate the critical velocity below which the wall becomes unstable. We consider the electroweak phase transition for several extensions of the Standard Model. We also estimate the wall velocity for each model in order to study the stability as a function of the parameters. In Sec. 4 we discuss on the possible consequences of the instability. Finally, in Sec. 5 we summarize our conclusions. Details of the calculation of the phase transition dynamics are contained in App. A. Further discussion on the critical velocity as well as a fit can be found in App. B.
2 Stationary wall propagation and hydrodynamic stability 2.1 First-order electroweak phase transition The relevant quantity describing the phase transition is the free energy density or finitetemperature effective potential 2 F (φ, T ). At a given temperature T , the minima of F give the possible thermal expectation values of the Higgs field φ, which determine the different phases. For the electroweak theory, we have a phase transition from the symmetric phase to the broken-symmetry phase at a critical temperature T c ∼ 100GeV. In the minimal SM, the electroweak phase transition is just a smooth crossover. In Sec. 3 we shall consider several extensions of the SM for which the electroweak phase transition is first-order. For a first-order phase transition, there is a range of temperatures around T c for which the effective potential has two minima separated by a barrier. For T > T c the symmetric minimum φ = 0 is the absolute minimum of F (φ, T ), while for T < T c the absolute minimum has a nonvanishing value φ b (T ). We shall use subindexes u and b for the unbrokenand broken-symmetry phase, respectively. These phases are thus characterized by the free energy densities F u (T ) = F (0, T ) and F b (T ) = F (φ b (T ), T ). The critical temperature is given by the equation F u (T c ) = F b (T c ). The energy density and pressure for each phase are obtained from the free energy density through ρ(T ) = F (T )−T F ′ (T ), p(T ) = −F (T ), where a prime indicates a derivative with respect to the temperature.
A first-order phase transition occurs via the nucleation and expansion of bubbles. As we shall see, the relevant parameters for our calculation will be the latent heat L defined as L ≡ ρ u (T c ) − ρ b (T c ), the enthalpy density before the phase transition, w u (T c ) = ρ u (T c ) + p u (T c ), and the nucleation temperature T n . The latter is the temperature at which bubbles effectively begin to nucleate. We describe its calculation in App. A. The enthalpy density and latent heat are given by After nucleating, bubbles expand due to the higher pressure of the stable phase. In most cases, the bubble walls quickly reach a terminal velocity due to the friction with the plasma. We shall concentrate in such a case. From the time a bubble nucleates to the time their walls collide with other bubbles, the temperature of the plasma varies due to the adiabatic expansion of the universe and due to the release of latent heat. As a consequence, the wall velocity will vary too. We shall use T = T n as a representative value for the temperature of the phase transition.

Microphysics and hydrodynamics
The motion of a bubble wall can be derived from the equation for the Higgs field in the plasma (see, e.g., [19,20,21,22,23,25,26,27,28]), with E i = p 2 + m 2 i , where m i are the φ-dependent particle masses, δf i are the deviations of particle densities from equilibrium, and the sum runs over all particle species. The last term gives the friction with the plasma. In order to transform Eq. (4) into an equation for the bubble wall, the general procedure is to use some approximation or ansatz for the field profile (which is static in the reference frame of the wall) and integrate across the wall. Since the deviations δf i depend on the wall velocity, the last term in (4) gives the friction force F fr , while the second term gives the driving force F dr . Thus, the steady state velocity of a bubble wall is given by the force balance F dr = F fr (for a more detailed explanation, see e.g. [29]).
The driving force is relatively easy to calculate. In particular, if the temperature remains constant across the wall, we have F dr = p b (T ) − p u (T ). On the other hand, the friction is proportional to the departure of the plasma particles from their equilibrium distributions. This departure from equilibrium depends not only on the interaction of the particles with the Higgs field at the wall, but also on the interactions of plasma particles away from the wall. Thus, the calculation involves solving a system of Boltzmann equations for the population densities of the relevant species. The Boltzmann equations include collision terms which must be computed by calculating the scattering rates for all the relevant processes. Such a calculation is often referred to as the microphysics calculation.
As already mentioned in Sec. 1, the microphysics calculation is very difficult. In particular, the computation of the collision terms was carried out only for a couple of models and in the non-relativistic (NR) case [21,22,23]. The result is a friction force of the form F fr = η NR v w . The friction coefficient η NR is very model dependent, and its calculation involves the use of several approximations. The ultra-relativistic (UR) limit has also been considered [30]. It turns out that this limit is even simpler than the NR case. The result is that the friction saturates for v w → 1, i.e., the friction force reaches a velocityindependent value F fr = η UR . Intermediate cases are much more difficult to treat. In a recent treatment [31], the friction was considered beyond the small wall velocity regime. However, the deviations from equilibrium were still considered to be small. In particular, the friction force calculated in Ref. [31] does not match the UR results of Ref. [30].
In order to overcome the difficulties of the microphysics calculation, a phenomenological approach has often been used, which consists in replacing the last term in Eq. (4) with an effective damping term of the form u µ ∂ µ φ, where u µ = (γ, γv) is the four velocity of the fluid (see, e.g., [12,32,33]). If we ignore hydrodynamics, this approach gives a friction force of the form F fr = ηv w in the NR limit, where the coefficient η is a free parameter coming from the phenomenological damping term. Hence, setting η to the value η NR from the microphysics calculation gives the correct friction force in this limit. This phenomenological model extrapolates the NR behavior ηv w to a force of the form ηγv for relativistic velocities, where v is the velocity of the fluid relative to the wall (which will vary across the wall, see below). However, this simple model does not give the friction saturation in the UR limit. In order to reproduce the saturating behavior, in Ref. [34] a phenomenological model which gives a friction force of the form ηv was considered. However, such a model with a single free parameter η can hardly match both the NR and UR forces η NR v and η UR . In Ref. [29] a phenomenological interpolation between the two regimes was considered.
For the subsonic velocities we are interested in, the exact dependence of the friction on the velocity does not introduce significant effects [29]. Therefore, we shall use the phenomenological scaling ηvγ, which will allow us to use the results of the stability analysis [18]. We shall discuss the implications of this choice in Sec. 4. In order to estimate the wall velocity for specific models, in Sec. 3 we shall set the parameter η to the value η NR coming from microphysics calculations.
Using this phenomenological approach to the friction, it is not difficult to include the hydrodynamics, i.e., to take into account the change of fluid variables across the wall. The fluid variables in each phase are related by matching conditions at the phase discontinuity. In the rest frame of the wall, we have [7] w u γ 2 where v is the component of the fluid velocity along the wall motion, v ⊥ is the velocity in the transverse direction, and γ = 1/ √ 1 − v 2 . By symmetry, we set v ⊥ = 0 for the stationary motion, but this component must be taken into account in the stability analysis. Furthermore, we define the incoming and outgoing flow velocities by −v u and −v b , respectively, so that we deal with positive values of the variables v u , v b . Using suitable approximations (see App. A), one obtains a friction force for a quantity f defined on each side of the wall. The thermodynamical quantities w(T ), p(T ) in Eqs. (5) are related by the equation of state. The treatment of hydrodynamics is considerably simplified by considering a simple approximation for the equation of state. This is particularly important for the stability analysis. In order to use the analytical results of Ref. [18], we shall consider the bag EOS. This is the simplest EOS which can describe a phase transition 3 . Due to the difficulty of the stability analysis, it is not trivial to generalize the results of Ref. [18] beyond these 3 One limitation of this model is the fact that it gives a constant speed of sound c s = 1/ √ 3 in both phases, while the actual value of c s may depart from this value [35].
approximations. With these approximations, the driving force becomes Several of the quantities in Eqs. (6)(7) are related through Eqs. (5) and boundary conditions (see App. A). As a result, the wall velocity v w depends only on the temperature T u , the coefficient η, and the parameter Notice that in general the driving force does not coincide with the pressure difference p b − p u , as one might expect, since it is nontrivially affected by hydrodynamics. Nevertheless, the fact that the pressure difference changes sign at T = T c is reflected in F dr , which is very sensitive to the departure from the critical temperature. In particular, the reheating of the fluid causes a decrease of F dr , acting as an effective friction [33,36].
There are in general several solutions for v w . The only solution for which the wall velocity is subsonic is the weak deflagration, and we shall only be interested in this case. Thus, in the bag approximation, we have v w < 1/ √ 3. The deflagration solution exists for large enough friction and small enough supercooling. For a deflagration, the fluid is reheated in front of the wall. Therefore, the temperature T u is higher than the nucleation temperature T n . The relation between T u and T n introduces an equation to solve together with that corresponding to the equilibrium of the forces (6-7). We write down all these equations in App. A.

Stability of the deflagration
The possible hydrodynamic instabilities of cosmological phase-transition fronts have been investigated in Refs. [15,16,17,18,37,38,39]. The linear stability analysis of the wallfluid configuration consists of considering small perturbations of the fluid variables on both sides of the wall, together with small deformations of the planar and infinitely thin wall. Below we briefly sketch the generalities of the calculation for the deflagration case. For detailed and more general analysis, see [18,39].
There are in principle seven variables, namely, the wall deformation ζ, the pressure fluctuations δp u , δp b , the variations of the fluid velocity along the propagation direction δv u , δv b , and the transverse velocities v ⊥ u , v ⊥ b . These perturbations depend on space and time. The three fluid fluctuations on a given side of the wall are related by the fluid equations. Linearizing these equations and looking for solutions of the form exp(ik · x ⊥ + qz +Ωt), one obtains dispersion relations for q(k,Ω), as well as algebraic equations relating the amplitudes of the different fluctuations. For the stability analysis one is interested in the unstable modes, which correspond to Ω > 0 and qz < 0 [7]. For weak deflagrations, we have one unstable mode in front of the wall and two unstable modes behind it. We are thus left with four unknowns, namely, the amplitudes of the three unstable fluid modes and that of the wall perturbation. Finally, the fluid perturbations on the two sides of the wall are related by junction conditions at the wall, which are the counterparts of Eqs. (5). There is also an equation for the perturbations of the surface, which is the counterpart of the equation F dr = F fr . As a consequence, one obtains a system of four algebraic equations for the four unknowns. The weak deflagration is linearly unstable if a nontrivial solution exists for this system.
Looking for linear instability is thus equivalent to demanding a 4 × 4 determinant to vanish. This gives an equation for the growth rate Ω as a function of the perturbation wavenumber k ≡ |k|. For a complete treatment and analytic expressions we refer the reader to Ref. [18]. The general result is that the deflagration is unstable for wavenumbers below a certain value k c which depends on the thermodynamic parameters and the wall velocity. Nevertheless, depending on the parameters the critical wavenumber k c may become negative, in which case all wavenumbers are stable. Demanding this to be the case (i.e., k c ≤ 0) one obtains a condition for the wall velocity, namely, v w ≥ v crit for stability of the deflagration. The critical velocity v crit corresponds to k c = 0, which is equivalent to the equation 4 where γ 2 For the bag model, the quantity β is given by For v w < v crit , perturbations with wavenumbers 0 < k < k c grow exponentially. As discussed in Sec. 4 below, for the electroweak phase transition the characteristic time for the development of the instabilities is generally much shorter than the time scales associated to bubble growth, even if v w is very close to v crit . It is worth remarking that the critical value of the wall velocity does not depend on the friction coefficient η. The parameter η determines the actual value of v w . The critical velocity v crit can thus be associated to a critical friction η crit (which could then be compared with the actual value of η in order to determine the stability). However, the form of the friction force (6) is implicit in the equations above. This is because, in the stability analysis, the coefficient η is eliminated by writing it as a function of the velocity and the driving force [16,18].
The critical velocity depends only on the dimensionless parametersL and T n /T c . Both parameters range between 0 and 1 and quantify the amounts of latent heat and supercooling, respectively. In Fig. 1 we show the curves of constant v crit in the space of these two parameters [we considered the parameter (T c − T n )/T c since in many physical cases the temperature is very close to T c ]. These curves are model independent. We also show the points in this plane corresponding to some of the specific models considered below (those which span larger regions of the plane). The limitL → 0, T c − T n → 0 corresponds to a second-order phase transition. The critical velocity increases with the amount of supercooling but is rather insensitive to the latent heat. The curves of constant v crit accumulate at the weak deflagration limit v crit = c s ≃ 0.577. For values of the parameters above this curve, there is no critical velocity and any subsonic velocity is unstable. In App. B we consider in more detail the dependence of v crit on these parameters and we give a simple fit for v crit (L, T n /T c ). In Fig. 1 it seems that, for physical models, none of the two parametersL or (T c − T n )/T c can reach the upper value 1. In fact, the supercooling parameter may in principle take values arbitrarily close to 1, since the nucleation temperature T n can be very small for models with barriers in the zero-temperature effective potential. The upper limits of (T c − T n )/T c in these curves are due to the break-down of our calculations for such strong phase transitions. On the other hand, it is true that the latent heat tends to take relatively small valuesL 0.1 for physical models, as seen in the figure. This is because, in general, the entropy released in the phase transition is only a fraction of the total entropy. This fraction is proportional to the fraction of degrees of freedom (d.o.f.) which are strongly coupled to the Higgs.
As can be seen in Fig. 1, the wide range of possible amounts of supercooling implies a wide range of possible values of the critical velocity v crit . On the other hand, the actual value of the wall velocity, v w , also grows with the amount of supercooling, and it is not easy to guess, without specific calculations, in which cases the deflagration will be unstable. We perform such calculations in the next section.

Electroweak phase transition models
Several extensions of the SM have been considered in the literature in order to increase the strength of the phase transition and obtain a sizeable electroweak baryogenesis. Constraints from the recent LHC results threaten the viability of baryogenesis in some models (see, e.g., [40,41,42]). Since our aim is a general investigation of the possible instability, we shall not discuss the implications of experimental constraints. Our results will give in principle an additional constraint for the baryogenesis scenario. Furthermore, we shall not limit ourselves to the case of strongly first-order phase transitions, since the instability of the stationary wall propagation may have cosmological consequences even for weakly first-order phase transitions.
A classification of models for the electroweak phase transition was recently given in Ref. [40]. Besides these model classes, we shall consider a model with TeV fermions introduced in Ref. [43] as well as two-loop effects. For simplicity, we shall consider a single background field φ. For several extensions of the SM, φ corresponds to the SM Higgs, H 0 ≡ φ/ √ 2. In extensions in which more than one scalar develop a vacuum expectation value (VEV), it is sometimes a good approximation to consider a single field φ, corresponding to a certain trajectory in field space.

Free energy and friction
We shall consider the one-loop finite-temperature effective potential given by where V 0 is the tree-level potential and V 1 , F 1 are the zero-temperature and finitetemperature parts of the one-loop correction. These corrections receive contributions from the SM particles and from beyond-SM particles. We shall consider a spontaneous symmetry-breaking potential of the form as well as some tree-level modifications. Here, the parameters m and λ are related to the Higgs mass and VEV by v = 2/λm = 246GeV, m H = √ 2λv 2 = 125GeV, and the constant term in Eq. (12) was added so that the potential vanishes at the minimum, i.e., V (v) = 0. We shall assume Higgs-dependent masses of the form and we shall consider the renormalized one-loop zero-temperature correction where the upper and lower signs correspond to bosons and fermions, respectively, and g i is the number of d.o.f. of particle species i. This expression corresponds to the renormalization conditions that the tree-level values of the minimum and Higgs mass are not shifted by radiative corrections, i.e., V ′ 1 (v) = 0, V ′′ 1 (v) = 0 (where a prime indicates a derivative with respect to φ). We have added a term m 4 i (v)/2 to the well known expression, so that the true vacuum energy density is not shifted either, i.e., V 1 (v) = 0. Thus, in the symmetric phase we will have a false vacuum energy density given by ρ vac = V 0 (0) + V 1 (0), which contributes to the Hubble rate during the phase transition. Finally, the one-loop finite-temperature correction, including the resummed daisy diagrams, is given by [44] The last term in Eq. (15) receives contributions from all the bosonic species. We have are the thermal masses. For the transverse polarizations of the gauge bosons we have Π(T ) = 0 and the last term in (15) vanishes.
In general, Eqs. (11)(12)(13)(14)(15) give a phase transition at a certain temperature T c from the high-temperature symmetric phase to the low-temperature broken-symmetry phase. The dynamics of the phase transition depends mostly on the difference Depending on the particle content of the model, the phase transition may be first order, i.e., the effective potential may have two minima separated by a barrier. This is better appreciated in the high-temperature approximation, i.e., when Eq. (15) admits an expansion in a power series of m/T . For masses of the form (13) we obtain the well known simple form [45] V with coefficients given approximately by withc i = 1 for bosons and 1/2 for fermions 5 . We have neglected for simplicity corrections to these coefficients which are suppressed by factors of order g i /(32π 2 ). The terms m 4 log m 2 in (14) cancel out with similar terms in the expansion of the thermal integrals in Eq. (15), and only a dependence of the form φ 4 log T 2 remains. This gives a soft dependence of the coefficient λ T on T , which was neglected in the approximation (18). For the coefficient E the sum runs only over transverse gauge bosons. Indeed, only the bosonic thermal integral has a cubic term ∼ (m/T ) 3 in its power expansion. The last term in Eq. (15) replaces m 3 with M 3 and we actually have For gauge bosons we have µ = 0, and for their transverse polarizations we also have Π = 0. Hence, only in this case we obtain a contribution to the term −ET φ 3 . It is well known that, for an effective potential of the form (17), it is precisely the cubic term the one which allows a first-order phase transition, while for E = 0 the transition is second order. Indeed, the value of the order parameter at the critical temperature is given by φ b (T c )/T c = 2E/λ Tc . For the minimal SM, the constant E is very small and Eq. (17) gives a very weakly first-order phase transition 6 . The relevant SM contributions to the one-loop effective potential come from the Z and W bosons, the top quark, and the Higgs and Goldstone bosons. The Higgs sector is usually ignored in the one-loop radiative corrections. This should be a good approximation in extensions of the SM which include particles with strong couplings to φ. We shall use this simplification. Therefore, the φdependent masses are of the form h i φ, with h i = m i /v, where m i are the physical masses at zero temperature. Thus, in the case of the SM we have Since only bosons contribute to the parameter E, extensions of the SM containing extra bosons are often considered in the literature. However, a cubic term in the effective potential is difficult to obtain due to the thermal mass Π ∼ h 2 T 2 in Eq. (19). This has lead to the investigation of scenarios in which a negative squared mass cancels the thermal mass, µ 2 ≃ −Π(T ). For high values of the couplings h i the expansion (17) breaks down. In this case, the one-loop terms φ 4 log φ in the zero-temperature effective potential (14) may strengthen the phase transition by causing a barrier at T = 0. Two loop contributions and tree-level modifications to the effective potential have also been considered in order to increase the strength of the electroweak phase transition. In the present paper we shall consider all these extensions of the SM.
A simple approximation for the general dependence of the friction coefficient η on the parameters of the model was derived in Refs. [46,47]. We shall use this approximation to estimate the wall velocity (for a different approach see [48]). In this approximation, η receives contributions from particles which obey the Boltzmann equation as well as from infrared excitations of bosonic fields, which are treated classically. The contribution of Boltzmann particles is given by where φ c = φ b (T c ), the function c 1 is given by andΓ is an average interaction rate arising from the collision terms of the Boltzmann equations. The friction decreases with this parameter since the deviations from the equilibrium distributions in front of the wall will be smaller if the processes are quicker. For the electroweak phase transition,Γ is typically ∼ 10 −2 T . The infrared bosons contribution is given by where m D is the Debye mass, given by m 2 D = (11/6)g 2 T 2 for the W and Z bosons of the SM, and m 2 D = h 2 T 2 /3 for a scalar singlet. The integral in (23) has an infrared cut-off φ 0 for small µ i , given by φ 0 = L −2 w − µ 2 i /h i for µ i < L −1 w , and φ 0 = 0 for µ i > L −1 w , where L w is the wall width. In the thin wall approximation, L w can be estimated as The two contributions dominate in different parameter regions, and we have η = η B + η ir . Depending on the model, we shall use either the full one-loop effective potential (14)(15) or its high-temperature expansion (17). In the former case, we shall use the friction coefficients (21) and (23), while in the latter case we shall use a similar high-temperature approximation for the friction coefficients, and L w ≈ φ 2 /σ. Here, χ i = 2 for fermions and χ i = m i (φ) /T for bosons, and σ is the surface tension of the bubble wall.

The SM with a low Higgs mass
For comparison with previous results, we consider first the unrealistic case of the SM with a light Higgs. We also consider larger Higgs masses (although for large m H the perturbative expansion breaks down) in order to analyze the dependence of the velocity on the strength of the phase transition. For this model we use the approximation (17) for the effective potential as well as the approximations (25)(26) for the friction. The result is shown in Fig 2. The critical velocity for this model (solid line) agrees in order of magnitude, but not exactly, with Ref. [16]. For instance, for m H = 60GeV, they obtain v crit ≃ 0.035, while we obtain v crit ≃ 0.022. As already mentioned, as a function of the thermodynamic parameters, our result for v crit is in good agreement with Ref. [16]. The present discrepancy is due to the rough estimation of the temperature T u for this model in Ref. [16]. The dotted line in Fig. 2 corresponds to the wall velocity obtained from the calculation of the friction using only the Boltzmann equations. We remark that our approximation contains a single, effective rateΓ instead of the several collision-term parameters coming from the different interactions. SettingΓ = 10 −2 T , the values of v w , as well as the dependence on m H , agree with the results of Ref. [21] (and are also close to those of [31]). We stress that the calculations of Refs. [21,31] omit the contribution of infrared boson excitations to the friction, which can make the wall velocity significantly smaller [23]. In the case of the SM, the infrared term receives contributions from the gauge bosons. The result for the complete friction is shown in a dashed line in Fig. 2. This result agrees with the estimations of Ref. [23].
We see that the infrared contribution to the friction dominates for weakly first-order phase transitions, and causes the deflagration to become unstable (i.e., v w < v crit ).

Negative squared mass (thermal cubic term)
The simplest extension of the SM consists of adding one or more gauge-singlet scalars S i (see, e.g., [49,50]). In many models, these bosons constitute a hidden sector which couples only to the SM Higgs doublet through the so called Higgs portal operator h 2 s H † H S 2 i (assuming, for simplicity, real fields and universal couplings h i = h s ). The scalars may have SU(2) × U(1)-invariant mass terms µ 2 s S 2 as well as quartic terms λ s S 4 . For the moment we shall not consider the possibility that S develops a VEV.
If h s φ b (T )/T is not too large, the free energy is of the form (17), with the cubic term replaced by the term (19). The latter is not as effective as a cubic term in strengthening the phase transition. The thermal mass is given by Π = (h 2 s + λ s )T 2 /3 [50]. A negative value of µ 2 s may enhance the strength of the phase transition, since for µ 2 s ≃ −Π(T ) the term (19) is effectively of the form −T φ 3 . This fact is exploited in the case of the MSSM in the light-stop scenario [51]. Notice that negative values of µ 2 s may induce a nonvanishing expectation value of the extra scalar. For the case of top squarks, this introduces the danger of color breaking minima. We shall not take into account this issue here. For the moment we will just consider the SM plus g s extra bosonic degrees of freedom. We thus have an effective potential of the form with D = D SM + g s h 2 s /24, E = E SM , T 2 0 = (m 2 h /4)/D, and λ = m 2 H /(2v 2 ). For definiteness, we consider the case g s = 6 and λ s = 0.5 (below we consider a different value of λ s ). In Fig. 3   We have checked that for different values of g s or h s the behavior of the various curves is qualitatively similar. For weakly first-order phase transitions with small values of φ b /T we have a small v w as well as a small v crit , while for more strongly first-order transitions both v w and v crit are higher. On the other hand, their values generally cross at a certain point (as can be appreciated in the left panel of Fig. 3). Thus, we have v w < v crit for weak enough phase transitions (φ b /T < 0.5) and v w > v crit for stronger phase transitions. In particular, for φ b /T ≥ 1 we obtain stable deflagrations. Notice, however, that the curves of v w and v crit begin to approach each other again as the strength of the phase transition continues increasing. As we shall see below, for very strong phase transitions we will have v w < v crit again.
For g s = 6 and h s ≃ 0.7 this model can be regarded as a toy model for the light stop scenario of the MSSM, consisting of the SM plus the light right-handed top squark (stop). At the one-loop order and disregarding the possibility that the extra scalar develops a VEV, the main quantitative difference with the realistic case appears in the resummed daisy diagrams. The main contribution to the thermal mass of the stop is of the form 4g 2 str T 2 /9, where g str is the strong gauge coupling [51]. In our numerical calculation, we obtain the same effect by setting the parameter λ s in Eq. (27) to the value 4g 2 str /3. This increases considerably the value of Π(T ) and, hence, the value of the negative squared mass needed to compensate this thermal mass. The collision terms in the Boltzmann equations are also different due to a different particle content. This effect may be relevant if the wall velocity is close to the critical value. For the MSSM we can choose the value of our effective rateΓ in order to match the result of the detailed microphysics calculation [22]. Turning off the infrared contribution to the friction and setting 7 m H = 110GeV as in [22], the wall velocity should vary around v w = 0.1 (the friction is higher than in the SM due to the coupling of the extra boson with the Higgs). We obtain the correct v w for The lower, blue curves of Fig. 4 show the estimated wall velocity (including the infrared part of the friction) together with the critical value for this case. The result is plotted as a function of µ 2 s as well as of the mass of the extra scalar, which is given by m s = h 2 s v 2 + µ 2 s . We considered values of µ 2 s from −Π(T 0 ). For this value of µ 2 s , the infrared contribution to the friction lowers the wall velocity from v w ≃ 0.1 to v w ≃ 0.05, and the effect is stronger for weaker phase transitions (dashed blue curve). We see that we have v w < v crit , except for µ 2 s very close to −Π(T 0 ). In any case, we obtained values of φ b (T c )/T c smaller than 0.7 even in this limit. Indeed, in this scenario a phase transition which is strong enough for baryogenesis is obtained only for unrealistic values of the Higgs mass. The situation improves when two-loop corrections are considered.

Two-loop effects (the MSSM)
Although the light stop scenario does not give a strong enough phase transition at oneloop order, two-loop corrections can make the phase transition strongly first-order even without requiring negative µ 2 s [52,53]. The most important two-loop corrections are of the form φ 2 log φ, We consider such a contribution by adding to the potential (27) a term of the form (28), with a coefficient C = 6h 4 s − 8g 2 strong h 2 s coming from the MSSM stop and gluon loops.
The results for the wall and critical velocities with this modification are shown in Fig. 4 (upper, black curves). We obtain a higher critical velocity v crit ≃ 0.1-0.15 (black solid line) as well as a higher wall velocity v w ≃ 0.1 (black dashed line), but we have v w < v crit . Notice that the strength of the phase transition is essentially due to the presence of the term (28) in the effective potential. As a consequence, the result has a very soft dependence on m s . The significant increase of the wall velocity with respect to the one-loop value is due to the fact that the two-loop term does not change the particle content of the model (hence the friction cannot increase significantly) but does increase the strength of the phase transition. In this approximation we obtain φ b /T 1 for all the range of µ 2 s considered. Since the wall velocity is so close to the critical value, the uncertainties in the calculation of v w become relevant. In this case the Boltzmann part of the friction dominates. This can be seen by turning off the infrared contribution, which results in the dash-dotted curve of Fig. 4. Since the friction is dominated by the Boltzmann term, the O(1) errors in the estimation of collision terms are consequential for the stability of the deflagration. To see the sensitivity to the parameterΓ, let us consider again the SM-like valueΓ = 10 −2 T instead of the MSSM-like valueΓ = 5 × 10 −3 T . This gives the dotted curve 8 , which is safe from the instability. On the other hand, a different approximation for the friction [48] (which gives similar values for the one-and two-loop effective potentials) gives v w ∼ 0.05. Such deflagrations would be clearly unstable, since the velocity is a factor of 2 below the critical value.
Therefore, a more accurate determination of v w would be important for this model. It is also important to remark that in baryogenesis calculations the wall velocity is often assumed to be v w 0.1, which is just below the lower bound v crit for this model. Moreover, this model is severely constrained by experimental data (see, e.g., [54]). We stress that we have considered a simplified version of the light stop scenario, which has several parameters we just have not taken into account. In spite of this, we do not expect a significant difference for the critical velocity in the realistic case, although we do expect O(1) factors in the wall velocity, as already discussed.

Tree-level effects
Real gauge-singlets allow cubic terms of the form S 3 or H † HS. Such corrections to the Lagrangian arise in extensions of the SM with a singlet (see, e.g., [55,56,57,58,59,60]) as well as in extensions of the MSSM (see,e.g., [61,62,63,64]). The presence of cubic terms in the tree-level potential makes it easier to get a strongly first-order electroweak phase transition. Indeed, the strength of the transition is dominated by such cubic terms, which provide a barrier already at zero temperature. In order to study this effect, one should consider the effective potential for the condensates of the two fields H and S. However, our numerical computation of the nucleation temperature is based on a singlevariable potential. In order to incorporate this kind of model into our generic analysis, we shall assume that the thermal tunneling occurs through a trajectory in configuration space which can be parameterized with a single field φ(x), and that along this trajectory the zero-temperature effective potential has a cubic term [40]. This is equivalent to considering a toy model which consists of adding a term −Aφ 3 to the tree-level potential (12), where A is a free parameter with mass dimensions [47]. In this model the parameters of the potential are related to the physical Higgs VEV and mass by 2m 2 = λv 2 − 3Av, We thus consider the high-temperature effective potential (27) plus a term −Aφ 3 . We consider values of the parameters as in the left panel of Fig. 3 (g s = 6, h s = 0.7, λ s = 0.5,Γ/T = 0.01) with µ s = 0. We show the result in the left panel of Fig. 5. We have considered values of the parameter A for which h s φ b (T n )/T n 1, so that the high-temperature expansion is valid. The strength of the phase transition grows quickly with A, and values φ b /T > 1 are reached for values of A which are much smaller than the scale v. In this parameter range we obtain subsonic velocities. These deflagrations are stable on almost the entire considered range.  Another possible tree-level modification is the introduction of non-renormalizable operators. In particular, a dimension-six term of the form (H † H −v 2 /2) 3 /Λ 2 gives a negative quartic coupling if the cutoff is low enough [65,66,67,68,69]. We thus consider the SM plus a term of this form, which gives a term (φ 2 − v 2 ) 3 /(8Λ 2 ) in the effective potential. Adding this term to the tree-level potential (12) does not shift the value of the minimum nor the Higgs mass, i.e., we have v 2 = 2m 2 /λ, m 2 H = 2λv 2 . The sextic operator introduces the terms which modify the quadratic and quartic terms of the effective potential and, in particular, may change their sign. Thus, for Λ < 3/2v 2 /m H ≃ 600GeV the quadratic term becomes positive (already at zero temperature). Nevertheless, for Λ < √ 3v 2 /m H ≃ 840GeV the quartic term becomes negative, causing a barrier in the zero-temperature effective potential, which is stabilized by the term +φ 6 /(8Λ 2 ). For this calculation we use the full one-loop correction (15) instead of its high-temperature approximation.
The results for this model are shown in the right panel of Fig. 5 as functions of the cutoff. Below Λ ≃ 850GeV we have φ b (T )/T > 1. Since the particle content is the same as in the SM, we set againΓ = 0.01. The values of the wall velocity (dashed line) are in good agreement with those of Ref. [31] for stronger phase transitions (for instance, for Λ = 700GeV we obtain v w = 0.45 while in [31] the result is v w = 0.46). For weaker phase transitions we obtain smaller values of v w (e.g., v w ≃ 0.17 for Λ = 900GeV, in contrast to v w = 0.27 in [31]). This is to be expected since, for weaker phase transitions, the infrared part of the friction becomes noticeable. Notice that v w is subsonic in all the range of Λ we considered, but it approaches the speed of sound for Λ ≃ 600GeV. Therefore, an O(1) variation of the friction coefficient may give supersonic solutions. To illustrate this we consider the valueΓ/T = 0.02 (dash-dotted line). In this case, detonation solutions appear below Λ ≃ 675GeV. Notice also that the curves for the two values ofΓ coincide for weaker phase transitions, where the infrared contribution dominates the friction.
Regarding the stability, we see that the wall velocity becomes smaller than the critical value (solid line) only for weak phase transitions with φ b /T < 0.3, corresponding to wall velocities v w 0.015.

Strong coupling
Let us consider again the extension of the SM with a singlet scalar. This time however, instead of considering a negative µ 2 s in order to change the strength of the phase transition, we shall set µ s = 0 and vary the coupling to the Higgs. Therefore, the φ-dependent mass is given by m 2 s (φ) = h 2 s φ 2 . Besides the cubic term, for high enough h s , the zero-temperature term m 4 s log m 2 s in Eq. (14) becomes relevant. A barrier in the effective potential appears at zero temperature, and the size of this barrier grows as h s is increased. This increases the strength of the phase transition and, as a consequence, the φ 4 log φ term becomes even more important. In this case we may have large values of the order parameter φ b /T [47,70], and we shall not use the high-temperature expansion.
For simplicity, we consider λ s = 0. The results are qualitatively similar in the general case (we have checked this) since λ s only affects the value of the thermal mass. We show the values of the velocity and the order parameter in the left panel of Fig. 6 as functions of h s for g s = 2 bosonic degrees of freedom (lower, blue curves) and g s = 12 (upper, black curves). The results are qualitatively very similar for the two cases, only that for smaller d.o.f. a higher coupling is needed to increase the strength of the transition. In both cases the deflagration is unstable for weak phase transitions (φ b /T 0.4), there is a range of stable deflagrations, and for stronger phase transitions (φ b /T 1.3) the velocity is again smaller than v crit . This tendency is already seen in previous cases. There is also a range of values of h s for which we have both φ b /T > 1 and stable deflagrations. In this range we have wall velocities v w ≃ 0.05-0.1, which is good for baryogenesis. We obtained subsonic velocities, but varying the friction by an O(1) factor detonations may appear (see, e.g., Ref. [71]).
Extra fermions strongly coupled to the Higgs field can also make the phase transition strongly first-order [43]. However, strongly coupled fermions may make the vacuum unstable due to the minus sign in front of the term m 4 f log m 2 f in (14). This problem can be solved by adding heavy bosons with the same couplings but with a large φ-independent which is the maximum value consistent with stability. For simplicity, Π b = 0 is assumed. In the right panel of Fig. 6 we show the results for this model as a function of the coupling h f . The velocities, as well as the order parameter, are smaller than in the case of extra scalars. However, the behavior is qualitatively similar. In particular, there is a range of values of h f for which the deflagration is stable (i.e., v w > v crit ). In this range we have v w ∼ 10 −2 . Values of the order parameter φ b /T ≃ 1 occur at the upper limit of this range and, thus, may be compatible with stability.

Consequences of the instability
So far we have only calculated the critical velocity below which the deflagration is unstable, and compared it with the actual value of the wall velocity. We now wish to discuss the possible consequences of the case v w < v crit .
Before doing so, it is worth discussing on the use of the phenomenological form ηvγ for the friction, explicit in the calculation of v w and implicit in the calculation of v crit (or, equivalently, of η crit ). Important deviations from this scaling may occur [31]. In particular, the friction should saturate in the ultra-relativistic regime [30], while this model grows unboundedly as v w → 1 due to the γ factor. We thus may expect significant errors for wall velocities close to the speed of sound or higher. Nevertheless, for small wall velocities (particularly those of interest for baryogenesis) our phenomenological model should give a good description. Indeed, notice that for the two models which admit a direct comparison with the results of [31] (namely, the SM with a low Higgs mass and the SM with a low cutoff, for which the particle content is that of the SM), our results for v w (before taking into account the infrared part of the friction) are in good agreement with those of [31]. In any case, since the scaling ηvγ generally overestimates the friction force, we may expect both v w and v crit to be actually higher than our estimations. Therefore, since deflagrations with v w < v crit are unstable, our estimation of v crit gives a conservative lower bound for the velocity of a stable deflagration.
Let us now assume that we are in the case v w < v crit . Then, perturbations of the wall-fluid system in a range of wavenumbers 0 < k < k c are unstable. In the first place, it is important to determine whether these instabilities will grow in a time comparable to the duration of the phase transition. There is a characteristic scale in the problem (see [18] for details), which is given by d = σ/F dr , where σ is the surface tension and F dr is the driving force. Thus, the deflagration wall is unstable for wavelengths λ which are higher than a critical value λ c ∝ d. The instabilities develop once the bubble reaches a size R b ∼ λ c . The quantity d is generally of the order of the initial bubble size. In contrast, the final bubble size R f is in general much higher than d. As a consequence, there will be instabilities in the range λ c < λ < R f .
To be more precise, in the limit of small supercooling, small latent heat, and small velocity we have λ c ∼ d/[L (1 − v w /v crit )]. Therefore, unless v w is very close to v crit , we have λ c ∼ d/L. AlthoughL may be quite small and, hence, λ c may be a few orders of magnitude larger than d, the final bubble size R f will still be many orders of magnitude larger [18]. This is because R f is related to the cosmological time H −1 (where H is the Hubble rate), while d is given by thermodynamical variables. Roughly, we have R f /d ∼ M P /T , where M P is the Planck mass. For the electroweak phase transition we have M P /T ∼ 10 17 .
Hence, the instabilities generally begin to grow when bubbles are still very small. On the other hand, the instabilities need time to develop. The growth rate Ω in the linear regime is proportional to λ−λ c and to 1−v w /v crit . A simple dynamical analysis shows that (in the approximation T n /T c ≃ 1,L ≪ 1, v w ≪ 1) the instabilities become important when the bubble size reaches the value 18]. Again, even thoughL may be quite small we will have R inst b ≪ R f unless v w is extremely close to v crit . Therefore, the instabilities become important very early in the development of the phase transition.
As already discussed, the instability of the deflagration may spoil the mechanism of electroweak baryogenesis if the wall accelerates and reaches supersonic velocities. In fact, the linear stability analysis breaks down as the perturbations grow. Therefore, Ω only indicates the initial acceleration rate of the unstable wall. It is in principle possible that nonlinear effects stabilize the propagation of the phase transition front. In such a case, one may expect that the wall reaches a new stationary regime, perhaps with a higher velocity. Notice, however, that this new regime will not correspond to a stable weak deflagration, unless the conditions which determine v w change. All other known stationary solutions (namely, Jouguet deflagrations and weak detonations) are supersonic.
There are (at least) two possible alternatives to this situation. Since the instabilities corrugate the wall, the growth of the bubble may be of dendritic type [72]. The hydrodynamics in this case may differ significantly from the case of planar (or spherical) walls, and there may be new stationary propagation modes. Another possibility is that, by the time the instabilities begin to become noticeable, the shock fronts coming from other bubbles already hit the wall. This may happen for small enough v w , since the shock velocity is supersonic. In such a case, the plasma outside the bubble will be reheated from the initial temperature T n to a higher temperature T r . This reheating will tend to decrease both v w and v crit , but the overall effect may be stabilizing. Indeed, the release of latent heat may reheat the plasma up to T c . In such a case the phase transition necessarily slows down [73,74,75]. The effects of such a reheating on electroweak baryogenesis have been investigated in Refs. [76,77,78].
Leaving aside these possibilities, the critical velocity sets a lower bound for the velocity of stationary phase transition fronts. This lower bound for the wall velocity may be used to constrain models of EWB in two ways. On the one hand, if the wall velocity for a given model is calculated and turns out to be below v crit , then in principle the wall will accelerate to velocities which are too high for EWB. This would rule out this model as a baryogenesis scenario. On the other hand, since an accurate calculation of v w is too difficult and model-dependent, using the value of v crit as a lower bound may be useful. This lower bound may constrain the baryogenesis mechanism if v crit turns out to be too large.
Interestingly, besides the potentially negative consequences for baryogenesis, the instability of slow walls may give rise to new cosmological consequences for weak phase transitions. Indeed, the instabilities lead to acceleration of the walls, anisotropic growth of bubbles and turbulence of the plasma. These effects may give origin to magnetic fields [24] or gravitational waves [18].

Conclusions
It is known that the stationary motion of phase transition fronts may be unstable for subsonic velocities. Specifically, the deflagration front is hydrodynamically unstable below a critical velocity v crit . In order to discuss the implications of this fact for the electroweak phase transition, we have computed both the wall velocity and the critical velocity for several extensions of the Standard Model. We have considered a significant variety of models and a wide range of parameters, including those which are favorable for baryogenesis.
The general result is that the deflagration tends to be unstable for very weak or very strong phase transitions, while for phase transitions with an order parameter φ b /T ≈ 1 the deflagrations are generally stable. In general, for any model there is a range of parameters for which we have stable deflagrations.
The stability condition v w > v crit constitutes in principle a restriction for electroweak baryogenesis (assuming that, if it is not fulfilled, the wall will accelerate to supersonic velocities). This condition, combined with that of avoiding the washout of the generated BAU (i.e., φ b /T 1), the requirement of enough CP violation, and experimental bounds, may restrict significantly some models. We remark that, even if we have a subsonic wall with v w > v crit , a departure of the wall velocity from the range 10 −2 < v w < 10 −1 , where the BAU has its maximum, may also prevent a quantitatively successful EWB. In all the models we considered we found a region of parameter space for which the deflagration is stable and the above conditions for φ b /T and v w are fulfilled.
We stress that the friction is very model-dependent, and current calculations have large errors which propagate to the wall velocity. For velocities close to the speed of sound, such O(1) factors may determine whether the wall propagates as a deflagration or a detonation. Similarly, for wall velocities close to the critical value, O(1) variations of the friction will determine the stability or instability of the deflagration. The calculation of the critical velocity is not easy either. We have used the analytic results of Ref. [18], which were obtained using several approximations, such as the bag EOS and the assumption of planar walls. Nevertheless, the value of v crit is less dependent on details of the specific model (i.e., it depends only on thermodynamical parameters). In the appendix we give a simple fit for v crit as a function of the parameters L/w + (T c ) and T n /T c .
We also remark that we have used for our calculations a simple phenomenological model for the friction force. We have argued that this approximation should be good at least for the case of small wall velocities which are required for baryogenesis. The actual value of v crit is possibly higher than our result. Hence, our calculation gives only a conservative lower bound for the velocity of a stable deflagration.
The nucleation rate vanishes at T = T c and grows rapidly as T descends below T c . At a temperature T < T c , the probability of finding a bubble in a causal volume V c ∼ (2/H) 3 is given by The time-temperature relation is given by dT /dt = −HT , where the expansion rate is given by H = 8πGρ u (T )/3, with G the Newton's constant. We define as usual the nucleation temperature T n at which nucleation effectively begins by the condition P (T n ) = 1. We compute T n numerically from Eqs. (31 -33).

A.2 Hydrodynamics
The equation for the wall can be obtained from the phenomenological equation for φ, which is similar to Eq. (4), with the last term replaced by a phenomenological term proportional to u µ ∂ µ φ. The usual procedure is to consider the equation in the reference frame of the planar wall, multiply by ∂φ/∂z, where z is the coordinate perpendicular to the wall, and then integrate across the wall, taking into account the variation of the fluid variables. To perform the integration, either an ansatz for the wall profile is used or some approximations are needed. We shall use the result of [18], which uses linear approximations for the variations of quantities inside the wall. We have We also use the bag EOS, where L is the latent heat and a is an effective radiation constant depending on the number of relativistic d.o.f. For the bag EOS, Eq. (34) becomes The temperature T b behind the wall is related to the temperature T u in front of it through the equations for the discontinuity, Eqs. (5). For the bag EOS we have whereL = L/(4aT 4 c /3) = L/w u (T c ). The velocities v u and v b are related by [81] v u = 1 with α u = L/ (4aT 4 u ) = (L/3)(T c /T u ) 4 .
We see that we have two kinds of solutions. Indeed, the + sign in front of the square root in Eq. (38) corresponds to detonations while the − sign corresponds to deflagrations. For detonations we have v b < v u while for deflagrations we have v b > v u . These solutions are further classified into weak, strong or Jouguet, depending on the relation of the outgoing velocity v b with the speed of sound c s = 1/ √ 3. For weak deflagrations we have v b < c s .
The fluid profiles away from the wall are quite simple for a planar interface. For a weak deflagration we can only have constant fluid velocities. In the reference frame of the bubble center, the velocity profile is the following. Inside the bubble the fluid is at rest (in the reference frame of the wall, this gives the condition v b = v w ). In front of the wall we have a constant fluid velocityṽ u > 0 (in the wall frame, this inequality transforms into the deflagration relation v u < v b ). The boundary condition of an unperturbed fluid far in front of the wall is fulfilled by a jump of the fluid velocity fromṽ u to 0 at a certain point. Such a discontinuity in the fluid profile is called a shock front.
The shock discontinuity is also determined by applying Eqs. (5) to the shock front, where the EOS is now the same in both sides of the interface. It turns out that the shock front moves supersonically. The region between the bubble wall and the shock front is reheated, i.e., the temperature T u is higher than the nucleation temperature T n . We have Using the relations (37)(38), the two equations (36) and (39) can be readily solved to obtain v w as a function of η and T n .

B The critical velocity
Using the relations (37)(38)(39) in Eq. (9), we obtain the equation for the critical velocity v crit as a function of T n /T c andL. The solution for v crit is plotted in Fig. 7 (solid lines) as a function ofL for different values of T n /T c . Technically, solving Eq. (9) becomes difficult in some cases due to the divergence of the quantity β as the product T u T b approaches T 2 c . This happens for instance for high values of the latent heat, since the reheating in front of the wall may cause T u to exceed T c . As can be seen from Eq. (36), the exact limit cannot be reached and we always have T u T b < T 2 c . As a consequence, there is always a critical velocity.
In the case of small supercooling and small latent heat, the wall velocity is generally small, since we have ηv w ≃ L(1 − T u /T c ). In this limit we also have T b ≃ T u and v b ≃ v u , and we have a simple expression for v crit as a function of T u [18], On the other hand, the temperature T u of the reheated plasma in front of the wall depends on v w . In the present approximation we have T u /T c ≃ T n /T c + (L/ √ 3v w ). Hence, for This simple approximation underestimates the value of the critical velocity. It becomes exact in the limit T n → T c , but only gives a good estimation for small amounts of supercooling. For instance, for T n /T c > 0.98 the error is less than a 5%, while for T n /T c = 0.9 the actual value of the critical velocity is about a 15% higher than the approximation (41). Nevertheless, a slight modification of Eq. (41) provides a good fit, v crit ≃L (1 +L) 2 √ 3 1 + 12 with ε = 0.5(1 − T n /T c ) and a = 1 − 2.7(1 − T n /T c ). The error of this approximation is less than a 5% for 0.75 < T n /T c < 1 and in all the range 0 <L < 1 (see Fig. 7).