Nonperturbative analysis of the gravitational waves from a first-order electroweak phase transition
Oliver Gould Jonathan Kozaczuk Lauri Niemi Michael J. Ramsey-Musolf Tuomas V. I. Tenkanen David J. Weir















  




oliver.gould@helsinki.fi



kozaczuk@illinois.edu



lauri.b.niemi@helsinki.fi



mjrm@physics.umass.edu



tenkanen@itp.unibe.ch




david.weir@helsinki.fi






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.



HIP-2019-5/TH
ACFI T19-04
Nonperturbative analysis of the gravitational waves from a

first-order electroweak phase transition
    David J. Weir
    December 30, 2023
====================================================================================================





§ 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 <cit.>, significantly affected the abundance of cosmic relics <cit.>, or produced a stochastic
background of gravitational waves (GWs) <cit.>. 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
<cit.>.
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 <cit.>.

In light of the successes of the LIGO and VIRGO
collaborations <cit.>, as
well as the LISA Pathfinder spacecraft <cit.>, a number
of space-based gravitational wave missions are now planned or in
preparation. Most notable of these is LISA <cit.>, due to launch before 2034;
others, such as DECIGO <cit.> and BBO <cit.>, 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 first-order 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) <cit.>, these include proposals for e^+e^- colliders such as
the ILC <cit.>, FCC-ee <cit.>, CEPC <cit.> and CLIC <cit.>, as well
as higher-energy proton-proton colliders such as the High-Energy LHC (HE-LHC) <cit.>,
FCC-hh <cit.> and SPPC <cit.>. 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 <cit.>.  Furthermore, perturbative studies are often plagued by issues related to gauge dependence <cit.>.
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. <cit.> 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 <cit.>. One has to tread very carefully in order to avoid double counting fluctuations <cit.>, to self-consistently account for the spacetime dependence of the bubble <cit.> and to maintain gauge invariance <cit.>.
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 <cit.>. 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 <cit.>. In the paper at hand, we perform such a benchmark comparison.

In practice, the most convenient approach to bypassing the infrared
problems inherent in perturbative calculations  is provided by dimensional reduction
(DR) <cit.>,
a technique for constructing an effective three-dimensional (3-d)
theory for the most infrared-sensitive light field modes. The
effective theory can then be studied on the lattice in an
infrared-safe
manner <cit.>. In the past, dimensional reduction has been applied to the
SM <cit.>, the
MSSM <cit.>
the Two-Higgs-Doublet model
(2HDM) <cit.>, 
as well as the real scalar singlet <cit.> (xSM) and triplet <cit.> (ΣSM) extensions of the SM (see also Ref. <cit.>). 
Perturbative analyses of the EWPT in
these models – usually based on the daisy-resummed one-loop effective
potential – have been widely studied in e.g. Refs. <cit.>
for the xSM, and for a selection of other models in e.g. Refs. <cit.>.

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
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. <cit.>; this strategy was first
adopted in Ref. <cit.>. Although computationally inexpensive, this approach
comes with its own limitations. One cannot study
interesting multistep transitions <cit.> 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. <cit.>). 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 model-independent 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 slow 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) <cit.>,
for which the required dimensional reduction has already been done in
Ref. <cit.>, but a systematic analysis of the
5-dimensional 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. <cit.> to deduce the
phase diagram of this model, comparing with the results of a traditional perturbative treatment as in Ref. <cit.>. 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 Section <ref> we
review the setup for recycling pre-existing nonperturbative results
for determining the phase structure of a given BSM scenario. In
Section <ref> 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 Section <ref> 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. Sec. <ref> discusses collider and gravitational wave probes
of the phase diagram in this theory, applying the results of Sections <ref>-<ref>. Finally, we conclude in Section
<ref>.



§ 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 first
order phase transition, in a theoretically sound manner.  The basic idea
of dimensional reduction is that at high temperatures, long-distance
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 non-zero 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. <cit.> and
<cit.>.

To illustrate the process of dimensional reduction in an extended scalar sector of the SM, we start with a four-dimensional (4-d) Euclidean Lagrangian in the imaginary-time formalism

    ℒ = ℒ (ϕ, A_μ, ψ, S, s),

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 ∼ g T). 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
non-zero 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

    ℒ_3 = ℒ_3 (ϕ_3, A_i, A_0, s_3),

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 non-zero 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 four-dimensional
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 g
T, 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
theory

    ℒ̅_3 = ℒ̅_3 (ϕ̅_3, A̅_i),

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 <ref>.  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 <cit.> (see also Section 6 of Ref. <cit.> 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. (<ref>) and a 3-d theory
containing the spatial SU(2) gauge field and a doublet
Higgs field,

    ℒ̅_3,scalar = D_i ϕ̅_3^† D_i ϕ̅_3 +
      μ̅_h,3^2 ϕ̅_3^†ϕ̅_3 + λ̅_h,3 (ϕ̅_3^†ϕ̅_3)^2,

with gauge coupling g̅_3. Note that nonperturbative effects of
the U(1)_Y gauge group were omitted in the original study
<cit.>, and we follow this same approach. Eq. (<ref>) is the same 3-d theory to which the
Standard Model was mapped in Ref. <cit.>, and then
studied nonperturbatively in
Refs. <cit.>. Furthermore, the
nucleation rate was studied nonperturbatively in
Ref. <cit.>. 

Remarkably, the theory described by Eqs. (<ref>)-(<ref>), 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 degrees of freedom 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. <cit.>), 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̅_3^2  and  
    y ≡μ̅_h,3^2/g̅_3^4,

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 <cit.>.
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. <cit.> showed that for 0 < x < 0.11 the
transition is of first-order, and becomes a crossover for x > 0.11
<cit.>.  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:


  * Dimensionally reduce the full 4-d theory to the infrared 3-d theory.

  * 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.

  * Identify regions for which 0<x<0.11 and y≈ 0. 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 non-dynamical.
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. (<ref>) and (<ref>) 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.



§ 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
<cit.>. 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_*
<cit.>. 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.




 §.§ Dimensional reduction for bubble nucleation
 

Given a theory that has been dimensionally reduced as in Sec. <ref> 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. (<ref>). Equilibrium properties of the 3-d
dimensionally reduced theory have been studied nonperturbatively in
Ref. <cit.> and in Ref. <cit.> 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 Table <ref>.





Following Ref. <cit.> 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

    α(T) ≡30/g(T) π^2L(T)/T^4,

where L(T) is the latent heat and g(T) the number of relativistic
degrees of freedom. When BSM scalar eigenstates are sufficiently
heavy, they are not relativistic degrees of freedom at the EWPT and we
have g(T ∼ T_*) = 106.75. 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. <cit.> and it reads

    L(T_c)/T^4_c = g^6_3/T_c^3(η_y(T_c) - η_x(T_c) dy_c/dx_c) Δ l_3,

where we have introduced the functions

    x(T)≡dx/dlog T   and  y(T)≡dy/dlog T .

For a given 4-d theory, these η-functions determine how the couplings
of the 3-d effective theory run with 
temperature. 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 
𝒪(g^-2),
the thermal mass corrections, whereas for x this leads to thermal corrections
first at 𝒪(g^2).
Thus the temperature dependence of x
is significantly weaker than that of y,

    x(T) ≪y(T),

whereas dy_c/dx_c=𝒪(1), hence

    L(T_c)/T^4_c≈g^6_3/T_c^3y(T_c)Δ l_3.


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 <cit.>,

    Γ = T/σ_el(g^2/4π)^5 T^4 Γ_3-d(x,y),

where σ_el∼ T/log(1/g) is the non-Abelian
“color" charge conductivity, arising in the effective description of
the real-time dynamics of the long wavelength degrees of freedom
<cit.>. 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. (<ref>)).

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 <cit.>

    Γ = 1/8π v_w^3(H dlogΓ/dlog T)^4,

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 <cit.>, 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

    β/H_*    = -dlogΓ/dlog T|_T=T_*
       ≈ -(x(T_*)∂/∂ x +
    y(T_*)∂/∂
      y)logΓ_3-d(x,y) 
       ≈ -y(T_*)∂/∂
      ylogΓ_3-d(x,y),

where all quantities are evaluated at the transition temperature[Note that in Ref. <cit.> the factor of η_y was mistakenly dropped, leading to their result for β/H_* being a factor of ∼ 4 too small.]. On the last line we have used Eq. (<ref>). 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. <cit.>, the quantity Γ_3-d(x,y) was calculated for x=0.036 and is reproduced here in Table <ref>. It is given as a function of (y-y_c), corresponding to δ m^2/g^4T^2 in the notation of Ref. <cit.>. Fitting a straight line to their results, one finds

    -logΓ_3-d(0.036,y)≈ 240 + 1.72× 10^4 (y-y_c),

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,

    β/H_*   ≈ 1.72× 10^4 η_y(T_*),
    α(T_c)    ≈ 2.2× 10^-3 η_y(T_c).

Details of the specific 4-d theory only enter these expressions through η_y. Interestingly, by combining Eqs. (<ref>) and (<ref>), one finds the relation

    β/H_*≈ 7.8× 10^6α(T_c),

 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 short-range
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 does not contain a tree-level
barrier. Hence one must integrate out some light degrees of freedom of
the 3-d theory to give the energy barrier between the vacua. The
degrees of freedom 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.



 §.§ 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. (<ref>)-(<ref>). 

The general recipe for this is
described as follows, assuming Eq. (<ref>) 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. (<ref>) one can find the phase transition temperature,
T_*, and finally evaluate Eqs. (<ref>) and
(<ref>) to give β/H_* and α. One can then
use the fit formula of Ref. <cit.> 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. <ref>. The gray shaded region in Fig. <ref> and the light blue region Fig. <ref>  can be reliably described by the infrared 3-d SM-like EFT. In Fig. <ref> we specifically state why this EFT description fails in different regions, as detailed below. The black line Fig. <ref>
corresponds to α and β/H_* values accessible by existing nonperturbative calculations. For reference, we also show the expected sensitivity of LISA <cit.> to this parameter space.
In Fig. <ref>, the region predicting a signal-to-noise ratio (SNR) of 10 at LISA is indicated in the bottom right corner. Sensitivity curves assuming a five year mission duration 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. <cit.> for a more detailed discussion of these distinct cases). In Fig. <ref>, 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. <ref> are determined as follows. First let us consider the range of allowed  η_y. 
This quantity, together with x_c, determines α(T_c) through Eq. (<ref>) and β/H_* through Eq. (<ref>) (for η(T_*)≈η(T_c)).
The value of η_y(T_c) is generically of order
𝒪(g^-2), absent new physics with large couplings. Its SM value for a 125GeV Higgs, η_y^SM(T_c)≈ 4.4,
receives corrections from BSM physics,

    η_y(T_c) = η_y^SM(T_c) + Δη_y^BSM(T_c),

where Δη_y^BSM(T_c) contains corrections which 
go to zero in the SM limit. In principle Δη_y^BSM(T_c) 
can be positive or negative, though all leading-order terms in 
η_y^SM(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 Δη_y^BSM(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)≪η_y^SM(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 fine-tuned 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. <ref>,
lead to unobservably, weak transitions from 
the perspective of their gravitational wave signal. 

Conversely, 
if Δη_y^BSM(T_c) is positive, one finds stronger 
transitions with shorter phase transition durations. The size of 
Δη_y^BSM(T_c) is nevertheless constrained by 
perturbativity. The perturbative loop expansion breaks down when couplings 
become large, limiting Δη_y^BSM(T_c) to be formally 𝒪(1/g^4) or smaller.
This implies η_y(T_c) ≲ 10 as a reasonable upper bound from perturbativity.
One possible caveat would be to instead add a large number
of weakly coupled degrees 
of freedom to increase Δη_y^BSM(T_c).
Though it is 
more subtle in this case, a sufficiently large number of degrees 
of freedom will also lead to a breakdown of perturbativity.
In either case, if Δη_y^BSM(T_c) is large,
the phase transition
duration becomes very short, shown as the top right arrow in 
Fig. <ref>. 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. (<ref>) and (<ref>) trace out the thick
black line shown in Fig. <ref>.

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. <cit.>. In the calculation of
Γ_3-d, we follow the perturbative semiclassical analysis
outlined in Ref. <cit.>, 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 <cit.>.

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. <ref>.
 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. <ref>). 
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 𝒪(y_t^6/(4π)^3), where y_t is the top Yukawa coupling
(see Appendix <ref>). 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,

    λ̅_h,3⟨ (ϕ̅_3^†ϕ̅_3)^2⟩_c ≫ 
     c_6⟨ (ϕ̅_3^†ϕ̅_3)^3⟩_c,

where the expectation values refer to the volume averaged 
condensates, evaluated at the critical point (see, for example, Ref. <cit.>).
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. (<ref>)-
(<ref>), 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. <ref> and  <ref>, repsectively.
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
SM-like EFT and gives a reasonably strong first-order phase transition
will lie in the light blue region of Fig. <ref>. The dark blue region in Fig. <ref> is for the specific case
of the singlet extended Standard Model (xSM), which we discuss further in Section <ref>.

From Fig. <ref> 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. <ref> and <ref>, 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. <ref>). 




§ 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 often-studied
real singlet extension of the Standard Model (xSM)
 <cit.>.  
Using the dimensional reduction of
Ref. <cit.>, in this section we first investigate the 5-dimensional
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. <ref>, comparing
with perturbative results, and discuss the prospects for exploring the phase diagram of this scenario
at collider experiments.



 §.§ Parameterization of the model



The most general scalar Lagrangian involving an additional real singlet scalar field can be written as <cit.>

    ℒ_scalar = (D_μϕ)^† (D_μϕ) - μ^2
      ϕ^†ϕ + λ (ϕ^†ϕ)^2 
    
      + 1/2
      (∂_μσ)^2 + 1/2μ_σ^2 σ^2 + b_1
      σ + 1/3 b_3 σ^3 + 1/4b_4
      σ^4 
    
      + 1/2 a_1 σϕ^†ϕ + 1/2
      a_2 σ^2 ϕ^†ϕ,

where D_μ is the usual covariant derivative for the Higgs doublet,
and parameters are renormalized in the scheme. In order to
relate these 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 ⟨σ⟩ =0, and that the Higgs
vev is ⟨ϕ^†ϕ⟩ = v^2/2 . This leads to the
following relations in the electroweak vacuum:

    μ^2 = v^2 λ,     b_1 = -a_1 v^2/4.

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 <cit.>.
 We treat the singlet self-couplings
b_3 and b_4, given at a fixed 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. <cit.>. 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 <ref>.





 §.§ 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. <cit.> and the results of Sections <ref>-<ref> can be applied. The singlet field in this case is non-dynamical 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 <cit.>. This can be seen explicitly in the expression in Ref. <cit.> for the 3-d Higgs quartic coupling after integrating out the superheavy degrees of freedom, λ_h,3, in terms of the model parameters:

    λ_h,3/T=λ - a_1^2/8μ^2_σ - a_1^2 b_3 b_1/4μ^6_σ + a_1 a_2 b_1/2μ^4_σ +  loop   effects.

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 𝒪(g^4) in the power counting scheme used in the dimensional reduction of Ref. <cit.>.

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 self-interactions. Our nonperturbative results are shown in Fig. <ref>. 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 singlet-induced dimension-6 operator coefficient, c_6, becomes larger than 10 times the corresponding SM value, c_6,  SM, assuming T≈ 140 GeV. An approximate expression for c_6 can be found in Eq. (<ref>). As discussed in Appendix <ref>, for c_6 ≳  10 c_6, SM the neglected higher dimension operators are expected to have 𝒪(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. <ref>, 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 higher-loop zero-T effects become non-negligible.  Below the μ_σ=200 GeV and c_6 ≳  10 c_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. <cit.>. Note also that similar results for the phase diagram in the real triplet-extended SM were obtained in Ref. <cit.>.



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. <ref> 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=10 c_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. <ref>. In these slices, portions of the gray shaded regions are excluded by 1-loop absolute vacuum stability, as discussed in Ref. <cit.> (see also Ref. <cit.>).





As can be seen from Figs. <ref>-<ref>, the neglected singlet-induced dimension-6 operator can become important in the first-order transition regions, especially for small |sinθ| and in the Z_2 limit. This is expected from simple power counting arguments. Comparing Eq. (<ref>) and Eq. (<ref>) for the approximate dimension-6 operator coefficient, we see that generally, for superheavy singlets such that μ_σ∼π T, the effect of the singlet on the dimension-6 operator is suppressed by ∼ a_2/π^2 relative to the leading effect on λ_3 (and hence the strength of the phase transition). However, for |sinθ| ≈ 0, the effect on λ_h,3 starts at 𝒪(a_2^2/16 π^2), rather than at  𝒪(g^2/8 π^2) for a_1∼ g T, and larger values of a_2 are required to achieve a given value of x. Thus, for small |sinθ|, the a_2/π^2 suppression of the dimension-6 operator is reduced, rendering these contributions more important. Our results motivate new nonperturbative studies including these dimension-6 operators, which will be required to more accurately study first-order transitions in the small |sinθ| and Z_2 regions. As discussed in Sec. <ref> below, these are precisely the regions requiring future colliders to experimentally probe.



 §.§ 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

    V_ eff=    V_0 + ∑_i ± n_i/64 π^2 m_i^4(h,σ)[log(m_i^2(h,σ)/μ_R^2) - c_i ] 
       + T^4/2π^2∑_i ± n_i J_±(m_i^2(h,σ)/T^2).

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 degrees of freedom for the species i, c_i are scheme-dependent constants, μ_R is the renormalization scale, and the J_± encode the thermal corrections (see e.g. Ref. <cit.> 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 self-energies in this step, although there has been recent progress in going beyond this approximation <cit.>. In what follows, we use the high-T thermal masses, the expressions for which can be found in e.g. Refs. <cit.>, as this is the most conventional approach. The effective potential is computed in Landau gauge utilizing the scheme.

With V_ eff defined above, we search for first-order electroweak phase transitions in the same way as that detailed in Ref. <cit.>, tracing the minima of the finite-temperature 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. Ref. <cit.> 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 gauge-invariant) 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 first-order 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 corresponding to ζ≈ 0.6-2 (see e.g. Refs. <cit.> 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 parameterized 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. <ref>-<ref> 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].

 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. <ref> 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. <ref>, 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 orange shaded regions of Figs. <ref>-<ref> show the parameter space for which v_c/T_c ∈ [0.3, 0.6] in the non-Z_2 case, where the transition is strengthened mainly by the tree-level contributions from the singlet, cf. Eq. (<ref>). These regions line up quite well with the first-order EWPT parameter space predicted by our nonperturbative analysis, especially where the (ϕ^†ϕ)^3 effects are expected to be small. 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.




§ 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. <ref>. 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.
 


 §.§ 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 Section <ref>. 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 <ref>.
Also shown are the same quantities calculated using the usual 4-d perturbative approach (4-d PT) – see Section <ref> – as well as our two-loop, fully-perturbative approach utilising the DR (3-d PT) – see Section <ref>.

A comparison of the nonperturbative and perturbative results in Table <ref> 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 <cit.>. 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 <ref> 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 <cit.>. The conclusion of the analysis in Appendix <ref> is that we expect our uncertainties in T_c, α(T_c) and β/H_* to be of order 1%, 40% and 60% 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 ≈ 16 c_6,  SM 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 <ref> 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, ∼ 50% in α(T_c), and ∼ 80% in β/H_*. Taking this into account we conclude that, except for T_c, 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 <ref>.

From Table <ref> and the fit formula of Ref. <cit.> 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 <cit.> and BBO <cit.> though the strength of the signal is too weak to be detected. We plot this point in Fig. <ref>, 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 xSM was performed in Ref. <cit.>, 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 in Fig. <ref>, 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. <cit.>).



 §.§ 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 <cit.>. 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 <cit.>, and the shift in the Higgs couplings from their SM predicted values <cit.>, most notably the h_1ZZ coupling <cit.>. Complementary information about the scalar potential can be inferred from di-Higgs processes <cit.> as well as direct production of the singlet and subsequent decays to SM gauge bosons or fermions <cit.> 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. <ref>. 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 <cit.>. 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. <cit.>. 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 p p → h_2 h_2 j j at a 100 TeV collider with 30 ab^-1 integrated luminosity. These projections are taken from Ref. <cit.>, 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 red dashed contour features a deviation in the Higgs self-coupling, δλ_111≡ |λ_111 - λ_111,  SM|/ λ_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 <cit.>. 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. <ref>, where we show current constraints on the parameter space: the red shaded regions are excluded either by electroweak precision measurements discussed in Ref. <cit.>, 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. <cit.> and the narrow width approximation, with the h_2 production cross-section taken from Refs. <cit.> 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. <cit.>. The parameter space outside the 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. <ref>. 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. <ref> 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). Ref. <cit.> 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. <cit.>). The 100 TeV collider reach is quite impressive, extending down to very small mixing angles in the parameter space shown.

From Fig. <ref> 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 |sinθ|. 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. <ref> 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.




§ 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 pre-existing lattice results in the dimensionally reduced effective theory for equilibrium <cit.> and non-equilibrium thermodynamic quantities <cit.>. These relatively few quantities are then used as input parameters for the analysis of numerical simulations of the coupled
fluid-scalar field system <cit.> 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 SM-like 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 non-dynamical 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 𝒪(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.




§ ACKNOWLEDGMENTS


The authors would like to thank Daniel Cutting, Patrick Draper, Mark
Hindmarsh, Stephan Huber, Venus Keus, Mikko Laine, Jonathan Manuel,
Jose M. No, Hiren Patel, Kari Rummukainen, Anders Tranberg and Aleksi
Vuorinen for discussions.  The work of JK was supported by NSF grant
PHY-1719642. MJRM was supported in part under U.S. Department of
Energy contract DE-SC0011095. LN was supported by Academy of Finland
grant no. 308791 and the Jenny and Antti Wihuri Foundation.
TT has been supported by the Vilho, Yrjö and
Kalle Väisälä Foundation, and by the European Research
Council grant no. 725369.  In addition, TT was supported by the Swiss
National Science Foundation (SNF) under grant 200020-168988.  DJW
(ORCID ID 0000-0001-6986-0517) was supported by Academy of Finland
grant no. 286769. OG and DJW are supported by the Research Funds of
the University of Helsinki. The work of JK, MJRM and DJW was performed
in part at the Aspen Center for Physics, which is supported by
National Science Foundation grant PHY-1607611.





§ DETAILS OF RELATING RENORMALIZED PARAMETERS TO PHYSICAL OBSERVABLES



As described briefly in Section <ref>, the input parameters
for our scans of the singlet-extended SM are the pole masses M_1, M_2 (which we have been writing as m_1, m_2, but in this section are denoted by capital letters to distinguish them from the corresponding masses) of the mass eigenstates
h_1 and h_2, as well as the mixing angle θ and the couplings
b_3,b_4,λ_221, which are assumed to be inputed directly
at a fixed scale. We parameterize the scalar doublet as

    ϕ = [                 ω^+; 1/√(2)(v + h + i z) ],

with ω^±, z being the Goldstone modes.

After rotating the scalar potential (<ref>) to a diagonal basis via

    h_1    = hcosθ - σsinθ, 
    
    h_2    = hsinθ + σcosθ,

the mass-eigenvalue relations can be inverted to express the parameters in the form

    μ^2    = 1/4( m^2_1 + m^2_2 + (m^2_1 - m^2_2)cos 2θ), 
    μ^2_σ   = 1/2( m^2_1 + m^2_2  + (m^2_2 - m^2_1)cos 2θ - v^2 a_2 ), 
    
    b_1    = -1/4 (m^2_2 - m^2_1) v sin 2θ, 
    
    a_1    = m^2_2 - m^2_1/vsin 2θ, 
    λ   = 1/4v^2( m^2_1 + m^2_2 + (m^2_1 - m^2_2 )cos 2θ).

Masses in the diagonal basis have been denoted by m_1 and m_2 to distinguish them from the pole masses M_1, M_2. The portal coupling a_2 can be solved in favor of the trilinear (h_2,h_2,h_1) coupling of the diagonalized theory via

    a_2    = (m^2_1 + 2m^2_2)(cos 2θ - 1) + 2 v λ_221θ + 2 v b_3 sin 2θ/v^2 (3cos 2θ - 1).

Inverted relations for the top Yukawa y_t and gauge couplings g,g' remain unchanged from the Standard Model case.

We follow Refs. <cit.> in relating the parameters appearing in the Lagrangian to physical pole masses. In practice this procedure consists of solving for the running masses from loop-corrected pole conditions of the form

    [p^2 - m^2 (Λ) + Π(p^2, Λ) ]_p^2 = M^2 = 0,

where Π denotes (the real part of) a self-energy calculated at one-loop level. In order to reduce the effect of logarithmic corrections, we calculate the self energies with the RG scale Λ set equal to (M_1+M_2)/2. The pole-mass conditions fix the renormalized masses of h_1,h_2 and gauge bosons W,Z at the input scale. The numerical values used for the pole masses are M_1 = 125.09 GeV,corresponding to the observed Higgs boson, M_W = 80.385 GeV and M_Z = 91.188 GeV. For fermions, the pole condition depends on the scalar and vector parts of the self-energy, and the relevant equation reads

    m_f^2 = M_f^2 (1 - 2 Σ_s(M_f^2,Λ) - 2 Σ_v(M_f^2,Λ)).

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 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 <cit.> α̂^-1(M_Z) ≈ 128, which we input into our parameterization.

The self-energies are renormalized by the counterterms that were previously used in Ref. <cit.>. 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 self-energies in the singlet-extended SM is straightforward, but since the actual formulae 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 non-linear system of pole-mass equations. The solutions are conveniently found iteratively, and inserting the loop-corrected masses and α̂ into the tree-level relations (<ref>)-(<ref>) gives the potential parameters at the input scale Λ=(M_1 + M_2)/2 – and similarly for y_t,g and g'.

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 towards smaller couplings and therefore have a non-negligible numerical impact on our EWPT analysis.




§ ACCURACY OF THE DIMENSIONAL REDUCTION



In Ref. <cit.> 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)
<cit.> or the real-triplet extension of the SM
<cit.>. 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 <cit.> in this work.


By varying the renormalisation 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.05 T – corresponding
to the average momentum of integration of the superheavy Matsubara modes
<cit.> – 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 <cit.>, also in
Ref. <cit.> 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. <cit.> (see Section 3.4.2 and
 Appendix C.7 therein for details).

Using the tools of Ref. <cit.>, 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. <cit.>, we
 obtain V_6,0 (φ^†φ)^3 with

    V_6,0≃(d/16 g^6 + 40 λ^3 )  I^4b_3 - 2 y^6_t I^4f_3 + 1/6 a^3_2 J̃^4b_3(μ_σ) 
    
    = 8 ζ(3)/3 (4π)^4 T^2( 3/64 g^6 - 21/2y^6_t + 30 λ^3 ) + a^3_2/6(1/32π^2 μ^2_σ + J_3(μ_σ) ),

where J_3(μ_σ)≡ 1/(2π^2 T^2) J_B”'(μ^2_σ/T^2), where J_B, called J_+ in our analysis, is given e.g. by Eq. (B1) of Ref. <cit.>. The other relevant master integrals can be found in Appendix B of Ref. <cit.> and in Ref. <cit.>.
The first three terms above are the SM contributions and can be compared to Eqs. (196) and (197) in <cit.>. We emphasize that whereas the sumintegral I^4b_3 contains only contributions from non-zero Matsubara modes, the sumintegral J̃^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

    1/24( a^2_1 a_2/μ^4_σ - a^3_1 b_3/μ^6_σ),

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 g^6 - 21/2y^6_t + 30 λ^3 ) + T^2 a^3_2/6(1/32π^2 μ^2_σ + J_3(μ_σ) ) + T^2/24( a^2_1 a_2/μ^4_σ - a^3_1 b_3/μ^6_σ).

With this result, and using the Landau gauge effective potential of
Ref. <cit.>, we observe that for the point in
Table <ref> 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 -40% and -60%,
 respectively. Overall, these estimates suggest that the uncertainty resulting from
 integrating out the singlet zero mode could be several tens of percent for points featuring relatively strong first-order transitions, which
 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. <cit.>, 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 <cit.>, and so we expect the change in the Higgs VEV to be 𝒪(c_6/c_6,  SM× 1%) for relatively weak transitions. On Figs. <ref>-<ref> we show where c_6=10 c_6,  SM, assuming T∼ 140 GeV. In the first-order transition regions for which this ratio is smaller than 10 or so, the corresponding shift in the VEV is expected be of order a few percent, and our nonperturbative analysis should yield results with accuracy comparable to that achieved in the original SM dimensional reduction and nonperturbative studies.  








