Collider and Gravitational Wave Complementarity in Exploring the Singlet Extension of the Standard Model

We present a dedicated complementarity study of gravitational wave and collider measurements of the simplest extension of the Higgs sector: the singlet scalar augmented Standard Model. We study the following issues: (i) the electroweak phase transition patterns admitted by the model, and the proportion of parameter space for each pattern; (ii) the regions of parameter space that give detectable gravitational waves at future space-based detectors; and (iii) the current and future collider measurements of di-Higgs production, as well as searches for a heavy weak diboson resonance, and how these searches interplay with regions of parameter space that exhibit strong gravitational wave signals. We carefully investigate the behavior of the normalized energy released during the phase transition as a function of the model parameters, address subtle issues pertaining to the bubble wall velocity, and provide a description of different fluid velocity profiles. On the collider side, we identify the subset of points that are most promising in terms of di-Higgs and weak diboson production studies while also giving detectable signals at LISA, setting the stage for future benchmark points that can be used by both communities.


Introduction
Since the first direct detection of gravitational waves (GWs) by the LIGO and Virgo collaborations [1], a new interface has arrived in particle physics -its intersection with GW astronomy. While ground based GW detectors have their best sensitivity at frequencies ∼ O(100) Hertz and their main targets are black hole and neutron star binaries, there is now growing interest in building space-based interferometer detectors for milli-Hertz or deci-Hertz frequencies. Many detectors have been proposed, such as the Laser Interferometer Space Antenna (LISA) [2], the Big Bang Observer (BBO), the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [3], Taiji [4] and Tianqin [5]. The physical sources of GWs in this frequency band include supermassive black hole binaries [6], extreme mass ratio inspirals [7] and the stochastic background of primordial GWs produced during first order cosmological phase transitions [8].
This offers tremendous opportunities for theorists, as a new window to the early Universe opens up. Aspects of dark sector physics and baryon asymmetry can now be framed fruitfully in a language that lends itself to data from the GW frontier. The key connection is phase transitions, which on the one hand are a primary target of future GW experiments, and on the other are important features of scalar potentials and hence have historically been the target of collider physics.
The purpose of our work is to explore the complementarity of future GW detectors and future particle colliders in probing phase transitions in the early Universe -in the simplest particle physics setting possible, but also with great attention to detail within such a setting. The natural choice is the electroweak phase transition (EWPT) [9] with the simplest extension of the Higgs sector: the singlet scalar augmented Standard Model or the xSM 1 . This model is capable of providing a strongly first order EWPT through a tree level barrier and is the simplest model in Class IIA of the tree level renormalizable operators described in [19] (see Ref.  for related studies). It has been extensively investigated in phenomenological studies [46][47][48][49], studies of EWPT [46,47,[50][51][52] and di-Higgs analyses [53] guided by the requirements of EWPT [48], and electroweak baryogenesis (EWBG).
We perform a detailed scan of this model, shedding light on the following issues: (i) the EWPT patterns admitted by the model, and the proportion of parameter space for each pattern; (ii) the regions of parameter space that give detectable GWs at future space-based detectors; (iii) the current and future collider measurements of di-Higgs production, as well as searches for a heavy weak diboson resonance, and how these searches interplay with regions of parameter space that exhibit strong GW signals; and (iv) the complementarity of collider and GW searches in probing this model.
We first carefully work out and incorporate all phenomenological constraints: boundedness of the Higgs potential from below, electroweak vacuum stability at zero temperature, perturbativity, perturbative unitarity, Higgs signal strength measurements and electroweak precision observables. Then, we identify the regions of parameter space which give large signal-to-noise-ratio (SNR) at LISA. We carefully address subtle issues pertaining to the bubble wall velocity v w , making a distinction between v w , which enters GW calculations, and the velocity v + that is used in EWBG calculations. The relation between these two velocities is determined from a hydrodynamic analysis by solving the velocity profile surrounding the bubble wall. We provide a description of different fluid velocity profiles and investigate the behavior of the normalized energy released during the phase transition, α, which primarily determines the SNR, as a function of the model parameters. On the collider side, we identify the subset of points with large SNR at LISA that are most promising in terms of di-Higgs and weak diboson production studies, setting the stage for future benchmark points.
Much remains to be understood about the Higgs sector. On the collider side, measuring the Higgs cubic and quartic couplings through double or triple Higgs production, both non-resonant as well as resonant, is an extremely difficult but central goal of future experiments (see e.g., [54][55][56][57][58][59][60]). While any deviation of the shape of the Higgs potential from what is expected within the Standard Model (SM) would hint to new physics, the sensitivities of such collider studies are found to be rather low. The detection of GWs from EWPT in future experiments can offer a complementary method of probing the currently largely unknown Higgs potential. Our work is a step in that direction.
The paper is structured as follows. In Sec. 2, we define the Higgs potential and set the notations. The standard phenomenological analysis is discussed in the following Sec. 3. The next Sec. 4 discuss the details of the EWPT and GW calculations, after which the results and discussions from the full scan is presented in Sec. 5 and we summarize in Sec. 6.

The Model
In this section, we fix our notation by defining the potential for the gauge singlet extended SM, known as the"xSM". This model is defined with the following potential setup [46][47][48]: where H T = (G + , (v EW + h + iG 0 )/ √ 2) is the SM Higgs doublet and S = v s + s the real scalar gauge singlet. All the model parameters in the above equation are real. The parameters µ and b 2 can be solved from the two minimization conditions around the EW vacuum(≡ (v EW , v s )), and λ, a 1 , a 2 can be replaced by physical parameters θ, m h 1 and m h 2 from the mass matrix diagonalization 2 : with a mixing angle θ. We note that h 1 is identified as the SM Higgs while h 2 is a heavier scalar. The coupling of h 1 with the SM particles is reduced by a factor of c θ while the coupling of h 2 2 Here s θ ≡ sin θ and c θ ≡ cos θ.
with SM particles is (−s θ ) times the corresponding SM couplings and vanishes in the case of zero mixing angle. With choices of parameter transformations described above, the potential is fully specified by the following five parameters: The model defined here has several variants in the literature. For example, since the potential can be defined with a translation in the S direction S → S = S − v s , such that S = 0, the resulting potential will take the same form as Eq. 2.1 but with the addition of a non-zero tadpole term b 1 S [53]. The potential and physics remain the same but the parameters in the potential will transform accordingly. The transformation rules to and from this basis are given in Appendix B. There is also a variant where there is a spontaneously broken Z 2 symmetry S → −S; this corresponds to a subset of the parameter space here where a 1 = b 3 = 0.
We further note that we do not include CP-violation in this study since the magnitude of the CP-violation is typically very constrained by current electric dipole moment searches (e.g., [36,61,62] or the included CP-violation may be large but has little effect on EWPT [63].

Phenomenological Constraints
In this section, we briefly discuss the phenomenological constraints used in our analysis, following the standard treatments given in Refs. [49,53,64]. The phenomenological discussion includes boundedness of the Higgs potential from below, EW vacuum stability at zero temperature, perturbativity, perturbative unitarity, Higgs signal strength measurements and electroweak precision observables.
First, the potential needs to be bounded from below. Requiring this for arbitrary field directions gives us the condition [53], λ > 0, b 4 > 0, a 2 −2 λb 4 . (3.1) Next, the EW vaccum also needs to be stable at zero temperature. Using physical parameters as input will automatically guarantee that the EW vacuum is a minimum. To ensure that the above EW vacuum is stable, one should require that no deeper minimum exists in the potential. In our analysis, we find all the minima by firstly solving ∂V /∂φ i = 0(φ 1 ≡ h, φ 2 ≡ s) and subsequently calculating eigenvalues of the Hessian matrix {∂ 2 V /∂φ i ∂φ j } to determine the nature of the extrema for each set of parameter input. Next, Higgs signal strength measurements in various channels require the couplings of h 1 to be not far from the SM Higgs couplings. In the xSM, the couplings of h 1 to SM particles are reduced by a factor of cos θ, therefore the Higgs signal strength is given by µ H = cos 2 θ. Experimentally, the most recent ATLAS and CMS combined fit of this value is µ H = 1.09 +0. 11 −0.10 [65] and a χ 2 analysis shows that | sin θ| > 0.33 are excluded at 95% CL [66].
Moreover, unitarity puts constraints on the high energy behavior of particle scatterings. Requiring further the perturbativity of these scatterings at high energy will lead to constraints on the model. This tree level perturbativity requirement is quantified as the condition that the partial wave amplitude a l (s) for all 2 → 2 processes satisfies |Re a l (s)| 1/2 for √ s → ∞. We consider all channels of scalar/vector boson 2 → 2 scatterings at the leading order in the high energy expansion, with details of the S-matrix given in Appendix. A.
Electroweak precision measurements, which mainly include the W boson mass measurement [67] and the oblique EW corrections [68,69], put further constraints on the model. The W boson mass m W can be calculated given experimentally measured values of G F , m Z and the fine structure constant at zero momentum transfer α(0) [67]. The function relating m W and these three parameters depends on the loop corrections of the vector boson self-energies. Comparing this calculated m W with the experimental measurement m exp W = 80.385 ± 0.015GeV [70][71][72] highly constrains the modification of the loop corrections by new physics effects. In this model, the modified loop corrections result from reduced Higgs couplings and from the presence of the heavier scalar h 2 and are only dependent on (θ, m h 2 ) at one-loop level. The same parameter dependence enters the oblique S, T, U parameters and it turns out that the W -mass constraint is much more stringent than that from the oblique corrections [49,67].
To give the reader a flavor of the above phenomenological constraints, we fix m h 2 = 300 GeV, θ = 0.2, b 4 = 4 and show the various bounds on the remaining two parameters (v s /v EW , b 3 /v EW ) in Fig. 1. This choice of m h 2 and θ evades the constraints from the W -mass as well as the oblique EW corrections and regions outside the color-shaded regions are excluded by the remaining constraints. It can be seen from this figure that the least constraining condition comes from the pertur-bative unitarity requirement for this parameter choice. The bounded-from-below condition is more restrictive and also separates the plane into two disconnected regions while the stability of the EW vacuum at zero temperature shrinks the allowed parameter space even more. We also overlaid on this plot the points which pass the various EWPT requirements and give GW signals with varying SNR. More details are given in the caption and in the following section.

Effective Potential
EWPT is an essential step in generating the observed baryon asymmetry in the universe by providing an out-of-equilibrium environment, one of the three Sakharov conditions [73], in the framework of electroweak baryogenesis (see [74] for a recent review). Augmented with the rapid baryon number violating Sphaleron process outside the electroweak bubbles and the CP-violating particle scatterings on the bubble walls, a net baryon number can be produced inside the bubbles. Aside from the particle interactions, which are used in EWBG calculations, the cosmological context that characterizes the dynamics of the EWPT can be calculated from the finite temperature effective potential. The standard procedure of calculating it includes adding the tree level effective potential, the Coleman-Weinberg term [75] and its finite temperature counterpart [76] as well as the daisy resummation [77,78]. Since the EWPT in this model is mainly driven by the cubic terms in the potential and out of concern of a gauge parameter dependence [79] of the effective potential calculated in the above standard procedure, we take here the high temperature expansion approximation, which is gauge invariant, in line with previous analyses of this model [46-48, 50, 80]. This effective potential is then given by 3 where Π h (T ) and Π s are the thermal masses of the fields, where the gauge and Yukawa couplings have been written in terms of the physical masses of W , Z and the t-quark. With this effective potential, the thermal history of the EW symmetry breaking can be analyzed. It depends mainly on the following key parameters: Here T c is the critical temperature at which the metastable vacuum and the stable one are degenerate. Below T c , the phase at the origin in the field space becomes metastable and the new phase becomes energetically preferable. The rate at which the tunneling happens is given by [81] Γ ∼ A(T )e −S 3 /T , (4.4) where S 3 is the 3-dimensional Euclidean action of the critical bubble, which minimizes the action and satisfies the bounce boundary conditions Here φ out denotes the two components vev of the fields outside the bubble, which is not necessarily the origin for two-step EWPT. The prefactor A(T ) ∝ T 4 on dimensional grounds. Its precise determination needs integrating out fluctuations around the above static bounce solution (see e.g., [82,83] for detailed calculations or [84] for a pedagogical introduction). For the EWPT to complete, a sufficiently large bubble nucleation rate is required to overcome the expansion rate. This is quantified as the condition that the probability for a single bubble to be nucleated within one horizon volume is O(1) at a certain temperature [85]: where V H (t) is the Horizon volume, M Pl is the Planck mass and ζ ∼ 3×10 −2 . From this equation, it follows that S 3 (T )/T ≈ 140 [86] and the temperature thus solved is defined as the nucleation temperature T n . Expanding the rate at T n , one can define the duration of the EWPT in terms of the inverse of the third parameter β [86]: where H n is the Hubble rate at T n . Next, α is the vacuum energy released from the EWPT normalized by the total radiation energy density (≡ ρ R ) at T n [87]: where ρ R = g * π 2 T 4 n /30 with g * ≈ 100 and φ b denotes the two components vev of the broken phase. In this expression, the first term is the free energy from the effective potential and the second term denotes the entropy production. Finally, v w is the bubble wall velocity.
Given that a first order EWPT can proceed and complete, the baryon asymmetry is generated outside the bubbles and then captured by the expanding bubble walls. When the EWPT finishes, the universe would be in the EW broken phase with non-zero baryon asymmetry. To ensure that these baryons would not be washed out, the Sphaleron rate needs to be sufficiently quenched inside the bubbles. This condition is known as the strongly first order EWPT (SFOEWPT) criterion [74,88]: The conventional choice of the temperature at which the above condition is evaluated is T c , but a more precise timing is the nucleation temperature T n , which we use here. Since generally T n < T c and v h (T n ) > v h (T c ), it might seem at first glance that the above condition is weaker when implemented at T n than at T c . However the implicit assumption associated with the former requires the capability of the EWPT to successfully nucleate, i.e., the condition Eq. 4.7 should be satisfied in the first place, which is typically a more stringent requirement of the potential. The presence of two scalar fields gives a richer pattern of EWPT and makes it possible to complete the EWPT with more than one step [85,89,90]. One can immediately imagine mainly the following EWPT types: where the last vacuum configuration (v H = 0, v S = 0) in each case would eventually evolve to the EW vacuum at T = 0 4 . Here pattern (A) is a one step EWPT from the origin in field space to the EW symmetry breaking vacuum directly, due mainly to the negative cubic term in the effective potential. This one step phase transition results in a typical GW spectrum as shown in the left panel of Fig. 3. Quite differently, patterns (B) and (C) are two-step EWPT, which differ only in how the vacuum transits for these two steps. For example, in case (B), the universe first goes to a vacuum which has non-zero vev for the singlet field and then transits to the would-be EW vacuum at high temperature. Case (C) is different in that it breaks the EW vacuum first and then further goes to the would-be vacuum in a subsequent step of phase transition. For each transit of the vacuum, it can be either first or second order, depending on whether there is a barrier separating the two vacua. We note that for case (C), baryon production generally needs to occur in the first step, otherwise, the exponentially reduced Sphaleron rate would greatly suppress the baryon number violating process in the second step as the EW symmetry is already broken outside the bubbles. Therefore the SFOEWPT criterion is imposed in the first step for this case.
We note that with the aid of the analytical methods presented in Ref. [51,85], it is possible to locate the region of the parameter space that gives exactly one specific type of EWPT by imposing various conditions on the input parameters. However, our task here is to reveal the overall behavior of the parameter space concerning EWPT and GW. Therefore we adopt here a scan-based analysis which covers the entire parameter space and for each scanned parameter space point, we determine its pattern of EWPT and calculate GW properties. This way, we can determine the most probable pattern of EWPT admitted by this model.

Hydrodynamics
Successful EWBG usually requires a subsonic v w to give sufficient time for chiral asymmetry propagation ahead of the wall and for conversion to baryon asymmetry through the Sphaleron process. On the other hand, a larger v w is better for GW production. Therefore a tension may arise between successful EWBG and a loud GW signal production. This problem can potentially be solved when the hydrodynamic properties of the fluid are taken into account [91]. This is because the expanding wall stirs the fluid surrounding the bubble wall and a non-zero velocity profile exists for the plasma ahead of the wall (see Ref. [92] for a recent combined analysis). In the bubble wall frame, this means the plasma outside the bubble will head towards the bubble wall with a velocity (≡ v + ) that can be different from v w . Therefore it is v + rather than v w that should be used in EWBG calculations. While the above argument still needs to be scrutinized taking into account the particle transport behavior around the bubble wall in the process of EWBG, we assume tentatively that this is true in this work. This hydrodynamic treatment hinges on solving the fluid velocity profile v(r, t) around the bubble wall given inputs of (α, v w ), where r is the distance from the bubble center and t is counted from the onset of the EWPT. Due to the properties of the problem here, v is a function solely of r/t ≡ ξ. The differential equation governing the velocity profile is derived from the conservation of the energy momentum tensor describing the fluid and scalar field [92]: where c s = 1/ √ 3 is the speed of sound in the plasma and µ(ξ, v) = (ξ − v)/(1 − ξv) is a Lorentz boost transformation. Far outside the bubble and deep inside the bubble, the plasma will not be stirred, that is v → 0 serves as the boundary condition. At the phase boundary, the velocity of the plasma inside and outside the bubble wall are denoted as v − and v + in the bubble wall frame, both heading towards the bubble center. The same energy momentum conservation, when applied across the bubble wall, gives a continuity equation connecting v − with v + . Therefore the whole fluid velocity profile can be solved from the center of the bubble to far outside the bubble where the plasma is unstirred.
The solutions of the fluid profiles can be classified into three modes depending on the value of v w . A set of profiles v(ξ) are shown in Fig. 2 for α = 0.1. For v w < c s , a deflagration mode is obtained, in which case, the plasma ahead of the bubble wall flows outward while it remains static inside the bubble, corresponding to the profiles with blue-dashed lines. It can also be seen from this figure that as v w increases in this mode, a discontinuity in v(ξ) appears outside the bubble and v(ξ) jumps to zero. This is the location of the shock front, and beyond this point the solution of Eq. 4.11 is invalid and a shock front develops such that v(ξ) goes to zero consistently. When v w surpasses c s but is less than a certain threshold ξ J (α), a supersonic deflagration mode [93] appears (magenta solid profiles) where the plasma inside the bubble has a non-zero profile, while still taking the form of deflagration outside the bubble. Here ξ J (α), as a function of α, corresponds to the Jouguet detonation [94], used in earlier studies. It is also evident that in this mode, as v w increases, the shock front becomes closer to the bubble wall until it coincides with the bubble wall, where v w = ξ J (α) and the fluid enters the third, detonation mode (brown dotted profiles). In this mode, the plasma outside the bubble has zero velocity and therefore v + = v w . If a subsonic velocity is required in EWBG, we conclude that the deflagration mode will not work for EWBG. On the contrary, v + < v w in the deflagration and supersonic deflagration modes and a solution for the tension between EWBG and GW might be achieved.
Therefore, instead of treating v w as a free parameter in the GW calculations, we require, given a certain input of α, the corresponding v + to have subsonic value, taken to be 0.05 here, a choice usually used in EWBG calculations [63,[95][96][97][98]). The procedure of achieving the above goal is as follows: for each given α we iterate over v w and solve the whole fluid profile until v + = 0.05 is reached. The resulting v w is used in GW calculations 5 .
With v(ξ) obtained, one can also calculate the bulk kinetic energy normalized by the vacuum energy released during the EWPT [92]: where ω(ξ) is the enthalpy density, varying as function of ξ, and can be solved once v(ξ) is found. The remaining part 1 − κ v ≡ κ T gives the fraction of the vacuum energy going to heat the plasma. Therefore a reheating temperature can be defined as This leads to an increase in entropy density and thus a dilution of the generated baryon asymmetry [89]. Typically in EWBG calculations, the wall curvature is neglected and the transport equations depend on a single coordinatez in the bubble wall rest frame, wherez > 0 (< 0) corresponds to broken (unbroken) phase. The solved baryon asymmetry density n B is a constant inside the bubbles(see, e.g., [99]): where s(T ) = 2g * π 2 T 3 /45 is the entropy density, Γ ws ≈ 120α 5 w T is the weak Sphaleron rate in the EW symmetric phase [100], λ ± = (v + ± v 2 + + 15Γ ws D q )/(2D q ) with D q the diffusion constant for quarks [100] and n L is the chiral asymmetry of left-handed doublet fields which serves as a source term in baryon asymmetry generation. The determination of n L is a key part in EWBG calculations and is decoupled from the analysis of EWPT dynamics here. In above expression, we have replaced v w by v + , to take into account the distinction between these two velocities. If the temperature at which n B is calculated is T n , then after the bubbles have collided, the temperature of the plasma is given, to a good approximation, by T * rather than T n or T c , which are conventionally used. The diluted baryon asymmetry is then given by where ξ D ≡ (1 + κ T α) −3/4 captures the dilution effect of the generated baryon asymmetry by reheating of the plasma. We then need to make sure that ξ D does not become too small, since otherwise a stronger CP-violation will be needed, which might be excluded by the stringent limits from electric dipole moment searches [101,102].

Stochastic Gravitational Waves
During the EWPT, bubbles of EW broken phase expand and collide with each other, which destroys the spherical symmetry of a single bubble, thus leading to the emission of gravitational waves [87]. Due to the nature of this process and according to the central limit theorem, the generated gravitational wave amplitude is a random variable which is isotropic, unpolarized and follows a Gaussian distribution. This therefore allows the description of gravitational wave amplitude using its twopoint correlation function and is parametrized by the gravitational wave energy density spectrum Ω GW (f ), as a function of frequency f . A natural consequence is that the GWs produced during the EWPT, when redshifted to the present, give a peak frequency at around the mili-Hertz range [9], falling right within the band of future space-based gravitational wave detectors.
It is now well known that there are mainly three sources of gravitational wave production in this process: bubble wall collisions [103][104][105][106][107][108], sound waves in the plasma [109,110] and magnetohydrodynamic turbulence (MHD) [109,110]. The total energy density spectrum can be obtained approximately by adding these contributions: neglected [112]. We thus include only the contribution of sound waves and turbulence in the gravitational wave spectrum calculations.
The dominant contribution comes from sound waves. By evolving the scalar-field and fluid model on 3-dimensional lattice, the gravitational wave energy density spectrum can be extracted, with an analytical fit formula available [110]: (4.17) Here H * is the Hubble parameter at T * when the phase transition has completed. It has a value close to that evaluated at the nucleation temperature T n for sufficiently short EWPT [8]. We take T * to be the reheating temperature, defined earlier in Eq. 4.13. Moreover, f sw is the present peak frequency which is the redshifted value of the peak frequency at the time of EWPT (= 2β/( √ 3v w )): Hz, (4.18) where κ v is defined in Eq. 4.12 and can be calculated as a function of (α, v w ) by solving the velocity profiles described in Sec. 4 [92]. It should be noted that a more recent numerical simulation by the same group [113,114] shows a slightly enhanced Ω sw h 2 and reduced peak frequency f sw . We also note that the results from these simulations are currently limited to regions of small v w and α and therefore their validity for ultra-relativistic v w and large α (say α 1) remains unknown. In the absence of numerical simulations for these choices of parameters at present, we assume that the results shown here apply for these cases and remind the reader to keep the above caveats in mind.
The fully ionized plasma at the time of EWPT can result in the formation of MHD turbulence, which gives another source of gravitational waves. The resulting contribution can also be modelled similarly with a fit formula [115,116], where f turb is the peak frequency and is given by, Hz. (4.20) Here the factor κ turb describes the fraction of energy transferred to the MHD turbulence and is given roughly by κ turb ≈ κ v with ≈ 5 ∼ 10% [110]. We take = 0.1 in this study. In both Eq. 4.17 and 4.19, the value of v w is found by requiring that v + = 0.05 by solving the velocity profiles, as discussed in the previous section. For the two-step EWPT, as discussed in last section, if both steps in case (B) and (C) are first order, then there would be two subsequent GW generation at generally different peak frequencies and amplitudes, corresponding to the example shown in the right panel of Fig. 3.
The detectability of the GWs is quantified by the signal-to-noise ratio (SNR), whose definition is given in Ref. [8]: (4.21) Here h 2 Ω exp (f ) is the experimental sensitivity and corresponds to the lower boundaries of the color-shaded regions in Fig. 3 for the shown detectors 6 . T is the mission duration in years for each experiment, assumed to be 5 here. The factor δ comes from the number of independent channels for cross-correlated detectors, which equals 2 for BBO as well as UDECIGO and 1 for the others [117]. In our numerical analysis, we stick to the most mature LISA detector with the C1 configuration, defined in Ref. [8]. To qualify for detection, the SNR needs to be larger than a threshold value, which depends on the details of the detector configuration. For example, for a four-link LISA configuration, the suggested value is 50 while for a six-link configuration, this value can be much lower (SNR = 10), since in this case a special noise reduction technique is available based on the correlations of outputs from the independent sets of interferometers of one detector [8].
As an example, we scan over the EW vacuum stability regions in the plane Fig. 1 and found the regions which can give successful bubble nucleations, satisfy the SFOEWPT criterion and generate GWs. These regions are plotted with blue (SNR < 10), green (50 > SNR > 10) and red (SNR > 50). Here most of the points give type (A) EWPT with only several points for type (B) or (C), denoted by diamond shapes.

Results and Discussions
In this section, we perform a full scan of the parameter space to address the following questions: (a) What kind of EWPT patterns can this model admit and in what proportion of the parameter space for each pattern?
(b) What is the region of parameter space that can give strong detectable gravitational waves at future space-based gravitational wave detectors?
(c) Do current collider measurements of double Higgs production and searches for a heavy resonance decaying to weak boson pairs exclude the points that give strong gravitational waves and could future high luminosity LHC (HL-LHC) at 3ab −1 probe the parameter space giving strong gravitational waves?
(d) How will a future space-based gravitational wave experiment complement current and future searches for a heavy scalar resonance?
The full scan is performed using the input of the tadpole basis parameters with the following ranges for parameters: where the lower range of a 2 is determined by the requirement that the potential is bounded from below. The scan takes into account the previously discussed theoretical and phenomenological requirements. Points which pass these selection criteria are fed into CosmoTransitions [118] for calculating the thermal history and the parameters relevant for EWPT. Those which can give a successful EWPT by meeting the bubble nucleation criteria are further scrutinized for the EWPT type and SFOEWPT conditions. The final remaining points are used to calculate the gravitational wave spectra, the SNR and collider observables.

EWPT and GW
We first give the answer to question (a): what kind of EWPT patterns can this model admit and in what proportion of the parameter space for each pattern ?
We find, of the xSM parameter space where a successful EWPT can be obtained, about 99% gives type (A) EWPT and the remaining slightly less than 1% can give type (B) EWPT. We do not observe type (C) EWPT. For type (A), 22% (19%) gives SNR larger than 10 (50). So there is a sufficiently large parameter space which can give detectable GW production.
The strength of the stochastic GW background is mainly governed by the two parameters α and β/H n , where a larger α and a smaller β/H n gives stronger GW SNR, as shown in the left panel of Fig. 4, where the colors denote SNR < 10 (blue), 50 > SNR > 10 (green) and SNR > 50 (red). We observe that the points which give detectable GWs lie in the bottom right region of the population. Physically, α quantifies the amount of energy released during the EWPT and therefore a larger α gives stronger GW signals. In addition, for fixed v w , a larger α leads to a larger fraction of energy transformed into the plasma kinetic energy, quantified by κ v , and therefore a further gain in GW production. A further enhancement for larger α comes from the fact that since we fixed v + = 0.05, increasing α also increases v w . It should be noted, even without an explicit calculation, that for each fixed value of α, the allowed values of v w are limited to a certain range (see e.g., Fig. 1 in Ref. [80]). This comes from two considerations: (1) admitting consistent hydrodynamic solutions of the plasma imposes a lower limit on v w ; (2) v w larger than ξ J (α) gives a detonation mode of the velocity profile, in which case v w = v + > c s and therefore v + is too large for EWBG to work. We further note that for α 1 and v w ∼ 1, the calculations of the GW spectra may become unreliable for the following reasons: (i) While the study of Ref. [112] suggests that the energy stored in the scalar field kinetic energy is negligible, a very large α might lead to a non-negligible contribution from the bubble collisions. Therefore a better understanding of the energy budget for this region is needed; (ii) the numerical simulations are all performed for relatively small α as well as v w and thus the use of these results for large α and v w may not be applicable; (iii) The universe is no longer radiation dominated at the EWPT but rather vacuum energy dominated. This has the consequence that bubbles might never meet to finish the EWPT and the universe would be trapped in the metastable phase (see Ref. [119] for a recent analysis). Despite these issues, we find 49% of points with SNR > 10 have α < 1 and removing the points with α > 1 does not change the main findings of our work.
We now turn to the parameter β/H n , which roughly characterizes the inverse time duration of the EWPT. A smaller β/H n or equivalently a longer EWPT generates stronger GW signals. This is due to the particular feature of the GWs coming from the sound waves in the plasma. As was found in the original papers on the importance of sound waves in generating the GWs [109,110], one enhancement comes from 1/(β/H n ) compared with the conventional bubble collision contribution. As long as the mean square fluid velocity of the plasma is non-negligible, GWs will continue being generated and the energy density of the GW is thus proportional to the duration fo the EWPT. It should be noted that β/H n also determines the peak frequency of the GW spectra. The bubble wall velocity v w also plays an important role here and the dependence of the SNR on v w is shown in the middle panel of Fig. 4, where the vertical axis is chosen to be T n . It is clear that points with larger SNR have larger v w since, for fixed v + , a larger α implies a larger v w . It can also be seen from this plot that the SNR increases as T n decreases. This is easily understood, since a smaller T n typically implies a larger amount of supercooling and therefore a larger α. The supercooling can be quantified by the fraction of the first term(≡ ∆ρ V ) of Eq. 4.9 in the total released vacuum energy, which we plot in the right panel. We can see from this figure that larger SNR indeed implies larger amount of supercooling. However the amount of supercooling as quantified by ∆ρ V /∆ρ is less than 0.6 for most of the parameter space. The remaining part comes from the second term of the definition of α.
The entropy production, if sizeable, can pose a problem for baryon asymmetry generation, as it will effectively dilute the baryon asymmetry n B /s by increasing s. In Sec. 4.2, we encode this effect in a dilution factor ξ D . Here since κ T is a function of v w and α while v w is also a function of α when v + is fixed, we find ξ D is solely a function of α. This functional relation is shown as the magenta line in the left panel of Fig. 5 and all points from the scan fall on this line. The message from this figure is that most of the points have ξ D 0.65 and those with a smaller α have a dilution factor closer to 1. In particular, the points with α 1 for which GW can be reliably calculated, the dilution effect is rather small as ξ D 0.8. Given the current relatively large uncertainties in the EWBG calculations, the dilution effect poses no real problem for the baryon asymmetry generation. Note that previous studies [89] used a different quantification of the dilution factor, with the definition: where s is the entropy density at T n and ∆s is calculated from the second term in the definition of D gives an overestimation of the dilution effect while ξ D firstly increases a little bit before slowly dropping.
Since the dilution factor we use here is based on a faithful hydrodynamic analysis, it gives a more precise description of the dilution effect. We also show ξ D calculated for all the points versus T n as a scatter plot in the right panel of Fig. 5, from which we find a larger dilution effect appears for typically smaller T n and those with α 1 fall in the high T n region.
The two-step EWPT, for which type (B) is the only observed here, constitutes about one percent of all the surviving parameter space. Of this tiny parameter space, more than half the points give detectable GWs.

Parameter Space Giving Detectable GWs
With a summary of the points described in previous section, we give in this section the answer to question (b), which, we recall, was: What is the region of parameter space that can give strong detectable gravitational waves at future space-based gravitational wave detectors?
The results are shown in terms of the three plots in Fig. 6. As was discussed in the previous section, a large α and small β/H n leads to loud GW signals. Even though the relation between (α, β/H n ) and the physical input parameters is not transparent as many numerical details are involved, it can still be revealed by the plots in Fig. 6. From the left panel in Fig. 6, we can see that the majority of the points are concentrated in two regions of parameter space where v s is rather small. In particular, we find 20 GeV |v s | 50 GeV for most points, with a peak distribution at around 20 GeV. The appearance of two regions comes from the bounded-from-below requirement of the potential, similar to Fig.1. While phenomenological constraints have the effect of shrinking both the regions, the appearance of points far outside the two regions indeed shows that the main cause of the narrow regions comes from the requirements of EWPT and GWs. Therefore it is fair to say that the region that gives detectable GWs from a type (A) EWPT mainly comes from the parameter space with smaller v s . On the other hand, the regions which provide type (B) EWPT are dramatically different from these regions, since most of the diamonds lie beyond the two narrow regions, as can be seen from the figure.
The middle figure shows these regions in the (m h 2 , sin θ) plane. It is clear that the points are concentrated around the region with larger m h 2 . For smaller m h 2 , the density of points becomes much smaller. To have a better understanding of the role of m h 2 in GW production, we show in the right panel its role in determining (α, β/H n ), denoted by the colors. In this figure, the points are separated into different bands characterized by the value of m h 2 . For fixed β/H n , a larger m h 2 gives a larger α, thus larger SNR. This explains the concentration of the points in the m h 2 direction in the middle figure. In the sin θ direction, the value of θ is more constrained for larger m h 2 . The outer boundary comes mainly from the W -mass constraint. The requirements from EWPT and larger GW signals also show their effects in this plot. For example, very small values of θ give rarer points. We also overlaid on this plot the various sensitivity projections from colliders in probing the value of θ, which includes HL-LHC, ILC with two configurations (ILC-1: 250GeV, 250fb −1 , ILC-3: 1TeV,1ab −1 ) and future circular e + e − colliders (240GeV, 1ab −1 ), all taken from Ref. [47]. We see that HL-LHC can barely probe any points; ILC-1 can probe a fraction of the small m h 2 points as well as a few large m h 2 points; ILC-3 can probe about a half of both light and heavy h 2 points; the future circular colliders can probe even more of the parameter space. We also can see that most of the points coming from the two-step EWPT lie at the very small θ region, even though a few do have larger θ. Therefore GW detections serve as a complementary probe of this region. We also note that for very small values of θ and m h 2 , the search for long lived particles can be used to probe this region (eg., the MATHUSLA detector) [120].

Correlation with Double Higgs Production Searches
Exploring possible deviations from the expected SM value of the cubic Higgs coupling through di-Higgs production is an important target of the HL-LHC. New physics scenarios, especially those designed for providing a SFOEWPT for baryon asymmetry generation, typically modify this coupling. Therefore di-Higgs production is correlated with EWPT and thus GW production. Future GW and collider experiments can then operate in a way that complement each other in exploring new physics scenarios. With the parameter space giving detectable GW identified in the previous section, we can find the correlation by calculating the corresponding di-Higgs cross sections and compare it with present di-Higgs measurements and with future projections. Figure 7. Representative resonant (left) and non-resonant (middle and right) Feynman diagrams contributing to di-Higgs production.
The leading order Feynman diagrams for double Higgs production occur at one-loop and consist of both the resonant and non-resonant channels, as shown in Fig. 7. The non-resonant channel includes the box diagrams and a triangle diagram involving the vertex h 1 h 1 h 1 . The resonant chan- nel is the production of a on-shell h 2 which subsequently decays into two Higgs, thus including the h 2 h 1 h 1 vertex. The amplitude at leading order was given in the early papers [121,122] with the result expressed in terms of Passarino-Veltman scalar integrals. This result has also been implemented into MadGraph [123] taking into account the presence of a heavier SM-like scalar 7 , which we use for calculating the corresponding cross sections for each point shown here. This takes as input the modified Higgs top Yukawa coupling, the Higgs trilinear coupling, the heavy scalar top coupling, the h 2 h 1 h 1 coupling and the mass as well as the decay width of h 2 . Since h 2 decays into SM particles with reduced coupling (− sin θ) as compared with the SM Higgs and also decays to a pair of h 1 , the total width is simply given by: where Γ SM (h 2 → X SM ) denotes an exact SM Higgs-like h 2 decaying into the SM particles. For the di-Higgs production, if the resonant production of h 1 h 1 via the h 2 resonance dominates the cross section, then the cross section can be written in the narrow width approximation as In reality, interference effects between the resonant and non-resonant diagrams may be important and lead to constructive or destructive effect on the final full cross section [66]. We thus compare, for each scanned point, the obtained cross section for both the full calculation and the above approximation from the purely resonant production. This is shown in the left and middle plots of Fig. 8 Figure 9. The upper limits on di-Higgs resonant production cross section from ATLAS and CMS combined searches, shown as solid green and brown lines for ATLAS and CMS, respectively. The dashed lines denote the corresponding future projections for 3ab −1 of data at the HL-LHC (13TeV). As in the other plots, we distinguish those points which give SNR > 50 (red) and those of 50 > SNR > 10 (green).
and drops sharply as m h 2 is increased (left panel). Since, as we have seen in previous sections, the points with large SNR are concentrated around the region with larger m h 2 , most of the points with detectable GWs turn out to give small di-Higgs production and even negligible resonant production. The colors in the left panels make it clear that most of the points which have larger m h 2 (and larger SNR) tend to give very small di-Higgs production, with a cross section of O(10)fb, while smaller m h 2 gives O(100)fb. Moreover, there is a sharp drop of the resonant production cross section. From the middle panel, we can see that the color of decreasing branching ratio h 2 → h 1 h 1 coincides partly with increasing m h 2 for the very large m h 2 points. The small branching ratio is found for a majority of points and is due to the smallness of λ 211 . This can be seen from the right panel, where this correlation is shown with the color denoting m h 2 . It is found that a majority of points which have large m h 2 give small branching ratio. This can partly explain the cause of the drop of the resonant production. On the experimental side, both the ATLAS and CMS collaborations have recently published their search results for non-resonant and resonant di-Higgs productions using the data collected in 2016 at 13 TeV, with nearly the same integrated luminosity. The CMS search result is based on the 35.9fb −1 data, in the di-Higgs decay channels bbγγ [124], bbτ + τ − [125], bbbb [126][127][128][129] and bbW W/ZZ [130], with a recent combination given in [131]. ATLAS used 36.1fb −1 data and searched in channels γγbb [132], bbτ + τ − [133], bbbb [134], W W ( * ) W W ( * ) [135] and bbW W * [136], with also a combination of the first three channels [137]. We use the ATLAS and CMS combined limits in the resonant production channels and show them with green and brown solid lines respectively in Fig. 9. For the points giving detectable GWs, we calculate the resonant cross sections from gluon fusion at NNLO+NNLL using the available result in Ref. [138]. We can see that none of the points with detectable GW gives cross section above this limit. With the anticipation of HL-LHC at a luminosity of 3ab −1 (13TeV), we can get the future projections of this limit by a simple rescaling and obtain the two dashed lines. For this projection, the region with lower m h 2 550GeV can be partly explored by CMS and a little bit higher for ATLAS, while the high mass region remains out of reach for di-Higgs searches. Yet, Some points of the scanned parameters space with observable SNR show a promising di-Higgs production cross section of 50 fb or more at the LHC which, in principle, can be probed with 3 ab −1 . Therefore GW measurements can complement collider searches by revealing the high m h 2 region of the xSM model.

Higgs Cubic and Quartic Couplings
Future precise measurements of the Higgs cubic and quartic self-couplings can be used to reconstruct the Higgs potential to confirm ultimately the mechanism of EW symmetry breaking 8 and shed light on the nature of the EWPT. The measurements of above double Higgs production can be used to determine the cubic coupling and there have been extensive studies on this topic [56]. The best sensitivities obtained for these future colliders is typically at O(1). Despite the more formidable challenges with the quartic coupling measurement, there is now growing interest in it. Several different methods have been proposed and studied: through triple Higgs production measurement [57], through double Higgs production at hadron colliders where the quartic coupling enters gg → hh at two-loop [59] or renormalizes the cubic coupling, and at lepton colliders(via Z-associated production e + e − → Zhh and VBF production e + e − → ννhh), where the quartic coupling is involved in the V V hh coupling at one loop [60]. For example, Ref. [60] found a precision of measurement of ∼ ±25 for (500GeV, 4ab −1 + 1 TeV, 2.5ab −1 ) and ∼ ±20 for (500GeV, 4ab −1 + 1 TeV, 8ab −1 ) at 1σC.L., when the cubic coupling is marginalized in their χ 2 analysis.
In the xSM, both the Higgs cubic and quartic couplings are modified compared with their SM counterparts: In the absence of mixing of the scalars(θ = 0), these couplings reduce to the corresponding SM val- When θ = 0, we parametrize the deviations of these couplings from the SM values as: and show in Fig. 10 Figure 10. The Higgs cubic and quartic couplings (∆κ 3 , ∆κ 4 ) for parameter space points giving detectable GW. Here the green points give SNR > 10 and the red gives SNR > 50. The bars denote the sensitivity of ∆κ 3 from a global analysis of future colliders in Ref. [56], for various detector scenarios shown on the right side of the figures. The brown solid and blue dashed lines are the 1σ contours for two different ILC scenarios taken from Ref. [60]. The bottom panel is a zoomed-in version of the top one. δκ 3 ∈ (0, 1) and δκ 4 ∈ (0, 4). (3) a correlation exists δκ 4 ≡ ηδκ 3 , with η ≈ 2.8 for δκ 3 0.4 and most points fall within η ∈ (2, 4). To understand these, we note, since phenomenological constraints requires a small θ, we expect the second feature to follow naturally. The other features can be understood by Taylor expanding the couplings for small θ and we find: In the above square brackets, the terms proportional to m 2 h 2 /m 2 h 1 dominate for the majority of the points since v s is concentrated at small values; b 3 is at most ∼ 10v EW , b 4 5 from the scan and m h 2 500GeV generally holds. Then the above approximations show positive δκ 3 and δκ 4 and give δκ 4 /δκ 3 ≈ 2.5, which is fairly close to η = 2.8. For relatively large θ, high order corrections need to be taken into account and above linear correlation would be changed.
To compare with the direct measurements of these couplings at future e + e − colliders and the HL-LHC, we added in Fig. 10 the precisions of these measurements from studies in the literature. The two elliptical 68%CL closed contours are taken from Ref. [60] which focuses on the quartic coupling, for two possible scenarios of the ILC. The bars are the precisions that can be reached from various considerations of future colliders, labelled on the right of the figure, taken from Ref. [56](for other studies, see e.g. [59,139,[139][140][141][142]). Here the inner and outer bar regions denote the 68%CL and 95%CL results. We can see, it is generically very hard for colliders to probe the cubic coupling at a precision that can reveal the points giving detectable GWs with high confidence level(say 95%) 9 . The most precise comes from the ILC when all possible runs at different luminosities are combined and with the data of HL-ILC included, which gives 0.4 ∼ 0.5 uncertainty on the measurement of δκ 3 at 95%CL. While the analysis in Ref. [56] does not include the quartic coupling, the contours from Ref. [60] do give a hint on its measurement and show that it is infeasible for the colliders to probe the parameter space giving detectable GWs. For the trilinear and quartic coupling deviations that we found, the impact on the triple Higgs cross section is mild for hadron colliders even for a future pp collider at 100 TeV [57,58], however, resonant contributions in xSM might enhance the cross section up to a factor of O(10) [143].
Therefore we expect future GW measurements can make a valuable complementary role in determining the Higgs self-couplings, especially the quartic coupling. While we do not have a statistical analysis here, Fig. 10 does tell us that δκ 4 is equally important as δκ 3 on GW signal generation since η is at most 4. Thus we expect a full statistical analysis would yield roughly the same precision on the determination of δκ 3 and δκ 4 , which is well improved compared with the situation at colliders.

Diboson Resonance Search Limits at Colliders
The W W and ZZ branching ratios become sizeable in parts of the parameter space where the trilinear coupling λ 211 is relatively small, as one can see from the rightmost panel of Fig. 8. In 9 It should be noted that both studies used some versions of the effective field theory approach to quantify the modification of the SM couplings due to possible new physics effects. Therefore the precisions overlaid in Fig. 10 might not be what the colliders can achieve if the xSM model was used in their studies. However we expect the two contours, taken from Ref. [60], to be largely unaffected since the heavier scalar contribution in their framework is suppressed by extra powers of s θ . We also expect that the bar regions, taken from Ref. [56], would get tighter since the set of parameters used in their study are highly correlated here.  Fig. 11, we show the branching ratios of the h 2 → W W, ZZ and h 2 → h 1 h 1 channels. We see that the W W, ZZ channels can be as big as 90% for a large range of h 2 masses which could show up at searches for weak diboson resonances. Combined, W W, ZZ and h 1 h 1 correspond to nearly all the decays of h 2 , which make them the best search channels for h 2 resonances at colliders. Besides the di-Higgs production measurements, which can be used to extract the Higgs cubic and quartic couplings, there also exist generic scalar resonance searches at the LHC. In particular, ATLAS and CMS have performed extensive analyses in the searches for a heavier SM-like scalar resonance in V V and V H decay channels of the heavy scalar (V = W/Z). ATLAS gives a recent combination of all previous analyses in bosonic and leptonic final states at √ s = 13TeV with 36fb −1 data collected in 2015 and 2016 [144]. The limits are drawn for h 2 production cross section in gluon fusion and vector boson fusion production channels. These two limits are shown in the left and right panels, respectively, in Fig. 12 with green solid lines, together with the detectable GW points. For cross section calculations, we use the set of result calculated to NNLO precision for VBF and for gluon fusion, we use NNLO+NNLL, as also used before in Fig. 9.
It is evident that the current limits from diboson searches are rather loose as most points fall under this line, with gluon fusion limit being able to touch a fraction of the lighter h 2 point. For the HL-LHC with ∼ 3ab −1 , we obtain estimates of future projections by a simple scaling factor and obtain the dashed lines for ∼ 3ab −1 at 13TeV (while HL-LHC would probably run at 14TeV). We can see in all cases that the HL-LHC will probe a larger fraction of the parameter space for both ggH and VBF channels. For ggH, this region covers a range from low to high masses. For VBF, it can cover a region of relatively heavy h 2 . Both channels are sensitive to h 1 h 1 cross section times branching ratio down to ∼ 1 fb in some favorable points of the parameters space. The points that can be probed by HL-LHC serve as promising targets for both colliders and GW detectors but a majority of the parameter space will probably be left to GW detectors.

Summary
In this paper, we embarked on a study of the singlet-extended SM Higgs sector. A detailed scan of the parameter space of this model was performed, incorporating all relevant phenomenological constraints, and regions with large SNR at LISA were identified. Subtle issues pertaining to the bubble wall velocity were discussed, and a range of velocity profiles described. Our main findings are the following. For the parameter space that satisfies all phenomenological constraints, gives successful EWPT and generates GWs, 99% leads to a one-step EWPT with the remaining to two-step EWPT and 22% generates detectable GWs(SNR > 10) at LISA. The main features of the parameter space that gives detectable GWs is: 20GeV |v s | 50GeV, where v s is the vev of the singlet field; it is more concentrated in the large m h 2 region, where m h 2 is the mass of the heavier scalar h 2 ; θ 0.2 for the majority of the space. Di-Higgs searches at both ATLAS and CMS are currently unable to probe this parameter space, but HL-LHC will be able to probe the lighter h 2 region while the heavier h 2 region will remain elusive. Weak diboson resonance searches cannot constrain xSM much either but the HL-LHC will be able to probe a large fraction of its parameters space in this channel. The Higgs cubic and quartic couplings are at O(1) deviations from the SM values and obey a relation δκ 4 ≈ (2 − 4)δκ 3 , where δκ 4 and δκ 3 are the relative deviations of the quartic and cubic couplings from their SM counterparts respectively.
Our results broadly indicate that high energy colliders and GW detectors are going to play complementary roles in probing the parameter space of scalar sectors. Several future directions can be contemplated. It would be interesting to understand how this complementarity plays out in two Higgs doublet models, as well as other scalar sector extensions classified in [19]. It would also be interesting to investigate the complementarity of GW and collider probes for phase transitions in the dark sector. We leave these questions for future study.
So the scalar couplings as well as their masses and mixing angles wont be affected by this translation. For easy comparison between these two representations, we show here the transformation rules between these two bases. Given potential parameters in the non-tadpole basis in Eq. 2.1, the parameters in the basis where b 1 = 0(denoted with a prime) can be obtained: a 2 v s ), while a 2 , λ, b 4 remains unchanged. On the other hand, given parameters in the tadpole basis where v s = 0 and b 1 = 0, the parameter set in the basis used in this work can be found: where x is to be solved from the cubic equation which might give more than one solutions. In the basis v s = 0, the degree of freedom carried by v s in the basis v s = 0 is transformed to a different parameter. For example, one can choose it to be a 2 and then the full set of independent parameters can be chosen as We note further there are also studies of this model where a Z 2 symmetry in the S fields are imposed and are spontaneously broken [49,64,66]. This specific model correspond to a special limit of the potential here.