Probing electroweak phase transition with multi-TeV muon colliders and gravitational waves

We study the complementarity of the proposed multi-TeV muon colliders and the near-future gravitational wave (GW) detectors to the first order electroweak phase transition (FOEWPT), taking the real scalar extended Standard Model as the representative model. A detailed collider simulation shows the FOEWPT parameter space can be greatly probed via the the vector boson fusion production of the singlet, and its subsequent decay to the di-Higgs or di-boson channels. Especially, almost all the parameter space yielding detectable GW signals can be probed by the muon colliders. Therefore, if we could detect stochastic GWs in the future, a muon collider could provide a hopeful crosscheck to identify their origin. On the other hand, there is considerable parameter space that escapes GW detections but is within the reach of the muon colliders. The precision measurements of Higgs couplings could also probe the FOEWPT parameter space efficiently.

The FOEWPT can manifest itself at two different kinds of experiments: the gravitational wave (GW) detectors and the high energy particle colliders. In the former case, the stochastic GWs generated during the FOEWPT are expected to be detectable at a few near-future space-based laser interferometers such as LISA [50], BBO [51], TianQin [52,53], Taiji [54,55] and DECIGO [56,57]. While in the latter case, the BSM physics related to FOEWPT might be probed at the colliders [58], such as the CERN LHC and future protonproton colliders including HE-LHC [59], SppC [60] and FCC-hh [61], or future electronpositron colliders such as CEPC [62], ILC [63] and FCC-ee [64]. Generally speaking, the hadron colliders have high energy reach but suffer from the huge QCD backgrounds, while the electron-positron colliders are very accurate but limited by the relatively low collision energy due to the large synchrotron radiation.
A muon collider might be able to offer both high collision energy and clean environment to probe the FOEWPT. On one hand, thanks to the suppressed synchrotron radiation compared to the electron, the energy of a muon collider can reach O(10) TeV. What's more, the entire muon collision energy can be used to probe the short-distance reactions (hard processes). In contrast, at a pp collider such as LHC, only a small fraction of the proton collision energy is available for the hard processes. On the other hand, due to the small QCD backgrounds, the muon collider is rather clean, allowing very precise measurements. The physics potential of a high energy muon collider has been discussed since the 1990s [65,66], while it receives a renewed interest recently [67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84].
In this work, we investigate the possibility of probing FOEWPT at a multi-TeV muon collider and the complementarity with the GW experiments, taking the xSM as the benchmark model. Although the xSM is simple, it has captured the most important features of the FOEWPT induced by tree level barrier via renormalizable operators [85], and can serve as the prototype of many BSM models that trigger the FOEWPT. For the muon collider setup, we follow Ref. [78] to consider collision energies of 3, 6, 10 and 30 TeV, with integrated luminosities of 1, 4, 10 and 90 ab −1 , respectively. This paper is organized as follows. We first introduce the xSM and derive its parameter space for FOEWPT in Section 2, where we also discuss the GW signals and their detectability at the future LISA detector. The phenomenology at high energy muon colliders is studied in Section 3, where both the direct (i.e. resonant production of the real singlet) and indirect (i.e. the Higgs coupling measurements) searches are considered. The complementarity between collider and GW experiments is also discussed. Finally, we conclude in Section 4.

The model
Up to renormalizable level, the scalar potential of xSM can be generally written as which has eight input parameters. However, one degree of freedom is unphysical due to the shift invariance of the potential under S → S + σ; in addition, the measured Higgs mass M h = 125.09 GeV and vacuum expectation value (VEV) v = 246 GeV put another two constraints, leaving us only five free physical input parameters.
and then the mass term of the two neutral scalars reads Diagonalizing M 2 s yields the mass eigenstates h 1 , h 2 and the mixing angle θ between them, Here we assume the lighter state h 1 is the SM Higgs-like boson.
The requirement that (v, v s ) is an extremum of Eq. (2.1) yields two relations [12] where the coefficients λ, a 1 and a 2 can be further expressed in terms of M h 1 , M h 2 and θ, with c θ and s θ being short for cos θ and sin θ, respectively. Fixing M h 1 = M h = 125.09 GeV and v = 246 GeV, we can use the following five parameters as input, and derive other parameters such as µ 2 , λ via Eq. (2.5) and Eq. (2.6). We use the strategy described in Appendix A to obtain the parameter space that satisfies the SM constraints. The dataset is stored in form of a list of the five input parameters in Eq. (2.7), and then used for the calculation of FOEWPT and GWs in the following subsection.

FOEWPT and GWs
The scalar potential V in Eq. (2.1) receives thermal corrections at finite temperature, becoming where we only keep the gauge invariant T 2 -order terms [86,87], and (2.9) In our convention, b 1 = 0, the tadpole term for s only arises at finite temperature. This term is found to be suppressed in most of the parameter space [7,11], however for completeness we also include it in numerical study. As we will see very soon, the tadpole has a nonnegligible impact on the FOEWPT pattern. Thermal corrections change the vacuum structure of the scalar potential. In suitable parameter space, there exists a critical temperature T c at which the potential V T in Eq. (2.8) has two degenerate vacua, one with h = 0 (EW-symmetric) and the other with h = 0 (EW-broken). Initially, the Universe stays in the EW-symmetric vacuum. As the Universe expands and the temperature falls below T c , the h = 0 vacuum is energetically preferred and the Universe acquires a probability of decaying to it. The decay rate per unit volume is [88] which we take as the criterion for a FOEWPT.
For each data point derived in last subsection, we calculate T n by solving Eq. (2.11) with the Python package CosmoTransitions [90]. Around 10% of the data can trigger a FOEWPT. The left panel of Fig. 1 shows the collection of FOEWPT data points by plotting the initial and final states of the vacuum decay (0, v i s ) → (v f , v f s ). For a successful EWBG, the phase transition should be strong [91,92] v f /T n 1, (2.12) such that the EW sphaleron process in the EW-broken vacuum is suppressed. Hereafter we will focus on data points satisfying Eq. (2.12). We found that many data points yield a decay pattern of (0, 0) → (v f , v f s ), but there are also considerable fraction of data that , then the fraction of data points falling in this pattern is around 8.8%, while in Ref. [12] the corresponding fraction is 99%. We have checked that the difference comes from the treatment of the thermal tadpole term in Eq. (2.8): we keep this term, while Ref. [12] drops it. Therefore, the tadpole term actually has a considerable impact on the FOEWPT pattern.
A FOEWPT generates stochastic GWs mainly through three sources: bubble collisions, sound waves in the plasma and the magneto-hydrodynamics turbulence [93]. After cosmological redshift, those GWs today typically peak at f ∼ mHz [94], which is the sensitive region of a few next-generation space-based interferometers mentioned in the introuction. To obtain the GW spectrum today, we derive the following two parameters for each FOEWPT data point is effective potential difference between the true and false vacua, and g * ∼ 100 is the number of relativistic degrees of freedom. In other words, α is the transition latent heat over the radiation energy, while β/H is the Universe expansion time scale over the phase transition duration. Refs. [94][95][96] point out that the GW spectrum Ω GW (f ) of a FOEWPT can be expressed as numerical functions of (α, β/H, v b ), where v b is the bubble expansion velocity. Taking v b = 0.6 as a benchmark, we are now able to calculate Ω GW (f ) for each FOEWPT data point. 2 The suppression factor coming from the short duration of the sound wave period has been taken into account [102,103].
The signal-to-noise ratio (SNR) characterizes the detectability of GWs signals at an interferometer. Taking the LISA detector as an example, we calculate the SNR as where Ω LISA is the sensitivity curve of the LISA detector [95], and T = 9.46 × 10 7 s the data-taking duration (around four years) [96]. According to Ref. [95], we use SNR > 10 (50) as the detection threshold for a six-link (four-link) configuration LISA. In the right panel of Fig. 1 we plot the (α, β/H) distribution as well as the SNRs of our data points. As shown 2 Note that v b is the bubble velocity with respect to the plasma at finite distance, while the velocity relevant for the EWBG calculation is actually vw, which is defined as the relative velocity to the plasma just in front of the wall. The relation between v b and vw can be solved using hydrodynamics [97,98], and it is possible to have a high v b (good for GW signals) and low vw (good for EWBG) simultaneously [12,[98][99][100][101].
in the figure, the data with large α (which means larger energy released in the transition) and smaller β/H (means longer duration of the transition) have larger SNRs. 3

Phenomenology at high energy muon colliders
Besides the GWs, the FOEWPT parameter space of the xSM can lead to signals of a resonantly produced heavy scalar (direct search), and corrections to the SM Higgs couplings (indirect search) at the colliders. The corresponding phenomenology has been studied at the LHC and the proposed pp or e + e − colliders [11,12,[99][100][101][109][110][111][112][113][114][115][116]. Typically, the direct search is implemented at pp colliders due to their high energy reach, while the indirect approach is preferred by the e + e − colliders because of the high accuracy. In this section we will demonstrate that, with the clean background and sufficient collision energy, a multi-TeV muon collider is able to perform both the direct and indirect searches, exhibiting a great potential to test the FOEWPT.

Production and decays of the heavy scalar
The heavy scalar h 2 can be produced at a lepton collider via the Zh 2 associated production or the vector boson fusion (VBF) process. At a collider with center-of-mass energy as high as a few TeV, the dominant channel is VBF 4 and the production rate is where the SM-like production rate σ SM h 2 is the SM Higgs VBF production cross section evaluated by replacing the Higgs mass with M h 2 . This is because the coupling of h 2 to the SM gauge bosons comes from mixing, see Eq. (2.4). In the left panel of Fig. 2, we have explicitly shown σ SM h 2 at the muon collider with different benchmark collision energies √ s. 5 From the figure, it is clear that we can easily obtain hundreds of fb of cross section for h 2 with O(TeV) mass. For a given √ s, the W + W − fusion contributes ∼ 90% of the total cross section. 3 There are some data points with α 1, implying a strong supercooling. In this case, it is suggested that it is the percolation temperature Tp rather than nucleation temperature Tn that should be used to calculate α and β/H [102,[104][105][106][107]. Since most of our data lie in the α 1 region, we adopt the approximation Tn ≈ Tp, and leave a more detailed treatment for the future work. Also note that we are using the traditional approach to derive FOEWPT profiles and calculate the GWs. It is shown that the alternative dimensional reduction approach can reduce the theoretical uncertainties significantly [108]. 4 The γh2 associated production (so-called "radiative return") can be comparable to the Zh2 production at a low energy muon collider [68]. 5 The cross sections are calculated by the MadGraph5aMC@NLO-v2.8.2 event generator [117] with the model file written using FeynRules [118]. The collider simulations in the next two subsections are also implemented with those two packages. The produced h 2 will subsequently decays to multiple final states, such as di-Higgs (h 1 h 1 ), di-boson (W + W − and ZZ) and di-fermion (e.g. tt). For the di-boson and difermion channels, where X denotes the SM vector boson or fermion, and Γ SM h 2 →XX is the decay width of the SM Higgs calculated at Higgs mass equal to M h 2 . For the di-Higgs channel, where the h 2 h 1 h 1 coupling is defined by and at tree level [112] The branching ratios are The branching ratios of the FOEWPT data points are projected to the Br(h 2 → h 1 h 1 )-Br(h 2 → V V ) plane in the right panel of Fig. 2, where V = W ± , Z. We see that the and Br(h 2 → tt) 20%. In the following two subsections, we choose the h 2 → h 1 h 1 → bbbb and h 2 → ZZ → + − + − as two complementary channels for collider simulations.

Direct search: the h
Directly characterizing the portal coupling between the singlet and the Higgs boson, the h 2 h 1 h 1 coupling is of our primary interests. The signal of such a coupling is a resonant di-Higgs production at the muon collider, VBF → h 2 → h 1 h 1 . As the SM Higgs dominantly decays into bb pairs, the major final state of the signal consists of four b-jets which can be reconstructed into two h 1 's and then one h 2 . For the ZZ fusion production channel, the final state contains two additional forward muons. Here we focus on the so-called inclusive channel by including both the W + W − fusion and ZZ fusion events without detecting the additional muons. In this case, the main backgrounds are the SM VBF h 1 h 1 and ZZ production, with h 1 → bb and Z → bb. 6 The signal and background events are generated at parton-level. We smear the jet four-momentum according to a jet energy resolution of ∆E/E = 10%, and assume a conservative b-tagging efficiency rate of 70%. The events are required to have exactly four b-jets satisfying the following basic acceptance cuts, where the pseudo-rapidity cut is based on a detector angular coverage of 10 • < θ < 170 • , and the recoil mass is defined as Next, we pair the four b-jets by minimizing 12) Figure 3. Left: after the basic acceptance cuts, the invariant mass distributions of the jet pairs and four-jet system for the signal and main backgrounds at the 10 TeV muon collider. Here we select M h2 = 600 GeV as the signal benchmark. Right: the expected probe limits on s 2 θ × Br(h 2 → h 1 h 1 ) for different muon collider setups. The scatter points are the FOEWPT data, in which red, green and blue colors represent SNR ∈ [50, +∞), [10,50) and [0, 10), respectively. The limit from ATLAS at the 13 TeV LHC with L = 36.1 fb −1 [119] and its extrapolation to the HL-LHC [12] are also shown for comparison.
The pairs (j 1 , j 2 ) and (j 3 , j 4 ) are then identified as the Higgs candidates, in which the harder pair is defined as (j 1 , j 2 ). As shown in blue in the left panel of Fig. 3, M j 1 j 2 and M j 3 j 4 peak around M h for the signal, while peak around M Z = 91.188 GeV for the ZZ background. Therefore, a invariant mass cut can significantly remove the ZZ background. While most SM h 1 h 1 events survive this cut, this background can be removed greatly by the cut on the four-jet system, as illustrated in orange in the left panel of Fig. 3. The cut flows for three chosen signal benchmarks at a 10 TeV muon collider are shown in Table 1, indicating Cut III is fairly powerful to improve the signal over background factor. Given the collision energy √ s and the integrated luminosity L, the signal and background event numbers are and make the interpolation to derive the s 2 θ ×Br(h 2 → h 1 h 1 ) reach as a function of M h 2 . The sensitivity of the muon collider to FOEWPT can be obtained by projecting the FOEWPT parameter space to such 2-dimension plane. This is done in the right panel of Fig. 3, in which the reach of different collider setups are plotted as different colored solid lines, and the FOEWPT data points lying above a specific line can be probed by the corresponding muon collider. Note that our projections are derived with a rather conservative b-tagging efficiency of 70%. A more optimistic efficiency such as 90% can improve the results by a factor of 2, while an analysis without b-tagging will weaken the limits by a factor of 2 or 3, as in this case the non-b jets (such as W ± /Z → jj) also contribute to the backgrounds. However, either case has little visual effects in the log coordinate.
The right panel of Fig. 3 demonstrates that the FOEWPT parameter space can be greatly probed by the muon colliders, and higher energy colliders (with also higher integrated luminosities) give better reach. The current and projected LHC reach is shown in black lines for comparison. Because of the high accuracy in the multi-jet final state, even a 3 TeV muon collider (1 ab −1 ) has a sensitivity more than one order of magnitude better than the HL-LHC (13 TeV, 3 ab −1 ), and a 30 TeV muon collider (90 ab −1 ) is able to probe s 2 θ × Br(h 2 → h 1 h 1 ) up to 10 −5 , covering almost all of the FOEWPT parameter space.
To manifest the complementarity with the GW experiments, we use different colors to mark the FOEWPT points with different SNRs: red, green and blue for SNR ∈ [50, +∞), [10,50) and [0, 10), respectively. Treating SNR = 10 as the detectable threshold, we see that those points which can be detected by LISA mostly lie in the reach of the muon colliders, especially for the √ s 6 TeV setups. This is a great opportunity to identify the origin of the stochastic GWs, if they were detected in the future. On the other hand, the muon colliders also have significant sensitivity to the blue data points which are not detectable at the LISA.
For muon colliders with √ s 10 TeV, there are still appreciable number of points which can not be reached, due to the tiny Br(h 2 → h 1 h 1 ) in those points. Hence, in the next subsection, we will change our strategy by looking for a complementary decay channel, h 2 → ZZ, to finally cover those points.

Direct search: the
The FOEWPT data points with tiny Br(h 2 → h 1 h 1 ) might be potentially probed via the h 2 → W + W − or h 2 → ZZ channels. To get a better accuracy we focus on the leptonic decay of the gauge bosons. Although the W + W − channel has a larger branching ratio, the neutrinos in the final state make this channel more challenging. In this subsection we would like to focus on the ZZ channel with Z → + − , where = e, µ, leading to a four-lepton final state.
Similar to the treatment in the di-Higgs channel in the previous subsection, three cuts are applied to the events. We first require the events to have exactly four charged leptons Here we select M h2 = 600 GeV as the signal benchmark. Right: the expected probe limits on s 2 θ × Br(h 2 → ZZ) for different muon collider setups. The scatter points are the FOEWPT data, in which red, green and blue colors represent SNR ∈ [50, +∞), [10,50) and [0, 10), respectively. The limit from ATLAS at the 13 TeV LHC with L = 36.1 fb −1 [120] and its extrapolation to the HL-LHC [12] are also shown for comparison.. with total zero charge, and satisfy the acceptance cuts p T > 30 GeV, |η | < 2.43, M recoil > 200 GeV, (Cut I) (3.17) where M recoil is defined similarly to Eq. (3.11). We then we pair the opposite-sign leptons by minimizing and put a second cut to the events. Note that this cut will select both the signal and background ZZ events. Finally, a cut for the four-lepton system can cut away the background significantly, as illustrated in the left panel of Fig. 4. The collider reach in ZZ channel can be shown as s 2 θ × Br(h 2 → ZZ) limits in the right panel of Fig. 4, in which different colored lines denote different muon collider setups, and the LHC current and future limits are also shown for comparison. Different from the case of h 1 h 1 channel, the reach of the HL-LHC is comparable to a 6 TeV muon collider (4 ab −1 ), thanks to the cleanness of the four-lepton final state. Nevertheless, better sensitivities can still be obtained for muon colliders with √ s 10 TeV, especially for a 30 TeV (90 ab −1 ) muon collider at which almost all the data points can be probed. We have checked that the FOEWPT data points escaping the h 1 h 1 channel search can be generally reached in the ZZ channel; as a result, even for a 3 TeV muon collider with an integrated luminosity of 1 ab −1 , the combination of h 1 h 1 and ZZ channels can cover almost the entire FOEWPT parameter space. The complementarity with GW experiments are also shown in the figure.

Indirect search: the Higgs coupling deviations
Besides the direct detection of the heavy scalar h 2 , measuring the couplings of the Higgslike boson h 1 also give hints of the FOEWPT, as those couplings usually deviate their corresponding SM values in the FOEWPT parameter space. For example, expanding the xSM Lagrangian as at tree level we obtain κ V = κ 3 = 1 for the SM, while for the xSM. Defining the deviations as we project the FOEWPT data points into the δκ 3 -δκ V plane in Fig. 5. One finds that δκ 3 is always positive (and 0.8). This can be understood by expanding the deviation at small mixing angle [12] where the M 2 h 2 /M 2 h term dominates the terms in the bracket, implying an enhanced Higgs triple coupling. Since we set θ 0.15 when scanning over the parameter space (see Appendix A), the δκ V distribution has a sharp edge at around 0.15 2 /2 ≈ 0.01.
Also shown in Fig. 5 are the projections of the reach for different setups of muon colliders. The corresponding probe limits are adopted from Ref. [78], which uses the VBF single Higgs production to study the h 1 V V coupling and the vector boson scattering di-Higgs production to study the triple Higgs coupling. It is clear that the FOEWPT parameter space can be probed very efficiently using via such indirect approach. A 3 TeV muon collider is already able to cover most of the data points, and a 30 TeV muon collider could test almost the whole parameter space.

Conclusion
FOEWPT is an important BSM phenomenon that might exist in the early Universe, and shed light on today's collider and GW experiments. In this article, we perform a complementarity study of the proposed high energy muon colliders and the near-future space-based GW detectors to the FOEWPT. Choosing the xSM as the benchmark model, we first derive the FOEWPT parameter space, and then test the possibility of detecting it via GW signals or muon collider experiments. In the calculation of FOEWPT, we have included the thermal tadpole term, which is dropped in a few previous references. It is shown that the inclusion of tadpole term reduces the fraction of "(0, 0) → (v f , v f s ) pattern" FOEWPT data points from 99% to 8.8%.
A considerable fraction of FOEWPT parameter space yields GW signals with SNR 10 at the LISA detector, thus might be probed. Since the TianQin and Taiji projects in China have projected sensitivities similar to the LISA, we expect the FOEWPT parameter space could be probed and crosschecked by those detectors in the near future.
For the muon colliders, we consider center of mass energies √ s = 3, 6, 10 and 30 TeV, with corresponding integrated luminosities L = 1, 4, 10 and 90 ab −1 , respectively. A detailed parton-level collider simulation is performed for the VBF production of the heavy real singlet and its subsequent decay to the di-Higgs (h 2 → h 1 h 1 → bbbb) and di-boson (h 2 → ZZ → + − + − ) channels. The results in Figs. 3 and 4 show that muon colliders offer great opportunity to probe the FOEWPT parameter space. In the di-Higgs channel, all muon collider setups have much higher reach compared to the HL-LHC; while in the di-boson channel, the reach of HL-LHC is comparable with a 6 TeV muon collider, however the muon colliders with √ s 10 TeV still work much better. Combining the di-Higgs and di-boson channels allows us to probe the whole parameter space. In addition, given the high accuracy of the muon colliders, the precision measurements of the Higgs gauge and triple couplings also help to test the FOEWPT.
As for the complementarity with the GW experiments, a remarkable result is that almost all parameter space yielding detectable GWs is within the reach of the muon colliders. This implies if in the future we really detected some signals at the GW detector, a muon collider could provide very useful crosscheck to locate their origin. We also find that there is large parameter space that is not detectable via GWs, but can be probed at the muon colliders.

Acknowledgments
We are grateful to Huai-Ke Guo and Ligong Bian for the very useful discussions and sharing the codes. KPX is supported by the Grant Korea NRF-2019R1C1C1010050.

A Deriving the phenomenologically allowed potential
This appendix demonstrates how to derive the parameter space of a xSM scalar potential that satisfies current phenomenological bounds. In summary, our strategy contains two steps: 1. Construct a potential in Eq. (2.1), and make sure it has a VEV in (h, s) = (v, 0), and a Higgs mass M h 1 = M h . In this case, generally b 1 = 0.
2. Shift the s field such that b 1 = 0, and match the new coefficients to the ones described in Eq. (2.7) of Section 2.1. In this case, generally v s = 0.
For the first step, an extremum at (v, 0) requires [121]  where the upper limits of a 2 , |b 3 | and b 4 come from the unitarity bound [115,121], while the lower limit of a 2 is required by a bounded below potential [121]. Note that Eq. (A.1) and Eq. (A.2) only guarantee (v, 0) is a local minimum, and there might be another deeper minimum. For a given set of Eq. (A.3), one needs to check whether (v, 0) is the vacuum (i.e. global minimum) by hand. We found that ∼ 38% of the sampling points yield (v, 0) as the vacuum. Such points are then phenomenologically allowed. For the second step, we shift s → s + σ to get the following redefinitions of the coefficients in Eq. (2.1) [8],