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

TL;DR

This study explores how gravitational wave detectors and collider experiments can jointly investigate the singlet scalar extension of the Standard Model, focusing on phase transition patterns, gravitational wave signals, and collider signatures.

Contribution

It provides a comprehensive analysis of the complementarity between gravitational wave signals and collider measurements in probing the singlet scalar extension of the Standard Model.

Findings

Identification of parameter regions with detectable gravitational waves.

Correlation between collider signals and gravitational wave signals.

Insights into phase transition dynamics and bubble wall velocities.

Abstract

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.

Figures

Equations

$$V(H,S)$$

$${\mu }^2=\lambda v_{\text{EW}}^2+\frac{1}{2}{v}_s(a_1+a_2{v}_s),$$

$$b_2=-\frac{1}{4v_s}[v_{\text{EW}}^2(a_1+2a_2{v}_s)+4v_s^2(b_3+b_4{v}_s)],$$

$$\lambda =\frac{m_{h_1}^2{c}_{\theta }^2+m_{h_2}^2{s}_{\theta }^2}{2v_{\text{EW}}^2},$$

$$a_1=\frac{2v_s}{v_{\text{EW}}^2}[2v_s^2(2b_4+{\tilde{b}}_3)-m_{h_1}^2-m_{h_2}^2+c_{2\theta }(m_{h_1}^2-m_{h_2}^2)],$$

$$a_2=\frac{-1}{2v_{\text{EW}}^2{v}_s}[-2v_s(m_{h_1}^2+m_{h_2}^2-4b_4{v}_s^2)$$

$$+(m_{h_1}^2-m_{h_2}^2)(2c_{2\theta }{v}_s-v_{\text{EW}}{s}_{2\theta })+4{\tilde{b}}_3{v}_s^3],$$

$$h_1=c_{\theta }h+s_{\theta }s,h_2=-s_{\theta }h+c_{\theta }s,$$

$$\text{centering}{v}_s,m_{h_2},\theta ,b_3,b_4.\text{@add@centering}$$

$$\lambda >0,b_4>0,a_2⩾-2\sqrt{\lambda b_4}.$$

$$V(h,s,T)=-\frac{1}{2}[{\mu }^2-{\Pi }_h(T)]h^2-\frac{1}{2}[-b_2-{\Pi }_s(T)]s^2$$

$$+\frac{1}{4}\lambda h^4+\frac{1}{4}{a}_1{h}^{2}s+\frac{1}{4}{a}_2{h}^2{s}^2+\frac{b_3}{3}{s}^3+\frac{b_4}{4}{s}^4,$$

$${\Pi }_h(T)=\left(\frac{2m_W^2+m_Z^2+2m_t^2}{4v^2}+\frac{\lambda }{2}+\frac{a_2}{24}\right)T^2,$$

$${\Pi }_s(T)=\left(\frac{a_2}{6}+\frac{b_4}{4}\right)T^2,$$

$$\text{centering}{T}_c,T_n,\alpha ,\beta ,v_w.\text{@add@centering}$$

$$\Gamma \sim \mathcal{A}(T)e^{-S_3/T},$$

$$S_3(\vec{\varphi },T)=4\pi \int r^{2}dr\left[\frac{1}{2}{\left(\frac{d\vec{\varphi }(r)}{dr}\right)}^2+V(\vec{\varphi },T)\right],$$

$$\displaystyle\frac{d\vec{\phi}(r)}{dr}\Big{|}_{r=0}=0,\quad\quad\vec{\phi}(r=\infty)=\vec{\phi}_{\text{out}}.$$

$$\displaystyle\int_{0}^{t_{n}}{\Gamma V_{H}(t)}dt=\int_{T_{n}}^{\infty}{dT\over T}\left(2\zeta M_{\rm Pl}\over T\right)^{4}e^{{-S_{3}/T}}={\cal O}(1),\$$

$$\beta \equiv H_n{T}_n{\frac{d(S_3/T)}{dT}|}_{T_n},$$

$$\displaystyle\alpha=\frac{\Delta\rho}{\rho_{R}}=\frac{1}{\rho_{R}}\left[-V(\vec{\phi}_{b},T)+T\frac{\partial V(\vec{\phi}_{b},T)}{\partial T}\right]\Bigg{|}_{T=T_{n}},\$$

$$\frac{v_{H}(T)}{T}\Big{|}_{T=T_{n}}\gtrsim 1.$$

$$2\frac{v}{\xi }=\frac{1-v\xi }{1-v^2}\left[\frac{{\mu }^2}{c_s^2}-1\right]{\partial }_{\xi }v,$$

$${\kappa }_v=\frac{3}{\Delta \rho v_w^3}\int \omega (\xi )\frac{v^2}{1-v^2}{\xi }^{2}d\xi ,$$

$$T_{\ast }=T_n(1+{\kappa }_T\alpha {)}^{1/4}.$$

$$n_B=\frac{3{\Gamma }_{\mathrm{ws}}}{D_q{\lambda }_{+}}{\int }_0^{-\infty }{n}_L(\bar{z})e^{-{\lambda }_{-}\bar{z}}d\bar{z},$$

$$\frac{n_B}{s}{|}_{T=T_{\ast }}={\xi }_D\frac{n_B}{s}{|}_{T=T_n},$$

$${\Omega }_{\text{GW}}{h}^2\simeq {\Omega }_{\text{col}}{h}^2+{\Omega }_{\text{sw}}{h}^2+{\Omega }_{\text{turb}}{h}^2.$$

$${\Omega }_{\text{sw}}{h}^2=2.65\times 10^{-6}\left(\frac{H_{\ast }}{\beta }\right){\left(\frac{{\kappa }_v\alpha }{1+\alpha }\right)}^2{\left(\frac{100}{g_{\ast }}\right)}^{1/3}$$

$$\times v_w{\left(\frac{f}{f_{\text{sw}}}\right)}^3{\left(\frac{7}{4+3(f/f_{\text{sw}}{)}^2}\right)}^{7/2}.$$

$$f_{\text{sw}}=1.9\times 10^{-5}\frac{1}{v_w}\left(\frac{\beta }{H_{\ast }}\right)\left(\frac{T_{\ast }}{100\text{GeV}}\right){\left(\frac{g_{\ast }}{100}\right)}^{1/6}\text{Hz},$$

$${\Omega }_{\text{turb}}{h}^2=3.35\times 10^{-4}\left(\frac{H_{\ast }}{\beta }\right){\left(\frac{{\kappa }_{\text{turb}}\alpha }{1+\alpha }\right)}^{3/2}{\left(\frac{100}{g_{\ast }}\right)}^{1/3}$$

$$\times v_w\frac{(f/f_{\text{turb}}{)}^3}{[1+(f/f_{\text{turb}}){]}^{11/3}(1+8\pi f/h_{\ast })},$$

$$f_{\text{turb}}=2.7\times 10^{-5}\frac{1}{v_w}\left(\frac{\beta }{H_{\ast }}\right)\left(\frac{T_{\ast }}{100\text{GeV}}\right){\left(\frac{g_{\ast }}{100}\right)}^{1/6}\text{Hz}.$$

Full text

1 Introduction

Since the first direct detection of gravitational waves (GWs) by the LIGO and Virgo collaborations Abbott:2016blz , a new interface has arrived in particle physics – its intersection with GW astronomy. While ground based GW detectors have their best sensitivity at frequencies $\sim\mathcal{O}(100)\,\text{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) Audley:2017drz , the Big Bang Observer (BBO), the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) Yagi:2011wg , Taiji Gong:2014mca and Tianqin Luo:2015ght . The physical sources of GWs in this frequency band include supermassive black hole binaries Klein:2015hvg , extreme mass ratio inspirals Babak:2017tow and the stochastic background of primordial GWs produced during first order cosmological phase transitions Caprini:2015zlo .

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 details within such a setting. The natural choice is the electroweak phase transition (EWPT)  Grojean:2006bp with the simplest extension of the Higgs sector: the singlet scalar augmented Standard Model or the xSM111Hidden sector phase transitions are also being actively investigated Schwaller:2015tja ; Alanne:2014bra ; Breitbach:2018ddu ; Aoki:2017aws ; Bai:2018dxf , and exploring complementarity in such settings is an interesting future direction. We refer to Ref. Caprini:2015zlo ; Cai:2017cbj ; Weir:2017wfa ; Caprini:2018mtu ; Mazumdar:2018dfl for recent work on these topics.. 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 Chung:2012vg (see Ref. Beniwal:2018ygo ; Croon:2018kqn ; Hashino:2018wee ; Marzo:2018nov ; Beniwal:2017eik ; Hashino:2018zsi ; Addazi:2017gpt ; Addazi:2017nmg ; Chiang:2017zbz ; Wan:2018udw ; Chen:2017cyc ; Vieu:2018nfq ; Basler:2017uxn ; Jinno:2017ixd ; Chao:2017ilw ; Huang:2014ifa ; Tsumura:2017knk ; Bian:2017wfv ; Hektor:2018esx ; Huang:2017rzf ; Ghorbani:2017lyk ; Kang:2017mkl ; Addazi:2017oge ; Chen:2017qcz ; Marzola:2017jzl ; Kobakhidze:2017mru ; Zhou:2018zli ; Beniwal:2018hyi ; Dev:2016feu ; Balazs:2016tbi ; Ahriche:2018rao ; Shajiee:2018jdq ; Bian:2018bxr ; Blinov:2015sna ; Inoue:2015pza ; Vaskonen:2016yiu ; No:2018fev ; Chiang:2018gsn ; Chala:2018opy ; Cai:2017tmh ; Alanne:2018brf ; Kannike:2019wsn ; Ashoorioon:2009nf ; Cheng:2018ajh ; Bi:2015qva for related studies on EWPT and GW). It has been extensively investigated in phenomenological studies Profumo:2007wc ; Profumo:2014opa ; Huang:2017jws ; Robens:2015gla , studies of EWPT Profumo:2014opa ; Profumo:2007wc ; Kotwal:2016tex ; Espinosa:2011ax ; Kozaczuk:2015owa and di-Higgs analyses Lewis:2017dme guided by the requirements of EWPT Huang:2017jws , 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, $\alpha$, 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., Azatov:2015oxa ; Alves:2017ued ; DiVita:2017vrr ; Plehn:2005nk ; Binoth:2006ym ; Bizon:2018syu ; Liu:2018peg ). 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.

2 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 Profumo:2007wc ; Profumo:2014opa ; Huang:2017jws :

[TABLE]

where $H^{\text{T}}=(G^{+},(v_{\text{EW}}+h+iG^{0})/\sqrt{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 $\mu$ and $b_{2}$ can be solved from the two minimization conditions around the EW vacuum($\equiv(v_{\text{EW}},v_{s})$),

[TABLE]

and $\lambda,a_{1},a_{2}$ can be replaced by physical parameters $\theta$, $m_{h_{1}}$ and $m_{h_{2}}$ from the mass matrix diagonalization 222Here $s_{\theta}\equiv\sin\theta$ and $c_{\theta}\equiv\cos\theta$.:

[TABLE]

where $\tilde{b}_{3}\equiv b_{3}/v_{s}$ and we have defined the physical fields $h_{1}$ and $h_{2}$ as

[TABLE]

with a mixing angle $\theta$. 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_{\theta}$ while the coupling of $h_{2}$ with SM particles is $(-s_{\theta})$ 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:

[TABLE]

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\rightarrow S^{\prime}=S-v_{s}$, such that $\langle S\rangle=0$, the resulting potential will take the same form as Eq. 1 but with the addition of a non-zero tadpole term $b_{1}S$ Lewis:2017dme . 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\rightarrow-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., Inoue:2014nva ; Chen:2017com ; Bian:2017wfv or the included CP-violation may be large but has little effect on EWPT Guo:2016ixx ).

3 Phenomenological Constraints

In this section, we briefly discuss the phenomenological constraints used in our analysis, following the standard treatments given in Refs. Pruna:2013bma ; Robens:2015gla ; Lewis:2017dme . 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 Lewis:2017dme  333 Note that these tree level relations will change when loop corrections are taken into account. However due to the way of calculating the effective potential in Eq.(4.1), these relations suffice to guarantee that the potential is bounded from below when $T\rightarrow 0$ in Eq. 7. ,

[TABLE]

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 $\partial V/\partial\phi_{i}=0$($\phi_{1}\equiv h$, $\phi_{2}\equiv s$) and subsequently calculating eigenvalues of the Hessian matrix $\{\partial^{2}V/\partial\phi_{i}\partial\phi_{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\theta$, therefore the Higgs signal strength is given by $\mu_{H}=\cos^{2}\theta$. Experimentally, the most recent ATLAS and CMS combined fit of this value is $\mu_{H}=1.09^{+0.11}_{-0.10}$ Khachatryan:2016vau and a $\chi^{2}$ analysis shows that $|\sin\theta|>0.33$ are excluded at $95\%$ CL Carena:2018vpt .

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\rightarrow 2$ processes satisfies $|\text{Re}\,a_{l}(s)|\lesssim 1/2$ for $\sqrt{s}\rightarrow\infty$. We consider all channels of scalar/vector boson $2\rightarrow 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 Lopez-Val:2014jva and the oblique EW corrections Peskin:1991sw ; Hagiwara:1994pw , 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 $\alpha(0)$ Lopez-Val:2014jva . 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_{W}^{\text{exp}}=80.385\pm 0.015\text{GeV}$ Alcaraz:2006mx ; Aaltonen:2012bp ; D0:2013jba 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 $(\theta,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 Lopez-Val:2014jva ; Robens:2015gla .

To give the reader a flavor of the above phenomenological constraints, we fix $m_{h_{2}}=300\,\text{GeV}$, $\theta=0.2$, $b_{4}=4$ and show the various bounds on the remaining two parameters $(v_{s}/v_{\text{EW}},b_{3}/v_{\text{EW}})$ in Fig. 1. This choice of $m_{h_{2}}$ and $\theta$ 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 perturbative 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.

4 EWPT and Gravitational Waves

4.1 Effective Potential

EWPT is an essential step 444 Other mechanisms generally do not need EWPT to generate the baryon asymmetry. For example, in leptogenesis, the out-of-equilibrium requirement is provided by the expanison of the universe and the lepton asymmetry is converted to the baryon asymmetry through the weak Sphaleron process. in generating the observed baryon asymmetry in the universe by providing an out-of-equilibrium environment, one of the three Sakharov conditions Sakharov:1967dj , in the framework of electroweak baryogenesis (see Morrissey:2012db 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 Coleman:1973jx and its finite temperature counterpart Quiros:1999jp as well as the daisy resummation Parwani:1991gq ; Gross:1980br . 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 Patel:2011th 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 Profumo:2007wc ; Profumo:2014opa ; Kotwal:2016tex ; Huang:2017jws ; Alves:2018oct . This effective potential is then given by 555 We also note that we have neglected a tadpole term proportional to $T^{2}s$, which originates from the $a_{1}$ and $b_{3}$ terms in the potential in Eq. 1, since it comes with a factor $v_{s}/v_{\text{EW}}$ and is suppressed for most of the parameter space giving detectable GWs, to be presented in later sections. Indeed its effect has been found to be numerically negligible from previous studies Profumo:2007wc ; Profumo:2014opa .

[TABLE]

where $\Pi_{h}(T)$ and $\Pi_{s}$ are the thermal masses of the fields,

[TABLE]

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:

[TABLE]

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 Turner:1992tz

[TABLE]

where $S_{3}$ is the 3-dimensional Euclidean action of the critical bubble, which minimizes the action

[TABLE]

and satisfies the bounce boundary conditions

[TABLE]

Here $\vec{\phi}_{\text{out}}$ denotes the two components vev of the fields outside the bubble, which is not necessarily the origin for two-step EWPT. The prefactor $\mathcal{A}(T)\propto T^{4}$ on dimensional grounds. Its precise determination needs integrating out fluctuations around the above static bounce solution (see e.g., Dunne:2005rt ; Andreassen:2016cvx for detailed calculations or Weinberg:1996kr 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 ${\cal O}(1)$ at a certain temperature Chao:2017vrq :

[TABLE]

where $V_{H}(t)$ is the Horizon volume, $M_{\rm Pl}$ is the Planck mass and $\zeta\sim 3\times 10^{-2}$. From this equation, it follows that $S_{3}(T)/T\approx 140$ Apreda:2001us 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 $\beta$ Apreda:2001us :

[TABLE]

where $H_{n}$ is the Hubble rate at $T_{n}$.

Next, $\alpha$ is the vacuum energy released from the EWPT normalized by the total radiation energy density ($\equiv\rho_{R}$) at $T_{n}$ Kamionkowski:1993fg :

[TABLE]

where $\rho_{R}=g_{\ast}\pi^{2}T_{n}^{4}/30$ with $g_{\ast}\approx 100$ and $\vec{\phi}_{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 Cline:2006ts ; Morrissey:2012db :

[TABLE]

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}$(see e.g., Ref. Grojean:2004xa ; Grojean:2006bp ; Chala:2018ari ), 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. 13 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 Patel:2012pi ; Ramsey-Musolf:2017tgh ; Chao:2017vrq . One can immediately imagine mainly the following EWPT types:

• (A):

$(0,0)\rightarrow(v_{H}\neq 0,v_{S}\neq 0)$

• (B):

$(0,0)\rightarrow(v_{H}=0,v_{S}\neq 0)\rightarrow(v_{H}\neq 0,v_{S}\neq 0)$

• (C):

$(0,0)\rightarrow(v_{H}\neq 0,v_{S}=0)\rightarrow(v_{H}\neq 0,v_{S}\neq 0)$

where the last vacuum configuration $(v_{H}\neq 0,v_{S}\neq 0)$ in each case would eventually evolve to the EW vacuum at $T=0$ 666More exotic patterns might appear but should be of negligible parameter space. For an example, see Ref. Angelescu:2018dkk .. 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. Espinosa:2011ax ; Chao:2017vrq , 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.

4.2 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}$ generally leads to more energy being released to the kinetic energy of the plasma and therefore a stronger 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 No:2011fi . 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. Espinosa:2010hh 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 ($\equiv 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 $(\alpha,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\equiv\xi$. The differential equation governing the velocity profile is derived from the conservation of the energy momentum tensor describing the fluid and scalar field Espinosa:2010hh :

[TABLE]

where $c_{s}=1/\sqrt{3}$ is the speed of sound in the plasma and $\mu(\xi,v)=(\xi-v)/(1-\xi v)$ is a Lorentz boost transformation. Far outside the bubble and deep inside the bubble, the plasma will not be stirred, that is $v\rightarrow 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(\xi)$ are shown in Fig. 2 for $\alpha=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(\xi)$ appears outside the bubble and $v(\xi)$ jumps to zero. This is the location of the shock front, and beyond this point the solution of Eq. 17 is invalid and a shock front develops such that $v(\xi)$ goes to zero consistently. When $v_{w}$ surpasses $c_{s}$ but is less than a certain threshold $\xi_{J}(\alpha)$, a supersonic deflagration mode KurkiSuonio:1995pp 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 $\xi_{J}(\alpha)$, as a function of $\alpha$, corresponds to the Jouguet detonation Steinhardt:1981ct , 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}=\xi_{J}(\alpha)$ 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 $\alpha$, the corresponding $v_{+}$ to have subsonic value, taken to be $0.05$ here, a choice usually used in EWBG calculations John:2000zq ; Cirigliano:2006dg ; Chung:2009qs ; Chao:2014dpa ; Guo:2016ixx ). The procedure of achieving the above goal is as follows: for each given $\alpha$ 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 777 For two-step EWPT, a small $v_{+}$ is not necessarily required for both steps of EWPT. However since $v_{w}$ is otherwise an almost free parameter, we stick to the choice $v_{+}=0.05$ for both steps. .

With $v(\xi)$ obtained, one can also calculate the bulk kinetic energy normalized by the vacuum energy released during the EWPT Espinosa:2010hh :

[TABLE]

where $\omega(\xi)$ is the enthalpy density, varying as function of $\xi$, and can be solved once $v(\xi)$ is found. The remaining part $1-\kappa_{v}\equiv\kappa_{T}$ gives the fraction of the vacuum energy going to heat the plasma. Therefore a reheating temperature can be defined as

[TABLE]

This leads to an increase in entropy density and thus a dilution of the generated baryon asymmetry Patel:2012pi . Typically in EWBG calculations, the wall curvature is neglected and the transport equations depend on a single coordinate $\bar{z}$ in the bubble wall rest frame, where $\bar{z}>0$ ($<0$) corresponds to broken (unbroken) phase. The solved baryon asymmetry density $n_{B}$ is a constant inside the bubbles(see, e.g., Lee:2004we ):

[TABLE]

where $s(T)=2g_{\ast}\pi^{2}T^{3}/45$ is the entropy density, $\Gamma_{\text{ws}}\approx 120\alpha_{w}^{5}T$ is the weak Sphaleron rate in the EW symmetric phase Bodeker:1999gx , $\lambda_{\pm}=(v_{+}\pm\sqrt{v_{+}^{2}+15\Gamma_{\text{ws}}D_{q}})/(2D_{q})$ with $D_{q}$ the diffusion constant for quarks Bodeker:1999gx 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_{\ast}$ rather than $T_{n}$ or $T_{c}$, which are conventionally used. The diluted baryon asymmetry is then given by

[TABLE]

where $\xi_{D}\equiv(1+\kappa_{T}\alpha)^{-3/4}$ captures the dilution effect of the generated baryon asymmetry by reheating of the plasma. We then need to make sure that $\xi_{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 Engel:2013lsa ; Chupp:2017rkp .

4.3 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 Kamionkowski:1993fg . 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 two-point correlation function and is parametrized by the gravitational wave energy density spectrum $\Omega_{\text{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 milli-Hertz range Grojean:2006bp , 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 Kosowsky:1991ua ; Kosowsky:1992rz ; Kosowsky:1992vn ; Huber:2008hg ; Jinno:2016vai ; Jinno:2017fby , sound waves in the plasma Hindmarsh:2013xza ; Hindmarsh:2015qta and magneto-hydrodynamic turbulence (MHD) Hindmarsh:2013xza ; Hindmarsh:2015qta . The total energy density spectrum can be obtained approximately by adding these contributions:

[TABLE]

Recent studies suggest that the energy deposited in the bubble walls is negligible, despite the possibility that the bubble walls can run away in some circumstances Bodeker:2009qy . Therefore while a bubble wall can reach relativistic speed, its contribution to gravitational waves can generally be neglected Bodeker:2017cim . 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 Hindmarsh:2015qta :

[TABLE]

Here $H_{\ast}$ is the Hubble parameter at $T_{\ast}$ when the phase transition has completed. It has a value close to that evaluated at the nucleation temperature $T_{n}$ for sufficiently short EWPT Caprini:2015zlo . We take $T_{\ast}$ to be the reheating temperature, defined earlier in Eq. 19. Moreover, $f_{\text{sw}}$ is the present peak frequency which is the redshifted value of the peak frequency at the time of EWPT ($=2\beta/(\sqrt{3}v_{w})$):

[TABLE]

where $\kappa_{v}$ is defined in Eq. 18 and can be calculated as a function of ($\alpha$, $v_{w}$) by solving the velocity profiles described in Sec. 4 Espinosa:2010hh . It should be noted that a more recent numerical simulation by the same group Hindmarsh:2017gnf ; Hindmarsh:2016lnk shows a slightly enhanced $\Omega_{\text{sw}}h^{2}$ and reduced peak frequency $f_{\text{sw}}$. We also note that the results from these simulations are currently limited to regions of small $v_{w}$ and $\alpha$ and therefore their validity for ultra-relativistic $v_{w}$ and large $\alpha$ (say $\alpha\gtrsim 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 Caprini:2009yp ; Binetruy:2012ze ,

[TABLE]

where $f_{\text{turb}}$ is the peak frequency and is given by,

[TABLE]

Here the factor $\kappa_{\text{turb}}$ describes the fraction of energy transferred to the MHD turbulence and is given roughly by $\kappa_{\text{turb}}\approx\epsilon\kappa_{v}$ with $\epsilon\approx 5\sim 10\%$ Hindmarsh:2015qta . We take $\epsilon=0.1$ in this study.

In both Eq. 23 and  25, 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. Caprini:2015zlo :

[TABLE]

Here $h^{2}\Omega_{\text{exp}}(f)$ is the experimental sensitivity and corresponds to the lower boundaries of the color-shaded regions in Fig. 3 for the shown detectors 888 There are possible astrophysical foregrounds coming from, e.g., the superposition of unresolved (i.e., low SNR) gravitational wave signals of the white dwarf binaries in our Galaxy Klein:2015hvg . Including these will slightly reduce the SNR calculated here. . $\mathcal{T}$ is the mission duration in years for each experiment, assumed to be $5$ here. The factor $\delta$ 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 Thrane:2013oya . In our numerical analysis, we stick to the most mature LISA detector with the C1 configuration, defined in Ref. Caprini:2015zlo . 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 ($\text{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 Caprini:2015zlo .

As an example, we scan over the EW vacuum stability regions in the plane $(v_{s}/v_{\text{EW}},b_{3}/v_{\text{EW}})$ of 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 ($\text{SNR}<10$), green ($50>\text{SNR}>10$) and red ($\text{SNR}>50$). Here most of the points give type (A) EWPT with only several points for type (B) or (C), denoted by diamond shapes.

5 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 $3\text{ab}^{-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?

We find that it is more efficient to cover the parameter space in the scan using the tadpole basis parameters. So we start the random number generating in the tadpole basis. Once a parameter space point is obtained from the sampling, it is converted to the non-tadpole basis parameters and for all subsequent phenomenological checks and calculations. The ranges of the tadpole basis parameters where we do the random sampling are the following:

[TABLE]

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 a modified version of CosmoTransitions Wainwright:2011kj for calculating the thermal history and the parameters relevant for EWPT 999 Other packages include Bubbleprofiler Athron:2019nbd and AnyBubble Masoumi:2016wot . It should also be noted that there generically exists a difficulty for solving bounce solutions in very thin-walled cases, the discussion of which can be found in above paper of Bubbleprofiler and also in Ref. Piscopo:2019txs where the neural networks is introduced to solve the bounce solutions. We have verified that the majority of the points we used are not very thin-walled. . 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.

5.1 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 $\alpha$ and $\beta/H_{n}$, where a larger $\alpha$ and a smaller $\beta/H_{n}$ gives stronger GW SNR, as shown in the left panel of Fig. 4, where the colors denote $\text{SNR}<10$ (blue), $50>\text{SNR}>10$ (green) and $\text{SNR}>50$ (red). We observe that the points which give detectable GWs lie in the bottom right region of the population.

Physically, $\alpha$ quantifies the amount of energy released during the EWPT and therefore a larger $\alpha$ gives stronger GW signals. In addition, for fixed $v_{w}$, a larger $\alpha$ leads to a larger fraction of energy transformed into the plasma kinetic energy, quantified by $\kappa_{v}$, and therefore a further gain in GW production. A further enhancement for larger $\alpha$ comes from the fact that since we fixed $v_{+}=0.05$, increasing $\alpha$ also increases $v_{w}$. It should be noted, even without an explicit calculation, that for each fixed value of $\alpha$, the allowed values of $v_{w}$ are limited to a certain range (see e.g., Fig. 1 in Ref. Alves:2018oct ). 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 $\xi_{J}(\alpha)$ 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 $\alpha\gtrsim 1$ and $v_{w}\sim 1$, the calculations of the GW spectra may become unreliable for the following reasons: ($i$) While the study of Ref. Bodeker:2017cim suggests that the energy stored in the scalar field kinetic energy is negligible, a very large $\alpha$ 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 $\alpha$ as well as $v_{w}$ and thus the use of these results for large $\alpha$ 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. Ellis:2018mja for a recent analysis). Despite these issues, we find $49\%$ of points with $\text{SNR}>10$ have $\alpha<1$ and removing the points with $\alpha>1$ does not change the main findings of our work.

We now turn to the parameter $\beta/H_{n}$, which roughly characterizes the inverse time duration of the EWPT. A smaller $\beta/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 Hindmarsh:2013xza ; Hindmarsh:2015qta , one enhancement comes from $1/(\beta/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 of the EWPT. It should be noted that $\beta/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 $\alpha$ 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 $\alpha$. The supercooling can be quantified by the fraction of the first term($\equiv\Delta\rho_{V}$) of Eq. 15 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 $\Delta\rho_{V}/\Delta\rho$ is less than $0.6$ for most of the parameter space. The remaining part comes from the second term of the definition of $\alpha$.

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 $\xi_{D}$. Here since $\kappa_{T}$ is a function of $v_{w}$ and $\alpha$ while $v_{w}$ is also a function of $\alpha$ when $v_{+}$ is fixed, we find $\xi_{D}$ is solely a function of $\alpha$. 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 $\xi_{D}\gtrsim 0.65$ and those with a smaller $\alpha$ have a dilution factor closer to $1$. In particular, the points with $\alpha\lesssim 1$ for which GW can be reliably calculated, the dilution effect is rather small as $\xi_{D}\gtrsim 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 Patel:2012pi used a different quantification of the dilution factor, with the definition:

[TABLE]

where $s$ is the entropy density at $T_{n}$ and $\Delta s$ is calculated from the second term in the definition of $\alpha$ in Eq. 15. To compare with the factor $\xi_{D}$, what we use here, we show values of this factor in the same plot of $\xi_{D}$ for every point that gives detectable GWs. It is evident from this figure that these two factors are roughly the same and both decrease linearly for $\alpha\lesssim 0.4$. For $\alpha\gtrsim 0.4$, $\xi_{D}^{(2)}$ gives an overestimation of the dilution effect while $\xi_{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 $\xi_{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 $\alpha\lesssim 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.

5.2 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 $\alpha$ and small $\beta/H_{n}$ leads to loud GW signals. Even though the relation between $(\alpha,\beta/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\,\text{GeV}\lesssim|v_{s}|\lesssim 50\,\text{GeV}$ for most points, with a peak distribution at around $20\,\text{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\theta)$ 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 $(\alpha,\beta/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 $\beta/H_{n}$, a larger $m_{h_{2}}$ gives a larger $\alpha$, thus larger SNR. This explains the concentration of the points in the $m_{h_{2}}$ direction in the middle figure. In the $\sin\theta$ direction, the value of $\theta$ 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 $\theta$ give rarer points. We also overlaid on this plot the various sensitivity projections from colliders in probing the value of $\theta$, which includes HL-LHC, ILC with two configurations (ILC-1: $250\text{GeV}$, $250\text{fb}^{-1}$, ILC-3: $1\text{TeV}$,$1\text{ab}^{-1}$) and future circular $e^{+}e^{-}$ colliders ($240\text{GeV},1\text{ab}^{-1}$), all taken from Ref. Profumo:2014opa . 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 $\theta$ region, even though a few do have larger $\theta$. Therefore GW detections serve as a complementary probe of this region. We also note that for very small values of $\theta$ and $m_{h_{2}}$, the search for long lived particles can be used to probe this region (eg., the MATHUSLA detector) Curtin:2018mvb .

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

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 channel 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 Eboli:1987dy ; Plehn:1996wb with the result expressed in terms of Passarino-Veltman scalar integrals. This result has also been implemented into MadGraph Alwall:2014hca taking into account the presence of a heavier SM-like scalar 101010https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/HiggsPairProduction, 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\theta)$ as compared with the SM Higgs and also decays to a pair of $h_{1}$, the total width is simply given by:

[TABLE]

where $\Gamma_{\text{SM}}(h_{2}\rightarrow 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

[TABLE]

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 Carena:2018vpt . 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 for $\sigma(pp\rightarrow h_{1}h_{1})$ versus $\sigma(pp\rightarrow h_{2}\rightarrow h_{1}h_{1})$ for all the points which give detectable GW signals, that is, those with $\text{SNR}>10$. These cross sections are both calculated at leading order but we have added a common K-factor of $2.27$ deFlorian:2013uza to take into account of higher order corrections. The colors in the left panel denote the values of $m_{h_{2}}$ and those in the middle denote $\text{BR}(h_{2}\rightarrow h_{1}h_{1})$. It is clear from these figures that the resonant cross section is always less than the full one-loop result 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 $\mathcal{O}(10)\text{fb}$, while smaller $m_{h_{2}}$ gives $\mathcal{O}(100)\text{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}\rightarrow 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 $\lambda_{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.9\text{fb}^{-1}$ data, in the di-Higgs decay channels $b\bar{b}\gamma\gamma$ Sirunyan:2018iwt , $b\bar{b}\tau^{+}\tau^{-}$ Sirunyan:2017djm , $b\bar{b}b\bar{b}$ Sirunyan:2017isc ; Sirunyan:2018zkk ; Sirunyan:2018qca ; Sirunyan:2018tki and $b\bar{b}WW/ZZ$ Sirunyan:2017guj , with a recent combination given in Sirunyan:2018two . ATLAS used $36.1\text{fb}^{-1}$ data and searched in channels $\gamma\gamma b\bar{b}$ Aaboud:2018ftw , $b\bar{b}\tau^{+}\tau^{-}$ Aaboud:2018sfw , $b\bar{b}b\bar{b}$ Aaboud:2018knk , $WW^{(\ast)}WW^{(\ast)}$ Aaboud:2018ksn and $b\bar{b}WW^{\ast}$ Aaboud:2018zhh , with also a combination of the first three channels ATLAS:2018otd . 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. deFlorian:2016spz . 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 $3\text{ab}^{-1}$ ($13\text{TeV}$), 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}}\lesssim 550\text{GeV}$ 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.

5.4 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 111111The Lorentz structure of $hWW$ coupling already gave us some insight about the nature of EW symmetry breaking at the leading order. 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 DiVita:2017vrr ; Adhikary:2017jtu ; Banerjee:2018yxy . The best sensitivities obtained for these future colliders is typically at $\mathcal{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 Plehn:2005nk , through double Higgs production at hadron colliders where the quartic coupling enters $gg\rightarrow hh$ at two-loop Bizon:2018syu or renormalizes the cubic coupling, and at lepton colliders(via Z-associated production $e^{+}e^{-}\rightarrow Zhh$ and VBF production $e^{+}e^{-}\rightarrow\nu\nu hh$), where the quartic coupling is involved in the $VVhh$ coupling at one loop Liu:2018peg . For example, Ref. Liu:2018peg found a precision of measurement of $\sim\pm 25$ for ($500\text{GeV},4\text{ab}^{-1}$ + 1 TeV, $2.5\text{ab}^{-1}$) and $\sim\pm 20$ for ($500\text{GeV},4\text{ab}^{-1}$ + 1 TeV, $8\text{ab}^{-1}$) at $1\sigma\text{C.L.}$, when the cubic coupling is marginalized in their $\chi^{2}$ analysis.

In the xSM, both the Higgs cubic and quartic couplings are modified compared with their SM counterparts:

[TABLE]

In the absence of mixing of the scalars($\theta=0$), these couplings reduce to the corresponding SM values $i\lambda_{h_{1}h_{1}h_{1}}=3m_{h_{1}}^{2}/v$ and $i\lambda_{h_{1}h_{1}h_{1}h_{1}}=3m_{h_{1}}^{2}/v^{2}$. When $\theta\neq 0$, we parametrize the deviations of these couplings from the SM values as:

[TABLE]

and show in Fig. 10 these values for the points that give detectable GWs. The features that we can read from this figure are:(1) both $\delta\kappa_{3}$ and $\delta\kappa_{4}$ are positive; (2) both variations are $\mathcal{O}(1)$ as $\delta\kappa_{3}\in(0,1)$ and $\delta\kappa_{4}\in(0,4)$; (3) a correlation exists $\delta\kappa_{4}\equiv\eta\delta\kappa_{3}$, with $\eta\approx 2.8$ for $\delta\kappa_{3}\lesssim 0.4$ and most points fall within $\eta\in(2,4)$. To understand these, we note, since phenomenological constraints requires a small $\theta$, we expect the second feature to follow naturally. The other features can be understood by Taylor expanding the couplings for small $\theta$ and we find:

[TABLE]

In the above square brackets, the terms proportional to $m_{h_{2}}^{2}/m_{h_{1}}^{2}$ dominate for the majority of the points since $v_{s}$ is concentrated at small values; $b_{3}$ is at most $\sim 10v_{\text{EW}}$, $b_{4}\lesssim 5$ from the scan and $m_{h_{2}}\gtrsim 500\text{GeV}$ generally holds. Then the above approximations show positive $\delta\kappa_{3}$ and $\delta\kappa_{4}$ and give $\delta\kappa_{4}/\delta\kappa_{3}\approx 2.5$, which is fairly close to $\eta=2.8$. For relatively large $\theta$, 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. Liu:2018peg 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. DiVita:2017vrr (for other studies, see e.g. Borowka:2018pxx ; Bizon:2018syu ; Kilian:2018bhs ; Jurciukonis:2018skr ; Borowka:2018pxx ; Maltoni:2018ttu ; Adhikary:2017jtu ; Banerjee:2018yxy ). 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%) 121212 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. Liu:2018peg , to be largely unaffected since the heavier scalar contribution in their framework is suppressed by extra powers of $s_{\theta}$. We also expect that the bar regions, taken from Ref. DiVita:2017vrr , would get tighter since the set of parameters used in their study are highly correlated here and the resonant contribution was not included in their analyses. . 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\sim 0.5$ uncertainty on the measurement of $\delta\kappa_{3}$ at $95\%$CL. While the analysis in Ref. DiVita:2017vrr does not include the quartic coupling, the contours from Ref. Liu:2018peg 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 Plehn:2005nk ; Binoth:2006ym , however, resonant contributions in xSM might enhance the cross section up to a factor of ${\cal O}(10)$ Chen:2015gva .

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 $\delta\kappa_{4}$ is equally important as $\delta\kappa_{3}$ on GW signal generation since $\eta$ is at most $4$. Thus we expect a full statistical analysis would yield roughly the same precision on the determination of $\delta\kappa_{3}$ and $\delta\kappa_{4}$, which is well improved compared with the situation at colliders.

5.5 Diboson Resonance Search Limits at Colliders

The $WW$ and $ZZ$ branching ratios become sizeable in parts of the parameter space where the trilinear coupling $\lambda_{211}$ is relatively small, as one can see from the rightmost panel of Fig. 8. In Fig. 11, we show the branching ratios of the $h_{2}\rightarrow WW,ZZ$ and $h_{2}\rightarrow h_{1}h_{1}$ channels. We see that the $WW,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, $WW,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 $VV$ and $VH$ 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 $\sqrt{s}=13\text{TeV}$ with $36\text{fb}^{-1}$ data collected in 2015 and 2016 Aaboud:2018bun . 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 $\sim 3\text{ab}^{-1}$, we obtain estimates of future projections by a simple scaling factor and obtain the dashed lines for $\sim 3\text{ab}^{-1}$ at $13\text{TeV}$ (while HL-LHC would probably run at $14\text{TeV}$). 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 $\sim 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.

6 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($\text{SNR}>10$) at LISA. The main features of the parameter space that gives detectable GWs is: $20\text{GeV}\lesssim|v_{s}|\lesssim 50\text{GeV}$, 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}$; $\theta\lesssim 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 $\mathcal{O}(1)$ deviations from the SM values and obey a relation $\delta\kappa_{4}\approx(2-4)\delta\kappa_{3}$, where $\delta\kappa_{4}$ and $\delta\kappa_{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 Chung:2012vg . 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.

7 Acknowledgments

A. Alves thanks Conselho Nacional de Desenvolvimento Científico (CNPq) for its financial support, grant 307265/2017-0. K. Sinha is supported by the U. S. Department of Energy grant desc0009956. T. Ghosh is supported by U. S. Department of Energy grant de-sc0010504. We would like to thank David Curtin, Ian Lewis, Hao-Lin Li and Ligong Bian for helpful discussions. We also thank Yi-Ming Hu and Zheng-Cheng Liang of the Tianqin program for sending us the Tianqin sensitivity curve and for valuable discussions.

Appendix A Perturbative Unitarity S Matrix

We consider a total of eleven $2\rightarrow 2$ channels of scalars and longitudinal gauge bosons scatterings. These are grouped into seven charge neutral channels $(h_{1}h_{1},h_{2}h_{2},h_{1}h_{2},h_{1}Z,h_{2}Z,ZZ,W^{+}W^{-})$, three charge-1 channels $(h_{1}W^{+},h_{2}W^{+},ZW^{+})$ and one charge-2 channel $(W^{+}W^{-})$. The leading partial wave amplitudes of these scatterings are given collectively by a symmetric matrix, which itself is a direct sum of the matrices from these three groups: $\mathcal{S}=\mathcal{S}_{0}\bigoplus\mathcal{S}_{1}\bigoplus\mathcal{S}_{2}$. The tree level perturbative unitarity requires that the absolute value of each eigenvalue of this matrix is less than $(1/2\times 16\pi)$. The non-zero elements of the $7\times 7$ matrix $\mathcal{S}_{0}$ is listed as follows(see e.g., Ref. Kanemura:2015ska for a detailed calculation):

[TABLE]

For charge-1 channels, we have:

[TABLE]

For the charge-2 channel with only one process, the matrix is simply given by $\mathcal{S}_{2}=(-2\lambda)$.

Appendix B Connection with Potential where $v_{s}=0$

The potential in Eq. 1 can be written into a different form by translating the coordinate system of $(H,S)$ such that the EW vacuum has $\langle S\rangle=0$ (see e.g., Lewis:2017dme ). In this basis, there will generally be an additional tadpole term ($b_{1}S$). Making this translation of field variables leads to the same potential being represented with different potential parameters, without changing the physics Espinosa:2011ax . 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. 1, the parameters in the basis where $b_{1}\neq 0$(denoted with a prime) can be obtained:

[TABLE]

while $a_{2},\lambda,b_{4}$ remains unchanged. On the other hand, given parameters in the tadpole basis where $v_{s}=0$ and $b_{1}\neq 0$, the parameter set in the basis used in this work can be found:

[TABLE]

where $x$ is to be solved from the cubic equation

[TABLE]

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}\neq 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

[TABLE]

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 Pruna:2013bma ; Robens:2015gla ; Carena:2018vpt . This specific model correspond to a special limit of the potential here.