Baryogenesis and gravity waves from a UV-completed electroweak phase transition

We study gravity wave production and baryogenesis at the electroweak phase transition, in a real singlet scalar extension of the Standard Model, including vector-like top partners to generate the CP violation needed for electroweak baryogenesis (EWBG). The singlet makes the phase transition strongly first-order through its coupling to the Higgs boson, and it spontaneously breaks CP invariance through a dimension-5 contribution to the top quark mass term, generated by integrating out the heavy top quark partners. We improve on previous studies by incorporating updated transport equations, compatible with large bubble wall velocities. The wall speed and thickness are computed directly from the microphysical parameters rather than treating them as free parameters, allowing for a first-principles computation of the baryon asymmetry. The size of the CP-violating dimension-5 operator needed for EWBG is constrained by collider, electroweak precision, and renormalization group running constraints. We identify regions of parameter space that can produce the observed baryon asymmetry or observable gravitational (GW) wave signals. Contrary to standard lore, we find that for strong deflagrations, the efficiencies of large baryon asymmetry production and strong GW signals can be positively correlated. However we find the overall likelihood of observably large GW signals to be smaller than estimated in previous studies. In particular, only detonation-type transitions are predicted to produce observably large gravitational waves.

Phase transitions in the early universe provide an opportunity for probing physics at high scales through cosmological observables, in particular, if the transition is first order. In that case, it may be possible to explain the origin of baryonic matter through electroweak baryogenesis (EWBG) [1][2][3][4] or variants thereof [5]. Such transitions can also produce relic gravitational waves (GWs) that may be detectable by future experiments like LISA [6,7], BBO [8], DECIGO [9,10] and AEDGE [11].
In the present work we continue the investigation started in Ref. [62], which determined v w and L w over much of the model parameter space, but did not try to predict the baryon asymmetry or GW production. Moreover, that study was limited to subsonic wall speeds, due to a breakdown of the fluid equations that determine the friction on the wall. Recently a set of improved fluid equations was postulated in Refs. [63,64], that do not suffer from the subsonic limitation. We use these in the present work in order to fully explore the parameter space, where high v w can be favorable to observable GWs, and also compatible with EWBG. It will be shown that for strong deflagrations, the fluid velocity in front of the wall saturates and even decreases with increasing wall velocity v w . Since the walls become thinner at the same time, the baryon asymmetry is enhanced at larger wall velocities for these transitions, becoming positively correlated with a strong GW signal. Despite this positive correlation, we find that producing the observed baryon asymmetry together with a GW signal detectable in next generation observations is not possible, in contrast to previous estimates [34,43]. The difference comes from several factors working in the same direction. For example, we find larger wall velocities and thicknesses than Ref. [43], which suppress the baryon asymmetry. Moreover, our GW fits include a recently derived suppression factor due to shock reheating [65,66], which leads to a much weaker GW signal for strong deflagrations.
A further improvement in this work is to present an ultraviolet completion of the effective coupling that gives rise to the CP-violation needed for EWBG. We introduce heavy vectorlike top partners which when integrated out induce a CP-violating coupling of the singlet scalar s to top quarks, giving the source term for EWBG. 1 Although the effective operator description of this term is quite adequate for quantitatively understanding EWBG [68,69], its resolution in terms of underlying physics is necessary for quantifying how large its coefficient can be, consistent with laboratory constraints. We present the details in section II, including comprehensive collider limits on the top partners and the subsequent constraints on the effective theory. The finite-temperature effective potential of the theory is also outlined there, along with a discussion of cosmological constraints on the small explicit breaking of the Z 2 symmetry, that is necessary for EWBG.
The paper continues in Sect. III with a brief description of our methodology for finding the high-temperature first-order phase transitions, and characterizing their strength. This is followed in Sect. IV by a detailed account of how the bubble wall speed and shape are determined. The techniques for computing the baryon asymmetry and GW production are described in Sect. V. We present the results of a Monte Carlo exploration of the model parameter space with respect to these observables in Sect. VI, with emphasis on the interplay between successful EWBG and potentially observable GWs. Conclusions are given in Sect. VII, followed by several appendices containing details about construction of the finite-temperature effective potential, solving junction conditions for the phase transition boundaries, and predicting GW production.

II. Z2-SYMMETRIC SINGLET MODEL
We study the Z 2 -symmetric singlet scalar extension of the SM with a real singlet s coupled to the Higgs doublet H. The scalar potential is We work in unitary gauge, which consists of taking H = h/ √ 2; the Goldstone bosons still contribute to the one-loop and thermal corrections, but they are set to zero in the tree-level potential. We assume µ 2 h < 0 and µ 2 s < 0, which implies that the potential has non-trivial minimums at v ≡ h = ±|µ h |/ √ λ h ≈ 246 GeV, s = 0 and h = 0, s = ±|µ s |/ √ λ s . The scalar fields' mass in the vacuum can then be written in terms of the parameters of the potential as m 2 h = −2µ 2 h ≈ (125 GeV) 2 and m 2 s = −λ hs µ 2 h /(2λ h ) + µ 2 s .
The other relevant interaction of s is a dimension-5 operator yielding an imaginary contribution to the top quark mass [70]: This term will be ignored during the discussion on the phase transition; however it is essential for generating the baryon asymmetry, since it gives the CP-violating source term when s temporarily gets a VEV in the bubble walls of the electroweak phase transition. In Eq.
(2) we have adopted a special limit of a more general model, in which the dimension-5 contribution is purely imaginary. This can be understood as a consequence of imposing CP in the effective Lagrangian, with s coupling like a pseudoscalar, s → −s. Hence it is consistent to omit terms odd in s in the scalar potential (1), even though Eq. (2) is odd in s. The CP symmetry prevents a VEV from being generated for s by loops.
The effective operator is generated by integrating out a heavy singlet vectorlike top quark partner T , whose mass term and couplings to the third generation quarks q L = (t L , b L ), Higgs and singlet fields are including also the SM q L -Higgs coupling. This is invariant under CP if s → −s. 2 Integrating out T leads to the effective operator in (2) with scale We consider experimental constraints on the scale Λ below.
In previous literature, thermal corrections were frequently approximated by including just the first term of the high-temperature expansion of the thermal functions presented in the Appendix B. However, this approximation fails at temperatures below the mass of particles strongly coupled to the Higgs, as can happen in models with a high degree of supercooling. Therefore, we employ the full one-loop thermal functions. This will be shown to have a large impact on the values of the tunneling action, and thus of the nucleation temperature. In addition to the tree-level potential and the thermal corrections, we also include the one-loop correction and the thermal mass Parwani resummation [71]. The complete effective potential then becomes The details are presented in Appendix A.

A. Laboratory constraints
It is important to determine how low the scale Λ of the dimension-5 operator in Eq. (4) can be, since it has a strong impact on the baryon asymmetry η b ; in the limit of large Λ, η b scales as 1/Λ. The relevant masses and couplings are constrained by direct searches for the top partner and precision electroweak studies. Moreover the properties of the singlet s are constrained by collider searches.
After electroweak symmetry breaking, a Dirac mass term (t L , For example, we consider a benchmark point with η 1 = 0.55 and a physical T mass M T = 800 GeV, which correspond to M = 794 GeV and mixing angles θ L = 0.126 and θ R = 0.027. The relations between y t and the physical top mass differ from the SM ones by less than 1%, which is allowed by current LHC constraints [72,73]. For sufficiently large η 2 , decays of T to ht/Zt/W b induced by mixing are highly subdominant to T → st, and searches for vector-like top partners that focus on the former channels are evaded. Near the Goldstone-equivalent limit (which should apply reasonably well for M T = 800 GeV and relatively small s masses, m s ∼ 100 GeV), the branching ratio for T → st is We roughly estimate from Refs. [74,75] that for M T = 800 GeV, vector-like quark searches that target SM final states are evaded provided B(T → st) > ∼ 90%, corresponding to η > ∼ 2.4 for our benchmark point. Ref. [76] (see Fig. 1 of contribution 5; also [77]) has reinterpreted collider bounds to constrain the parameter space (m s , M T ) for models in which T → st dominates, finding that top partner masses above ∼ 750 GeV are allowed in the case where s decays 100% into two gluons. This is true in our model, where the dominant s decays are induced by the loop diagrams shown in Fig. 1. One can estimate that the gluon final state dominates over that of b quarks by a factor of (g 2 s m s /g 2 w m b ) 2 10 3 , and over decays into photons by (g s /e) 4 ∼ 300. Precision electroweak data constrain the additional contributions to the oblique parameters, especially T , which is corrected by [78] where T sm = 1.19 is the SM value, c L = cos θ L , s L = sin θ L , and r = (M T /m t ) 2 ; the upper limit is from section 10 of [79]. The benchmark point chosen above almost saturates this constraint, giving ∆T 0.09.
There are also direct searches for resonant production of the singlet, by gluon-gluon fusion. The coupling of s to t in the mass eigenstate basis is y st = η 2 cos θ R sin θ L ∼ η 2 θ L , while that to T is y sT = −η 2 cos θ L sin θ R ∼ −η 2 θ R . The squared matrix element for the decays s → gg is [80]  where τ i = 4m 2 i /m 2 s . The parton-level production cross section for gg → s isσ = π|M| 2 δ(ŝ − m 2 s )/(256ŝ) where the 256 comes from averaging over gluon colors and spins. Integrating this over the gluon PDFs gives the hadron-level cross section in which dependence on m s drops out except in the parton luminosity factor L g . This production is probed via decays s → γγ, whose branching ratio is approximately B(s → γγ) = (8/9)α 2 /α 2 s [80]. For the dominant s → gg decay into gluons, in principle LHC dijet resonance searches could be constraining, but these exist only for m s 500 GeV which is beyond the range of interest for the present study. To a good approximation, σ(pp → s) is determined by m s and Λ. In Fig. 2(a) we show limits from ATLAS [81,82] and CMS [83] on σB(s → γγ) as a function of m s , along with the predictions for various Λ, and in Fig. 2(b) we show the associated lower bounds on Λ. In the low-mass region (65 GeV < m s < 110 GeV), lower bounds on Λ range roughly from 400 GeV to 650 GeV; in the intermediate-mass region (110 GeV < m s < 160 GeV), Λ is not yet constrained by diphoton resonance searches, and for much of the high-mass region (m s > 160 GeV), Λ is bounded to be above 1 TeV. For our subsequent scans of parameter space, we adopt a fixed reference value for Λ, which is large enough to be consistent with much of the low-m s region. Because Λ ref is well below the lower-bounds on Λ in the high-mass region, we confine our scans to m s < 160 GeV for consistency. 3 The constraints from precision electroweak data, diphoton resonance searches, and vector-like quark searches are shown in the η 1 -η 2 plane in Fig. 3, for M T = 800 GeV, where we approximate the T search constraints by the requirement B(T → st) > 0.9, and for M T = 1300 GeV, heavy enough to evade T searches for any B(T → st). For the chosen m s , it is apparent that the reference value Λ = 540 GeV is attainable for η 2 > ∼ 2.5 for M T = 800 GeV and η 2 > ∼ 3 for M T = 1300 GeV. For slightly heavier s in the window 110 GeV < m s < 160 GeV, diphoton resonance searches are evaded and the red contours disappear. In this case even lower values of Λ are allowed provided one is willing to consider larger values of η 2 . Since the baryon asymmetry η b scales roughly as 1/Λ, it is straightforward to reinterpret our final results for larger (or smaller) Λ. From the results of Section VI one can infer that a significant fraction of models remain viable for baryogenesis for Λ = 2Λ ref (or for even larger Λ), a scale consistent with more modest couplings, η 2 ∼ 1.5.
Allowing for very large values of η 2 could invalidate the effective theory above the heavy top partner threshold M at scales only slightly larger than M , which would require us to specify additional new physics in order to have a complete description. There are two principal challenges arising from the running of the couplings, where µ denotes the renormalization scale. The most serious problem is that for large values of η 2 , the selfcoupling λ s is quickly driven to zero, and the scalar potential becomes unstable. The second is that η 2 reaches a For selected T and s masses, constraints on η1 and η2 from precision electroweak data (green), diphoton resonance searches [82,83](red), and searches for vector-like quarks [74] (blue), along with contours of Λ in GeV. The allowed region is unshaded.
Landau pole at somewhat higher scales. The first problem could be ameliorated by coupling additional scalars to s, without impacting our results for EWBG or GWs. For this reason, we do not limit the scope of our investigation based on the running of λ s . Regarding the second problem, we note that even for η 2 = 3, the Landau pole is nearly an order of magnitude above M , which we consider to be an acceptably large range of validity for the effective theory.

B. Explicit breaking of Z2 symmetry
Since we are considering a scenario where the Z 2 symmetry s → −s is spontaneously broken during the early universe and restored at the EWPT, domain walls form before the EWPT, and the universe will consist of domains with random signs of the s condensate. The source term for EWBG that arises from Eq. (2) is linear in s, resulting in baryon asymmetries of opposite signs, that could average to zero after completion of the EWPT. To avoid this outcome, the Z 2 symmetry should be explicitly broken, by potential terms with small coefficients µ b , µ b . We have used the freedom of shifting s by a constant to remove a possible tadpole of s at the true vacuum (h, s) = (v, 0).
The presence of the biasing potential V b can prevent the baryon washout in several ways. First, if the transition to the broken-s phase is of second order, even a small tilt can suffice to make the lower-energy vacuum dominate. Second, in a first order transition, symmetry breaking terms can bias the bubble nucleation rates to prefer the lower-energy vacuum. Indeed, the number of bubbles nucleated during the transition is n ∼ t * tc dt Γ(t), where t * is the time when transition completes, and Γ(t) ∼ exp(−S 3 /T ). Writing the action as S 3± =S 3 ∓ δS in the two respective vacua, the relative number density of bubbles in each phase at the end of the transition becomes n + /n − ≈ exp(2 δS * /T * ). In general [86] S 3 ∝ E, where E is the coefficient of the cubic term in the potential. Using this scaling we may write δS * = (δE/E 0 )S * 3 , where typically S * 3 /T * ≈ 100. In our model E 0 ≈ (3λ s ) 3/2 T /12π, so taking V b = µ b s 3 , corresponding to δE = µ b , and T * ≈ 100 GeV, the condition for single-phase vacuum dominance becomes µ b > ∼ 0.1 λ 3/2 s GeV. Barring very large λ s , this condition is easily met with no limitations on our analysis.
Even if a domain wall network forms, the higher-energy domains will collapse due to pressure gradients, and we should ensure that this process completes before the EWPT. The collapse starts with the acceleration of a wall at relative position R according toR = −∆V /τ , where τ ∼ √ λ s w 3 is the surface tension (distinct from the tension σ used above in the nucleation estimate), ∆V ∼ V b (0, w) ∼ µ b w 3 is the difference in the vacuum energies, and w ∼ µ s / √ λ s is the singlet VEV. Using H = 1/2t and T ≈ 100 GeV, one finds that walls reach light speed in time which is practically instantaneous on the timescales of interest, for reasonable values of µ b . We note that global symmetries like Z 2 are expected to be broken by quantum gravity effects, so that it could be reasonable to anticipate µ b ∼ v 2 /M p ∼ 0.1 eV, which is large enough from the perspective of Eq. (15).
The higher energy domains subsequently collapse at the speed of light, since there is no appreciable friction. The time required for this process to complete is determined by R * = 2a(t 1 ) t2 t1 dt/a(t), where R * is the comoving size of the domain wall separation. By the Kibble mechanism one expects that R * = AH −1 * with A 1, leading to the ratio of domain wall collapse to formation times t 2 /t 1 = (1 + A/2) 2 . The temperature interval corresponding to this time interval is ∆T /T ≈ A, assuming that the growth phase also proceeded at the speed of light.
The temperature of the first phase transition, T 1 can be estimated as that when ∂ 2 V /∂s 2 becomes negative. In the approximation of neglecting V b , and keeping only leading terms in the high-T expansion, one finds T 2 c /c s where T c is the critical temperature of the EWPT, and c s = (3λ s + 2λ hs ) /12. Thus the temperature difference between transitions is of order ∆T 1c ∼ λ h w 2 /(c s T c ). Requiring that ∆T 1c /T c > A then gives Given that A ∼ (T * /S * 3 )(∆T /T ) * ∼ 10 −2 -10 −4 [50], this is a very weak constraint. We conclude that it is easy to avoid cosmological problems associated with the domain walls by small symmetry breaking terms, that do not affect the rest of our analysis.

III. PHASE TRANSITION AND BUBBLE NUCLEATION
In the examples of interest for this work, the phase transition in the Z 2 -symmetric singlet model proceeds in two steps: starting from the high-temperature global minimum h = s = 0, a transition first occurs to nonzero s, while the Higgs field remains at h = 0. This is followed by the EWPT, in which s returns to zero and h develops its VEV. The h 2 s 2 interaction provides the potential barrier to make this a first order transition.
As usual, the first order transition occurs at the bubble nucleation temperature T n , which is below the critical temperature T c , where the two potential minima become degenerate, Bubble nucleation occurs when the vacuum decay rate per unit volume Γ d becomes comparable to H 4 , the Hubble rate per Hubble volume. The decay rate is [87] where S 3 is the O(3) symmetric action, The precise criterion that we use for nucleation is which is satisfied when S 3 /T n ∼ = 140 [88]. We used the package CosmoTransitions [89] to calculate S 3 . The action obtained with the full potential can differ significantly from the commonly used thin wall approximation [90,91] or the approximation of evaluating it along the minimal integration path for the potential [43]. We compare the predictions for nucleation of these approximations to the full one-loop result, for several exemplary models, in Table III. The approximate methods tend to underestimate the action, giving a higher nucleation temperature; hence we use the values derived from the full one-loop action in the following.   There are two complementary parameters for characterizing the strength of the first order transition. One is the ratio of the Higgs VEV to the temperature at the time of nucleation, v n /T n , which is especially relevant for EWBG, as we will discuss in Sect. V B. The other, which is more important for GW production, is the ratio of released vacuum energy density to the radiation energy density [92,93]: where ρ γ = g * π 2 T 4 n /30, g * is the effective number of degrees of freedom in the plasma (we use g * = 106.75) and ∆ denotes the difference between the unbroken and broken phase. α quantifies the amount of supercooling that occurs prior to nucleation, which determines how much free energy is available for the production of GWs.

IV. WALL VELOCITY AND SHAPE
The derivation of the wall velocity and field profiles is a technically demanding problem [50], that was first addressed in the context of Higgs plus singlet models in Refs. [54,56,94], in various approximations. One must solve the equations of motion (EOM) for the scalar sector coupled to a perfect fluid, where z is the direction normal to the wall, that is to a good approximation planar by the time it has reached its terminal velocity. We use a sign convention where the wall is moving to the left, so that z > 0 corresponds to the broken phase. The sum is over all the relevant species coupled to h or s in the plasma, with N i and m i respectively denoting the number of degrees of freedom and the field-dependent mass of the corresponding species, and δf i the deviation from equilibrium of its distribution function. All the temperature-dependent quantities appearing in these equations are evaluated at T + , which is the plasma's temperature just in front of the wall. We calculate T + in Appendix B using the method described in Ref. [93], and δf i will be computed in Sect. IV A.
The terms in Eqs. (22) with δf i represent the friction 4 of the plasma on the wall, that leads to a terminal wall speed v w < 1, unless the friction is too small and the wall runs away to speeds close to that of light. Following previous work, we take the dominant sources of friction to be from the top quark (i = t) and electroweak gauge bosons (i = W ), neglecting the contributions to friction from the Higgs itself and from the singlet. This approximation is bolstered by the smaller number of degrees of freedom N h = N s = 1 compared to N t = 12 and N W = 9, as well as the smallness of the Higgs self-coupling λ h and the not-too-large values of the cross-coupling λ hs that will be favored in the subsequent analysis. Then the friction term for the s equation of motion vanishes, since s couples only to itself and to the Higgs, apart from its suppressed dimension-5 coupling to t. This allows for some simplification in the following procedure.
In Ref. [62], a similar study of the present model was done, where no a priori restriction of the wall shape was assumed, but it was found that the actual shapes conform to a very good approximation to the tanh profiles where h 0 and s 0 are respectively the vacuum expectation values (VEV) of the h and s fields in the broken and unbroken phases. Hence we adopt the ansatz (23), which allows the singlet and Higgs wall profiles to have different widths, and to be offset from each other by a distance L s δ. The s field's VEV is taken to be the usual one evaluated at T + , which solves the equation dV eff (0, s; T + )/ds s=s0 = 0. The situation is more complicated for the h field, for which the Higgs VEV should be evaluated at T − , the plasma's temperature behind the wall. Since we are fixing a constant temperature T + in the potential, the change in the effective action due to the shift in the background temperature must be accounted for by the perturbation in the broken phase. As a consequence we are choosing h 0 so that it solves the equation This choice guarantees that the Higgs EOM is satisfied far behind the wall. We will estimate the uncertainty of our results due to this approximation in Sect. VI D.
To approximately solve the Higgs EOM, one can define two independent moments M 1,2 of E h (z), and assume that they both vanish at the optimal values of v w and L h . A convenient choice is [56] These also have intuitive physical interpretations that naturally distinguish them as good predictors of the wall speed and thickness, respectively. M 1 is a measure of the net pressure on the wall, so that Eq. (25)   illustrated in Fig. 4, which shows the dependence of M 1 and M 2 on v w and L h .
We chose a different approach to determine the singlet wall parameters L s and δ. Instead of solving moment equations analogous to (25,26), one can determine their values by minimizing the s field action with respect to L s and δ. Here s * is a field configuration with arbitrary fixed parameters L * s and δ * , that we choose to be L * s = L h and δ * = 0. The second term is just a constant, but it allows for the convergence of the integral by canceling the contributions of V eff at z → ±∞. This method has the advantage that it does not depend on any arbitrary choice of moments, and it is more efficient to numerically minimize the function of two variables than to solve the system of equations for the moments of the EOMs.

A. Transport equations for fluid perturbations
The final step toward the complete determination of the velocity and the shape of the wall is to compute the distribution functions' deviations from equilibrium δf i , by solving the Boltzmann equation for each relevant species in the plasma. The method of approximating the full Boltzmann equation by a truncated set of coupled fluid equations was originally carried out in Ref. [50], for the regime of slowly-moving walls (see also Ref. [56]). This approach was recently improved in Ref. [64] in order to be able to treat wall speeds close to or exceeding the speed of sound consistently. We briefly summarize the formalism, which we use in the present study.
The out-of-equilibrium distribution function can be parametrized in the wall frame as where β = 1/T + and the ± is + for fermions and − for bosons. δτ and µ are the dimensionless temperature and chemical potential perturbations from equilibrium, and δf u is a velocity perturbation whose form is unspecified, but is constrained by d 3 p δf u = 0. By assuming that the perturbations are small, one can expand f to linear order in µ, δτ and the velocity perturbation δf u to obtain with To simplify the problem, one models the plasma as being made of three different species: the top quark, the W bosons (shorthand for W ± and Z) and a background fluid, which includes all the remaining degrees of freedom. It is convenient to write the velocity perturbation as u ∝ d 3 p (p z /E) δf u when constructing the moments of the linearized Boltzmann equation. By taking three such moments, using the weighting factors 1, E and p z /E, the perturbations are determined by transport equations where prime denotes d/dz, q i = (µ i , δτ i , u i ) , q = (q W , q t ) , the Γ matrices are collision terms, and S is the source term, whose definitions, as well as those of the the matrices A, Γ,Ã −1 bg , Γ bg,t , Γ bg,W , can be found in Ref. [64]. If A and Γ were independent of z, one could use the Green's function method to solve Eq. (31); however, A is a function of m i (z)/T . To deal with this dependence on z, we discretize space, z → z 0 + n∆z with n = 0, · · · , N − 1, and Fourier transform Eq. (31), where the tilde denotes the discrete Fourier transform. This is a linear system that is straightforward to numerically solve forq k . Onceq k is known, it can be transformed back and interpolated to obtain q(z). Eq. (32) can then be integrated using a Runge-Kutta algorithm.
Finally, one can substitute Eq. (29) into the Higgs EOM (22) to express the friction in terms of the fluid perturbations µ i , δτ i and u i . This leads to the result where the functions C m,n v and D m,n v can be found in Ref. [64].

V. COSMOLOGICAL SIGNATURES
We have now established the machinery needed to compute all the relevant properties of the first order phase transition bubbles, starting from the fundamental parameters of the microscopic Lagrangian. In this section we describe how to apply these results for the estimation of GW spectra and the baryon asymmetry.

A. Gravitational Waves
We follow the methodology of Refs. [7,65,66,93,95] to estimate future gravitational wave detectors' sensitivity to the GW signals that can be produced by a first-order electroweak phase transition in the models under consideration. The GW spectrum Ω gw (f ) is the contribution per frequency octave to the energy density in gravitational waves, i.e., Ω gw d ln f is the fraction of energy density compared to the critical density of the universe. The spectrum gets separate contributions from the scalar fields, sound waves in the plasma and magnetohydrodynamical turbulence created by the phase transition: Each of these contributions depends on the wall velocity v w , the supercooling parameter α (Eq. (21)), and the inverse duration of the phase transition, defined as Another useful quantity is the mean bubble separation, which can be written in terms of v w and β as [7] It has been shown in Ref. [57] that interactions with gauge bosons prevent the wall from running away indefinitely towards γ → ∞. In that case, the contribution from the scalar fields has been shown to be negligible. Furthermore, the estimates for the magnetohydrodynamical turbulence are very uncertain and sensitive to the details of the phase transition dynamics [96], and are expected to be much smaller than the contribution from sound waves. Hence, we consider only the effects from the latter, and set Ω m (f ) = Ω φ (f ) = 0. For convenience, we reproduce the numerical fits of the GW spectra derived in Refs. [7,65,66,93,95] in appendix C.
We will use these predictions with respect to four proposed space-based GW detectors: LISA [97], AEDGE [11], BBO [98] and DECIGO [9]. A successful GW detection depends upon having a large enough signal-to-noise ratio [99], where Ω sens (f ) denotes the sensitivity of the detector 5 and T is the duration of the mission. The sensitivity curves for the detector LISA, BBO and DECIGO were obtained from Ref. [100]. Whenever SNR is greater than a given threshold SNR thr , we conclude that the signal can be detected. In general, this threshold can depend upon the configuration of the detector. For all the experiments, we take SNR thr = 10 and T = 1.26×10 8 s. In the following, SNR max will designate the maximum signal-to-noise ratio detected by one of the detectors: While Ω sens (f ) can be obtained from the noise spectrum of a detector, it is not practical to compare it to the GW spectrum directly; one needs to compute the SNR to determine if a signal is detectable. A useful tool for visualizing the sensitivity of a detector is the peak-integrated sensivity curve (PISC) defined in Refs. [101][102][103], which is a generalization of the power-law sensitivity curve [104]. The main advantage of the former is that it does not assume a power-law spectrum, hence it conserves all the information about the SNR. In the simple case where one considers the contribution from only one GW source, the PISC can be obtained by factorizing the GW spectrum as where f p and Ω p = max[Ω gw (f )] are the peak frequency and GW amplitude and S is a function that parametrizes the spectrum's shape, with a maximum at f = f p and S(f p , f p ) = 1. One can then write the SNR as with the PISC By construction, any GW signal that peaks above the PISC has SNR > SNR thr and can therefore be detected.

B. Baryogenesis
The mechanism of electroweak baryogenesis is sensitive to the speed and shape of the bubble wall during the phase transition. In most previous studies, these quantities were treated as free parameters to be varied, but in this work we have already derived them, as was discussed in Section IV. An important requirement for EWBG is to avoid the washout, by baryon-violating sphaleron interactions, of the generated asymmetry inside the bubbles of broken phase, once they have formed. This leads to the well-known constraint [105] v n T n > 1.1 , which was derived within the SM for low Higgs masses where a first order EWPT was possible. The bound can be slightly higher (up to 1.2) in singlet-extended models [106], depending upon the parameters, due to the sphaleron energy being modified. Here we adopt the SM constraint (43); we checked that taking the more stringent bound 1.2 removes ∼ 5% of viable models in the scan over parameter space to be described below.
Near the bubble wall, CP-violating processes associated with the effective interaction in Eq. (2) give rise to perturbations of the plasma, that result in a local chemical potential µ B L for left-handed baryons, which by imposing the chemical equilibrium of strong-sphaleron interactions, is related to those of the t L , t c R and b L quarks by where the K a 1 functions were defined in [107] (K a 1 = D a 0 in the notation of [63]). The µ B L potential biases sphalerons, leading to baryon number violation, whose associated Boltzmann equation can be integrated to obtain the baryon to photon ratio 6 where f sph quantifies the diminution of the sphaleron rate in the broken phase [108,109]. The most challenging step for the computation of EWBG is in the determination of the chemical potentials µ t L , µ t c R and µ b L appearing in Eq. (44). They satisfy fluid equations resembling the network (31,32), except that the potentials relevant for EWBG are CP-odd, whereas those determining the wall profiles are CP-even.
The CP-odd transport equations have been discussed extensively in the literature, leading to two schools of thought as to how best to compute the source term for the CP asymmetries. These are commonly known as the VEV-insertion [110,111] or WKB (semiclassical) [112][113][114][115][116][117] methods, respectively. A detailed discussion and comparison of the two approaches was recently given in Ref. [63], which quantified the well-known fact that the VEV-insertion source tends to predict a larger baryon asymmetry than the WKB source, by a factor of ∼ 10. In the present work we adopt the WKB approach, which was updated in Ref. [63] to allow for consistently treating walls moving near or above the sound speed. In addition, that reference computed the source term arising from the same effective interaction (2) as in the present model, so we can directly adopt the CP-odd fluid equations studied there.

VI. MONTE CARLO RESULTS
To study the properties of the phase transition, we performed a scan over the parameter space of the models, imposing several constraints. We found that variations in λ s do not qualitatively change the results, prompting us to initially fix its value at λ s = 1, leaving λ hs and m s as the free scalar potential parameters. We will first discuss this slice of parameter space, and later consider the quantitative dependence on λ s . We also chose Λ = 540 GeV, which is conservative since there are no collider constraints on its value for singlet masses in the region m s = [110,160] GeV. Recall that Λ is important for the determination of the baryon asymmetry η b , which is expected to scale roughly as 1/Λ. Finally, in order to prevent Higgs invisible decays, we imposed m s > m h /2.
We used a Markov Chain Monte Carlo algorithm to efficiently explore the regions of parameter space having desired phase transition properties. Starting with an initial model satisfying the sphaleron bound (43), one generates a new trial model by randomly varying the parameters λ i by small increments δ i . The trial model is added to the chain using a conditional probability that favors models having strong first order phase transitions, and for which a solution to the nucleation condition (20) can be found. We adjust the δ i so that roughly half of the models are kept in successive trials, with larger values of δ i being more likely to result in a rejection. This procedure yielded 842 models with strong phase transitions, of which 712 were amenable to finding solutions for the moment equations (25)(26). Our analysis typically works for γ 10; for faster walls, the algorithm for determining the wall properties becomes numerically unstable and does not yield reliable results. This is due to the large (500 × 500) matrix (A −1 Γ) of eq. (33) becoming singular as v w → 1. It is therefore difficult to determine the type of solution of the 130 remaining models using our methodology alone: they could either stabilize at ultrarelativistic speeds, or (from a naive perspective-see below) run away indefinitely towards γ → ∞. The value of the baryon asymmetry should not be affected by this ambiguity since it is negligible for v w ≈ 1. The GW spectrum produced during the phase transition is sensitive to this distinction since runaway walls have a nonnegligible fraction of their energy stored in the wall, while for non-runaway walls, the energy gets dissipated into the plasma, so the fraction of energy in the wall becomes negligible. This ambiguity can be lifted using the result of Ref. [57], which found that in the limit γ → ∞, interactions between gauge bosons and the wall create a pressure proportional to γ, preventing it from running away. 7 We therefore assume that the 130 models without a solution to the moment equations (25)(26) correspond to non-runaway walls with v w ≈ 1. The results of this scan, showing the calculated wall velocity, signal-to-noise ratio of gravity waves observable by at least one of the proposed experiments (LISA, AEDGE, BBO or DECIGO), and the predicted baryon asymmetry (in units of the observed value) are presented in Fig. 5, in the plane of of λ hs versus m s .

A. Deflagration versus detonation solutions
A striking feature of these results is that all the detonation solutions have v w ≈ 1. 8 We have tested that this is not specific to the choice of fixed parameter values, but also holds for all models having 0.01 < λ s < 8 and Λ > 110 GeV; hence it seems to be a general property of phase transitions in the Z 2 -symmetric singlet framework. One can understand this behavior by considering the net pressure opposing the wall's expansion, M 1 (recall Eq. (25-26)), as a function of the wall velocity, as illustrated in Fig. 6. It shows how M 1 differs when evaluated with the appropriate quantities v + , T + rather than the incorrect ones v w , T n . Using the latter, we would find no solution to the equation M 1 = 0 for the exemplary model used in Fig. 6, and would then incorrectly conclude that it satisfies v w ≈ 1. The relevant quantities are those measured right in front of the wall, v + and T + . The speed v + is smaller than v w for v w < ξ J , which would lower the pressure against the wall (ξ J is the Jouguet velocity, defined as the smallest velocity a detonation solution can have). However, in the same region, the temperature T + is larger than T n , which causes the pressure to increase. The latter effect turns out to dominate over the 7 More recently, the authors of Ref. [58] have carried out an all-orders resummation at leading-log acuracy, finding that the pressure is in fact proportional to γ 2 for fast-moving walls. 8 Strictly speaking there are models with vw < 1 detonation solutions but these always have another solution at a lower velocity corresponding to a deflagration or hybrid wall. Then only the latter solution is physically relevant, since the bubble is created at vw = 0 and accelerates until it reaches the solution with the lowest velocity. former. Indeed, the actual pressure, represented by the solid blue line in Fig. 6, increases much more rapidly than M 1 (v w , T n ) close to the speed of sound. This qualitative difference allows for a solution to M 1 = 0, which would have been missed if we had used the naive quantities v w and T n .
We find that the previous statements apply quite generally: for all models, T + > T n when v w < ξ J , and this always leads to a much higher pressure on the wall, even if the difference between T + and T n is quite small; the pressure barrier at v w = ξ J is always greater than the maximum possible value for a detonation solution. Therefore, if the phase transition is strong enough to overcome the pressure barrier at ξ J , the solution becomes a detonation, but the pressure in the region v w > ξ J is never enough to prevent it from accelerating towards v w ≈ 1. If the phase transition is weaker, the pressure barrier is high enough to impede the detonation, and it becomes a deflagration or hybrid solution.
The wall thickness and speed for the models with deflagration 9 solutions are shown in Fig. 7, which demonstrates that the behaviors for subsonic (deflagration) and supersonic (hybrid) walls are qualitatively different. Subsonic walls generally have v + ≈ v w , which is expected since the fluid should not be strongly perturbed by a slowly moving wall. The wall width is not uniquely determined by v w , but there exists a clear correlation, with slower walls being thicker. For supersonic cases, the correlation between v + and v w gets inverted: higher wall velocity leads to lower v + . The wall width becomes uniquely determined by v w and the relation between these two variables is to a good approximation linear. One observes that stronger phase transitions, quantified by higher values of α, generally produce faster and thinner walls. Even for the strongest transitions our solutions still have wall thickness LT > ∼ 3. Since the semiclassical force mostly affects particles with momenta k z ∼ T , we find L k z > ∼ 3, so that the semiclassical approximation is still valid. In fact the semiclassical picture has been shown to remain valid for surprisingly narrow walls [118], working very well for L k z ≈ 4 and still reasonably for L k z ≈ 2. There is a linear correlation between the h and s wall widths, but the slope is not 1; in all cases, we find that L h > L s . The distribution of wall offset values δ is also indicated in Fig. 7(c).

B. Baryogenesis and gravity wave production
Of the 842 sampled models, 517 are able to generate the baryon asymmetry at a level large enough to agree with observations, and 20 detonation walls can produce observable gravitational waves. We found no detectable deflagration solutions. More detailed results are presented in Table II. The complementarity of the experiments considered here, with respect to the present model, can be appreciated by considering the relation between the maximum GW amplitude 10 max[Ω gw h 2 ] and the frequency of this peak amplitude f max , as shown in Fig. 8 (a). The peak frequency of the strongest detonation walls are positioned exactly in LISA's region of maximal sensitivity, while the peak frequency of the deflgration solutions are closer to the peak sensitivity of AEDGE, DECIGO and BBO. The complete spectrum's shape are also shown in Fig. 8 (b,c) for deflagration and detonation solutions respectively. We conclude that detonation walls could be probed by LISA, DECIGO and BBO, but not by AEDGE.
In previous studies, where the wall velocity was considered as a free parameter, there was an expectation that baryogenesis would be less efficient with increasing v w , whereas gravity waves would become more so. In the present study, where v w is not adjustable but is a derived parameter, we surprisingly find that rather than EWBG and stronger GWs being anticorrelated, instead they are positively correlated, as is illustrated in Fig. 9 (a). This can be understood from the fact (see Fig. 7 (b)) that L h is a decreasing function of v w , which enhances EWBG. Moreover, the relevant velocity for EWBG is v + , which is a decreasing function of v w for supersonic walls, and is bounded by v + < c s ; this effect also enhances EWBG for fast-moving walls. The actual relation between η b and v w is shown in Fig. 9 (b) and, at least for supersonic walls, there is a positive correlation between these two variables. Fig. 9 also indicates that the supercooling parameter α is positively correlated with both η b and SNR max : stronger phase transitions generally lead to both higher GW and baryon production. Detailed predictions for EWBG in the Z 2 symmetric model were previously made in Refs. [43] and [34], as opposed to merely requiring the sphaleron bound (43) to be satisfied. Comparisons with the present work are hindered by the fact that different source terms for the CP asymmetry were assumed. In Ref. [43], the dimension-6 coupling i(y t / √ 2)(s/Λ) 2h t L t R was used, rather than the dimension-5 coupling in Eq. (2). Moreover, a value v w = 0.2 was taken for the wall velocity, and an estimate L h = v n / √ 8V b was made for the wall width, where v n is the Higgs VEV at the nucleation temperature, and V b is the potential barrier between the two minima. For the same potential parameters (λ s = 0.1) as in [43], we find no values of v w below 0.43, and our determination of L h is two to three times larger than the estimate in [43]. Both of these discrepancies would lead to overestimating the efficiency of EWBG, helping to explain why Ref. [43] obtains a high frequency of successful models, despite the extra suppression that should result from using a dimension-6 source term.
In Ref. [34], the dimension-5 coupling to leptons rather than the top quark was studied, and a different formalism (the VEV insertion approximation) for computing the CP asymmetry was employed, which tends to give significantly larger estimates for the baryon asymmetry than the WKB method that we adopt [63]. For the parameters of the benchmark models taken in that paper, we find significantly higher wall velocities, v w ∼ 0.6-0.7 than the values v w 0.1 that were needed to match the observed baryon asymmetry there. This can be compensated by increasing the CP-violating phase φ = 0.02 assumed there by a factor of ∼ 10. We are reanalyzing this alternative source term within the EWBG formalism used in the present paper (work in progress).

C. Dependence on λs and Λ
To study the quantitative dependence on the singlet self-coupling λ s , we performed 3 other scans similar to the one previously described, taking λ s = 0.01, 0.1 and 8 (the largest value being near the limit of perturbative unitarity) and Λ = 540 GeV. The results of these scans are summarized in Table II. We find that EWBG remains efficient for λ s 0.1. Again, we found no deflagration walls producing detectable GW, and no models detectable by AEDGE. These results confirm that only detonation solutions, which are not good candidates for EWBG, could be probed by GW detectors. Increasing λ s generally leads to stronger phase transitions, resulting in more models with successful EWBG and detectable GWs.
The value of Λ (recall Eq. (4)) can in principle also have an effect on the strength of the phase transition, through the effective potential's dependence on the top quark mass. The leading thermal term added to the potential varies like h 2 s 2 T 2 /Λ 2 , which becomes negligible at high Λ, but could significantly modify the behavior of the phase transition for Λ ∼ T n , resulting in a larger baryon asymmetry and GW production. We have verified that this term is already subdominant when Λ = 540 GeV. However, for m s > 110 GeV, the weaker constraints allow for values of Λ as low as 300 GeV, which could have an important effect on the phase transition.
To test the sensitivity to lower values of Λ, we repeated the previous scans using Λ = Λ min (m s ), where Λ min is The results are shown in Table II 11 . As one could anticipate from the relation η b ∼ 1/Λ, EWBG is more efficient at lower values of Λ. One can also see that the number of detonation walls or walls generating detectable GW does not change substantially, which indicates that the lower values of Λ do not change the character of the phase transition.

D. Theoretical uncertainties
In Ref. [64], the integrals that determine the collision rates Γ appearing in the Boltzmann equation network (31-32) were reevaluated, and it was noticed that the leading log approximation that was used in their derivation leads to theoretical uncertainties of O(1) in the fractional error. To study the impact of these uncertainties on our results, we recomputed the wall velocity with uniformly rescaled collision rates, Γ → 2Γ and Γ → Γ/2. The ensuing variations of velocity ∆v and wall width ∆L are shown in Figs. 10 (a) and (b) respectively. The effect on v w can be significant for slow walls, leading to a ±40 % change when v w ∼ 0.2. On the other hand for nearly supersonic walls, v w c s , the wall speed is quite insensitive to Γ. The variation of L h is generally below 5%, much smaller than the corresponding variation in Γ.
This behavior is not surprising since, near the speed of sound, the pressure on the wall is mainly determined by the variation of T + , which does not depend on Γ. Likewise, the results for the baryon asymmetry and GW production turn out to be relatively robust against variations in Γ. This is demonstrated by the error intervals in the λ s = 1 row of Table II. The error on the ratio of models satisfying η b /η obs > 1 or SNR i > 10 is of order 10%, which is much smaller than the range of variation in Γ.
Another source of uncertainty is the discrepancy between the temperatures computed with the Boltzmann equation (see Section IV A) and the conservation of the energy-momentum tensor (see Appendix B). Ideally one should obtain T + = T BE (z → −∞) and T − = T BE (z → ∞), where T BE (z) = T + (1 + δτ bg (z)) is the local temperature calculated with the Boltzmann equation. The first condition is always satisfied since we impose the boundary condition δτ bg (−∞) = 0, but we fail to recover the second one due to the different approximations made in the two methods. The discrepancy becomes larger as v w approaches the Jouguet velocity ξ J , where T + increases compared to T − ≈ T n (see Fig. 6 (b)). On the other hand, δτ bg does not change significantly in the same region. Hence, we observe an error in the temperature of order ∆T = T − − T BE (∞) ≈ T − − T + . Since the temperature is not accurate in the broken phase, the Higgs EOM is not automatically satisfied asymptotically. To solve that problem, we shift the actual Higgs VEV h − evaluated in the broken phase by an amount −∆h, so that the adjusted VEV h 0 = h − − ∆h asymptotically solves the EOM (see Eq. (24)). This gives an additional source of uncertainty for v w and L h .
We estimate the errors induced on v w and L h by ∆T and ∆h, assuming they are small enough to justify keeping just the first order terms. Assuming that v w is completely determined by the solution of M 1 = 0 and L h by M 2 = 0, the error on these solutions can be obtained by expanding around the estimated values. For example, for the error in the wall velocity is estimated by where , and we integrate over the temperature variation because M 1 is a functional of T (z). Since v w is the solution of M 1 (v w , h 0 , T (z)) = 0, the absolute errors on v w and L h are estimated as where ∆ T M i = dz(δM i /δT (z))∆T (z). Notice that Eq. (49) overestimates the errors since ∆ T M i and ∆ h M i have opposite signs. From Eqs. (22,25,26), one can see that the functional derivative δM i /δT (z) can be approximated by d dT (∂V eff /∂h), so that where F 1 = h and F 2 = h (2h − h 0 ). We can simplify this integral with the approximation ∆T (z) ≈ (T − − T + )[1 + tanh(z/L h )]/2. Furthermore, we approximate d dT ∂V ∂h as being constant and half of its maximal value, occurring near z = 0. Then where C 1 = dzF 1 (z)[1 + tanh(z/L h )]/2 = h 0 /2 and C 2 = h 2 0 /6. Substituting this expression in Eq. (49), we finally obtain that the errors on v w and L h are given by The relative errors are presented in Fig. 10 (c) for the scan with λ s = 1 and Λ = 540 GeV. The error on v w is below 7% for 97% of the models, and exhibits no strong correlation with v w . This happens because ∆T = T − − T + and dM 1 /dv w are roughly proportional (see Fig. 6), and therefore cancel each others' contributions. The relative error on L h is small at low velocity (or large L h ), but becomes more significant near the speed of sound, however without ever exceeding 10%.

E. Comparison of the GW signal with previous studies
We end this section with a brief comparison with recent studies of the GW produced during a first-order electroweak phase transition. With the prospect of the upcoming LISA experiment, numerous forecasts of the GW spectrum have been made for various extensions of the Standard Model [47,49,[120][121][122]. Most of these find regions of model parameter space that would produce detectable GWs. Here we focus on studies of the singlet scalar extensions [34,43,44,48,123,124].
Our results agree qualitatively with the conclusions of previous work, in the prediction of GWs detectable by LISA, DECIGO and BBO. However there are distinctions stemming from differences in methodology. To compute the GW contribution from the sound waves, previous authors used the numerical fit presented in Ref. [6], while we used the updated formulas of Refs. [65,66]. This leads to a smaller peak frequency, decreasing the number of detectable models. Ref. [6] also does not include the factor 1 − (1 + 2HR/ √ K sw ) −1/2 in the GW amplitude (see Appendix C). We find that this factor is generally quite small (of order 10 −3 -10 −2 for deflagrations and 10 −2 -10 −1 for detonations); hence the predicted GW signals are considerably reduced.
Another significant difference arises from our determination of the wall velocity, which was treated as a free parameter in previous work, whereas we have computed it from the microphysics. The GW spectrum and hence signal-to-noise ratio and ultimately the detectability are strongly dependent on the wall speed. For example, Ref. [123] assumed v w = 0.95 for all models, which considerably enhanced GW production and led to more optimistic predictions. Moreover, using a fixed value for v w hides the discontinuous transition between the deflagration and detonation solutions shown in Fig. 8.

VII. CONCLUSION
In this work we have taken a first step toward making complete predictions for baryogenesis and gravity waves from a first order electroweak phase transition, starting from a renormalizable Lagrangian that gives rise to the effective operator needed for CP-violation. This is in contrast to previous studies in which quantities like the bubble wall velocity or thickness were treated as free parameters, instead of being derived from the microphysical input parameters as we have done here. This is a necessary step for properly assessing the chances of having successful EWBG and potentially observable GWs, since the two observables are correlated in a nontrivial way, when they are both computed from first principles.
We have incorporated improved fluid equations, both for the CP-even perturbations that determine the friction acting on the bubble wall [64], and for the CP-odd ones that are necessary for baryogenesis [63], that can properly account for wall speeds close to the sound barrier. Earlier versions of these equations were singular at the sound speed, making reliable predictions impossible for fast-moving walls. Contrary to previous lore, we find that EWBG can be more efficient for faster walls, due in part to the tendency for fast walls to be thinner.
The Z 2 -symmetric singlet model with vector-like top partners, analyzed in this work, was chosen for its simplicity, but the methods we used can be applied to other particle physics models that could enhance the EWPT. For example, singlets with no Z 2 symmetry have additional parameters, and would thus be likely to have more freedom to simultaneously yield large GW production and sufficient baryogenesis. It would be interesting to identify other UV-completed models with these properties. A limitation we identified with the Z 2 -symmetric model is that for the large values of the η 2 coupling that are desired for EWBG, the singlet self-coupling is rapidly driven toward zero by renormalization group running, above the top partner threshold.
For future work, some improvements could be made to the analysis presented here. The wall velocity might be more accurately determined at low v w by using collision rates for the fluid perturbation equations beyond leading-log accuracy, and by including the singlet and Higgs out-of-equilibrium (friction) contributions. Another limitation is that the current state-of-the-art for predicting the GW spectrum is subject to large systematic uncertainties for wall velocities close to the speed of sound. Since a large fraction of deflagration transitions have 0.5 < ∼ v w < ∼ ξ J , our analysis of the GW production could greatly benefit from more accurate fits in that range of wall speeds.
V CW is the Coleman-Weinberg potential in the MS renormalization scheme that incorporates the vacuum one-loop corrections and V T is the thermal potential: and the n i 's are the particle's number of degrees of freedom: n W T = 4, n W L = n Z T = 2, n Z L = n γ L = 1, n 1,2 = 1, n χ = 3, n t = −12.
We adopt the method developed by Parwani [71] The thermal masses for the longitudinal mode of the photon and Z boson are At low temperature (m 2 i /T 2 1), one would expect all the thermal effects to be Boltzmann suppressed, since the species i becomes essentially absent from the plasma. This is manifestly the case for V T , since the thermal integrals decay exponentially in the limit M 2 i /T 2 ≈ m 2 i /T 2 1. However, in the same limit, V CW would depend quadratically on T if we used the thermal masses defined above. This would spoil the potential's low-T behaviour. Therefore, we define a regulated thermal mass 12M 2 is a regulator chosen to recover the right behaviour in the low and high-T limit. In order to do so, it should be a smooth function satisfying R(x = 0) = 1 and R(x) ∼ e − √ |x| when |x| 1. We choose here the integrated Boltzmann number density function given by where K 2 is the modified Bessel function of the second kind and [x] = x tanh(x) is a smoothed absolute value.
The last term of Eq. (A1) contains the following counterterms: δV (h, s) = Ah 2 + Bh 4 + Cs 2 + D, which are fixed by requiring the renormalization conditions 0 = ∂V eff ∂h h=v,s=0,T =0 While the use of the resummed one-loop potential is a clear improvement over the leading thermal-mass-corrected approximation, one should keep in mind that higher loop corrections and even nonperturbative physics may be relevant, in particular for very strong transitions [125][126][127].
Appendix B: Relativistic fluid equation We here calculate the hydrodynamical properties of the plasma close to the wall using the method described in Ref. [93]. The quantities of interest are the temperatures T ± and the velocities of the plasma measured in the wall frame v ± . The subscript + and − indicate that the quantity is measured in front or behind the wall respectively.
By integrating the conservation of the energy-momentum tensor equation across the wall, one can show that the quantities T ± and v ± are related by the equations where α + and r are defined as These quantities are often approximated by the so-called bag equation of state, which is given in Ref. [93]. This approximation is expected to hold when the masses of the plasma's degrees of freedom are very different from T , which is not necessarily true in the broken phase. Therefore, we keep the full relations (B2) in our calculations.
Subsonic walls always come with a shock wave in front of the phase transition front. The Eqs. B1 can be used to relate T ± and v ± at the wall and the shock wave, but we need to understand how the temperature and fluid velocity evolve between these two regions. Assuming a spherical bubble and a thin wall, one can derive from the conservation of the energy-momentum tensor the following differential equations where v is the fluid velocity in the frame of the bubble's center and ξ = r/t is the independent variable, with r the distance from the bubble center t the time since the bubble nucleation. With that choice of coordinates, the wall is positioned at ξ = v w . µ is the Lorentz-transformed fluid velocity and c s is the speed of sound in the plasma The last approximation is valid for relativistic fluids, which models well the unbroken phase. In the broken phase, the particles get a mass that can be of the same order as the temperature, and it causes the speed of sound to become slightly smaller.
One can find three different types of solutions for the fluid's velocity profile: deflagration walls (v w < c − s ) have a shock wave propagating in front of the wall, detonation walls (v w > ξ J ) have a rarefaction wave behind it and hybrid walls (c − s < v w < ξ J ) have both shock and rarefaction waves. ξ J is the model-dependent Jouguet velocity, which is defined as the smallest velocity a detonation solution can have. Each type of wall have different boundary conditions that determine the characteristics of the solution. Detonation walls are supersonic solutions where the fluid in front of the wall is unperturbed. Therefore, it satisfies the boundary conditions v + = v w and T + = T n . For that type of solution, Eqs. (B1) can be solved directly for v − and T − . Subsonic walls always have a deflagration solution with a shock wave at a position ξ sh that solves the equation v − sh ξ sh = (c + s ) 2 , where v − sh is the fluid's velocity just behind the shock wave measured in the shock wave's frame. It satisfies the boundary conditions v − = v w and T + sh = T n . Because these boundary conditions are given at two different points, the solution of this system can be somewhat more involved than for the detonation case. Indeed, one has to use a shooting method which consists of choosing an arbitrary value for T − , solving Eqs. (B1) for T + and v + , integrating Eqs. (B3) with the initial values T (v w ) = T + and v(v w ) = µ(v w , v + ) until the equation µ(ξ, v(ξ))ξ = (c + s ) 2 gets satisfied. One can then restart this procedure with a different value of T − until the Eqs. (B1) are satisfied at the shock wave. Hybrid walls satisfy v + < c − s < v w and they have the boundary conditions v − = c − s and T + sh = T n , which make them very similar to the deflagration walls.
Numerical fits for the efficiency coefficient κ sw (the fractions of the available vacuum energy that go into kinetic energy) were presented in [93]. For non-runaway walls, these fits depend on the wall velocity and are given by We caution that while these fits, when used as input for a signal-to-noise estimate, are useful to get an overall estimate for the GW signal in a given model, their precise predictions should be interpreted with care. The fit for the sound wave production is reliable for relatively weak transitions α < 0.1, which is the range where most of our models fall. For stronger transitions the fit can overestimate the GW-signal by as much as a factor of thousand (strong deflagrations) [128]. In addition to the strength of the transition, fit parameters have also been shown to be sensitive to the shape of the effective potential [129] and the wall velocity [7,66]. As explained in Ref. [7] Eqs. (C1-C3) are not expected to be accurate for 0.5 < ∼ v w < ∼ ξ J , which includes a large fraction of the deflagration models found in this work. Thus, pending improvements in the theoretical predictions for GW spectra in this range of wall speeds, the results should not be regarded as conclusive.