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

TL;DR

This paper provides a nonperturbative analysis of gravitational waves from a first-order electroweak phase transition, showing such signals are generally too weak for detection unless higher-dimension operators or light BSM fields are considered.

Contribution

It introduces the first end-to-end nonperturbative framework for analyzing gravitational waves from EWPTs using dimensionally reduced effective field theory and compares results with perturbative methods.

Findings

EWPT gravitational wave signals are typically too weak for detection.

Higher-dimension operators or light BSM fields are needed for observable signals.

Collider experiments are more promising for exploring phase structure.

Abstract

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.

Tables

Table 1. Table 1: Nonperturbatively obtained 3-d phase transition parameters from Refs.  Kajantie  et al. ( 1996a ); Rummukainen  et al. ( 1998 ); Moore and Rummukainen ( 2001 ) (see also Ref.  Gurtler  et al. ( 1997 ) ). The entry marked with a dagger is interpolated based solely on values for other x c subscript 𝑥 𝑐 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 ) U 1 \mathrm{U}(1) field, which makes the phase transition slightly stronger  Kajantie  et al. ( 1997 ) . This moves the critical line endpoint to x c ≈ 0.11 subscript 𝑥 𝑐 0.11 x_{c}\approx 0.11   Rummukainen  et al. ( 1998 ) .

$x_c$	$y_c$	$\Delta {\ell }_3$	Source
$0.01830$	$0.05904\left(56\right)$	$4.07\left(13\right)$	Kajantie et al. (1996a)
$0.036$	$†0.030$	$1.265$	Moore and Rummukainen (2001)
$0.06444$	$-0.00146\left(35\right)$	$0.491\left(8\right)$	Kajantie et al. (1996a)
$0.08970$	$-0.01531\left(69\right)$	$0.302\left(18\right)$	Kajantie et al. (1996a)
$0.0983\left(15\right)$	$-0.018$	$0$	Rummukainen et al. (1998)

Table 2. Table 2: The dimensionless rate of bubble nucleation at x = 0.036 𝑥 0.036 x=0.036 as calculated nonperturbatively in Ref.  Moore and Rummukainen ( 2001 ) . In that reference − log ⁡ Γ 3 − d subscript Γ 3 d -\log\Gamma_{\mathrm{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 𝑦 subscript 𝑦 𝑐 𝛿 superscript 𝑚 2 superscript 𝑔 4 superscript 𝑇 2 (y-y_{c})\equiv\delta m^{2}/g^{4}T^{2} .

$x_c$	$\left(y-y_c\right)$	$-\log {\Gamma }_{3-d}$	$\mathrm{Source}$
$0.036$	$-0.00835$	$97.0\pm 0.7$	Moore and Rummukainen (2001)
$0.036$	$-0.00951$	$74.1\pm 0.6$	Moore and Rummukainen (2001)
$0.036$	$-0.01057$	$58.8\pm 0.7$	Moore and Rummukainen (2001)

Table 3. Table 3: Sample of thermodynamic quantities relevant for determination of gravitational wave power spectrum. These values correspond to the xSM input parameters m 2 = 400 subscript 𝑚 2 400 m_{2}=400 GeV, b 3 = 0 subscript 𝑏 3 0 b_{3}=0 GeV, b 4 = 0.25 subscript 𝑏 4 0.25 b_{4}=0.25 , sin ⁡ ( θ ) = 0.2 𝜃 0.2 \sin(\theta)=0.2 , λ 221 / v = 2.22 subscript 𝜆 221 𝑣 2.22 \lambda_{221}/v=2.22 . For these inputs the parameter λ 221 / v subscript 𝜆 221 𝑣 \lambda_{221}/v is tuned such that at the critical point x c = 0.036 subscript 𝑥 𝑐 0.036 x_{c}=0.036 . This allows us to repurpose the only pre-existing nonperturbative result for bubble nucleation. In Appendix  B we estimate that the uncertainty in the quantities in the first row to be a few tens of percent.

	$T_c$/GeV	$T_n$/GeV	$\alpha \left(T_c\right)$	$\beta /H_{\ast }$
NP	$140.4$	$140.2$	$0.011$	$8.20\times {10}^4$
3-d PT	$140.4$	$140.0$	$0.010$	$6.11\times {10}^4$
4-d PT	$131.0$	$130.7$	$0.004$	$5.59\times {10}^4$

Equations

$$\mathcal{L}=\mathcal{L}(\varphi ,A_{\mu },\psi ,S,s),$$

$${\mathcal{L}}_3={\mathcal{L}}_3({\varphi }_3,A_i,A_0,s_3),$$

$${\bar{\mathcal{L}}}_3={\bar{\mathcal{L}}}_3({\bar{\varphi }}_3,{\bar{A}}_i),$$

$${\bar{\mathcal{L}}}_{3,\text{scalar}}=D_i{\bar{\varphi }}_3^{†}{D}_i{\bar{\varphi }}_3+{\bar{\mu }}_{h,3}^2{\bar{\varphi }}_3^{†}{\bar{\varphi }}_3+{\bar{\lambda }}_{h,3}({\bar{\varphi }}_3^{†}{\bar{\varphi }}_3{)}^2,$$

$$x\equiv \frac{{\bar{\lambda }}_{h,3}}{{\bar{g}}_3^2}\text{and}y\equiv \frac{{\bar{\mu }}_{h,3}^2}{{\bar{g}}_3^4},$$

$$\alpha (T)\equiv \frac{30}{g(T){\pi }^2}\frac{L(T)}{T^4},$$

$$\frac{L(T_c)}{T_c^4}=\frac{g_3^6}{T_c^3}\left({\eta }_y(T_c)-{\eta }_x(T_c)\frac{d{y}_c}{d{x}_c}\right)\Delta l_3,$$

$${\eta }_x(T)\equiv \frac{dx}{d\log T}\text{and}{\eta }_y(T)\equiv \frac{dy}{d\log T}.$$

$${\eta }_x(T)\ll {\eta }_y(T),$$

$$\frac{L(T_c)}{T_c^4}\approx \frac{g_3^6}{T_c^3}{\eta }_y(T_c)\Delta l_3.$$

$$\Gamma =\frac{T}{{\sigma }_{\mathrm{el}}}{\left(\frac{g^2}{4\pi }\right)}^5{T}^4{\Gamma }_{3-d}(x,y),$$

$$\Gamma =\frac{1}{8\pi v_w^3}{\left(H\frac{d\log \Gamma }{d\log T}\right)}^4,$$

$$\frac{\beta }{H_{\ast }}$$

$$\approx -\left({\eta }_x(T_{\ast })\frac{\partial }{\partial x}+{\eta }_y(T_{\ast })\frac{\partial }{\partial y}\right)\log {\Gamma }_{3-d}(x,y)$$

$$\approx -{\eta }_y(T_{\ast })\frac{\partial }{\partial y}\log {\Gamma }_{3-d}(x,y),$$

$$-\log {\Gamma }_{3-d}(0.036,y)\approx 240+1.72\times 10^4(y-y_c),$$

$$\frac{\beta }{H_{\ast }}$$

$$\alpha (T_c)$$

$$\frac{\beta }{H_{\ast }}\approx 7.8\times 10^6\alpha (T_c),$$

$${\eta }_y(T_c)={\eta }_y^{\mathrm{SM}}(T_c)+\Delta {\eta }_y^{\mathrm{BSM}}(T_c),$$

$${\bar{\lambda }}_{h,3}⟨({\bar{\varphi }}_3^{†}{\bar{\varphi }}_3{)}^2{⟩}_c\gg c_6⟨({\bar{\varphi }}_3^{†}{\bar{\varphi }}_3{)}^3{⟩}_c,$$

$$\begin{gathered} {\mathcal{L}}_{\text{scalar}}=(D_{\mu }\varphi {)}^{†}(D_{\mu }\varphi )-{\mu }^2{\varphi }^{†}\varphi +\lambda ({\varphi }^{†}\varphi {)}^2 \\ +\frac{1}{2}({\partial }_{\mu }\sigma {)}^2+\frac{1}{2}{\mu }_{\sigma }^2{\sigma }^2+b_1\sigma +\frac{1}{3}{b}_3{\sigma }^3+\frac{1}{4}{b}_4{\sigma }^4 \\ +\frac{1}{2}{a}_1\sigma {\varphi }^{†}\varphi +\frac{1}{2}{a}_2{\sigma }^2{\varphi }^{†}\varphi , \end{gathered}$$

$${\mu }^2=v^2\lambda ,b_1=-a_1\frac{v^2}{4}.$$

$$\frac{{\lambda }_{h,3}}{T}=\lambda -\frac{a_1^2}{8{\mu }_{\sigma }^2}-\frac{a_1^2{b}_3{b}_1}{4{\mu }_{\sigma }^6}+\frac{a_1{a}_2{b}_1}{2{\mu }_{\sigma }^4}+\mathrm{loop}\mathrm{effects}.$$

$$V_{\mathrm{eff}}=$$

$$+\frac{T^4}{2{\pi }^2}\sum _i\pm n_i{J}_{\pm }\left(\frac{m_i^2(h,\sigma )}{T^2}\right).$$

$$\zeta \in [0.3,0.6].$$

$$\varphi =\left(\begin{array}{c}{\omega }^{+}\\ \frac{1}{\sqrt{2}}\left(v+h+iz\right)\end{array}\right),$$

$$h_1$$

$$h_2$$

$${\mu }^2$$

$${\mu }_{\sigma }^2$$

$$b_1$$

$$a_1$$

Full text

††preprint: HIP-2019-5/TH††preprint: ACFI T19-04

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 Morrissey and Ramsey-Musolf (2012), significantly affected the abundance of cosmic relics Wainwright and Profumo (2009), or produced a stochastic background of gravitational waves (GWs) Caprini et al. (2016). 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 Kajantie et al. (1996a, b, 1997); Csikor et al. (1999). 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 Caprini et al. (2016); Weir (2018).

In light of the successes of the LIGO and VIRGO collaborations Abbott et al. (2016, 2017), as well as the LISA Pathfinder spacecraft Armano et al. (2016), a number of space-based gravitational wave missions are now planned or in preparation. Most notable of these is LISA Audley et al. (2017), due to launch before 2034; others, such as DECIGO Kawamura et al. (2011) and BBO Harry et al. (2006), 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) Cepeda et al. (2019), these include proposals for $e^{+}e^{-}$ colliders such as the ILC Baer et al. (2013), FCC-ee Abada et al. (2019a), CEPC Dong and Li (2018) and CLIC Linssen et al. (2012), as well as higher-energy proton-proton colliders such as the High-Energy LHC (HE-LHC) Abada et al. (2019b), FCC-hh Abada et al. (2019c) and SPPC Ahmad et al. (2015). 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 Linde (1980); Gross et al. (1981). Furthermore, perturbative studies are often plagued by issues related to gauge dependence Patel and Ramsey-Musolf (2011); Garny and Konstandin (2012). 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. Chung et al. (2013) 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 Langer (1974); Berges et al. (1997); Litim (1997). One has to tread very carefully in order to avoid double counting fluctuations Weinberg (1993); Gleiser et al. (1993); Alford and March-Russell (1994), to self-consistently account for the spacetime dependence of the bubble Surig (1998); Garbrecht and Millington (2015) and to maintain gauge invariance Metaxas and Weinberg (1996); Baacke and Heitmann (1999); Patel and Ramsey-Musolf (2011); Garny and Konstandin (2012). 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 Strumia and Tetradis (1999). 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 Moore and Rummukainen (2001). 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) Kajantie et al. (1996c); Appelquist and Carazzone (1975); Braaten and Nieto (1995); Farakos et al. (1994), 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 Kajantie et al. (1996a, 1997); Farakos et al. (1996, 1995). In the past, dimensional reduction has been applied to the SM Kajantie et al. (1996a, b, 1997), the MSSM Laine (1996); Cline and Kainulainen (1996); Farrar and Losada (1997); Cline and Kainulainen (1998); Laine (1997); Laine and Rummukainen (2001); Laine et al. (2014) the Two-Higgs-Doublet model (2HDM) Losada (1997); Andersen (1999); Andersen et al. (2018); Gorda et al. (2018), as well as the real scalar singlet Brauner et al. (2017) (xSM) and triplet Niemi et al. (2018) ($\Sigma$SM) extensions of the SM (see also Ref. Tenkanen (5 03)). 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. Huber and Schmidt (2001); O’Connell et al. (2007); Ahriche (2007); Profumo et al. (2007); Espinosa et al. (2012); Cline and Kainulainen (2013); Damgaard et al. (2013); Profumo et al. (2015); Curtin et al. (2014); Kozaczuk (2015); Huang et al. (2016a); Damgaard et al. (2016); Curtin et al. (2018); Beniwal et al. (2017); Kurup and Perelstein (2017); Alves et al. (2018a, b); Beniwal et al. (2018); Chen et al. (2017) for the xSM, and for a selection of other models in e.g. Refs. Apreda et al. (2002); Dorsch et al. (2014, 2017); Basler et al. (2017, 2018); Patel and Ramsey-Musolf (2013).

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. Kajantie et al. (1996a, 1997); this strategy was first adopted in Ref. Cline and Kainulainen (1996). Although computationally inexpensive, this approach comes with its own limitations. One cannot study interesting multistep transitions Patel and Ramsey-Musolf (2013); Blinov et al. (2015) 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. Tenkanen (5 03)). 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) Barger et al. (2008); Ashoorioon and Konstandin (2009); Robens and Stefaniak (2015); Kanemura et al. (2017, 2016); Beniwal et al. (2017), for which the required dimensional reduction has already been done in Ref. Brauner et al. (2017), 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. Kajantie et al. (1996a) to deduce the phase diagram of this model, comparing with the results of a traditional perturbative treatment as in Ref. Chen et al. (2017). 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 II we review the setup for recycling pre-existing nonperturbative results for determining the phase structure of a given BSM scenario. In Section 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 Section 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. Sec. V discusses collider and gravitational wave probes of the phase diagram in this theory, applying the results of Sections II-III. Finally, we conclude in Section 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 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 $\pi 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. Kajantie et al. (1996c); Laine (1997); Laine and Vuorinen (2016) and Gorda et al. (2018); Niemi et al. (2018).

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

[TABLE]

where the 4-d fields depend on the imaginary time $\tau$ and spatial coordinates $\vec{x}$. The SM fields $\phi$, $A_{\mu}$ and $\psi$ – 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}\sim\pi T$) and the scalar $s$ is heavy, meaning its mass is suppressed in a power counting sense ($m_{s}\sim gT$). We emphasize that here $m_{S}$ and $m_{s}$ are running masses in e.g. the $\overline{\text{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

[TABLE]

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 $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 $\phi$ 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

[TABLE]

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 Farakos et al. (1994, 1995) (see also Section 6 of Ref. Laine and Vuorinen (2016) 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 $\mathrm{SU(2)}$ gauge field and a doublet Higgs field,

[TABLE]

with gauge coupling $\bar{g}_{3}$. Note that nonperturbative effects of the $\mathrm{U}(1)_{Y}$ gauge group were omitted in the original study Kajantie et al. (1996a), and we follow this same approach. Eq. (4) is the same 3-d theory to which the Standard Model was mapped in Ref. Kajantie et al. (1996c), and then studied nonperturbatively in Refs. Kajantie et al. (1996a, b). Furthermore, the nucleation rate was studied nonperturbatively in Ref. Moore and Rummukainen (2001).

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 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 $(\bar{\phi}_{3}^{\dagger}\bar{\phi}_{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. Andersen et al. (2018)), 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,

[TABLE]

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 Kajantie et al. (1996a). Furthermore, the condition $y\approx 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. Kajantie et al. (1996a) showed that for $0<x<0.11$ the transition is of first-order, and becomes a crossover for $x>0.11$ Rummukainen et al. (1998). 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. 2. 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. 3. 3.

Identify regions for which $0<x<0.11$ and $y\approx 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. (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 Hindmarsh et al. (2017). The key quantities that we would like to determine nonperturbatively are the strength of the phase transition, $\alpha$, the phase transition temperature, $T_{*}$, and the inverse phase transition duration, $\beta/H_{*}$ Caprini et al. (2016); Weir (2018). The latter two can be found once one knows the bubble nucleation rate, $\Gamma$. 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_{\mathrm{w}}$, for which no nonperturbative studies exist. In our study, we treat it as an input parameter.

III.1 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. Kajantie et al. (1996a) and in Ref. Moore and Rummukainen (2001) 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 1.

Following Ref. Weir (2018) we define the phase transition strength $\alpha$ as the ratio of the latent heat of the transition to the radiation energy density in the symmetric phase at the time of transition

[TABLE]

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\sim T_{*})=106.75$. Note that the gravitational wave spectrum is determined by $\alpha(T_{*})$, not $\alpha(T_{c})$. However, only $\alpha(T_{c})$ has been determined nonperturbatively, so we concentrate on this in what follows and approximate $\alpha(T_{*})\approx\alpha(T_{c})$ which is typically justified absent significant supercooling.

The relation of the latent heat at $T_{c}$ to its 3-d counterpart, $\Delta l_{3}$, can be found in Ref. Kajantie et al. (1996a) and it reads

[TABLE]

where we have introduced the functions

[TABLE]

For a given 4-d theory, these $\eta$-functions determine how the couplings of the 3-d effective theory run with temperature. They are analogous to the usual $\beta$-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 $\mathcal{O(}g^{-2})$, the thermal mass corrections, whereas for $x$ this leads to thermal corrections first at $\mathcal{O(}g^{2})$. Thus the temperature dependence of $x$ is significantly weaker than that of $y$,

[TABLE]

whereas $dy_{c}/dx_{c}=\mathcal{O(}1)$, hence

[TABLE]

In order to calculate the phase transition temperature and inverse phase transition duration, one needs to know the rate of bubble nucleation, $\Gamma$. The 4-d nucleation rate is related to its 3-d counterpart by Moore and Rummukainen (2001),

[TABLE]

where $\sigma_{\mathrm{el}}\sim 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 Bodeker (1998); Arnold et al. (1999). The powers of $g^{2}/(4\pi)$ in the prefactor arise due to the relation between dimensionally reduced lattice quantities and physical quantities. Finally, $\Gamma_{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/\mathrm{e}\approx 0.37$ of the universe remains in the symmetric phase. At this point the following equality holds Enqvist et al. (1992)

[TABLE]

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 Enqvist et al. (1992); Anderson and Hall (1992), yielding what is often termed the nucleation temperature $T_{\mathrm{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 $\beta/H_{*}$ and hence is important if one wants to make accurate predictions.

The inverse duration of the phase transition is then determined by

[TABLE]

where all quantities are evaluated at the transition temperature111Note that in Ref. Moore and Rummukainen (2001) the factor of $\eta_{y}$ was mistakenly dropped, leading to their result for $\beta/H_{*}$ being a factor of $\sim 4$ too small.. 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. Moore and Rummukainen (2001), the quantity $\Gamma_{3-d}(x,y)$ was calculated for $x=0.036$ and is reproduced here in Table 2. It is given as a function of $(y-y_{c})$, corresponding to $\delta m^{2}/g^{4}T^{2}$ in the notation of Ref. Moore and Rummukainen (2001). Fitting a straight line to their results, one finds

[TABLE]

for $(y-y_{c})\in[-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$,

[TABLE]

Details of the specific 4-d theory only enter these expressions through $\eta_{y}$. Interestingly, by combining Eqs. (15) and (16), one finds the relation

[TABLE]

where we have used $\eta_{y}(T_{*})\approx\eta_{y}(T_{c})$. This defines a line in the $(\alpha(T_{c}),\beta/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 $\eta_{y}(T_{c})$,

An analogous relation will hold for other values of $x$: the ratio of $\beta/H_{*}$ to $\alpha$ 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.

III.2 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 $\beta/H_{*}$ and $\alpha$. One can then use the fit formula of Ref. Hindmarsh et al. (2017) 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, $\alpha$, $\beta/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 and the light blue region Fig. 2 can be reliably described by the infrared 3-d SM-like EFT. In Fig. 2 we specifically state why this EFT description fails in different regions, as detailed below. The black line Fig. 2 corresponds to $\alpha$ and $\beta/H_{*}$ values accessible by existing nonperturbative calculations. For reference, we also show the expected sensitivity of LISA Audley et al. (2017) to this parameter space. In Fig. 2, 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. Ellis et al. (2018) for a more detailed discussion of these distinct cases). In Fig. 2, 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 $\eta_{y}$. This quantity, together with $x_{c}$, determines $\alpha(T_{c})$ through Eq. (7) and $\beta/H_{*}$ through Eq. (13) (for $\eta(T_{*})\approx\eta(T_{c})$). The value of $\eta_{y}(T_{c})$ is generically of order $\mathcal{O(}g^{-2})$, absent new physics with large couplings. Its SM value for a $125~{}$GeV Higgs, $\eta_{y}^{\mathrm{SM}}(T_{c})\approx 4.4$, receives corrections from BSM physics,

[TABLE]

where $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ contains corrections which go to zero in the SM limit. In principle $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ can be positive or negative, though all leading-order terms in $\eta_{y}^{\mathrm{SM}}(T_{c})$ are positive definite.

There is a well-defined range of $\eta_{y}(T_{c})$ expected in generic BSM theories mapping on to the 3-d SM-like EFT. First, note that If $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ is negative, $\eta_{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 $\eta_{y}(T_{c})\ll\eta_{y}^{\mathrm{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 $\eta_{y}(T_{c})\gtrsim(g^{2}/4\pi)\eta^{\mathrm{SM}}_{y}(T_{c})\approx 0.15$. Here we have conservatively taken the loop expansion parameter to be $\sim g^{2}/4\pi$: the presence of any larger couplings would strengthen this approximate lower bound. Furthermore, such small $\eta_{y}(T_{c})$, shown as the bottom left arrow in Fig. 2, lead to unobservably, weak transitions from the perspective of their gravitational wave signal.

Conversely, if $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ is positive, one finds stronger transitions with shorter phase transition durations. The size of $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ is nevertheless constrained by perturbativity. The perturbative loop expansion breaks down when couplings become large, limiting $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ to be formally $\mathcal{O}(1/g^{4})$ or smaller. This implies $\eta_{y}(T_{c})\lesssim 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 $\Delta\eta_{y}^{\mathrm{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 $\Delta\eta_{y}^{\mathrm{BSM}}(T_{c})$ is large, the phase transition duration becomes very short, shown as the top right arrow in Fig. 2. 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 $\eta_{y}(T_{c})\in[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.

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. Gorda et al. (2018). In the calculation of $\Gamma_{3-d}$, we follow the perturbative semiclassical analysis outlined in Ref. Moore and Rummukainen (2001), 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$ Weinberg (1993); Moore and Rummukainen (2001).

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. 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). 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}(\phi_{3}^{\dagger}\phi_{3})^{3}$. The leading SM contribution to $c_{6}$ is $\mathcal{O(}y_{t}^{6}/(4\pi)^{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,

[TABLE]

where the expectation values refer to the volume averaged condensates, evaluated at the critical point (see, for example, Ref. Farakos et al. (1995)). 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\gg 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 $\Delta l_{3}$ and $\Gamma_{3-d}$ allows us to trace out the entire region on the $(\alpha,\beta/H_{*})$ plane that can be described by the SM-like 3-d EFT in Eqs. (3)- (4), as argued above. Considering $x\in[0.01,0.06]$ and the aforementioned range of $\eta_{y}(T_{c})$, we arrive at the gray and light blue regions of Figs. 2 and 2, 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. 2. The dark blue region in Fig. 2 is for the specific case of the singlet extended Standard Model (xSM), which we discuss further in Section V.1.

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 and 2, 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 often-studied real singlet extension of the Standard Model (xSM)  Barger et al. (2008); Profumo et al. (2007); Ashoorioon and Konstandin (2009); Espinosa et al. (2012); Robens and Stefaniak (2015); Kanemura et al. (2017, 2016); Chen et al. (2017); Beniwal et al. (2017). Using the dimensional reduction of Ref. Brauner et al. (2017), 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. V, comparing with perturbative results, and discuss the prospects for exploring the phase diagram of this scenario at collider experiments.

IV.1 Parameterization of the model

The most general scalar Lagrangian involving an additional real singlet scalar field can be written as Barger et al. (2008); Profumo et al. (2007); Espinosa et al. (2012); Chen et al. (2015, 2017)

[TABLE]

where $D_{\mu}$ is the usual covariant derivative for the Higgs doublet, and parameters are renormalized in the $\overline{\text{MS}}$ scheme. In order to relate these $\overline{\text{MS}}$ parameters to physical quantities, we reparametrize the potential as follows. After electroweak symmetry breaking, we assume that the singlet field $\sigma$ has a zero-temperature vev $\langle\sigma\rangle=0$, and that the Higgs vev is $\langle\phi^{\dagger}\phi\rangle=v^{2}/2$ . This leads to the following relations in the electroweak vacuum:

[TABLE]

Furthermore, by diagonalizing the potential and solving for the mass eigenvalues – with the eigenstates denoted $h_{1}$ and $h_{2}$ – the parameters $\lambda,\mu$ and $a_{1}$ can be eliminated in favor of physical masses $m_{1},m_{2}$ and the mixing angle $\theta$. 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 $\lambda_{221}$ as an input parameter, which is one of the couplings relevant for scalar pair production processes at colliders Curtin et al. (2014); Craig et al. (2016); Chen et al. (2015, 2017). We treat the singlet self-couplings $b_{3}$ and $b_{4}$, given at a fixed $\overline{\text{MS}}$ scale, as input parameters. As a subset of the general model parameter space, we will also consider the case in which $\sigma\rightarrow-\sigma$ under a discrete symmetry so that $a_{1}$, $b_{1}$, $b_{3}$ and $\sin\theta$ 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. Kajantie et al. (1996c); Laine et al. (2017); Niemi et al. (2018). 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.

IV.2 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. Brauner et al. (2017) and the results of Sections II-III 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 Profumo et al. (2007). This can be seen explicitly in the expression in Ref. Brauner et al. (2017) for the 3-d Higgs quartic coupling after integrating out the superheavy degrees of freedom, $\lambda_{h,3}$, in terms of the model parameters:

[TABLE]

The presence of the $Z_{2}$-breaking couplings $a_{1},b_{1},b_{3}$ lead to a tree-level reduction of $\lambda_{h,3}$, whereas absent these terms the leading effect starts at $\mathcal{O}(g^{4})$ in the power counting scheme used in the dimensional reduction of Ref. Brauner et al. (2017).

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\theta=0$, $b_{3}=0$, and we can simply vary $m_{2}$ and $\lambda_{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. 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 singlet-induced dimension-6 operator coefficient, $c_{6}$, becomes larger than 10 times the corresponding SM value, $c_{6,{\rm SM}}$, assuming $T\approx 140$ GeV. An approximate expression for $c_{6}$ can be found in Eq. (38). As discussed in Appendix B, for $c_{6}\gtrsim 10c_{6,{\rm SM}}$ the neglected higher dimension operators are expected to have $\mathcal{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 $\lambda_{221}$ above this contour. In the gray shaded region, the singlet mass parameter $\mu^{2}_{\sigma}<0$ and the superheavy dimensional reduction completely breaks down. We also show a contour for which the singlet mass parameter $\mu_{\sigma}=200$ GeV, above which our assumption $\mu_{\sigma}\sim\pi 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 higher-loop zero-$T$ effects become non-negligible. Below the $\mu_{\sigma}=200$ GeV and $c_{6}\gtrsim 10c_{6,{\rm 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 $\lambda_{221}$. This validates the conclusions drawn from perturbation theory in e.g. Refs. Craig et al. (2016); Curtin et al. (2014). Note also that similar results for the phase diagram in the real triplet-extended SM were obtained in Ref. Niemi et al. (2018).

We now move on to the more general case without a $Z_{2}$ symmetry. We first fix $b_{3}=0$ and vary $\sin\theta$ and $\lambda_{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 $(\phi^{\dagger}\phi)^{3}$ operator in the 3-d EFT to explore one-step phase transitions. We also show where $c_{6}=10c_{6,{\rm 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 $\mu^{2}_{\sigma}<0$ for which the superheavy dimensional reduction completely breaks down. Also shown is a contour indicating where $\mu_{\sigma}=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\theta=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. Chen et al. (2017) (see also Ref. Gonderinger et al. (2010)).

As can be seen from Figs. 3-5, the neglected singlet-induced dimension-6 operator can become important in the first-order transition regions, especially for small $|\sin\theta|$ and in the $Z_{2}$ limit. This is expected from simple power counting arguments. Comparing Eq. (22) and Eq. (38) for the approximate dimension-6 operator coefficient, we see that generally, for superheavy singlets such that $\mu_{\sigma}\sim\pi T$, the effect of the singlet on the dimension-6 operator is suppressed by $\sim a_{2}/\pi^{2}$ relative to the leading effect on $\lambda_{3}$ (and hence the strength of the phase transition). However, for $|\sin\theta|\approx 0$, the effect on $\lambda_{h,3}$ starts at $\mathcal{O}(a_{2}^{2}/16\pi^{2})$, rather than at $\mathcal{O}(g^{2}/8\pi^{2})$ for $a_{1}\sim gT$, and larger values of $a_{2}$ are required to achieve a given value of $x$. Thus, for small $|\sin\theta|$, the $a_{2}/\pi^{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\theta|$ and $Z_{2}$ regions. As discussed in Sec. V below, these are precisely the regions requiring future colliders to experimentally probe.

IV.3 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_{\rm eff}$. To one loop order, and in a mass-independent renormalization scheme, it is given by

[TABLE]

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, $\mu_{R}$ is the renormalization scale, and the $J_{\pm}$ encode the thermal corrections (see e.g. Ref. Quiros (1999) 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 Curtin et al. (2018). In what follows, we use the high-$T$ thermal masses, the expressions for which can be found in e.g. Refs. Chen et al. (2017), as this is the most conventional approach. The effective potential is computed in Landau gauge utilizing the $\overline{\text{MS}}$ scheme.

With $V_{\rm eff}$ defined above, we search for first-order electroweak phase transitions in the same way as that detailed in Ref. Chen et al. (2017), 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. Patel and Ramsey-Musolf (2011) 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}\geq\zeta$. The value of $\zeta$ 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 $\zeta\approx 0.6-2$ (see e.g. Refs. Quiros (1999); Patel and Ramsey-Musolf (2011) for a more detailed discussion). From the standpoint of observable gravitational waves, $\zeta$ typically needs to be larger, although the strength of the phase transition in this case is usually parameterized by the more physical quantity $\alpha$. For the purpose of benchmarking perturbation theory, we consider values of $\zeta$ 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

[TABLE]

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}\in[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 orange shaded regions of Figs. 4-5 show the parameter space for which $v_{c}/T_{c}\in[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. (22). These regions line up quite well with the first-order EWPT parameter space predicted by our nonperturbative analysis, especially where the $(\phi^{\dagger}\phi)^{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}\gtrsim 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.

V.1 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 III.2. We show the values of $T_{c}$, $T_{n}$, $\alpha(T_{c})$ and $\beta/H_{*}$ calculated using our nonperturbative approach (NP) for a representative sample point in Table 3. Also shown are the same quantities calculated using the usual 4-d perturbative approach (4-d PT) – see Section IV.2 – as well as our two-loop, fully-perturbative approach utilising the DR (3-d PT) – see Section III.2.

A comparison of the nonperturbative and perturbative results in Table 3 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 Linde (1980); Gross et al. (1981); Weinberg (1993); Gleiser et al. (1993); Alford and March-Russell (1994); Surig (1998); Garbrecht and Millington (2015); Metaxas and Weinberg (1996); Baacke and Heitmann (1999); Patel and Ramsey-Musolf (2011); Garny and Konstandin (2012); Strumia and Tetradis (1999). There is nevertheless a 10% discrepancy in $T_{c}$, a factor of 3 discrepancy in $\alpha(T_{c})$ and a 30% discrepancy in $\beta/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 Niemi et al. (2018). The conclusion of the analysis in Appendix B is that we expect our uncertainties in $T_{c}$, $\alpha(T_{c})$ and $\beta/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}\approx 16c_{6,{\rm 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 3 the renormalization scale was set to $\mu_{R}=(m_{1}+m_{2})/2$. Varying $\mu_{R}$ between $m_{Z}$ and $2m_{2}$, we find uncertainties of $\sim 5\%$ in $T_{c}$, $T_{n}$, $\sim 50\%$ in $\alpha(T_{c})$, and $\sim 80\%$ in $\beta/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 3.

From Table 3 and the fit formula of Ref. Hindmarsh et al. (2017) 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 Kawamura et al. (2011) and BBO Harry et al. (2006) though the strength of the signal is too weak to be detected. We plot this point in Fig. 2, 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. Kozaczuk (2015), which found $v_{w}\sim 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 $\eta_{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. 2, indicating values of $\alpha$ and $\beta/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. Alves et al. (2018a)).

V.2 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 Profumo et al. (2007); Noble and Perelstein (2008); Assamagan et al. (2016). Some of these are directly related to $\lambda_{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 Craig et al. (2016); Curtin et al. (2014); Chen et al. (2017), and the shift in the Higgs couplings from their SM predicted values Profumo et al. (2007); Katz and Perelstein (2014); Profumo et al. (2015), most notably the $h_{1}ZZ$ coupling Curtin et al. (2014); Huang et al. (2016b). Complementary information about the scalar potential can be inferred from di-Higgs processes Noble and Perelstein (2008); Huang et al. (2016a); Kotwal et al. (2016); Huang et al. (2016b, 2017) as well as direct production of the singlet and subsequent decays to SM gauge bosons or fermions Buttazzo et al. (2015) 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 Craig et al. (2016); Curtin et al. (2014). An $e^{+}e^{-}$ collider operating as a Higgs factory is expected to achieve excellent precision in measuring the $Zh$ cross-section, $\sigma_{Zh}$ (where $h$ is the SM-like Higgs). The parameter space lying above the blue dashed contour feature $\delta\sigma_{Zh}>0.5\%$, computed at one-loop using the expressions found in Ref. Chen et al. (2017). 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\rightarrow h_{2}h_{2}jj$ at a 100 TeV collider with 30 ab*-1* integrated luminosity. These projections are taken from Ref. Curtin et al. (2018), 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, $\delta\lambda_{111}\equiv|\lambda_{111}-\lambda_{111,{\rm SM}}|/\lambda_{111,{\rm 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* Abada et al. (2019c). 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. Robens and Stefaniak (2015), 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. Aaboud et al. (2018) and the narrow width approximation, with the $h_{2}$ production cross-section taken from Refs. Dittmaier et al. (2011); Andersen et al. (2013) for a SM-like Higgs of the same mass, rescaled by $\sin^{2}\theta$. We also show projected sensitivities for the same searches ($pp\rightarrow h_{2}\rightarrow VV$) at the 13 TeV HL-LHC with 3 ab*-1* (orange dashed contours), obtained from a simple rescaling of the sensitivity in Ref. Aaboud et al. (2018). 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. 6. The blue dashed contours indicate regions for which $\delta\sigma_{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\rightarrow h_{2}\rightarrow h_{1}h_{1}$). Ref. Kotwal et al. (2016) suggests that a 100 TeV collider with 30 ab*-1* integrated luminosity will be able to probe $\sigma\times$BR$(pp\rightarrow h_{2}\rightarrow h_{1}h_{1})$ down to $\sim 100$ fb at 95% C.L. for $m_{2}\approx 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. Hig ). 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 $|\sin\theta|$. 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\rightarrow h_{2}\rightarrow VV,h_{1}h_{1}$ and $e^{+}e^{-}\rightarrow 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 pre-existing lattice results in the dimensionally reduced effective theory for equilibrium Kajantie et al. (1996a) and non-equilibrium thermodynamic quantities Moore and Rummukainen (2001). These relatively few quantities are then used as input parameters for the analysis of numerical simulations of the coupled fluid-scalar field system Hindmarsh et al. (2017) 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 $\mathcal{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.

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.

Appendix A Details of relating renormalized parameters to physical observables

As described briefly in Section IV.1, 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 $\overline{\text{MS}}$ masses) of the mass eigenstates $h_{1}$ and $h_{2}$, as well as the mixing angle $\theta$ and the couplings $b_{3},b_{4},\lambda_{221}$, which are assumed to be inputed directly at a fixed $\overline{\text{MS}}$ scale. We parameterize the scalar doublet as

[TABLE]

with $\omega^{\pm},z$ being the Goldstone modes.

After rotating the scalar potential (20) to a diagonal basis via

[TABLE]

the mass-eigenvalue relations can be inverted to express the $\overline{\text{MS}}$ parameters in the form

[TABLE]

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

[TABLE]

Inverted relations for the top Yukawa $y_{t}$ and gauge couplings ${g,g^{\prime}}$ remain unchanged from the Standard Model case.

We follow Refs. Kajantie et al. (1996c); Laine et al. (2017); Niemi et al. (2018) 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

[TABLE]

where $\Pi$ 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 $\Lambda$ 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

[TABLE]

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 $\hat{\alpha}$, in the $\overline{\text{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 Erler (1999) $\hat{\alpha}^{-1}(M_{Z})\approx 128$, which we input into our parameterization.

The self-energies are renormalized by the counterterms that were previously used in Ref. Brauner et al. (2017). 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 $\hat{\alpha}$ into the tree-level relations (28)-(33) gives the potential parameters at the input scale $\Lambda=(M_{1}+M_{2})/2$ – and similarly for $y_{t},g$ and $g^{\prime}$.

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 $\lambda$, which all play an important role in strengthening the phase transition. However, for masses $M_{2}\gtrsim 400$ GeV the loop correction to the singlet mass parameter $\mu^{2}_{\sigma}$ can be large enough to cause deviations of order $100\%$ from the tree-level value if additionally $\lambda_{221}/v\gg 1$ or $a_{1}/v\gg 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.

Appendix B Accuracy of the dimensional reduction

In Ref. Brauner et al. (2017) 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) Andersen et al. (2018); Gorda et al. (2018) or the real-triplet extension of the SM Niemi et al. (2018). 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 Chen et al. (2017) 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\pi e^{-\gamma}T\approx 7.05T$ – corresponding to the average momentum of integration of the superheavy Matsubara modes Farakos et al. (1994) – 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 Kajantie et al. (1996c), also in Ref. Brauner et al. (2017) 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}(\phi_{3}^{\dagger}\phi_{3})^{3}$ and analyze the effective potential with this operator in 3-d perturbation theory. A similar approach was recently used in Ref. Gorda et al. (2018) (see Section 3.4.2 and Appendix C.7 therein for details).

Using the tools of Ref. Brauner et al. (2017), we can estimate the matching coefficient for the $(\phi_{3}^{\dagger}\phi_{3})^{3}$ operator. In particular, utilizing the effective potential, Eqs. (3.27) and (3.29) in Ref. Brauner et al. (2017), we obtain $V_{6,0}(\varphi^{\dagger}\varphi)^{3}$ with

[TABLE]

where $J_{3}(\mu_{\sigma})\equiv 1/(2\pi^{2}T^{2})J_{B}^{\prime\prime\prime}(\mu^{2}_{\sigma}/T^{2})$, where $J_{B}$, called $J_{+}$ in our analysis, is given e.g. by Eq. (B1) of Ref. Niemi et al. (2018). The other relevant master integrals can be found in Appendix B of Ref. Brauner et al. (2017) and in Ref. Niemi et al. (2018). The first three terms above are the SM contributions and can be compared to Eqs. (196) and (197) in Kajantie et al. (1996c). We emphasize that whereas the sumintegral $I^{4b}_{3}$ contains only contributions from non-zero Matsubara modes, the sumintegral $\tilde{J}^{4b}_{3}(\mu_{\sigma})$ 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-$\sigma$-reducible diagrams. However, in the non-$Z_{2}$-symmetric case, there are already contributions at tree-level

[TABLE]

and we use these contributions to estimate the leading order effect. Therefore, our estimate of leading behavior reads

[TABLE]

With this result, and using the Landau gauge effective potential of Ref. Gorda et al. (2018), we observe that for the point in Table 3 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 $\alpha$ and $\beta/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. Kajantie et al. (1996c), the shift in the VEV scales approximately linearly with $c_{6}/\lambda_{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 Kajantie et al. (1996c), and so we expect the change in the Higgs VEV to be $\mathcal{O}(c_{6}/c_{6,{\rm SM}}\times 1\%)$ for relatively weak transitions. On Figs. 3-5 we show where $c_{6}=10c_{6,{\rm SM}}$, assuming $T\sim 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.