Complementarity of the future $e^+ e^-$ colliders and gravitational waves in the probe of complex singlet extension to the Standard Model

In this work, we study the future probes of the complex singlet extension to the Standard Model (cxSM). This model is possible to realize a strongly first-order electroweak phase transition (SFOEWPT). The cxSM naturally provides dark matter (DM) candidate, with or without an exact $\mathbb{Z}_2$ symmetry in the scalar sector. The benchmark models which can realize the SFOEWPT are selected, and passed to the current observational constraints to the DM candidates, including the relic densities and the direct detection limits set by the latest XENON1T results. We then calculate the one-loop corrections to the SM-like Higgs boson decays and the precision electroweak parameters due to the cxSM scalar sector. We perform a global fit to the benchmark models and study the extent to which they can be probed by the future high-energy $e^+ e^-$ colliders, such as CEPC and FCC-ee. Besides, the gravitational wave (GW) signals generated by the benchmark models are also evaluated. We further find that the future GW detector, such as LISA, is complementary in probing the benchmark models that are beyond the sensitivity of the future precision tests at the $e^+ e^-$ colliders.


Introduction
The observed baryon asymmetry of the Universe (BAU) and the nature of dark matter (DM) are two of the leading puzzles that motivate the new physics beyond the Standard Model (BSM). One compelling scenario to achieve the BAU is the electroweak baryogenesis (EWBG) [1][2][3][4][5][6]. To preserve the baryon asymmetry generated, a strongly first-order electroweak phase transition (SFOEWPT) is necessary. It is well-known that the Standard Model (SM) itself cannot realize the SFOEWPT, since the 125 GeV Higgs boson discovered at the LHC [7,8] is too heavy [9][10][11][12]. On the other hand, there is no viable DM candidate in the SM. To achieve the SFOEWPT and provide possible DM candidate, the SM should be extended.
The simplest realization of the SFOEWPT can be achieved through adding one real scalar singlet to the SM Higgs sector [13][14][15][16][17][18][19][20][21]. If we impose the Z 2 symmetry under which only the real scalar is odd, this extension can also provide a cold DM candidate since the discrete symmetry forbids the mixing between the neutral doublet and the real singlet. This scenario admits strongly first order and two-step phase transition in which the singlet scalar acquires a vacuum expectation value (vev) before the electroweak symmetry breaking. However, in this scenario, the deviations in the hZZ and hhh couplings are induced at loop level. Thus, no future Higgs factories have the required sensitivities to probe the evidence of such SFOEWPT [16,22]. Besides of the extension of one real scalar singlet, the SFOEWPT can also be realized in the complex scalar extension to the SM (cxSM), as discussed in Refs. [23][24][25][26][27][28]. DM candidates can naturally arise in the cxSM, in both Z 2 symmetric and Z 2 breaking scenarios 1 . Hence, the cxSM is appealing to address both the SFOEWPT issue and the DM candidate at one shot. The next question is whether the Z 2 symmetric and Z 2 breaking scenarios of the cxSM with SFOEWPT and DM candidate can be distinguished by the future experiments.
Direct searches for the extended scalar sector beyond the SM have been carried out in the Large Hadron Collider (LHC) experiments [34][35][36][37][38][39][40][41][42][43]. No signal has been reported so far. Due to the small mixing effects of the SU(2) L singlets, it is expected that the direct search for the scalars from the complex singlet is very challenging at the LHC [17,44,45] 2 . Complementary to the direct searches, the precision measurements of the Higgs boson properties could shed light on the underlying new physics. Several well-known proposals have been made to build the nextgenerational Higgs factory, such as the Circular Electron Positron Collider (CEPC) in China [47,48], the electron-positron stage of the Future Circular Collider (FCC-ee) at CERN [49], and the International Linear Collider (ILC) in Japan [50,51]. Each facility is proposed to run at √ s = 240 − 250 GeV to produce 10 5 − 10 6 SM-like Higgs bosons, aiming to reach sub-percentage precision measurement of its couplings. Besides, they will also run at the Z-pole to improve the precisions on the measurement of SM parameters by a factor of 20 − 200 over the results from Large Electron Positron (LEP) Collider [52]. With such incredible improvements in the precision measurements, a number of studies have been carried out to look for the BSM effects through both tree-level and one-loop corrections to the productions [16,[53][54][55] and decays [22,[56][57][58][59] of Higgs 2 The complex singlet extension to the SM

The Higgs potential and global symmetries
We consider the extended Higgs sector beyond the SM by introducing a complex scalar singlet S of the SU(2) L . The most general scalar potential in this extension is expressed as [23] V (Φ , S) = µ 2 |Φ| 2 + λ|Φ| 4 + δ 2 2 where Φ is the SU(2) L Higgs doublet breaking the electroweak symmetry. The parameters in the first line of Eq. (2.1) are real, and the other parameters in the second and third lines of Eq. (2.1) are generally complex. Two possible global symmetries can be imposed to the above Higgs potential: • A discrete Z 2 symmetry of S → −S can be imposed to eliminate all terms with odd powers of S, which include δ 1 , a 1 , c 1 ,2 .
• A global U(1) symmetry of S → e iα S eliminates all terms with complex coefficients (δ 1 ,3 , If the complex scalar field S does not obtain a zero-temperature vev, the discrete Z 2 symmetry has to be introduced to stabilize the scalar singlet and enable DM candidate(s). Under a further global U(1) symmetry, this cxSM model yields two degenerate stable DM particles (the two components in S). This case with only the terms in the first line of Eq. (2.1) is very similar to the real singlet model. By including one U(1) breaking term, for instance the b 1 term, the real and imaginary parts of S are still stable but not identical anymore. Below we refer this more general case as the Z 2 symmetric scenario with the following scalar potential On the other hand, if the S field receives a zero-temperature vev and thus the real component of S mixes with the neutral Higgs of Φ, the U(1) and Z 2 symmetries are both spontaneously broken by the singlet vev and the Goldstone boson from the imaginary part of S is stable but massless. To provide a viable DM candidate, a soft breaking of the global U(1) symmetry is introduced to generate a mass for it. The U(1) breaking requires that one or more terms in the second and third lines of Eq. (2.1) does not vanish. We demand b 1 = 0 here as well, and the U(1) symmetry is both spontaneously and softly broken. Now the spontaneously broken Z 2 symmetry may lead to the cosmological domain wall problem [76,77]. To solve this problem, one can further introduce one or more of δ 1 , a 1 , c 1,2 terms to explicitly break the Z 2 symmetry. We consider the following potential with a non-vanishing a 1 as in Ref. [23] V (Φ , S) We refer the above potential as the Z 2 breaking scenario below. One should keep in mind that, although we follow the choices of Ref. [23] in the rest of this paper, the scalar potential for achieving the above purposes is not unique.

The Z 2 symmetric scenario
To minimize the scalar potential, we represent the complex scalar singlet as S = 1 √ 2 (S + iA) and the Higgs doublet as Φ = (0 , h/ √ 2) T . In the Z 2 symmetric case, we only have the SM Higgs doublet developing a vev (v) and the a 1 term is vanishing. From Eq. (2.2), the field-dependent scalar potential at the tree level becomes By minimizing the potential, one arrives at the following condition The mass spectrum is obtained as follows With the exact Z 2 symmetry, h and S do not mix. Both S and A are stable and regarded as the DM candidates in our following discussions. Altogether, the parameters in the generic basis and the physical basis are generic basis : µ 2 , λ , δ 2 , |b 1 | , b 2 , d 2 ; (2.7a) with the fixed inputs as M h = 125 GeV and v ≈ 246 GeV. The ranges of remaining parameters we take for the scan are

The Z 2 breaking scenario
In the Z 2 breaking scenario, the field-dependent scalar potential at the tree level reads Both h and S obtain vevs in this case. The corresponding minimization conditions are The mass spectrum of the scalars for the Z 2 breaking scenario is obtained as follows The mass eigenstates after diagonalizing the CP-even scalars are with the masses of h 1 and h 2 being M 1 and M 2 , respectively. The CP-odd component A will not develop a vev and will be treated as the DM candidate for the later discussion.
In terms of the mass eigenstates, our parameter inputs can be traded into the CP-even scalar masses and the mixing angle as λ = 1 2v 2 cos 2 θM 2 1 + sin 2 θM 2 2 , (2.13a) Altogether, the parameters in the generic basis and the physical basis are generic basis : µ 2 , λ , δ 2 , a 1 ,

The Higgs self-couplings
In the Z 2 symmetric scenario, the relevant cubic and quartic Higgs self-couplings are listed below The cubic and quartic Higgs self-couplings of the SM-like Higgs bosons are the same as those in the SM case, while the other two Higgs couplings λ hSS and λ hAA are relevant for the Higgs boson self-energy corrections at the one loop level.
In the Z 2 breaking scenario, the relevant cubic Higgs and quartic self-couplings in the physical basis can be expressed as follows The cubic and quartic self-couplings recover the SM couplings when the mixing angle θ → 0. We define the deviations of the cubic and quartic self-couplings of the SM-like Higgs as δκ 3 ≡ λ 111 /λ SM hhh − 1 and δκ 4 ≡ λ 1111 /λ SM hhhh − 1. The correlation between them guarantees the treelevel driven SFOEWPT.

Unitarity and stability
In order to have a well-defined Higgs potential, a set of theoretical constraints should be taken into account. The Lee-Quigg-Thacker unitarity bound [78,79] should be imposed so that the quartic couplings are not too large. In both Z 2 symmetric and Z 2 breaking scenarios, the quartic terms of the Higgs potential are the following By taking the neutral states of |π The s-wave unitarity conditions are imposed such that |ã i 0 | ≤ 1, withã i 0 being all eigenvalues of matrix a 0 above. By using the relations in Eqs. (2.13), the perturbative unitarity condition can impose bounds to the Higgs boson masses and mixings. In addition, one should impose the following tree-level stability conditions so that the scalar potential is bounded from below at the large field values Here, the last term is necessary for δ 2 < 0.

The global minimum
In terms of the classical fields, there may be three different configurations for the symmetry breaking: for the Z 2 symmetric (breaking) cases. As the temperature cools down, the symmetry breaking may occur either by one step via O → B, or by two steps via O → A → B. The one-step phase transition occurs if the configuration-B is the only possible Higgs potential minimum, and the twostep phase transition occurs if both configure-A and configuration-B coexist as the Higgs potential minimum.
The EWSB vacuum solution of B should be the lowest one of the scalar potential, while the origin point of O should be the highest one. The vacuum configurations of A and B are obtained by solving the following cubic equations The numerical solutions are then fed into V 0 (A) and V 0 (B), and the global minimum condition V 0 (B) ≤ V 0 (A) will be imposed.

The constraints on the DM candidates
In the Z 2 symmetric scenario with v s = 0, both S and A are regarded as the DM candidates, which gives a two-component DM case. In the Z 2 breaking scenario, only the CP-odd scalar A becomes the DM candidate. The annihilation processes that contribute to the DM relic density in the two cases are shown in Fig. 1. The relic density typically exhibits one (two) dip(s) with the DM mass being around M h /2 (M 1 /2 or M 2 /2), due to the enhancement of the annihilation cross section near the h (h 1 , h 2 ) resonance(s) in the Z 2 symmetric (breaking) scenario. Several ongoing DD experiments are looking for DM scattering off atomic nuclei, including XENON1T [80] and PandaX-II [81]. No conclusive observation has been reported so far. For the DM mass range of O(10) − O(10 3 ) GeV, XENON1T has set the most stringent lower limit on the SI DM-nucleon scattering cross section as σ SI 10 −46 − 10 −44 cm 2 [80,82]. For the Z 2 symmetric case, the SI scattering processes are mediated only by the SM-like Higgs boson h; while for the Z 2 breaking case, the SI scattering processes are mediated by two CP-even scalars of h 1 ,2 .  Figure 2. The rescaled SI cross sections of the DM candidate for the Z 2 symmetric scenario (a) and the Z 2 breaking scenario (b). The grey points are those oversaturate the relic density. The blue points satisfy the relic density requirement but have been ruled out by the XENON1T limit (solid line). The red points are those satisfy both the relic density requirement and the current direct detection limit by the XENON1T, and pass the SFOEWPT criterion. The future projected limit by the XENONnT is also displayed (dashed line).
The corresponding cross sections are given by [24,83,84] with the nucleon form factors of f T s and f (p) T q . The possible cancellation between the h 1 and h 2 diagrams [30], as indicated by Eq. (3.6b), leads to further suppressed scattering cross section for the Z 2 breaking scenario compared with the Z 2 symmetric scenario.
In practice, we first produce the FeynRules [85] files by implementing the cxSM model parameters and interactions. The model files are then fed to MicrOMEGAs [86] to calculate the DM relic density for the cxSM model (denoted as Ω cxSM h 2 ) and the SI scattering cross section σ SI . The above quark/gluon-nucleon form factors are taken as the default values in MicrOMEGAs. The current measurements of the cold DM relic density are given as Ω DM h 2 = 0.1138 ± 0.0045 (WMAP) [87] or Ω DM h 2 = 0.1196 ± 0.0031 (Planck) [88]. After the scan of parameter spaces in Eqs. (2.8) and (2.15) by imposing the above theoretical constraints, the survived points that oversaturate the relic density are further rejected. For those points that undersatuarate the relic density, we rescale the SI cross section by  and compare with the latest limit from the XENON1T [80]. In Fig. 2, the rescaled SI cross sections of model points are evaluated and the current limit set by the XENON1T experiment is added as reference. The model points satisfying both the relic density constraint Ω cxSM h 2 < Ω DM h 2 and the current XENON1T DD limit and passing the SFOEWPT criterion are marked in red. We will apply this constraint on our benchmark points for later studies. We also display the future direct detection limit set by the XENONnT 4 . As seen from the plots, the red points are expected to be testable by the future facilities such as XENONnT, PandaX-4T [89] or LZ [90]. 4 The SFOEWPT, GW signals and the precision test at e + e − colliders 4.1 The finite-temperature effective potential In order to evaluate the EWPT in the cxSM, we follow the recipe of Ref. [25] by using the hightemperature expanded effective potential in order to avoid the gauge dependence problem. The EWPT is driven by the cubic terms in the effective potential. Thus, we take the following hightemperature expansion of where the finite temperature corrections are given by the thermal mass contributions Π h (T ) and Π S (T ).
The history of the phase transitions from high temperature to the vacuum today is displayed in Fig. 3 for the Z 2 symmetric scenario (left) and the Z 2 breaking scenario (right). One can see that the Universe follows a two-step symmetry breaking in both cases in the space of two order parameters for doublet and singlet scalars. The global minimum of both cases at high temperature happens at the spot "O" with restored electroweak symmetry. For the Z 2 symmetric scenario with zero |a 1 |, when the temperature of the Universe falls down to "A", we expect a first minimum with S = 0 and h = 0. Along with the further temperature decreasing, a second minimum of with S = 0 and h = 0 develops, which eventually becomes the present vacuum at the spot "B". The critical temperature T c is given when two minima of "A" and "B" are degenerate. As there is a barrier between these two minima, a first-order EWPT happens. For the Z 2 breaking scenario with a 1 = 0, the origin at high temperature is shifted along the S direction from spot O to O' due to the non-vanishing a 1 term. The electroweak symmetry then follows a two-step phase transition process as well.

The GW signals
The GW signals generated during the EWPT depend on the evaluation of the tunneling rate per unit time per unit volume, which is given by [91] Γ with S 3 being the Euclidean action of the critical bubble that minimizes the finite-temperature action of The bubble nucleation temperature T n is defined as the probability for a single bubble to be nucleated within one horizon volume being O(1), that is where M pl = 1.2 × 10 19 GeV is the Planck mass, and ζ 3 × 10 −2 . Numerically, this equation implies that S 3 (T n )/T n ≈ 140 [92]. Two other parameters that are directly relevant to the GW signal calculations are given by where ρ vac stands for the latent heat released during the EWPT, H n is the Hubble parameter at T n , and ρ * rad = g * π 2 T 4 n /30 with g * representing the relativistic degrees of freedom at T n . Typically, a relatively larger α accompanied with a small β/H n will trigger the SFOEWPT and a significant GW signal.
The observed GW signal is characterized by the energy spectrum Ω GW (f )h 2 [92] Ω The total energy spectrum here is dominated by the summation of two terms: (1) the sound waves after the bubble collisions; (2) MHD turbulence The GW signals from the sound wave contribution is given by where κ v represents the fraction of latent heat transferred into the bulk motion of the fluid, and was estimated in Ref. [67]. The peak frequency f sw is rescaled from its values at the phase transition by The MHD turbulence contribution to the GW energy spectrum is written as where κ tu ≈ 0.1κ v . The peak frequency from the MHD turbulence term is given by To compatible with EWBG, the wall velocity v w here is obtained as a function of α. [93][94][95][96] after taking into account Hydrodynamics. The discovery prospects of the GW signals are determined by the signal-to-noise ratio (SNR) [67] where Ω exp (f )h 2 stands for the experimental sensitivity for the proposed GW programs. T is the mission duration in years for each experiment, and we assume it to be five here. The factor δ counts the number of independent channels for cross-correlated detectors, which is taken to be 1 for the LISA program [97]. In practice, we evaluate the SNRs for each benchmark points that achieve the SFOEWPT. For the LISA program, we take the threshold SNR of 50 for discovery. This corresponds to the least sensitive configuration of C4 with four links [67].

4.3
The precision test at the future e + e − colliders

The one-loop corrections to the Higgs boson couplings
The 125 GeV SM-like Higgs boson can receive corrections from both SM sector as well as the extended scalar sector in the cxSM. The SM-like Higgs couplings normalized to its SM value, κ, is defined as [59]: where g SM(cxSM) tree and g

SM(cxSM) loop
are the couplings in the SM (cxSM) at tree and one-loop level, respectively. In the Z 2 breaking case, the couplings of SM-like Higgs boson h 1 to all SM fields are universally proportional to a factor of cos θ due to the singlet-doublet mixing. The new one-loop contributions 5 from cxSM to h 1 bb and h 1 ZZ are shown in Fig. 4. Note that, in the Z 2 symmetric case, although we have hSS and hAA coupling, the contributions in Fig. 4 are zero due to the vanishing S(A)V V and S(A)f f couplings. As a consequence, in the Z 2 symmetric case, the modification of the couplings mainly comes from the Higgs self energy corrections from S and A loops. Hence, the deviations of κ's in the Z 2 symmetric case are quite universal. The general vertices of the SM-like Higgs with a pair of gauge bosons hV V and SM fermions hff take the following forms where (q , p 1 , p 2 ) represent the momenta of the SM-like Higgs boson and two final-state particles. The κ i for each vertex is given by Γ 1 hV V and Γ S hff for the hV V and hff vertices as In practice, the one-loop corrections to the SM-like Higgs boson couplings are evaluated by adopting the on-shell renormalization scheme [98]. All counter terms, renormalization constants and renormalization conditions are implemented into model files of FeynArts [99], which is then used to generate all possible one-loop diagrams for corresponding couplings for cxSM. After that, FormCalc [100] is used to calculate the full loop level couplings. The numerical results are performed by LoopTools [101] 6 . However, the observables in each experiment are the signal strength µ i 's instead of κ's. Thus, for each channel, we will calculate the signal strength by where κ i,f are the normalized coupling relevant for production and decay, and κ width represents the ratio of the total width of the SM-like Higgs boson in cxSM to that in SM. To constrain the model parameters from the current and future Higgs boson precision measurements, a global fit to the observed signal strength is performed with the profile likelihood method. The χ 2 is defined as where we sum over all measurements available in experiments and neglect the correlations between different measurements. In our analyses, µ obs i are set to be the SM value µ obs i = 1, for the future colliders. The estimated errors σ µ i are listed in Tab. 1 for the future circular e + e − colliders (CEPC and FCC-ee), and also in Tab. 2

The electroweak precision tests
Beside of the SM-like Higgs boson couplings, the model will also change the electroweak observables. To take this into account, the Peskin-Takuechi parameters [102]  to represent the electroweak precision measurements. However, S and A have vanishing gauge couplings in the Z 2 symmetric case. Thus, they do not modify the S, T and U parameters. In the Z 2 breaking case, the expressions for the modifications are given by All these modifications are proportional to the CP-even mixing angle of s θ . The loop functions of B 0 and B 00 follow the convention in LoopTools [101]. We perform a global fit to the electroweak observables using Gfitter [103] with current electroweak precisions [52] and future prospects [47,49,104]. Unlike the case for Higgs signal strength, the χ 2 constructed from S, T and U also includes the correlations among them. The corresponding χ 2 is thus defined as with X i = (∆S , ∆T , ∆U ) being the contributions from the cxSM, andX i being the corresponding best-fit central values 7 . The σ 2 ij ≡ σ i ρ ij σ j are the error matrix with uncertainties σ i and correlation matrix ρ ij given in Tab. 3 for different experiments.

Numerical results
Practically, we implement the tree-level effective potential and the high-temperature expansion in Eq. (4.1) into the CosmoTransitions [105]. The temperature-dependent minima of "A" parametrized by ( h , S ) = (0 , ϕ A S ), and "B" parametrized by ( h , S ) = (ϕ B h , 0) for Z 2 symmetric scenario or ( h , S ) = (ϕ B h , ϕ B S ) for Z 2 breaking scenario are similarly evaluated by using Eq. (4.1). For the numerical presentation below, we take the data points that not only evade all theoretical constraints and DM constraints, but also achieve the SFOEWPT. The SFOEWPT is characterized by obeying the condition ϕ B h /T n ≡ v n /T n 1 based on the requirement of the baryon number preservation criterion [106][107][108]. The CosmoTransitions [105] is used for solving the nucleation temperatures T n , as well as the GW signal parameters of α and β/H n . The solutions of (T n , α , β/H n ) for each parameter points will be used for the SNRs of the GW signals according to Eq. (4.12). For the precision test of cxSM at future colliders, the results of χ 2 in Eq. (4.17) and Eq. (4.19) are linearly combined for both Higgs boson and the electroweak precision measurements.

SFOEWPT and GW
The GW spectrum is characterized by parameters of α and β/H n defined in Eq. (4.5), and their values can be fixed by the cxSM potential. As α and β/H n represent the latent heat released by EWPT and the reversed duration of the EWPT, respectively, significant GW observation typically prefers larger α and smaller β/H n values. In Fig. 5, we display the parameter dependences of (α , β/H n ) on the cxSM parameters for both Z 2 symmetric and Z 2 breaking scenarios. For the Z 2 symmetric case, the α and β/H n parameters are displayed in the (d 2 , δ 2 ) plane. Among the model points we generated, the values of α and β/H n are found to be uniformly distributed with respect to the d 2 in the cxSM potential; while α (β/H n ) values are enhanced (suppressed) with relatively larger inputs of δ 2 . This is because small |Φ| 2 |S| 2 coupling of δ 2 tends to reject the potential barrier and also the SFOEWPT. For the Z 2 breaking case, we show the values of α and β/H n in the plane of (a   along the ϕ S direction should not be rather sizable to break the discrete symmetry in this case. Although we present all date points satisfying the SFOEWPT criterion here, one should note that those points with β/H n 10 4 produce too large peak frequencies as from Eqs. (4.9) (4.11), and too small power spectrum as from Eqs. (4.8) (4.10). Thus, such points are impossible to be detected by the GW detectors which are mostly sensitive to milliHz frequencies. Figure 7. The GW spectrum for two benchmark points from Z 2 and / Z 2 cases.
In Fig. 6, we show the bubble nucleation temperature T n in the GW parameter plane of (α , β/H n ). The lower T n one obtains, the stronger the EWPT becomes. In principle, as a result, we can have increased α and decreased β/H n . The realistic situation of their relationship might be more complicated to achieve the bubble nucleation condition while comparing different specific models. For instance, for the Z 2 breaking case, the origin is shifted at high temperature. As seen in the plots, the Z 2 breaking scenario exhibits lower T n and larger α as a result of a two-step bubble nucleation. The values of β/H n parameter span a broader space and can be relatively larger for the Z 2 breaking case, compared with the Z 2 symmetric case with nearly the same temperature.  Table 4. Two benchmark points for the Z 2 symmetric and Z 2 breaking cases. The marks of × ( √ ) represent whether the benchmark point is without (within) the precision of the corresponding collider searches.
To compare the GW signal spectra, we list two benchmark points for Z 2 symmetric and Z 2 breaking cases in Tab. 4. The CP-odd scalar masses of these two benchmark points are close to each other, i.e. M A = 950 GeV (Z 2 symmetric) and M A = 963 GeV (Z 2 breaking) respectively. The Z 2 symmetric benchmark point cannot be searched for via the precision measurements of the Higgs boson decays at the future HL-LHC and e + e − colliders, while it yields a SNR of O(10 4 ) at the LISA. The Z 2 breaking benchmark point can be probed via both the precision measurements of the Higgs boson decays at the future e + e − colliders and the GW spectrum at the LISA, with an SNR of O(10 6 ). Their GW spectra Ωh 2 versus the frequency f are displayed in Fig. 7, together with the viable signal regions of different ongoing/upcoming GW detection programs. The benchmark point in the Z 2 symmetric case exhibits higher T n and corresponding decreased α. The value of β/H n in Z 2 breaking scenario is larger than that in Z 2 symmetric case, which leads to a higher peak of frequency in the Z 2 breaking case.

5.2
The precision tests at the colliders Figure 8. The normalized SM-like Higgs boson couplings of κ b (upper pannels) and κ Z (lower pannels) versus d 2 and δ 2 for the Z 2 symmetric scenario. The red points beyond search sensitivities of the HL-LHC and future e + e − colliders are denoted as "Nightmare" model points.
The precision tests are made by the combined χ 2 fit of the SM-like Higgs boson measurements and the electroweak precision measurements according to Eqs. (4.17) and (4.19). In Figs. 8 and 9, we display two couplings of κ b and κ Z for both Z 2 symmetric and Z 2 breaking cases. For the Z 2 symmetric case, the normalized Higgs boson couplings are displayed versus parameters of (d 2 , δ 2 ); while for the Z 2 breaking case, they are shown for the physical parameters of (a 1/3 1 , v s , θ). The model points in grey are within the sensitivities of both the HL-LHC and any of the future e + e − colliders (including CEPC, FCC-ee and ILC), while the blue points are only within the sensitivities of e + e − colliders. The red points are those beyond search sensitivities of the HL-LHC and future e + e − colliders. We denote them as "nightmare" model points [16,22]. We found for the Z 2 symmetric case, as shown in Fig. 8, all points give the deviation of the SM-like Higgs boson couplings from the SM values by at most 0.5 % and are thus beyond the sensitivities of the HL-LHC and any future e + e − colliders. This is due to the absence of tree-level correction and small oneloop corrections from the SM-like Higgs boson self-energy terms. In comparison, in Z 2 breaking scenario, there can be sizable one-loop corrections to the SM-like Higgs boson decays through the cxSM sector, as were previously shown in Fig. 4, besides the tree-level correction given by g cxSM tree /g SM tree = cos θ. Thus, as seen in Fig. 9, some points are within the search sensitivities of the HL-LHC while some other points with larger κ couplings can be probed by future e + e − colliders.
In Fig. 10, we combine the experimental sensitivities of the colliders and the GW signal probes via the LISA interferometer in the (α , β/H n ) plane. For the Z 2 symmetric case, the one-loop corrections from the cxSM sector are typically small. Correspondingly, the model points are beyond the sensitivity of precision tests at any current or planned colliders. The LISA interferometer is likely to probe the models points with relatively large values of α and small values of β/H n . For the Z 2 breaking case, the corrections of treel-level and one-loop effects from the cxSM sector become significant. Correspondingly, we found a majority of model points are within the search sensitivities of the HL-LHC and the future e + e − runs. A smaller fraction of model points (denoted in red) are beyond the search limits by either the Higgs measurements at the future colliders or the LISA interferometer. Nevertheless, these points are all within the sensitivity of future DM DD experiments.
Finally, in Fig. 11, we show the expected sensitivities of future colliders to δκ 3 in the Z 2 breaking case. To the right of the colored (grey) bars, the corresponding colliders are sensitive to the measurement of cubic Higgs coupling at 68% (95%) C.L. . Note that, in the Z 2 symmetric case, the cubic and quartic Higgs couplings are the same as those in SM. Thus, we do not have  Figure 11. The precision measurement of δκ 3 in the Z 2 breaking case from different collider experiments. The vertical position of the points is irrelevant. The colored and grey shaded regions correspond to 68% and 95% C.L. regions, respectively, for δκ 3 from Ref. [109]. any sensitivity in these measurements. While, in the Z 2 breaking case, the Higgs self-couplings do differ from the SM as shown in Eqs. (2.17). The sensitivity of future colliders to the cubic coupling, however, is much lower than the precision measurements of other couplings.

Conclusion
In this work, we study the future experimental tests of the cxSM. In both Z 2 symmetric and Z 2 breaking scenarios in this model, the realization of the SFOEWPT is achieved and scalar DM candidate(s) can be provided. Future experimental facilities at the high intensity/energy frontiers, such as GW detection and e + e − colliders, can test the visible parameter space of this complex scalar model and discriminate the scenarios with or without the discrete Z 2 symmetry. We apply theoretical constraints and DM constraints from relic density and the latest XENON1T limit on the parameter space of the cxSM, and also require they pass the SFOEWPT criterion. By combining the χ 2 fit of the SM-like Higgs boson measurements and the electroweak precision measurements, we estimate whether the model points can be accessible at the future e + e − colliders.
In the Z 2 symmetric scenario, the complex scalar singlet S does not develop a vev and a quadratic term of S is introduced to break a global U(1) symmetry. As a result, the real part of S does not mix with the SM Higgs and both the real and imaginary parts become the DM candidates. This scenario admits a two-step SFOEWPT in the way that the scalar singlet acquires a vev at high temperature prior to the electroweak symmetry breaking. In this scenairo, as S does not mix with the SM Higgs doublet, there is no tree-level correction to the SM-like Higgs couplings and the oneloop corrections from the Higgs boson self-energy terms are very small. We find that none of the generated model points are within the sensitivities of the HL-LHC and any future e + e − colliders.
In the Z 2 breaking scenario, the complex scalar singlet S develops a non-vanishing vev and an additional linear term of S is introduced to break the discrete Z 2 symmetry. Thus, besides the sizable loop corrections, the mixing of the complex singlet and the SM Higgs doublet induces a tree-level correction to the SM-like Higgs couplings which is the cosine of the mixing angle. The CP-odd component of the complex singlet serves as the DM candidate. This scenario also achieves a two-step SFOEWPT driven at tree-level. It turns out that a majority of model points can be covered by the precision Higgs measurements at future colliders.
We also find that, in both Z 2 symmetric and Z 2 breaking scenarios, some of the points without the sensitivities of future colliders are accompanied with sizable signal-to-noise ratio around f ∼ O(10 −4 ) − O(1) Hz for their GW signals. Future space-based GW interferometer, such as LISA, can thus probe such "nightmare" parameter space. In addition, all the model points realizing the SFOEWPT and satisfying DM constraints are within the sensitivity of future DM DD experiments.