Nonperturbative analysis of the gravitational waves from a first-order electroweak phase transition

We present the first end-to-end nonperturbative analysis of the gravitational wave power spectrum from a thermal first-order electroweak phase transition (EWPT), using the framework of dimensionally reduced effective field theory and pre-existing nonperturbative simulation results. We are able to show that a first-order EWPT in any beyond the Standard Model (BSM) scenario that can be described by a Standard Model-like effective theory at long distances will produce gravitational wave signatures too weak to be observed at existing and planned detectors. This implies that colliders are likely to provide the best chance of exploring the phase structure of such theories, while transitions strong enough to be detected at gravitational wave experiments require either previously neglected higher-dimension operators or light BSM fields to be included in the dimensionally reduced effective theory and therefore necessitate dedicated nonperturbative studies. As a concrete application, we analyze the real singlet-extended Standard Model and identify regions of parameter space with single-step first-order transitions, comparing our findings to those obtained using a fully perturbative method. We discuss the prospects for exploring the electroweak phase diagram in this model at collider and gravitational wave experiments in light of our nonperturbative results.


I. INTRODUCTION
The nature of electroweak symmetry breaking in the early Universe has important implications for cosmology and particle physics. An electroweak phase transition (EWPT) could have given rise to the observed matter/ antimatter asymmetry of the universe [1], significantly affected the abundance of cosmic relics [2], or produced a stochastic background of gravitational waves (GWs) [3]. In the pure Standard Model (SM), electroweak symmetry breaking is expected to have occurred at a crossover for the observed value of the Higgs mass [4][5][6][7]. However, many of its extensions predict much stronger transitions, and in particular first-order phase transitions that may have sourced gravitational waves that can be observed at present day or future experiments [3,8].
In light of the successes of the LIGO and VIRGO collaborations [9,10], as well as the LISA Pathfinder spacecraft [11], a number of space-based gravitational wave missions are now planned or in preparation. Most notable of these is LISA [12], due to launch before 2034; others, such as DECIGO [13] and BBO [14], may follow. These space-based interferometers will have arm lengths long enough to search for a stochastic GW background produced during the epoch of electroweak symmetry breaking in the early Universe, and therefore at a firstorder electroweak phase transition. There are also many future collider projects currently being discussed that have the potential of probing the nature of the electroweak phase transition. In addition to the high-luminosity LHC (HL-LHC) [15], these include proposals for e þ e − colliders such as the ILC [16], FCC-ee [17], CEPC [18] and CLIC [19], as well as higher-energy proton-proton colliders such as the High-Energy LHC (HE-LHC) [20], FCC-hh [21] and SPPC [22]. Various measurements at these machines will provide valuable insight into the nature of electroweak symmetry breaking and complement the observations at gravitational wave experiments.
Given the rich experimental program aimed at exploring the electroweak phase transition, it is crucial to accurately characterize and make reliable predictions for the nature of the EWPT in extensions of the Standard Model. Due to the infamous Linde problem, perturbative methods in thermal field theory are known to suffer from severe infrared divergences near phase transitions [23,24]. Furthermore, perturbative studies are often plagued by issues related to gauge dependence [25,26]. Establishing the existence of a strong first-order EWPT, and a precise determination of thermodynamic quantities related to it, therefore requires nonperturbative methods. Such an approach is crucial for providing robust phenomenological predictions for gravitational wave and collider experiments.
A particular example of the limitations associated with perturbative calculations is the evaluation of the bubble nucleation rate, especially when the energy barrier between the vacua is radiatively induced. In this case thermal fluctuations are responsible for creating the barrier between the symmetric and broken phases (see Ref. [27] for a discussion). For such transitions, the perturbative calculation of bubble nucleation requires a separation of scales between those fluctuations which are responsible for the energy barrier and those which describe the bubble nucleation [28][29][30]. One has to tread very carefully in order to avoid double counting fluctuations [31][32][33], to selfconsistently account for the spacetime dependence of the bubble [34,35] and to maintain gauge invariance [25,26,36,37]. Furthermore, explicit calculations taking into account the necessary separation of scales have shown a strong dependence on the precise matching scale, perhaps suggesting a breakdown of the usual semiclassical methods [38]. Beyond the leading order, calculating corrections to the bubble nucleation rate is a formidable task involving the evaluation of functional determinants. In practice, these issues are rarely addressed by bubble nucleation calculations in the literature. Instead practical shortcuts are adopted which, though not rigorously theoretically justified, are hoped to account for the relevant physics and provide the correct nucleation rate to within order of magnitude. For these reasons, it is important to know how well the usual approach performs quantitatively. This can only be answered by comparison to a fully self-consistent, nonperturbative approach [39]. In the paper at hand, we perform such a benchmark comparison.
Even in the effective 3-d description, a lattice study of the full parameter space in BSM scenarios is computationally challenging. An economic alternative to dynamical simulations becomes possible if the new fields are sufficiently heavy, so that their effect on the dynamics of the transition is suppressed. The BSM degrees of freedom (d.o.f.) can then be integrated out in the process of dimensional reduction to obtain a "SM-like effective" 3-d theory whose phase structure is well understood from the earlier lattice studies of Refs. [4,6]; this strategy was first adopted in Ref. [47]. Although computationally inexpensive, this approach comes with its own limitations. One cannot study interesting multistep transitions [84,85] where new BSM fields can be light near the transition, nor even one-step transitions in which the new light BSM field plays an active, dynamical role. Furthermore, integrating out the BSM fields can lead to reduced accuracy when studying the dynamics of the transition, since the effects of higher dimensional operators can be large but have been neglected in all nonperturbative studies of the 3-d infrared effective theory (cf. Ref. [59]). Nevertheless, before performing new simulations with additional dynamical BSM fields or higher dimension operators, it is informative to investigate how much the existing nonperturbative results can tell us about the theory.
One objective of this study is to extract as much modelindependent information from existing lattice results as possible and apply them to gravitational wave predictions. The 3-d infrared SM effective field theory (EFT) is particularly simple and we are able to determine the range of phase transition parameters relevant for computing the GW spectrum accessible by this particular effective theory, in which effects from higher-dimension operators and light BSM fields are negligible. We find that a first-order EWPT in any BSM theory that can be described by the minimal SM-like 3-d EFT upon integrating out the additional fields will lead to signatures that are undetectable at any planned gravitational wave detector. That is, one cannot get a sufficiently strong and long-lasting transition from any theory which looks SM-like in the infrared. As a result, collider experiments are likely to provide the most sensitive probe of the phase diagram in such scenarios. This also means that future lattice studies incorporating the effects of higher dimension operators or light BSM fields will be required in order to make theoretically sound predictions for gravitational wave experiments.
As a concrete application, the latter portion of this paper fleshes out these results in a specific BSM scenario, the real singlet-extended Standard Model (xSM) [73,[86][87][88][89][90], for which the required dimensional reduction has already been done in Ref. [57], but a systematic analysis of the 5dimensional parameter space was not performed. We analyze the phenomenologically interesting parameter space with a sufficiently heavy singlet, and make use of the nonperturbative results of Ref. [4] to deduce the phase diagram of this model, comparing with the results of a traditional perturbative treatment as in Ref. [78]. We present the first nonperturbative results for the gravitational wave power spectrum in the xSM with the level of rigor that can be achieved with existing lattice results. As expected, the corresponding transitions are too weak and cannot be realistically probed by gravitational wave experiments. Nevertheless, we use these results to perform a comparison between perturbative treatments and the nonperturbative predictions, finding reasonable agreement within systematic uncertainties for the various thermodynamic quantities relevant for the gravitational wave spectrum. While a future nonperturbative study of this model will be required to make firm predictions for gravitational wave experiments, our analysis of the phase diagram provides robust targets for collider searches, which will have a realistic opportunity to probe the first-order regions accessible by existing lattice results.
This article is organized as follows. In Sec. II we review the setup for recycling pre-existing nonperturbative results for determining the phase structure of a given BSM scenario. In Sec. III we discuss in detail how to use the existing nonperturbative results to predict gravitational wave signals from first-order electroweak phase transitions in general BSM scenarios that map on to the infrared 3-d SM-like EFT. In Sec. IV we focus on the real singlet extension of the Standard Model and present our nonperturbative results for the phase diagram of the model, comparing against predictions obtained from perturbation theory. Section V discusses collider and gravitational wave probes of the phase diagram in this theory, applying the results of Secs. II-III. Finally, we conclude in Sec. VI.

II. DIMENSIONAL REDUCTION AND NONPERTURBATIVE DETERMINATION OF THE PHASE DIAGRAM
The technique of dimensional reduction allows us to determine the phase structure of a given model, and in particular establish the existence of a first order phase transition, in a theoretically sound manner. The basic idea of dimensional reduction is that at high temperatures, longdistance physics decouples from the temporal direction, and the equilibrium properties of the system can then be described by a spatial 3-d theory. In the Matsubara formalism, this can be understood as a consequence of a thermal scale hierarchy generated by the heat bath. Namely, all nonzero Matsubara modes come with an effective mass correction proportional to πT-a scale dubbed superheavy-and their effect is thus suppressed compared to modes with zero Matsubara frequencies. For more details regarding dimensional reduction, see, e.g., Refs. [40,50,91] and [56,58].
To illustrate the process of dimensional reduction in an extended scalar sector of the SM, we start with a fourdimensional (4-d) Euclidean Lagrangian in the imaginarytime formalism where the 4-d fields depend on the imaginary time τ and spatial coordinates ⃗ x. The SM fields ϕ, A μ and ψrepresenting the Higgs doublet, the gauge fields and the fermions, respectively-are accompanied by new scalars that we divide into two categories: by assumption the field S is superheavy (m S ∼ πT) and the scalar s is heavy, meaning its mass is suppressed in a power counting sense (m s ∼ gT). We emphasize that here m S and m s are running masses in e.g., the MS scheme, and not physical masses. All fermionic modes, and bosonic hard modes, i.e., those with nonzero Matsubara frequency, have masses that are superheavy. Integrating out these superheavy field modes (including the zero-mode of the scalar S) produces a bosonic 3-d theory with the Lagrangian where the purely spatial 3-d fields (that correspond to zero modes of original fields) obtain masses that are generally smaller than the superheavy scale. The effects of the nonzero modes are captured by the parameters and field normalizations of the effective theory, which are functions of the temperature and physical quantities of the original 4-d theory. Note that since the heat bath breaks fourdimensional Lorentz invariance, the temporal components of the gauge fields A 0 are treated as Lorentz scalars in the effective 3-d theory. Thermal masses generated for these gauge field zero-modes are of order gT, i.e., at the heavy scale. In contrast, the spatial components of the non-Abelian gauge bosons remain massless. The scalar fields undergoing a phase transition can also become very light, with masses formally of order g 2 T, due to a cancellation between vacuum masses and thermal corrections. In the scenarios of interest below, it is assumed that only the Higgs doublet ϕ becomes light near the critical temperature. The EWPT is then driven by the Higgs field while the scalars S, s are spectators. In order to affect the transition dynamics, the new scalars need to have sufficiently strong interactions with Higgs. In practice, fields belonging to the heavy scale can be integrated out as well, in the conventional effective field theory sense. This results in a still simpler 3-d theorȳ where the physics of the heavy scale is captured in the parameters and field normalizations of the final light scale effective theory. The steps of dimensional reduction are illustrated in Fig. 1. In the end, one arrives at an effective theory solely for the infrared-sensitive zero modes. Therefore the procedure of constructing the effective theory, while perturbative in nature, is completely infrared-safe and can be understood as a sophisticated method for performing thermal resummation [43,45] (see also Sec. 6 of Ref. [91] for a more pedagogical discussion of this point). At the light scale, the physics of the phase transition becomes inherently nonperturbative due to infrared divergences in the effective expansion parameter. The theory can, however, be studied on the lattice. The end result of dimensional reduction for theories with only heavy and superheavy BSM physics is therefore a set of matching relations between the renormalized parameters of the theory in Eq. (1) and a 3-d theory containing the spatial SU(2) gauge field and a doublet Higgs field, with gauge couplingḡ 3 . Note that nonperturbative effects of the Uð1Þ Y gauge group were omitted in the original study [4], and we follow this same approach. Equation (4) is the same 3-d theory to which the Standard Model was mapped in Ref. [40], and then studied nonperturbatively in Refs. [4,5]. Furthermore, the nucleation rate was studied nonperturbatively in Ref. [39]. Remarkably, the theory described by Eqs. (3)-(4), which we subsequently refer to as the 3-d SM-like EFT, is universal in the sense that it can be used to describe all extensions of the SM as long as the new d.o.f. remain heavy near the transition after thermal mass corrections are taken into account. Furthermore, the effects of the neglected higher dimensional operators, such as ðφ † 3φ3 Þ 3 , should remain small. In many cases, BSM interactions with the Higgs doublet generate operators of dimension six and higher (counting dimensions as in 4-d) that can have a significant impact on the transition dynamics (cf. Ref. [55]), and so dropping these operators results in reduced accuracy. In our results below, we will specify where we expect these operators to become important.
In scenarios described by the 3-d SM-like EFT in the infrared, the dynamics are fully determined by only two parameters, x ≡λ h;3 g 2 3 and y ≡μ which are, as a result of dimensional reduction, functions of the temperature and couplings of the full four dimensional theory. The tree-level value for the critical temperature T c is given by the line y ¼ 0. The true, fully nonperturbative critical line can be obtained from a Monte Carlo simulation and turns out to be rather close to the tree-level approximation [4]. Furthermore, the condition y ≈ 0 requires the thermal correction to the doublet mass parameter to effectively cancel its vacuum mass. This therefore fixes the critical temperature in our analysis to be of the same order of magnitude as the physical Higgs mass (this also approximately sets the scale at which gravitational waves will be emitted). The strength of the transition is then controlled by the parameter x, and the lattice study of Ref. [4] showed that for 0 < x < 0.11 the transition is of first-order, and becomes a crossover for x > 0.11 [92]. Note that if, for a given choice of 4-d parameters, x becomes negative, the tree-level potential of the effective theory is unbounded from below and cannot be used for simulations. This generally signals that the dynamical effects from the BSM fields are too large to be neglected, or that the aforementioned higher-dimensional operators should be included in the effective theory. Summarizing, one can use existing nonperturbative results to study the electroweak phase diagram in BSM scenarios by performing the following steps: (1) Dimensionally reduce the full 4-d theory to the infrared 3-d theory.
(2) In the parameter space of the theory that can be accurately described by the 3-d SM-like EFT (i.e., where there are no additional light fields or important higher-dimension operator effects), compute x and y. These quantities fully determine the 3-d dynamics.
In this parameter space, nonperturbative simulations predict a first-order electroweak phase transition. We stress again that the procedure above applies to transitions driven by the Higgs field and for which the heavy BSM fields are nondynamical. Extending our results to dynamical BSM fields requires retaining these fields in the effective 3-d theory; this, in turn will require fresh nonperturbative simulations in the future.
The simple EFT described by Eqs. (3) and (4) can also be used to predict quantities relevant for the gravitational wave power spectrum arising from first-order transitions. In the following section, we describe how this can be done.

III. THE GRAVITATIONAL WAVE POWER SPECTRUM
At cosmological first-order phase transitions, excepting the case of very large supercooling, gravitational waves are primarily produced by the collisions of bubbles and the subsequent evolution of the resulting fluid sound waves. Numerical simulations of the coupled field-fluid system give the spectrum of gravitational waves as a function of various equilibrium and dynamical properties of the transition [93]. The key quantities that we would like to determine nonperturbatively are the strength of the phase transition, α, the phase transition temperature, T Ã , and the inverse phase transition duration, β=H Ã [3,8]. The latter two can be found once one knows the bubble nucleation rate, Γ. In this section we review how this can be done making use of existing 3-d lattice results. There is one additional relevant quantity: the bubble wall speed, v w , for which no nonperturbative studies exist. In our study, we treat it as an input parameter.

A. Dimensional reduction for bubble nucleation
Given a theory that has been dimensionally reduced as in Sec. II and predicting a first-order electroweak phase transition, we can also study bubble nucleation in the 3-d theory. The philosophy here is the usual one of dimensional reduction: 4-d quantities are computed by matching to the analogous quantities computed in the 3-d effective theory, Eq. (4). Equilibrium properties of the 3-d dimensionally reduced theory have been studied nonperturbatively in Ref. [4] and in Ref. [39] the rate of bubble nucleation was determined. This allows an end-to-end nonperturbative description of the phase transition. The relevant results are collected in Tables I and II. Following Ref. [8] we define the phase transition strength α as the ratio of the latent heat of the transition to the radiation energy density in the symmetric phase at the time of transition where LðTÞ is the latent heat and gðTÞ the number of relativistic d.o.f. When BSM scalar eigenstates are sufficiently heavy, they are not relativistic d.o.f. at the EWPT and we have gðT ∼ T Ã Þ ≈ 100. Note that the gravitational wave spectrum is determined by αðT Ã Þ, not αðT c Þ. However, only αðT c Þ has been determined nonperturbatively, so we concentrate on this in what follows and approximate αðT Ã Þ ≈ αðT c Þ which is typically justified absent significant supercooling. The relation of the latent heat at T c to its 3-d counterpart, Δl 3 , can be found in Ref. [4] and it reads where we have introduced the functions For a given 4-d theory, these η-functions determine how the couplings of the 3-d effective theory run with temperature.  [4,39,92] (see also Ref. [94]). The entry marked with a dagger is interpolated based solely on values for other x c . The final row shows the point at which the phase transition becomes second order, separating first-order transitions from crossovers. Note that the results in this table do not include the effect of the U(1) field, which makes the phase transition slightly stronger [6]. This moves the critical line endpoint to x c ≈ 0.11 [92].  [39]. In that reference −log Γ 3-d is given in the last column of Table IV. In terms of their notation, ðy − y c Þ ≡ δm 2 =g 4 T 2 . They are analogous to the usual β-functions of quantum field theory, in the sense that they describe the change in parameters with respect to energy scale variations, though their appearance is connected to the process of DR. For a generic BSM theory, thermal corrections to couplings arise first at one-loop. For y, this corresponds to thermal corrections at Oðg −2 Þ, the thermal mass corrections, whereas for x this leads to thermal corrections first at Oðg 2 Þ. Thus the temperature dependence of x is significantly weaker than that of y, whereas dy c =dx c ¼ Oð1Þ, hence In order to calculate the phase transition temperature and inverse phase transition duration, one needs to know the rate of bubble nucleation, Γ. The 4-d nucleation rate is related to its 3-d counterpart by [39], where σ el ∼ T=logð1=gÞ is the non-Abelian "color" charge conductivity, arising in the effective description of the realtime dynamics of the long wavelength d.o.f. [95,96]. The powers of g 2 =ð4πÞ in the prefactor arise due to the relation between dimensionally reduced lattice quantities and physical quantities. Finally, Γ 3−d ðx; yÞ is the (dimensionless) rate of bubble nucleation as calculated on the lattice in the 3-d effective theory. It depends only on the 3-d parameters x and y [see Eq. (5)]. The phase transition temperature, T Ã , may be defined to be the temperature at which only a fraction 1=e ≈ 0.37 of the universe remains in the symmetric phase. At this point the following equality holds [97] where H is the Hubble parameter. Note that this condition is different from the one-bubble-per-Hubble-volume condition, which determines the very beginning of bubble nucleation [97,98], yielding what is often termed the nucleation temperature T n . However, it is at the later time corresponding to T Ã , when the universe contains many bubbles, that bubble collisions occur, imprinting their length-scale on the fluid sound waves and the resulting gravitational wave spectrum. This distinction can lead to a factor of 2 difference in β=H Ã and hence is important if one wants to make accurate predictions.
The inverse duration of the phase transition is then determined by where all quantities are evaluated at the transition temperature. 1 On the last line we have used Eq. (9). From semiclassical calculations of the nucleation rate, one also expects the x dependence of the rate to be weaker than the y dependence, making the term we have kept even more dominant. In Ref. [39], the quantity Γ 3−d ðx; yÞ was calculated for x ¼ 0.036 and is reproduced here in Table II. It is given as a function of ðy − y c Þ, corresponding to δm 2 =g 4 T 2 in the notation of Ref. [39]. Fitting a straight line to their results, one finds for ðy − y c Þ ∈ ½−0.01057; −0.00835, a range which is sufficient to incorporate the value of y at the transition temperature, T Ã . Thus, for this value of x, Details of the specific 4-d theory only enter these expressions through η y . Interestingly, by combining Eqs. (15) and (16), one finds the relation where we have used η y ðT Ã Þ ≈ η y ðT c Þ. This defines a line in the ðαðT c Þ; β=H Ã Þ plane on which lie all 4-d theories which reduce to the SM-like 3-d EFT with x ¼ 0.036. A given 4-d theory's position along the line is determined solely by η y ðT c Þ, An analogous relation will hold for other values of x: the ratio of β=H Ã to α is strikingly independent of all shortrange 4-d physics. Physically though, the dominance of 3-d, long-range physics is no surprise. The tree-level 3-d theory is the result of integrating out the superheavy and heavy modes of the 4-d theory. However, the 3-d theory 1 Note that in Ref. [39] the factor of η y was mistakenly dropped, leading to their result for β=H Ã being a factor of ∼4 too small. does not contain a tree-level barrier. Hence one must integrate out some light d.o.f. of the 3-d theory to give the energy barrier between the vacua. The d.o.f. describing the bubble nucleation are then necessarily lighter still. Thus, if any 4-d theory reduces to this 3-d theory, the dynamics of its transition are largely independent of the original 4-d theory, excepting the value of x c to which it maps.
Note that we generically find rather large values for β=H. This is completely determined by the long-distance, 3-d physics, dictated by the large numerical coefficient of η y in Eq. (15), also noting that η y ¼ Oðg −2 Þ. One can trace the large values for β=H to the extremely narrow metastability range, Δy, of the symmetric phase. The lower spinodal decomposition point, at which point the symmetric phase becomes absolutely unstable, is at a value of y ¼ y c − Δy only slightly below the critical one, Δy ≪ 1. This was studied in Sec. 7 of Ref. [4], in which it was found that Δy ≈ 0.007, for the case x ¼ 0.06444 (note, in Ref. [4] the result is written in terms of the "temperature," T Ã ). If the universe were cooled suddenly to y ≤ y c − Δy, spinodal decomposition would take place: without a barrier the transition to the broken phase would take place essentially instantly. Thus the nucleation rate goes from zero at y ¼ y c to taking place in a microscopic timescale at y ¼ y c − Δy. Hence, written as a function of y=Δy, the logarithm of the rate goes from minus infinity to Oð1Þ as its argument changes by unity. Thus, is large, both multiplicative terms on the right-hand side being large. Perturbatively, we can go a step further and trace the narrowness of the metastable region to the radiatively induced nature of the transition. This implies that, near the critical point, the 3-d Higgs mass-squared parameter and radiative corrections are of the same order-ofmagnitude. As radiative corrections are suppressed by loop factors, 1=ð16πÞ, the transition takes place when the 3-d Higgs mass-squared parameter is very small, y c ≪ 1. Further, at leading perturbative order, the lower spinodal decomposition point is at y ¼ 0, when the 3-d Higgs mass-squared parameter is zero. Thus Δy ≪ 1, and the above argument goes through. Although this argument is perturbative, the nonperturbative results bear it out.

B. Implications for gravitational wave experiments
We are now in a position to make theoretically sound, nonperturbative predictions for the gravitational wave spectrum resulting from a first-order electroweak phase transition in any BSM theory that can be mapped into the EFT given by Eqs. (3)-(4).
The general recipe for this is described as follows, assuming Eq. (9) holds. First one solves the matching relations to find points in the 4-d parameter space that give x c ¼ 0.036, so that we can use the only existing nonperturbative result for the bubble nucleation rate. Then by solving Eq. (12) one can find the phase transition temperature, T Ã , and finally evaluate Eqs. (15) and (16) to give β=H Ã and α. One can then use the fit formula from the numerical simulations of Ref. [93] to determine the GW power spectrum today. Variation of the bubble wall speed affects the signal strength weakly, well within the range of other sources of uncertainty, and so we simply set it to 1 in what follows.
The aforementioned analysis, together with generic field theoretic arguments, yields a well-defined range of the gravitational wave parameters, α, β=H Ã , and T Ã accessible by theories matching on to the 3-d SM-like EFT. We present these results in Fig. 2. The gray shaded region in Fig. 2(a) and the light blue region Fig. 2(b) can be reliably described by the infrared 3-d SM-like EFT. In Fig. 2(a) we specifically state why this EFT description fails in different regions, as detailed below. The black line Fig. 2(b) corresponds to α and β=H Ã values accessible by existing nonperturbative calculations, i.e., x ¼ 0.036. For reference, we also show the expected sensitivity of LISA [12] to this parameter space. In Fig. 2(a), the region predicting a signalto-noise ratio (SNR) of 10 at LISA is indicated in the bottom right corner. Sensitivity curves assuming a ten year mission duration with 75% duty cycle are shown for two cases: one where the sound source is on for the shock formation time, and the other where the sound source lasts the full Hubble time (see Ref. [99] for a more detailed discussion of these distinct cases). In Fig. 2(b), we show additional LISA SNR curves assuming that the sound source lasts for a Hubble time. The dashed contours give the shock formation time, and only in the shaded region is the shock formation time longer than a Hubble time. In both panels we set T Ã ¼ 140 GeV and v w ¼ 1. Neither the sensitivity curves nor the region for which the SM-like EFT are valid are particularly sensitive to T Ã ; the scale is implicitly set by the physical Higgs mass.
The regions accessible by the 3-d SM-like EFT and existing nonperturbative results in Fig. 2 are determined as follows. First let us consider the range of allowed η y . This quantity, together with x c , determines αðT c Þ through Eq. (7) and β=H Ã through Eq. (13) (for ηðT Ã Þ ≈ ηðT c Þ). The value of η y ðT c Þ is generically of order Og −2 Þ, absent new physics with large couplings. Its SM value for a 125 GeV Higgs, η SM y ðT c Þ ≈ 4.4, receives corrections from BSM physics, where Δη BSM y ðT c Þ contains corrections which go to zero in the SM limit. In principle Δη BSM y ðT c Þ can be positive or negative, though all leading-order terms in η SM y ðT c Þ are positive definite.
There is a well-defined range of η y ðT c Þ expected in generic BSM theories mapping on to the 3-d SM-like EFT. First, note that If Δη BSM y ðT c Þ is negative, η y ðT c Þ is reduced with respect to its SM value, leading to weaker transitions with a longer phase transition duration. In special cases such that the SM and BSM contributions cancel, one can in principle get very weak transitions. However, to achieve η y ðT c Þ ≪ η SM y ðT c Þ requires order-by-order cancellations between SM and BSM physics and hence is fine-tuned. Any generic BSM physics which does not lead to such finetuned cancellations at more than the leading order in the loop expansion will satisfy η y ðT c Þ ≳ ðg 2 =4πÞη SM y ðT c Þ ≈ 0.15. Here we have conservatively taken the loop expansion parameter to be ∼g 2 =4π: the presence of any larger couplings would strengthen this approximate lower bound. Furthermore, such small η y ðT c Þ, shown as the bottom left arrow in Fig. 2 ðT c Þ is large, the phase transition duration becomes very short, shown as the top right arrow in Fig. 2(a). This pushes the peak of the gravitational wave spectrum to higher frequencies, which are unfavorable for gravitational wave experiments.
From this reasoning, we conclude that any generic BSM physics which couples perturbatively to the Higgs and does not lead to nontrivial cancellations at more than one loop order will satisfy η y ðT c Þ ∈ ½0.15; 10. Assuming this we find that, for x ¼ 0.036, Eqs. (15) and (16) trace out the thick black line shown in Fig. 2(b).
At other values of x, we do not have a nonperturbative determination of the phase transition parameters. In this case, we use a fully perturbative approach utilizing the DR (3-d PT). In the 3-d PT approach, our calculations are performed at two-loop order and in the Landau gauge, cf. Ref. [56]. In the calculation of Γ 3−d , we follow the perturbative semiclassical analysis outlined in Ref. [39], including the effects of wave function renormalization. This gives the tunneling action accurate to next-to-leading order in x. Note that the semiclassical prefactor contributes at next-to-next-to-leading order in x [31,39].
Larger values of x result in weaker first-order transitions and shorter phase transition durations, until x ¼ 0.11 where the transition is a crossover and hence there is no bubble nucleation. This corresponds to the upper left arrow in Fig. 2(a). Conversely, for smaller values of x, the phase transition becomes stronger and longer in duration, hence smaller x is the interesting region for gravitational wave production [bottom right arrow in Fig. 2(a)]. However, there is a limit to how small one can reliably go in x, as the dimension-6 operators that have been dropped in the DR can become comparable to, or larger than, the dimension-4 term, which is proportional to x. The dominant higher dimensional operator is c 6 ðϕ † 3 ϕ 3 Þ 3 . The leading SM contribution to c 6 is Oðy 6 t =ð4πÞ 3 Þ, where y t is the top Yukawa coupling (see Appendix B). There will also be contributions from any BSM physics that couples to the Higgs. For the effect of the dimension-6 operator to be small, the following strong inequality should hold, where the expectation values refer to the volume averaged condensates, evaluated at the critical point (see, for example, Ref. [45]). Calculating these condensates using the two-loop effective potential, which should be reasonably accurate at small x, and taking c 6 to be its Standard Model value, we arrive at x ≫ 0.01. However, it should be remembered that the precise condition will depend on BSM contributions to c 6 . In practice, such small values of x are likely overoptimistic and encourage us to investigate an extended 3-d EFT including higher dimensional operators in the future. Using the 3-d PT approach to calculate Δl 3 and Γ 3−d allows us to trace out the entire region on the ðα; β=H Ã Þ plane that can be described by the SM-like 3-d EFT in Eqs. (3)-(4), as argued above. Considering x ∈ ½0.01; 0.06 and the aforementioned range of η y ðT c Þ, we arrive at the gray and light blue regions of Figs. 2(a) and 2(b), respectively. Here we have chosen the upper value x ¼ 0.06 since larger values are not accurately described by our 3-d PT approach and our next-to-leading order calculation of the bubble nucleation rate becomes numerically challenging in this regime. Nevertheless values in the range 0.06 < x < 0.11 give only very weak first order phase transitions. Any BSM theory which reduces to the 3-d SMlike EFT and gives a reasonably strong first-order phase transition will lie in the light blue region of Fig. 2(b). The dark blue region in Fig. 2(b) is for the specific case of the singlet extended Standard Model (xSM), which we discuss further in Section VA.
From Fig. 2 we can reasonably conclude that, in the region where the effects of heavy BSM fields on the EWPT can be captured by the 3-d SM-like EFT, a given model does not predict a sufficiently strong transition to produce an observable gravitational wave signal at LISA. In fact, we have verified that no planned gravitational wave experiment is expected to be sensitive to the gray and light blue regions of Figs. 2(a) and 2(b), respectively. In order to explore the phase diagrams of theories mapping onto these regions, other experimental approaches, such as searches for BSM physics at colliders, will be required. We will see how this plays out in the specific case of the real singlet-extended SM below. For stronger transitions as required for an observable GW signal there are two options: either the BSM field must play an active, dynamical role or higher dimension operators in the 3-d EFT must be included. A nonperturbative description will therefore require new 3-d lattice simulations including the dynamical BSM field, or inclusion of higher-dimension operators involving the scalar doublet (see Fig. 2).

IV. APPLICATION TO THE STANDARD MODEL WITH REAL SINGLET SCALAR
Up until this point our analysis has been quite general. As an illustration of our methods, we apply our nonperturbative analysis to the specific case of the oftenstudied real singlet extension of the Standard Model (xSM) [63,64,73,78,[86][87][88][89][90]. Using the dimensional reduction of Ref. [57], in this section we first investigate the 5dimensional parameter space of the model to determine the phase structure in the regions where the singlet is sufficiently heavy to be integrated out. We will then compute the gravitational wave spectrum nonperturbatively in Sec. V, comparing with perturbative results, and discuss the prospects for exploring the phase diagram of this scenario at collider experiments.

A. Parametrization of the model
The most general scalar Lagrangian involving an additional real singlet scalar field can be written as [63,64,78,86,100] where D μ is the usual covariant derivative for the Higgs doublet, and parameters are renormalized in the MS scheme. In order to relate these MS parameters to physical quantities, we reparametrize the potential as follows. After electroweak symmetry breaking, we assume that the singlet field σ has a zero-temperature VEV hσi ¼ 0, and that the Higgs vev is hϕ † ϕi ¼ v 2 =2. 2 This leads to the following relations in the electroweak vacuum: Furthermore, by diagonalizing the potential and solving for the mass eigenvalues-with the eigenstates denoted h 1 and h 2 -the parameters λ, μ and a 1 can be eliminated in favor of physical masses m 1 , m 2 and the mixing angle θ. We assume that h 1 is the measured Higgs boson with mass m 1 ¼ 125 GeV. Finally, instead of the quartic portal coupling a 2 , we take the trilinear ðh2; h2; h1Þ coupling λ 221 as an input parameter, which is one of the couplings relevant for scalar pair production processes at colliders [68,78,100,101]. We treat the singlet self-couplings b 3 and b 4 , given at a fixed MS scale, as input parameters. As a subset of the general model parameter space, we will also consider the case in which σ → −σ under a discrete symmetry so that a 1 , b 1 , b 3 , and sin θ are all set to zero in a technically natural way. We refer to this as the Z 2 -symmetric xSM.
In scalar extensions of the SM, first-order phase transitions are often associated with large portal couplings to the Higgs doublet. We therefore expect zero-temperature effects from vacuum renormalization to have a significant effect on the properties of the EWPT. To account for these corrections, we employ a one-loop renormalization procedure in the T ¼ 0 vacuum similar to that of Refs. [40,58,102]. In short, the running parameters are related to observables at an initial scale, which we take to be an average of the scalar masses, by requiring that the physical masses match poles of the loop-corrected propagators. Technical details of this calculation are presented in Appendix A.

B. Identifying regions with a first-order transition
We are interested in parameter space across which the singlet can be reasonably approximated as superheavy so that the dimensional reduction of Ref. [57] and the results of Secs. II-III can be applied. The singlet field in this case is nondynamical at the transition and integrated out. Since this involves effectively replacing the singlet background field by a solution to its equation of motion (EOM), transitions between different singlet vacua will not be captured by this approach, and so in the superheavy regime one is restricted to considering one-step transitions across which the singlet VEV does not change significantly. This approach is therefore expected to capture the "weakest" transitions (i.e., those with the smallest latent heat or order parameter). As singlet-induced higher dimension operators are required to be negligible in order to match on to the 3-d SM-like EFT in the infrared, the dominant effect strengthening the electroweak phase transition will be the reduction of the effective Higgs quartic coupling, with the barrier separating the phases generated radiatively in the EFT. This corresponds to the same mechanism responsible for driving the EWPT first-order in the Standard Model by reducing the Higgs mass. The reduction of the effective quartic coupling is primarily due to tree-(loop-) level effects in the non-Z 2 (Z 2 ) cases, respectively [63]. This can be seen explicitly in the expression in Ref. [57] for the 3-d Higgs quartic coupling after integrating out the superheavy d.o.f., λ h;3 , in terms of the model parameters: The presence of the Z 2 -breaking couplings a 1 , b 1 , b 3 lead to a tree-level reduction of λ h;3 , whereas absent these terms the leading effect starts at Oðg 4 Þ in the power counting scheme used in the dimensional reduction of Ref. [57]. The terms on the right-hand side (RHS) of Eq. (23) can significantly reduce the 3-d quartic coupling near the transition temperature relative to the SM value, λ SM 3;h , with additional singlet contributions to the potential maintaining the correct physical Higgs mass at T ¼ 0. If we were to integrate the singlet out at zero temperature, these additional contributions to the potential not captured by Eq. (23) would manifest as higher-dimension operators involving h, and would necessarily be important in order to maintain m h ¼ 125 GeV with a reduced quartic. On the other hand, when integrating out the singlet at finite temperature, we require higher-dimension operator effects to be small [see Eq. (20)] in order to match onto existing lattice results (see Eq. (B3) for an explicit expression for the h 6 operator). Both requirements can in principle be satisfied simultaneously, since at high temperatures, the range of field values relevant for the transition decreases (the VEV in the broken phase becomes smaller as the temperature increases, since μ 2 h;3 decreases) and numerically the higher dimension operators become less important in determining quantities of interest at the transition. Nevertheless, for a fixed zero temperature Higgs mass, there is generally tension between requiring λ 3;h ≪ λ SM 3;h and the effects of higher-dimension operators to be small near the transition temperature. This tension will be evident in our results: in most cases of interest, the size of the dim-6 operator is at least an order of magnitude larger than its SM value. New nonperturbative analyses will be required in order to address this issue across the full parameter space of the model. In the interim, however, we proceed as far as we can using existing results, noting this persistent tension.
We would like to use our nonperturbative framework to ascertain where in singlet model parameter space the aforementioned effects are large enough to drive the EWPT first-order. Let us first consider the Z 2 -symmetric limit of the model. In this case, sin θ ¼ 0, b 3 ¼ 0, and we can simply vary m 2 and λ 221 . We fix b 4 ¼ 0.25, as the properties of the transition are not very sensitive to singlet selfinteractions. Our nonperturbative results are shown in Fig. 3. A first-order electroweak phase transition is predicted in the light green region. In the darker green region, the 3-d parameter x < 0, so that higher-dimensional operators are required to probe one-step transitions nonperturbatively. We also show where the neglected singletinduced dimension-6 operator coefficient, c 6 , becomes larger than 40 times the corresponding SM value, c 6;SM , assuming T ≈ 140 GeV. An approximate expression for c 6 can be found in Eq. (B3). As discussed in Appendix B, for c 6 ≳ 40c 6;SM the neglected higher dimension operators are expected to have Oð10%Þ or larger effects on the 3-d Higgs VEV in the first-order transition regions (in 3-d perturbation theory using Landau gauge), and so our nonperturbative analysis already becomes less reliable for λ 221 above this contour. In the gray shaded region, the singlet mass parameter μ 2 σ < 0 and the superheavy dimensional reduction completely breaks down. We also show a contour for which the singlet mass parameter μ σ ¼ 200 GeV, above which our assumption μ σ ∼ πT is already questionable.
As can be seen from Fig. 3, most of the first-order phase transition parameter space accessible to current nonperturbative studies lies in the region where the DR and neglect of higher dimension operators is not well justified. For larger masses the superheavy approximation improves, however the portal couplings required for a first-order transition becomes larger for heavier singlets and therefore higherloop zero-T effects become non-negligible. Below the μ σ ¼ 200 GeV and c 6 ≳ 40c 6;SM contours, the superheavy singlet DR and neglect of higher-dimension operators is justified, and a first-order electroweak phase transition is robustly excluded by our nonperturbative results for small λ 221 . This validates the conclusions drawn from perturbation theory in, e.g., Refs. [68,101]. Note also that similar results for the phase diagram in the real triplet-extended SM were obtained in Ref. [58].
We now move on to the more general case without a Z 2 symmetry. We first fix b 3 ¼ 0 and vary sin θ and λ 221 . Our results are shown in Fig. 4 for m 2 ¼ 240 GeV and m 2 ¼ 400 GeV. In the light green regions our analysis predicts a first-order electroweak phase transition, while the darker green region features x < 0, and therefore require the ðϕ † ϕÞ 3 operator in the 3-d EFT to explore one-step phase transitions. We also show where c 6 ¼ 40c 6;SM , assuming T ¼ 140 GeV. To the right of the dotted green contour, singlet-induced dimension-6 effects can already be non-negligible. The gray region features μ 2 σ < 0 for which the superheavy dimensional reduction completely breaks down. Also shown is a contour indicating where μ σ ¼ 200 GeV. In the m 2 ¼ 240 GeV case the superheavy approximation can lead to inaccuracies in determining the first-order phase transition parameter space. For m 2 ¼ 400 GeV the superheavy approximation is better justified, and we expect our nonperturbative analysis to be more accurate. Corresponding results for the parameter space with varying b 3 are shown for m 2 ¼ 240, 400 GeV and sin θ ¼ 0.05, 0.2 in Fig. 5. In these slices, portions of the gray shaded regions are excluded by 1-loop absolute vacuum stability, as discussed in Ref. [78] (see also Ref. [103]).

C. Comparison with perturbation theory
It is useful to compare the results above to those obtained by the conventional perturbative approach. In a perturbative setting, one typically analyzes the phase structure of the theory using the finite temperature effective potential, V eff . To one loop order, and in a mass-independent renormalization scheme, it is given by FIG. 3. Phase structure of the xSM in the Z 2 symmetric limit. Our nonperturbative approach predicts a first-order electroweak phase transition in the light green region. The darker green shaded region features x < 0 so that the higher dimension operator ðϕ † ϕÞ 3 must be kept in the dimensionally reduced theory to resolve one-step transitions. Furthermore, above the dotted green contour, the dimension-6 operator coefficient, c 6 , becomes larger than 40 times the corresponding SM value, 40jc 6;SM j ¼ 0.067. In Section B we find that in this case the neglected ðϕ † ϕÞ 3 term can already have significant effects. The gray shaded region corresponds to where the singlet mass parameter μ 2 σ < 0 so that the superheavy dimensional reduction completely breaks down. The parameter space in which 4-d perturbation theory predicts a first-order transition with v c =T c ¼ 0.3-0.6 is shaded orange. We also show contours indicating the size of the singlet MS mass parameter μ σ . When μ σ ⪅ gT ≈ 100 GeV, the superheavy singlet approximation becomes compromised. A first-order EWPT is robustly excluded for small values of λ 221 , where our nonperturbative treatment is well justified.
Here the sums run over all species coupled to the scalar fields, with upper (lower) signs for bosons (fermions), n i are the number of d.o.f. for the species i, c i are schemedependent constants, μ R is the renormalization scale, and the J AE encode the thermal corrections (see e.g., Ref. [104] for full expressions). To improve the perturbative expansion away from the symmetric phase, one can insert thermally corrected masses for the m i above. This procedure resums so-called daisy diagrams and delays the infrared breakdown of perturbation theory to smaller field values. Typically a high-temperature expansion is made for the thermal selfenergies in this step, although there has been recent progress in going beyond this approximation [72]. In what follows, we use the high-T thermal masses, the expressions for which can be found in, e.g., Refs. [78], as this is the most conventional approach. The effective potential is computed in Landau gauge utilizing the MS scheme.
With V eff defined above, we search for first-order electroweak phase transitions in the same way as that detailed in Ref. [78], tracing the minima of the finitetemperature effective potential up in temperature until the electroweak minimum is degenerate with a vacuum with restored electroweak symmetry. This defines the critical temperature, T c , and critical field value, v c . It is important to note that this approach yields gauge-dependent results for the critical temperature and order parameter of the phase transition. Reference [25] provides a systematic method for avoiding the spurious gauge dependence in determining T c . It would be interesting to perform a comparison of the results obtained by this method with our (manifestly gaugeinvariant) nonperturbative approach in the future. However, in our present study we restrict ourselves to the conventional (gauge-dependent) method described above.
In the perturbative approach, a thermal cubic term is always present in the finite temperature effective potential. This means that perturbation theory always predicts a firstorder transition. This is a spurious effect due to the IR breakdown of the perturbative expansion away from the broken phase and prevents one from resolving the endpoint of the transition. Therefore, to extract meaningful information about the phase structure of the theory in this approach, one must restrict themselves to considering sufficiently "strong" first-order transitions. The strength of the phase transition can be parametrized by the order parameter v c =T c . Conventionally, strong first-order transitions are defined by imposing a minimum value for the order parameter, v c =T c ≥ ζ. The value of ζ depends on the particular question of interest. In the setting of electroweak baryogenesis, it is determined by requiring sufficient sphaleron suppression in the broken phase, typically The darker shaded region features x < 0 so that the higher dimension operator ðϕ † ϕÞ 3 must be kept in the dimensionally reduced theory to resolve one-step electroweak transitions. Furthermore, to the right of the dotted green contour, the dimension-6 operator coefficient, c 6 , becomes larger than 40 times the corresponding SM value, c 6;SM , and the neglected ðϕ † ϕÞ 3 term can already have significant effects. We also show contours indicating the size of the singlet MS mass parameter μ σ . When μ σ ⪅ gT ≈ 100 GeV, the superheavy singlet approximation becomes compromised. The gray shaded regions correspond to where the singlet mass parameter μ 2 σ < 0 so that the superheavy dimensional reduction completely breaks down. The parameter space in which 4-d perturbation theory predicts a first-order transition with v c =T c ¼ 0.3-0.6 is shaded orange, and matches up well with the nonperturbative first-order regions especially away from the Z 2 limit, where the couplings required for a first-order EWPT are not as large.
corresponding to ζ ≈ 0.6-2 (see, e.g., Refs. [25,104] for a more detailed discussion). From the standpoint of observable gravitational waves, ζ typically needs to be larger, although the strength of the phase transition in this case is usually parametrized by the more physical quantity α. For the purpose of benchmarking perturbation theory, we consider values of ζ that most closely reproduce the first-order transition region accessible by our nonperturbative analysis and in which the various assumptions of Secs. II-III are justified. This will allow us to compare how the phase structure of the theory depends on the various model parameters in the perturbative and nonperturbative approaches. In practice this amounts to considering ζ ∈ ½0.3; 0.6: ð25Þ Such transitions are quite weak from the standpoint of sphaleron suppression and gravitational wave generation (as we will see below) but serve to delineate the regions of parameter space in which current nonperturbative studies concretely indicate a genuine first-order electroweak phase transition.
The region with v c =T c ∈ ½0.3; 0.6 in 1-loop 4-d perturbation theory is indicated by the shaded orange region in Fig. 3 for the Z 2 -symmetric case. The couplings required for a "strong" first-order EWPT in the 4-d approach tend to be somewhat larger than those corresponding to first-order EWPTs in our nonperturbative analysis. However, it is important to note that the superheavy approximation is not justified across most of the first-order PT region in Fig. 3, and that the neglected dimension-6 operators can be important here, as discussed above. Our 4-d treatment retains the singlet field in the effective potential and therefore accounts for the effects corresponding to the higher-dimension operators. Future nonperturbative studies including these effects will be required to conclusively benchmark the predictions from perturbation theory in this limit.
The especially away from the Z 2 limit, where the couplings required for a first-order EWPT are not as large. The parametric dependence of the strong first-order region predicted by 4-d perturbation theory reliably traces that of the first-order EWPT regions obtained by our nonperturbative analysis in portions of the parameter space where our nonperturbative treatment is well-justified. Our results suggest that, in the particular 4-d approach considered, a genuine first-order EWPT tends to correspond roughly to v c =T c ≳ 0.3. It should be emphasized that this conclusion is gauge-dependent and not necessarily true for other gauge choices.

V. EXPLORING THE PHASE DIAGRAM AT GRAVITATIONAL WAVE AND COLLIDER EXPERIMENTS
Having determined the singlet model parameter space where the electroweak phase transition goes from a thermodynamic cross-over to a first-order transition, we would like to consider the prospects for studying the nature of the electroweak phase transition experimentally in this scenario. Below, we discuss gravitational wave and collider probes of the first-order EWPT parameter space obtained from the nonperturbative methods described in Sec. III. While the transitions accessible by the repurposed nonperturbative results are quite weak and therefore challenging to detect at experiments like LISA, it is nevertheless informative to compare the results from the 3-d and 4-d approaches. This serves as a benchmark of perturbation theory, and indicates the levels of uncertainty inherent in perturbative estimates of the resulting gravitational wave signal.

A. Benchmarking for gravitational wave predictions
For a nonperturbative determination of the gravitational wave power spectrum in the xSM, we follow the general recipe outlined in Sec. III B. We show the values of T c , T n , αðT c Þ, and β=H Ã calculated using our nonperturbative approach (NP) for a representative sample point in Table III. Also shown are the same quantities calculated using the usual 4-d perturbative approach (4-d PT)-see Section IV B-as well as our two-loop, fully-perturbative approach utilising the DR (3-d PT)-see Section III B.
A comparison of the nonperturbative and perturbative results in Table III shows agreement on the order of magnitude of all quantities. This provides an important validation for the usual 4-d perturbative methods, at least in this region of parameter space, and a benchmark of the numerical importance of underlying uncontrolled approximations [23][24][25][26][31][32][33][34][35][36][37][38]. There is nevertheless a 10% discrepancy in T c , a factor of 3 discrepancy in αðT c Þ and a 30% discrepancy in β=H Ã . However, the ratio T n =T c agrees very well, to better than 1%. We have verified that similar trends apply to other points as well in the relevant parts of the xSM parameter space.
Due to the Linde problem (the breakdown of perturbation theory), it is not possible to make a fully reliable error estimate within perturbative theory. For our nonperturbative approach however, this is possible, and in Appendix B we discuss the various sources of uncertainty and their expected sizes. The main sources of uncertainty in the nonperturbative calculation are due to renormalization scale dependence and the neglect of dimension-6 operators in the dimensionally reduced theory [58]. The conclusion of the analysis in Appendix B is that we expect our uncertainties in T c , T n , αðT c Þ, and β=H Ã to be approximately 8%, 8%, 55%, and 80%, respectively. The dominant uncertainty is due to the neglect of dimension-6 operators, which get relatively large corrections from singlet-Higgs interactions, in particular from a 1 (for reference, c 6 ≈40jc 6;SM j for this benchmark point). From the 4-d perturbative standpoint, there are significant uncertainties due to the residual renormalization scale dependence of our one-loop approach. For the values reflected in Table III the renormalization scale was set to μ R ¼ ðm 1 þ m 2 Þ=2. Varying μ R between m Z and 2m 2 , we find uncertainties of ∼5% in T c , T n , ∼40% in αðT c Þ, and ∼80% in β=H Ã . Taking this into account we conclude that the perturbative and nonperturbative approaches agree within uncertainties. Note that at smaller values of x, corrections from dimension-6 operators will be comparatively larger, hence we expect that one cannot trust our DR for the xSM at x significantly smaller than x ¼ 0.036, reflected in Table III.
From Table III and the fit formula of Ref. [93] one can deduce that, for this benchmark point in the xSM, the peak frequency of the gravitational wave signal today is of order 1 Hz. This is near the peak sensitivity of both DECIGO [13] and BBO [14] though the strength of the signal is too weak to be detected. We plot this point in Fig. 2(b), along with the predictions from 3-d and 4-d perturbation theory. The sensitivity curves shown assume v w ¼ 1. A perturbative calculation of the wall velocity for similar points in the For these inputs the parameter λ 221 =v is tuned such that at the critical point x c ¼ 0.036. This allows us to repurpose the only preexisting nonperturbative result for bubble nucleation. In Appendix B we estimate that the uncertainty in the quantities in the first row. The error in T c and T n is only of order ten percent, whereas the error in αðT c Þ and β=H Ã is of order one. xSM was performed in Ref. [69], which found v w ∼ 0.1 for relatively weak transitions such as those we consider here. Assuming this value only weakens the gravitational wave signal, and so we have shown projections for v w ¼ 1 to remain conservative in our conclusions. Our scan over the xSM parameter space finds that η y ðT c Þ lies in the range [3.85,5.45] almost independently of x. This region is shaded dark blue on the right in Fig. 2(b), indicating values of α and β=H Ã in the parameter space where the singlet model is described by the 3-d SM-like EFT in the infrared. This region again lies outside that detectable by planned gravitational wave experiments. For a stronger transition there are two options: either the BSM field must play an active, dynamical role or the dimension-6 operators in the 3-d EFT must be included. In the case of the xSM, this is in agreement with previous perturbative studies (see, e.g., Ref. [75]).

B. Implications for colliders
While gravitational wave experiments such as LISA will likely be able to probe strong electroweak transitions in this model (for which nonperturbative studies do not exist), the transitional regions between a cross-over and a first-order EWPT can still be interesting cosmologically and will likely require other probes to access. Fortunately, colliders can be sensitive to precisely these regions and complement gravitational wave experiments in exploring the singlet model electroweak phase diagram. There has been a significant amount of work on the collider phenomenology of singlet models over the years, and a large number of possible experimental signatures of the new scalar have been proposed [63,105,106]. Some of these are directly related to λ 221 that is one of the couplings correlated with the nature of the electroweak phase transition in both the perturbative and nonperturbative approaches. Two such observables are the h 2 h 2 production cross section at hadron colliders [68,78,101], and the shift in the Higgs couplings from their SM predicted values [63,67,107], most notably the h 1 ZZ coupling [68,108]. Complementary information about the scalar potential can be inferred from di-Higgs processes [70,105,[108][109][110] as well as direct production of the singlet and subsequent decays to SM gauge bosons or fermions [111] at hadron colliders.
We demonstrate the interplay between some of the aforementioned experimental signatures on the parameter space predicting a first-order electroweak phase transition in Fig. 6. The left panel shows collider projections for the Z 2 -symmetric case. This model is very difficult to probe at the LHC but provides a compelling target for future colliders [68,101]. An e þ e − collider operating as a Higgs factory is expected to achieve excellent precision in measuring the Zh cross section, σ Zh (where h is the SM-like Higgs). The parameter space lying above the blue dashed contour feature δσ Zh > 0.5%, computed at one-loop using the expressions found in Ref. [78]. This corresponds approximately to the expected sensitivity of future circular Higgs factories such as the FCC-ee, or CEPC to this observable. Regions above the purple dashed contour could be probed at the 95% confidence level by searches for pp → h 2 h 2 jj at a 100 TeV collider with 30 ab −1 integrated luminosity. These projections are taken from Ref. [72], corrected to account for our slightly different choice of renormalization scale. While this search will likely not be sensitive to the EWPT parameter space accessible by our nonperturbative analysis, it should be able to probe stronger two-step phase transitions, which will require new lattice studies to explore nonperturbatively. The region above the FIG. 6. As in Fig. 3 (Z 2 ) and the right panel of 4 (non-Z 2 ) but including limits and projections for several collider observables, discussed in the text. The red shaded regions on the right are currently excluded either by the electroweak precision measurements discussed in Ref. [88], or direct searches for the singlet-like scalar in final states involving two gauge bosons [112]. The high luminosity LHC, a future Higgs factory and 100 TeV collider could together probe a significant amount of parameter space consistent with a firstorder electroweak phase transition as predicted by our nonperturbative analysis. red dashed contour features a deviation in the Higgs selfcoupling, δλ 111 ≡ jλ 111 − λ 111;SM j=λ 111;SM greater than 5%, which corresponds roughly to the level of sensitivity expected to be achieved at a future 100 TeV collider with 30 ab −1 [21]. Such measurements will be able to probe the region for which existing nonperturbative studies predict a first-order phase transition, although it is likely that including the relevant higher-dimension operators will be required to improve the nonperturbative predictions. Similar results for the non-Z 2 parameter space are shown for m 2 ¼ 400 GeV, b 3 ¼ 0 on the RHS of Fig. 6, where we show current constraints on the parameter space: the red shaded regions are excluded either by electroweak precision measurements discussed in Ref. [88], or direct searches for the singlet-like scalar in final states involving two gauge bosons. In deriving the latter, we use the 13 TeV ATLAS search in Ref. [112] and the narrow width approximation, with the h 2 production cross-section taken from Refs. [113,114] for a SM-like Higgs of the same mass, rescaled by sin 2 θ. We also show projected sensitivities for the same searches (pp → h 2 → VV) at the 13 TeV HL-LHC with 3 ab −1 (orange dashed contours), obtained from a simple rescaling of the sensitivity in Ref. [112]. The parameter space outside these sets of contours would be probed at the 95% confidence level and include a sizable portion of the first-order transition region.
We also show projections for future colliders in Fig. 6. The blue dashed contours indicate regions for which δσ Zh > 0.5% (here again h is understood to refer to the SM-like Higgs, h 1 in our case). The parameter space outside of these contours is expected to be probed by future Higgs factories. Also shown in Fig. 6 is the approximate sensitivity of a 100 TeV pp collider to the first-order EWPT parameter space through (resonant) di-Higgs production (pp → h 2 → h 1 h 1 ). Reference [109] suggests that a 100 TeV collider with 30 ab −1 integrated luminosity will be able to probe σ × BRðpp → h 2 → h 1 h 1 Þ down to ∼100 fb at 95% C.L. for m 2 ≈ 400 GeV. The projections shown assume this reach, and make use of the leading order gluon fusion h 2 production cross section in the narrow width approximation, multiplied by a k-factor of 2.3 to account for higher-order corrections (this k-factor was chosen to yield the gluon fusion cross section for a 125 SM-like Higgs at 100 TeV reported in Ref. [115]). The 100 TeV collider reach is quite impressive, extending down to very small mixing angles in the parameter space shown.
From Fig. 6 it is clear that the LHC and future colliders will be able to access much of the parameter space predicting a genuine first-order EWPT in this model. The most difficult to access region is that featuring very small jsin θj. However, this region is tuned (away from the Z 2 limit) and could in principle be probed if a future e þ e − collider were able to measure the Zh cross section to 0.25% precision. Given the different parametric dependence of the cross sections for pp → h 2 → VV, h 1 h 1 and e þ e − → Zh 1 , it will be important to perform all of these different measurements, if possible, in order to gain clear insight into the nature of the electroweak phase transition. Similar conclusions hold for different slices of the model parameter space. From the results shown in Fig. 6 it is clear that colliders will play an important role in exploring the phase structure in this model and complement the reach of gravitational wave experiments in this regard. The regions of parameter space most challenging to access for colliders will require new simulations including the effect of the higher dimension operators to more accurately study nonperturbatively.

VI. CONCLUSIONS
In this work, we have demonstrated the full pipeline for obtaining theoretically sound, nonperturbative results for the gravitational wave power spectrum arising from a first-order phase transition in the early Universe. Our proof of principle calculation collects and utilizes several tools developed in the literature for the SM with a light Higgs, and repurposes preexisting lattice results in the dimensionally reduced effective theory for equilibrium [4] and nonequilibrium thermodynamic quantities [39]. These relatively few quantities are then used as input parameters for the analysis of numerical simulations of the coupled fluid-scalar field system [93] in order to determine the prediction for the gravitational wave power spectrum from a first-order electroweak phase transition in general extensions of the Standard Model matching on to the infrared 3-d SM-like EFT.
The physics of the phase transition is largely determined by the long-distance 3-d effective field theory. Our analysis suggests that any SM extension which reduces to the SMlike 3-d EFT will have a rather weak EWPT and therefore be unobservable at gravitational wave experiments in the foreseeable future. To accommodate a strong transition producing an observably large gravitational wave signal, either new light fields or sizable dimension-6 operators must turn up in the IR effective theory. Otherwise, collider experiments are likely to provide the best chance of probing first-order electroweak phase transitions in these scenarios, which can still be interesting cosmologically despite the small corresponding gravitational wave signal.
We have applied these results to the real singlet extension of the Standard Model, identifying the regions where the SM-like electroweak crossover becomes a first-order transition. In the regions we have been able to analyze by repurposing pre-existing nonperturbative results-where the singlet field plays a nondynamical role and induces negligible effects through higher dimension operators-the gravitational wave signal is unfortunately far too weak to be detected with LISA or other near-future gravitational wave experiments, despite the fact that the electroweak phase transition can be of first-order. On the other hand, a combination of measurements at the high luminosity LHC and future colliders can in fact probe much of the first-order transition regions in the singlet model. Our nonperturbative results provide new targets for exploring the phase structure of the theory at colliders that cannot be delineated in a perturbative approach.
Observable gravitational wave signals in the xSM require the singlet field to play a dynamical role in the transition or induce non-negligible effects through higher dimension operators. This in turn makes the demand for new nonpertubative simulations even more urgent. Nevertheless, using exiting results we are able to benchmark the predictions of perturbation theory for the phase transition parameters relevant for computing the gravitational wave power spectrum. We find that perturbative results reproduce the latent heat and duration of the phase transition to within Oð1Þ factors, and the nucleation temperature to within 10% of the nonperturbative prediction, agreeing to within the expected uncertainties in our approach. It will be worthwhile and interesting to perform the analogous comparison for points within the LISA sensitivity band once new nonperturbative simulations including higher dimension operators or the dynamical singlet field become available.
In practice, we only solve the correction to the top quark mass using M t ¼ 173.1 GeV and neglect the procedure for other fermions as their Yukawa couplings have little effect on the phase transition. However, all fermions are included in the loop corrections. The correction to the Higgs VEV is calculated by relating it to the electromagnetic fine structure constantα, in the MS scheme, and evaluating the one-loop correction to on-shell Thomson scattering. Due to the absence of new charged particles in the singlet extension, this correction is equal to its SM value of [116]α −1 ðM Z Þ ≈ 128, which we input into our parametrization.
The self-energies are renormalized by the counterterms that were previously used in Ref. [57]. Note, however, that in this reference counterterms were listed in Landau gauge, whereas the vacuum renormalization is more convenient to carry out in the Feynman-t'Hooft gauge. Dependence on the gauge parameter is easily added to the counterterms by first evaluating the field renormalization counterterm and requiring the bare parameters to be gauge independent. A one-loop evaluation and renormalization of the selfenergies in the singlet-extended SM is straightforward, but since the actual formulas are fairly long, we do not list them here. The self-energies are functions of the renormalized masses and couplings, so for a self-consistent calculation we have to solve a nonlinear system of polemass equations. The solutions are conveniently found iteratively, and inserting the loop-corrected masses andα into the tree-level relations (A4)-(A9) gives the potential parameters at the input scale Λ ¼ ðM 1 þ M 2 Þ=2-and similarly for y t , g and g 0 .
As a result of this vacuum renormalization procedure, the renormalized parameters that enter the phase transition calculations are shifted by roughly 5%-30% relative to the values one would obtain by matching to physical observables directly at tree level. In the parameter scans, the largest corrections occur for the portal couplings a 1 , a 2 as well as for λ, which all play an important role in strengthening the phase transition. However, for masses M 2 ≳ 400 GeV the loop correction to the singlet mass parameter μ 2 σ can be large enough to cause deviations of order 100% from the tree-level value if additionally λ 221 =v ≫ 1 or a 1 =v ≫ 1, which may hint at bad convergence of perturbation theory even without considering infrared effects at high temperature. Overall, the effect of these zero-temperature corrections is to shift the phase transition toward smaller couplings and therefore have a non-negligible numerical impact on our EWPT analysis.

APPENDIX B: ACCURACY OF THE DIMENSIONAL REDUCTION
In Ref. [57] dimensional reduction was performed at one-loop level, leaving out the two-loop contribution to the mass parameter of the doublet field. This means that our determination of the critical temperature in the xSM is not as accurate as in the similar studies in the two Higgs doublet model (2HDM) [55,56] or the real-triplet extension of the SM [58]. However, we expect that the character of transition can still be determined with satisfactory accuracy and regions of first-order transition can be identified. In fact, the one-loop accuracy of the dimensional reduction matches the accuracy used in the perturbative four-dimensional finite-temperature effective potential [78] in this work.
By varying the renormalization scale of the 4-dimensional theory in the DR matching relations, we can estimate the systematic uncertainty in our analysis arising from scale dependence. By varying scale from T to 4πe −γ T ≈ 7.05Tcorresponding to the average momentum of integration of the superheavy Matsubara modes [43]-we find that we have about 10 to 20 percent uncertainty in determination of T c , x c , and the location of FOPT regions.
Furthermore, as in the original study of the SM [40], also in Ref. [57] for the xSM, the dimension-6 operators were dropped from the effective theory. While it is difficult to estimate the effect of the dimension-6 operators comprehensively, we can include leading order contributions to the operator c 6 ðϕ † 3 ϕ 3 Þ 3 and analyze the effective potential with this operator in 3-d perturbation theory. A similar approach was recently used in Ref. [56] (see Sec. III D 2 and Appendix C.7 therein for details).
Using the tools of Ref. [57], we can estimate the matching coefficient for the ðϕ † 3 ϕ 3 Þ 3 operator. In particular, utilizing the effective potential, Eqs. (3.27) and (3.29) in Ref. [57], we obtain V 6;0 ðφ † φÞ 3 with where J 3 ðμ σ Þ ≡ 1=ð2π 2 T 2 ÞJ 000 B ðμ 2 σ =T 2 Þ, where J B , called J þ in our analysis, is given, e.g., by Eq. (B1) of Ref. [58]. The other relevant master integrals can be found in Appendix B of Ref. [57] and in Ref. [58]. The first three terms above are the SM contributions and can be compared to Eqs. (196) and (197) in [40]. We emphasize that whereas the sumintegral I 4b 3 contains only contributions from nonzero Matsubara modes, the sumintegralJ 4b 3 ðμ σ Þ contains the zero mode as well, and the mass is now assumed to be superheavy. In the above expression, we have imposed a Z 2 -symmetry for simplicity, i.e., we set b 1 ¼ b 3 ¼ a 1 ¼ 0. Without imposing a Z 2 -symmetry, singlet contributions to the doublet 6-point function at one-loop become more complicated, as in addition to contributions from the effective potential, one needs to include all one-σ-reducible diagrams. However, in the non-Z 2 -symmetric case, there are already contributions at tree-level and we use these contributions to estimate the leading order effect. Therefore, our estimate of leading behavior reads c 6 ≃ 8ζð3Þ 3ð4πÞ 4 3 64 With this result, and using the Landau gauge effective potential of Ref. [56], we observe that for the point in Table III the critical temperature changes much less than one percent when including this operator, while the scalar VEV changes about −20%. In 3-d perturbation theory, values for α and β=H Ã change about −55% and 80%, respectively. Overall, these estimates show that for such a relatively strong first-order phase transition, integrating out the singlet zero mode and truncating the theory at the minimal 3-d EFT leads to large uncertainties. This uncertainty ultimately propagates to the nonperturbative analysis of the 3-d theory.
For weaker first-order transitions, these effects are expected to be less significant. From the arguments of Ref. [40], the shift in the VEV scales approximately linearly with c 6 =λ h;3 . The SM dimension-6 operators cause a shift of around 1% for values of x corresponding to first-order transitions near the crossover region [40], and so we expect the change in the Higgs VEV to be Oðc 6 =c 6;SM × 1%Þ for relatively weak transitions. On Figs. 3-5 we show where c 6 ¼ 40c 6;SM , assuming T ∼ 140 GeV. In the first-order transition regions near the c 6 ¼ 40c 6;SM contours, the dimension 6 effects should be large, as we found at our benchmark point. However, for c 6 =c 6;SM closer to one, our nonperturbative analysis should yield results with accuracy comparable to that achieved in the original SM dimensional reduction and nonperturbative studies.