Wall speed and shape in singlet-assisted strong electroweak phase transitions

Models with singlet fields coupling to the Higgs can enable a strongly first order electroweak phase transition, of interest for baryogenesis and gravity waves. We improve on previous attempts to self-consistently solve for the bubble wall properties--wall speed $v_w$ and shape--in a highly predictive class of models with $Z_2$-symmetric singlet potentials. A new algorithm is implemented to determine $v_w$ and the wall profiles throughout the singlet parameter space in the case of subsonic walls, focusing on models with strong enough phase transitions to satisfy the sphaleron washout constraint for electroweak baryogenesis. We find speeds as low as $v_w \cong 0.22$ in our scan over parameter space, and the singlet must be relatively light to have a subsonic wall, $m_s \lesssim 110$ GeV.


I. INTRODUCTION
The electroweak phase transition (EWPT) in the early universe has been intensively studied as a possible source for the cosmic baryon symmetry and gravitational waves. Within the standard model (SM) neither of these interesting outcomes are possible, given the known mass of the Higgs boson, because the phase transition is a smooth crossover [1,2], whereas a first-order EWPT is required for electroweak baryogenesis (EWBG) and production of observable gravity waves (for reviews, see for example refs. [3,4]). New physics, typically in the form of scalar fields coupling to the Higgs boson, can however lead to a first-order transition, with consequent nucleation of bubbles of the true (electroweak symmetry broken) vacuum, at the onset of the transition.
In order to make quantitative predictions for either baryogenesis or gravitational wave production in a given model, it is necessary to understand the detailed properties of the phase transition bubbles, especially the shape of the bubble walls (typically modeled as a tanh with some thickness L w ) and the terminal velocity v w attained by them, once the forces of internal pressure and external friction from the plasma have balanced each other. This calculation, first carried out for the SM in refs. [5][6][7] (assuming a light Higgs boson), and in the Minimal Supersymmetric Standard Model, where the phase transition is enhanced by light stops, in refs. [8,9], * avi.friedlander@queensu.ca † banta@physics.ucsb.edu ‡ jcline@physics.mcgill.ca § dtuckers@williams.edu turns out to be quite challenging because the frictional force, which requires solving the Boltzmann equations for the perturbations of the plasma caused by the wall, depends on the same wall properties that one is trying to determine. A self-consistent procedure to solve this system is numerically expensive, and for this reason many studies of EWBG or gravitational wave production leave L w and v w as phenomenological parameters that can be freely varied, or in a somewhat better approximation, calculated by modeling the friction in a phenomenological way [10][11][12]. However to assess the prospects for a specfic model to yield interesting results, one must eventually carry out the actual computation of L w and v w with the actual friction term derived from the fluid perturbations. Accurate estimates of these parameters are needed to make quantitative predictions for baryogenesis or gravity wave production.
The procedure becomes even more laborious in the case where an extra singlet field couples to the Higgs, in order to facilitate the first order transition, and also gets a vacuum expectation value (VEV) in the bubble wall [13,14]. In that case one must solve for both field profiles, which has been attempted in refs. [15,16], subject to some limiting approximations. In particular, these previous works assumed that the bubble wall shapes are described by tanh profiles. In reality, the Higgs field and the singlet can have shapes that differ from such an assumption, and it is not obvious how strongly this affects the determination of v w . One of our main purposes is to overcome this limitation by developing an algorithm to determine the actual wall profiles along with v w .
Moreover previous studies of singlet-assisted strong EWPTs have focused on a few benchmark models. In the present work we make a comprehen-sive scan of the parameter space for a class of models, where the singlet potential has the Z 2 symmetry s → −s, and the singlet VEV disappears at low temperataures. This choice has the virtue of simplicity, being characterized by three parameters, the singlet mass m s , its cross coupling λ hs to the Higgs, and the VEV w 0 of s in the false vacuum where h = 0. The barrier between the true and false vacua provided by the λ hs h 2 s 2 interaction is already present at tree level, and is what enables the phase transition to be strongly first order [13,14]. Moreover with s = 0 at T = 0, the new sources of CP violation needed for EWBG are not overly constrained by experimental limits on electric dipole moments.
The paper is organized as follows. Section II describes the singlet scalar model used throughout the paper. Section III outlines the main features of the electroweak phase transition dynamics that will be studied in detail in the following. In section IV the methodology for determining the wall dynamics, including its velocity, are described; the results of those calculations are presented in section V. Conclusions are given in section VI. Appendices contain details concerning the finite-temperature effective potential (appendix A), diffusion equations used to determine the fluid perturbations (appendix B), and exceptional models that have peculiar features in their potentials (appendix C).

II. THE MODEL
A simple extension of the SM is the addition of a scalar singlet s that couples only to the Higgs field, and has the Z 2 symmetry s → −s. Its zerotemperature tree level potential is given by where H is the Standard Model Higgs doublet, and λ h , v 0 are the Higgs self-coupling and VEV respectively. There are three new parameters λ s , w 0 , and λ hs , that describe the singlet's self coupling, its VEV when in the false minimum where H = 0, and the coupling between H and s. There is no loss of generality by omitting a separate m 2 0 s 2 mass term. The physical singlet mass in the electroweak broken vacuum is given by We restrict the parameters so that m 2 s > 0, implying that s = 0 in the true vacuum. The Higgs doublet components are where h denotes the background Higgs field, and the χ's are the Goldstone bosons. The full effective potential takes into account oneloop corrections and temperature effects, where V 0 is the tree-level potential (1), V 1 is the one-loop correction, V CT contains the counterterms associated with V 1 , and V T is the thermal contribution, including ring resummation of thermal masses. These expressions are standard, and we have described them in detail in appendix A. V eff is determined by the measured SM parameters and the three new ones, that we henceforth take to be w 0 , λ hs and the singlet mass m s by trading λ s for m s through eq. (2).

III. PHASE TRANSITION
Because of our assumption that s = 0 at low temperatures, the λ hs h 2 s 2 coupling creates a barrier in field space between the false and true vacua, and gives rise to a two-step phase transition. At temperatures T w, the minimum of the potential is at the origin, where both electroweak symmetry and the Z 2 symmetry are restored. At some temperature T , the first transition occurs, where s → w (see fig. 1). As T decreases, s and the corresponding minimum of the potential becomes metastable. At the critical temperature T c , the two minima become degenerate, and at a slightly lower temperature T n nucleation of bubbles begins, signaling the second transition where electroweak symmetry is broken while the Z 2 symmetry is restored. These two transitions are summarized as It is the second transition that is important for baryogenesis. We note that while domain walls could form during the first transition, the restoration of the Z 2 symmetry during the second transition will cause them to annihilate. This occurs long before they can dominate the energy density of the universe, hence we expect that no cosmological problems will arise from this brief appearance of domain walls [13].
The dynamics of the phase transition depend strongly on T n , the temperature at which the probability of a bubble nucleating within one Hubble volume per Hubble time is O(1) [17]. The tunneling probability goes as exp(−S 3 /T ), where S 3 is the three-dimensional Euclidean action. We used Cos-moTransitions [18] to find phase transition candidates and to determine T n 1 . The criterion for the nucleation temperature is taken to be S 3 /T n = 140 [17].
Since our investigation is motivated by electroweak baryogenesis, we focus attention on first order transitions that are strong enough to preserve the baryon asymmetry from washout by residual sphaleron interaction inside the bubbles. This requires the Higgs VEV at the nucleation temperature to satisfy [19] v n T n 1.1 .

IV. WALL DYNAMICS
The bubble-wall dynamics are determined by the interactions of the Higgs and singlet fields with a thermal fluid consisting of top quarks, electroweak gauge bosons, and any other particles to which the scalars couple significantly. After the bubble nucleates, it expands due to the outward pressure caused by the potential difference between the phases of the scalar fields on either side. The interactions of the wall with the surrounding fluid counteract the expansion by a friction force that depends on the speed and shape of the wall. 1 CosmoTransitions occasionally fails to find transitions when they should exist; in such cases, changing the value of λ hs by O(10 −6 ) can overcome the problem. Moreover, Cos-moTransitions often reports more phase transitions than expected for a given model; we find that defining the EWPT as the most recent first order transition where the Higgs' VEV in the unbroken phase is smaller than the nucleation temperature gives correct results. Of those, only phase transitions that ended with no singlet VEV were studied. We found it a useful cross-check to require the transition identified by CosmoTransitions to have a critical temperature that matched our own calculations for the parameter point in question.
If the friction is strong enough, the bubbles reach a steady-state velocity whose value is relevant for gravitational waves and baryogenesis. The terminal v w depends upon the field profiles that solve the equations of motion. These in turn depend upon the temperature T w of the wall and the v w -dependent friction exerted by the plasma on the wall, leading to T w > T n , due to heating by the fluid. A self-consistent solution thus requires simultaneously solving for v w and the scalar field profiles in the wall.

IV.1. Deflagration profiles
The wall temperature, and the rest of the dynamics of the bubble, depend on whether the phase transition proceeds through deflagrations or detonations (hybrids of these two are also possible). For subsonic bubble walls, with v w < 1/ √ 3, the bubbles grow via deflagrations [10], in which the wall is preceded by a shock front that moves through the fluid, perturbing it, increasing the temperature from T n to T s and causing the wall to move (see fig. 2). The fluid velocity decreases until the point where the wall passes it, so that the fluid behind the wall is at rest relative to that preceding the shock front. In this work we limit our investigation to the case of deflagrations, hence subsonic walls, deferring the study of supersonic walls to the future [20].
Because of the heating, the bubble wall dynamics are not determined at temperature T n , but rather the temperature of the fluid near the wall. The calculation is performed at times sufficiently long after nucleation that the bubble has reached a steadystate velocity, and the profiles of the fluid perturbations vary on scales much larger than the wall thickness. Therefore one approximates the wall as a discontinuity, such that the fluid temperature is T + (T − ) just in front of (behind) the wall. The fluid velocity likewise is discontinuous there.
Since it is often convenient to switch between reference frames, we adopt the notation v xy for the velocity of x in the reference frame y. In this context, x and y either refer to the wall (w), shock front (s), or the fluid at position 1 (in front of the wall), 2 (behind the shock front), or u (the unperturbed "universe" frame, in front of the shock front or behind the wall). The wall velocity, which is measured with respect to the fluid directly in front of the wall, is v w1 in this notation. We note that for any x and y, v xy = −v yx . A diagram depicting the geometry and labels is shown in fig. 2.
The relationships between the various fluid velocities and temperatures are found by integrating the stress tensor T µν across either of the two interfaces shown in fig. 2. Approximating the fluid as perfect, these depend only on the fluid density and pressure. The equations of state can be expressed as [16] where ρ ± is the fluid density on either side of the bubble wall and p ± is the pressure. a ± and ± are given by F ± is the free energy of the fluid evaluated at the respective VEVs outside and inside wall and T = T n . It is given by Here g * = 107.75 − 24.5 = 83.25 is the effective number of degrees of freedom, apart from the t, W/Z, h, χ and s, whose contributions are already included in V eff . In general, a ± (T ) and ± (T ) in (6,7) should be evaluated at the temperatures T ± , but for transitions typically of interest for baryogenesis, where there is a limited degree of supercooling, the T -dependence of a(T ) and (T ) is insignificant.
Integrating T µν across the bubble wall provides the relations [10] v Integrating across the shock front, the temperature changes, but not the the field values, leading to Fluid velocities in the wall and shock wave frame are related by Lorentz transforming to the u frame using The relationship between the fluid velocity and temperature behind the shock front and in front of the wall can be approximated by linearizing the stress-energy equations with respect to small fluid velocities in the universe frame [11], to obtain For a given guess of the wall velocity, v w ≡ v w1 , and an associated nucleation temperature, eqs. (11)(12)(13)(14)(15) can be used to solve for the eight remaining variables: , v 2u , and v s2 . An example of the solution to these equations for sample set of parameters is shown in fig. 3.
3. An example of how the wall temperature changes as a function of vw for a sample model with ms = 63 GeV, w0 = 130 GeV, and λ hs = 0.9

IV.2. Equations of Motion
The equation of motion of a scalar field coupled to a perfect fluid has been derived by enforcing the conservation of the stress-energy tensor in the WKB approximation [7], or starting from the Kadanoff-Baym equations [15]. Both methods lead to is the zero-temperature effective potential, the sum is over all particles that couple to φ, n i is the number of degrees of freedom of particle i, m i is its field-dependent mass, and f i ( p, x) is its phase space distribution. By separating f i = f 0,i + δf i into equilibrium and out-of-equilibrium components, the equation of motion takes a more useful form. The integral over f 0,i is equivalent to accounting for the T -dependence of the effective potential, giving (17) describes the friction force that comes from the dissipative interactions between the scalar field and the surrounding fluid.
In the following, we assume that the dominant sources of friction are the top quark and electroweak gauge bosons. The lighter fermions, gluons and photons can be safely ignored because of their negligible couplings to the Higgs field. The self-couplings and mixing of the scalar fields are assumed to be subdominant due to their fewer degrees of freedom relative to vector bosons or quarks. It is possible that including these contributions would lead to moder-ately slower walls, which could be advantageous for baryogenesis, but we defer this issue to future study.
For the electroweak phase transition considered here, there are two relevant scalar fields, each with its own equation of motion, that must be simultaneously solved. Eq. (17) can be further simplified by accounting for the spherical symmetry of the wall and going to the planar limit, which reduces the system to one spatial dimension, and by considering only the steady-state regime. Therefore, the equations of motion for a bubble wall traveling in the negative z direction become where primes denote derivatives with respect to z. Strictly speaking, the existence of a steady-state solution to the equations of motion does not guarantee that those solutions will in fact be realized in the physical setting, but this issue is beyond the scope of this paper.

IV.3. Friction
The friction experienced by the wall depends on δf i ( p, z), the deviation from equilibrium of W/Z and t. We adopt the fluid approximation framework developed in [7], in which the friction is fully described by three fluids: that of the top quark, the massive gauge bosons, and the other particles, denoted as the "background." We label the gauge boson contribution by W although it also includes Z. For simplicity W and Z are grouped together due to their similar couplings, and assigned a mass-squared that is the weighted average of m 2 W and m 2 Z . The background fluid encompasses all the fields that are assumed to contribute negligible friction, but which nevertheless play an important role in the wall dynamics. We consider friction only from fluid excitations with large momentum, such that the wavelength is shorter than the width of the wall. It has been shown that IR excitations in the massive gauge boson fluid can be important [21], but we have checked numerically that these are subdominant for parameters of interest in the present study, using the same approximations to evaluate the IR contributions as in ref. [16], but taking care to impose the perturbative cutoff m W (z) > g 2 T [21].
The phase space distribution for the t and W flu-ids can be parametrized as where the +/− is for fermions/bosons and accounts for perturbations in the fluids. δµ i (z), δτ i (z), δv i (z) are respectively the perturbations in the chemical potential, the relative temperature and the velocity. The subscript bg denotes the background fluid. All the perturbations are relative to the fluid directly in front of the wall where µ = 0 and T = T + as described in section IV.1. Deviations from equilibrium in the fluids are governed by the Boltzmann equation Rather than solving the full Boltzmann equation, one linearizes it in δ i (z), and converts it to a system of ordinary differential equations by taking three moments: 3 . A detailed derivation is provided in appendix B, leading to the coupled matrix equations while the source terms are The coefficients c i j and d i j denote the integrals and where f 0,i is the equilibrium distribution function for particle i. In previous literature (with the exception of ref. [22]), these coefficients were evaluated in the massless approximation, where d i j = c i j , setting m i /T = 0. But for some phase transitions satisfying the sphaleron bound (5), m t /T > 1 in the broken phase; hence we we derived the full mass dependence for this work. For the background fluid, one can show that where the background fluids are approximated as massless [15]. The Γ i matrices in eqs. (23,24) quantify the interactions of the fluids. Γ t and Γ W take the form where the matrix elements are numerical constants calculated in [7] to be 2 Using energy-momentum conservation, the background fluid collision terms are given by where n t = 12 and n W = 9 are the number of degrees of freedom in the respective components. The background fluid is assumed to be in chemical equilibrium, implying that δµ bg = 0. This assumption removes the top row of eq. (25). The remaining two rows determine δτ bg and δv bg in terms of q W and q t , where A −1 bg23 denotes the matrix where the bottom right block of A bg is inverted and the rest of the matrix elements are zero.
Equations (23) and (24) can then be expressed in 6 × 6 matrix form, in the rest frame of the bubble wall, with and The W and t fluid perturbations are determined by solving eq. (35) using the relaxation method as described in ref. [23], since shooting tends to be unstable. The background fluid perturbations are found by integrating eq. (34). One can carry out this procedure for given values of the wall velocity and shape, and from the ensuing perturbations compute the friction term in the Higgs field equation of motion (18) using An example of the solutions for the perturbations is shown in fig. 4.

IV.4. Solving the Equations of Motion
With the friction calculated in eq. (38), the equations of motion that must be solved to determine v w and the shape of the wall are Deep into the bubble interior, eq. (39) is not exactly satisfied once we adopt our approximation schemes for calculating the effective potential and the perturbations. The Higgs' VEV is unchanging there, so the kinetic term is zero. Similarly the perturbations in the W and t fluids go to zero on both sides of the wall. This implies that the terms proportional to δτ bg must exactly cancel out the potential term. When the perturbations are determined as described above and the potential term is calculated with the Higgs VEV that minimizes the potential inside the bubble, the two terms do not cancel as they should. This is due to differences in the derivation of the friction terms in comparison to the effective potential. Firstly, the fluid perturbations are only determined to linear order whereas the temperatures that go into the effective potential, T + and T − , were calculated including non-linearities in the fluid equations. This means that while in theory T + − T − = T + δτ bg , their relationship is only approximate. The other cause is that the scalar fields were treated as massless background fields in the friction calculation but their full contribution was included in the effective potential. There are three ways to account for this inconsistency: the Higgs VEV inside the bubble can be chosen not to minimize the potential but instead to cancel the friction term, the entire friction can be scaled to cancel the potential term but maintaining the friction shape in z, or just the background perturbation contribution to the friction can be scaled to cancel the potential term. We adopt the last option, which we found to be the most conservative choice (leading to slightly larger wall velocities). The equations of motion that we actually use to determine the wall dynamics then become where y is an O(1) parameter chosen so that the equations are satisfied for larger positive values of z.
For a given value of v w , the relaxation method can be used to find the shapes of h(z) and s(z) that come closest to solving the equations of motion. One must then vary v w and find a complete solution to the equations, by iterating this procedure. A reasonable initial guess for both v w and the wall shape is required, leading us to solve the equations in two stages. The first part is to guess v w and the wall shape using the tanh ansatz employed in previous studies of wall velocities [7], [15], [16]. The second uses these as a starting point to numerically determine v w and the wall shapes.
The tanh ansatz in the first stage assumes that the Higgs profile has the form where v(T − ) is the Higgs VEV at temperature T − and L w is the width of the wall. The friction and shape of the singlet profile are independent of each other, so there is no need to impose a tanh ansatz for s; rather its profile is found by numerically solving its equation of motion. This reduces the problem to finding values of v w and L w that come closest to solving the Higgs equation of motion 3 . No choice of v w and L w will exactly solve eq. (41), since the true shape is not a tanh function. Instead we follow ref. [15] by calculating two moments of E h in eq. (41) and finding the values of v w and L w that make them vanish. The two moments are taken to be since with this choice the Jacobian matrix ∂(E 1 , E 2 )/∂(v w , L w ) is always far from being singular.
The first stage of the algorithm can then be summarized as: In the second stage, we aimed to relax the tanh profile assumption for h(z) and to determine its shape more exactly. Using the values of v w and L w from the first part as new initial guesses, we solved both h and s equations of motion simultaneously, using relaxation. A challenge here is that the friction on the wall, which is expensive to compute, depends on the background h(z) solution. To speed up the algorithm, we recomputed the friction only after several relaxation steps. This procedure leads to eventual convergence, unless the initial guess for v w is too poor. Convergence was tested by seeing how closely the two equations of motion were satisfied, using the squared error statistic The best value of v w was determined by varying v w in the region of the guess from step 1 as to minimize E tot . An example of the wall shapes that solve the equations of motion is given in fig. 5. It demonstrates that the actual profiles can differ significantly from the tanh ansatz.

V. RESULTS AND DISCUSSION
A scan of the parameter space of the scalar singlet model was performed in the ranges We imposed the lower bound m s > m h /2 so that collider constraints from invisible Higgs decay (h → ss) do not apply [24]. This is a mild restriction, since for m s < m h /2 and not too close to the upper limit, these constraints imply the invisible branching ratio is 25%, hence λ hs 0.01, which is too small to give rise to a strong phase transition.

V.1. Wall Velocity Results
Our determinations of the wall speed over the full parameter space are illustrated in in fig. 6, showing contours of v w in the plane of m s versus λ hs , for a series of w 0 values. The grey hatched regions indicate parameters for which no transitions with subsonic walls were found. One can see that models with heavier singlets and larger λ hs couplings tend to produce faster-moving walls. Generally we find a minimum value for v w , which depends on w 0 and is smallest for w 0 ∼ 120 GeV, where the lowest speed v w ∼ = 0.22 is found. The parameters specifying a few benchmark models and their resulting phase transition properties are shown in table I.
Since it is numerically expensive to compute v w for a given model from first principles, it is useful to look for relations between it and other quantities characterizing the strength of the phase transition, that are easier to compute. In fact we observe a strong correlation between v w and the double ratio where v/T is evaluated respectively at the nucleation and the critical temperatures. This is a measure of the degree of supercooling, and its correlation with v w is plotted in fig. 7, showing that v w increases rapidly with r − 1. We find an analytic fit v w ∼ = 0.53 (1 − 0.77 r −34 ), with deviations of order ±0.03. The maximum value of r found for subsonic bubble walls was r ∼ = 1.15. It remains close to unity even for strong transitions, validating the assumption made in section IV.1 that the equations of state at the nucleation and critical temperatures do not differ significantly from each other. The fact that a cutoff on r exists, above which it is unlikely to produce subsonic walls, can be seen in fig.  8, which shows all the models tested, including those  found not to have slow bubble walls. It clearly shows that for r 1.1, few transitions produce subsonic walls, whereas below that point all the models tested were found to do so. Fig. 6 shows that subsonic walls require the singlet to be relatively light, m s 110 GeV, often with a relatively large coupling to the Higgs, λ hs ∼ 1. If s is long-lived enough to escape detection within a collider, Ref. [25] suggests that a singlet with these properties may be a realistic target for the highluminosity LHC (in the MET plus forward jets channel, from vector-boson fusion production of an offshell Higgs). On the other hand, if we take the model at face value, as a complete model with a standard thermal history, Ref. [25] also finds that the LUX direct detection experiment [26] rules out m s 120 GeV even though s would make a subdominant contribution to the dark matter. Of course, additional model ingredients can easily make s unstable on cosmological time scales without affecting our phasetransition and wall-velocity results.

V.2. Wall Shape Results
Although Fig. 5 shows that the wall shapes deviate from a tanh profile, it is nevertheless a useful approximation for concisely encoding information about the wall shapes. We have accordingly analyzed our results from the fully numerical algorithm to find the best-fit tanh profiles, including a possible offset δ z between the Higgs and the singlet profiles: where we have allowed for independent widths L h and L s of the Higgs and singlet profiles.
To display results for the wall thicknesses, we have opted to use dimensionless combinations like L h T + , where T + is the temperature of the wall. If one wants to translate these into absolute thicknesses, it can be done using the strong correlation between L h T + and L h in GeV −1 units, shown in Fig. 9. Since all models with subsonic walls have nucleation temperatures in the range 80 GeV ≤ T n ≤ 140 GeV (see Fig. 8), and for slow walls the wall temperature does not deviate much from the nucleation temperature, the relationship between these two ways of characterizing L h is linear with relatively little scatter: L h T + ∼ = L h 120 GeV. This reflects the fact that the deviations of wall temperature from the mean value T + = 120 GeV are relatively small.
Contour plots of L h T + , L s T + , and δ z T + similar to those for v w are presented in Figs. 10 and 11. We find that faster walls tend to be thinner and have smaller offsets. These relationships are plotted in Figs. 12-14, which show strong correlations, especially in the case of L h . With rare exceptions, L s < L h , with L s typically smaller than L h by 20-30%.
As discussed in Appendix C, a small number of points were found to have extra potential minima or plateaus in the middle of the wall. The results for these points are likely to change with a more careful treatment of the effective potential. However, these points are rare, and including them does not change any of our conclusions.

VI. CONCLUSION
This work has laid out a more quantitative methodology than has been previously used, for calculating the wall velocity of bubbles during the electroweak phase transition with an additional scalar field. We improved on previous similar studies by solving for the actual profiles of the scalar fields, rather than just parametrizing them using a tanh ansatz. Other improvements made here include use of the one-loop Coleman-Weinberg contributions to the potential including the effect of thermal masses, accounting for the sphericity of the bubbles, not expanding fluid perturbations to first order in m/T , and performing a scan over the three-dimensional parameter space.
Scanning over the parameter space reveals that the scalar singlet model is able to produce slow bubble walls that are preferable for electroweak baryogensis to occur, down to a minimum wall velocity of v w ∼ = 0. 22. These examples of slow-moving walls only occur in phase transitions with small amounts of supercooling.
There are a few ways in which this study can be extended by future work. The precision of the wall velocity calculation could be improved by including additional sources of friction such as from the scalar fields and IR gauge boson modes. matter anti-matter asymmetry that would be produced. Lastly, extending this analysis to apply for faster walls or even supersonic walls could be of interest for studying other effects of the phase transition such as gravitational waves. We are currently studying these issues [20]. ACKNOWLEDGMENTS We thank B. Laurent for useful discussions and comments on the draft. JC and DTS thank the Aspen Center for Physics for providing a stimulating environment where this work was initiated. The computations in this work were run using equipment funded by the Canada Foundation for Innovation and supported by the Centre for Advanced Computing at Queen's University. Besides the previously mentioned CosmoTransitions, the code used for the calculations utilised Eigen [27] and the GNU Scientific Library [28]. AF and JC are supported by NSERC (Natural Sciences and Engineering Research Council, Canada).

Appendix A: Effective potential
The one-loop contribution to the potential can be approximated as where n i is the number of degrees of freedom of each particle. For the scalar fields, longitudinal W/Z and top quark c i = 3/2 but for the transverse gauge bosons c i = 1/2, in the MS scheme. The top quark is the only fermion included in the sum since the contributions from lighter fermions are suppressed by their small Yukawa couplings. χ stands for the Goldstone boson contributions.
The one-loop contribution acquires a temperature dependence through the thermal masses of the particles, in this method of carrying out the ring resummation [29]. It has been shown that for sufficiently strong phase transitions, a more careful treatment of thermal masses can be important [30].
The scalar masses in eq. (A1) are given by the eigenvalues of the mass matrix: where φ i and φ j are the five scalar fields summed over in eq. (A1)and The three mass eigenvalues associated with the Goldstone bosons vanish in the vacuum state making those terms in eq. (A1) formally divergent. This is properly dealt with by introducing a scale coinciding with the Higgs mass, m h , to cut off the IR divergence [31].
The masses associated with the longitudinal modes of the gauge bosons in eq. (A1) are given by the eigenvalues of the mass matrix: The rest of the field-dependent masses in eq. (A1) are given by: The counterterm contribution to the potential can be parameterized as The five counterterms were chosen to ensure that the full effective potential at T = 0 maintains its treelevel values for the scalar masses, potential minima, and scalar mixing. This is done by imposing the following conditions at T = 0: ∂V ∂h h=v0,s=0 = ∂V ∂s h=0,s=w0 = 0 (A9) and where m s = 1 2 λ hs v 2 0 − λ s w 2 0 is the mass of the scalar singlet in the true vacuum.
The resulting counterterm parameters are found to be and Lastly, the temperature dependence of the potential is given by where J F and J B are functions which describe fermions and bosons temperature-dependent contribution to the one-loop potential. The functions are calculated from and These equations fully describe the one-loop potential of the scalar fields.

Appendix B: Linearized Boltzmann Equations
The following derivation of the linearized moments to the Boltzmann equation, which are used to to determine the friction of the equation of motion, follows closely to that originally expressed in [7]. The difference between that derivation and the one here is that the full dependence of m/T is included here instead of expanding to lowest order. This allows for stronger phase transitions to be quantitatively studied.
As noted in eqs. (20)(21)(22) the fluids are described by the distribution function where the +/− is for fermions/bosons and The background fluid is in chemical equillibrium so for the rest of the derivation δµ bg = 0 Deviations from equilibrium in the fluids are governed by the Boltzmann equation The left side of eq. (B3) can be expanded as In the fluid's reference frame Starting with the last term As will be shown the perturbations are sourced by a term proportional to This may raise the concern that if m i /T is not small and δ i ∝ , does the linear approximation break down? The tanh ansatz can be used to set a rough condition on the relation between v n /T n and LT under which taking the linear order is valid. That will be derived at the end of this section.
Next one observes that ∂ t = v w ∂ z in the fluid's reference frame, so to linear order in the perturbations Going back to eq. (B4), the term independent of δ i acts as the source term in the perturbations equations.
Therefore the Boltzmann equation becomes  This condition is easily met by all the wall found to have subsonic walls in this paper therefore indicating that the linear order approximation is valid.

Appendix C: Questionable Transitions
Some of the models included in this analysis have features in their potential that are likely to be sensitive to higher-order corrections. The two types of features seen are additional local minima and plateaus in the potential in the middle of the bubble wall. An example of a model with such a minimum can be seen in fig. 15 which shows a potential with an additional local minimum in the h-s plane and spatially along the wall respectively. Similarly, an example of a wall profile with a potential plateau can be seen in fig. 16. These features appear in the potential when the bosonic squared masses become negative. This phenomenon, was first discussed in [31], where it was pointed out that this is sensitive to IR corrections to the potential. Including the thermal masses in the one-loop calculation of the effective potential reduced the frequency that models with these features appear. Due to the fact they are sensitive to IR corrections, these models are more likely to change from a more complete treatments of thermal resummation as discussed in [30] and including higher order loops in the effective potential. While that may make these specific models less robust, removing them does not change any of the trends or conclusions discussed in this paper and therefore even if more complete calculations change these models, the overall results should hold. A final observation relating to these points is that they tend to have lower wall velocities and shapes that deviate more strongly from a tanh profile, as seen in fig. 17. If these minima are indeed physical, or if another physics model can produce similar physical effects, then removing the tanh ansatz is even more important for capturing the true wall behaviour.