Neutral scalar pair production via W-boson fusion at multi–TeV muon colliders

Measurements of multi-scalar Higgs production constitute a central goal of upcoming collider experiments, including the High-Luminosity Large Hadron Collider (HL-LHC), prospective $ e^{+}e^{-} $ lepton colliders, and future multi–TeV muon colliders. Such precise measurements provide a promising opportunity to probe the structure of the Higgs scalar potential and to gain deeper insights into the underlying electroweak symmetry-breaking (EWSB) mechanism within the SM and its extensions. The ATLAS Collaboration has performed extensive searches for SM-like Higgs boson pair production in the $ b\bar{b}\gamma\gamma $ final state at a CoM energy of $ \sqrt{s}=13\; \text{TeV} $ in proton–proton ($ pp $) collisions at the LHC [1, 2]. Additional measurements of Higgs boson pair production in association with a vector boson (W or Z) at $ \sqrt{s}=13\; \text{TeV} $ have also been reported in Ref. [3]. The CMS Collaboration has investigated nonresonant Higgs boson pair production in the four-bottom-quark final state in $ pp $ collisions [4, 5]. Complementary studies of nonresonant Higgs pair production in the $ b\bar{b}\tau^-\tau^+ $ [6] and four-bottom-quark [7] final states have been carried out by ATLAS. A combined analysis across the $ b\bar{b}\gamma\gamma $, $ b\bar{b}\tau^-\tau^+ $, and $ b\bar{b}b\bar{b} $ channels is presented in Ref. [8]. Further studies of Higgs pair production have been reported in Refs. [9−15].

Phenomenological studies of Higgs boson pair production at colliders within the SM and its extensions have been extensively performed. In particular, di-Higgs production in the 2HDM has been investigated at the LHC in parameter regions consistent with the diphoton excess and the muon $ (g-2) $ anomaly [16]. Theoretical computations for Higgs boson pair production in the $ b\bar{b}\mu^+\mu^- $ final state at the LHC have been discussed in Ref. [17]. Moreover, probing the electroweak sector via the process $ W_L W_L \to n\times h $ has been examined within various effective field theory (EFT) frameworks [18]. Theoretical predictions for di-Higgs production in the $ b\bar{b}\ell^-\ell^+ $ final state with missing transverse momentum at the HL-LHC have been presented in Ref. [19]. Furthermore, double Higgs production via W-boson fusion at TeV-scale $ e^+e^- $ colliders has been analyzed within EFT frameworks, including studies of sensitivity to BSM Higgs couplings [20], and multi-Higgs measurements in the SMEFT have been examined in [21]. Additional analyses of Higgs pair production have been considered in EFTs [22] as well as in extended Higgs sector models [23]. Di-Higgs production has also been investigated in the Higgs singlet extension [24] and the inert doublet model with two active doublets [25]. Further studies of di-Higgs production in various BSM frameworks, such as the SMEFT [26], heavy neutral lepton extensions of the SM at future lepton colliders [27], and scenarios involving axion-like particles [28], have also been carried out. The potential of probing new physics through di-Higgs production at the LHC has been explored in Ref. [29]. Higher-order corrections to Higgs boson pair production have been widely investigated. In particular, QCD corrections to Higgs boson pair production, including decays into the $ b\bar{b}\tau^+\tau^- $ final state, have been reported in Ref. [30]. In addition, double-logarithmic enhancements in Higgs boson pair production in the high-energy limit have been studied in Ref. [31]. Next-to-leading-order (NLO) electroweak corrections have been presented in Refs. [32, 33], while higher-order contributions within the Higgs Effective Field Theory (HEFT) framework have been investigated in Refs. [34, 35], including top-Yukawa- and light-quark-induced electroweak effects [36]. Furthermore, NLO QCD and electroweak corrections to double Higgs production have been computed in Refs. [37−50].

Future lepton colliders (LCs) provide a cleaner environment than hadron colliders, which are subject to large QCD backgrounds, enabling higher-precision measurements. Proposed multi–TeV muon colliders further access a higher-energy regime, allowing direct probes of new physics. Moreover, since scalar particles couple strongly to vector bosons, vector-boson fusion processes provide excellent prospects for probing double scalar production. Phenomenological studies of vector-boson-fusion processes within the SM at the LHC and at future multi–TeV muon colliders can be found in Refs. [67−71]. To date, no calculations have included full one-loop electroweak corrections for double scalar production via W-boson fusion. In this work, we present the first evaluation of full one-loop electroweak radiative corrections to the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ in the SM. We then compute neutral scalar pair production via W-boson fusion at multi–TeV muon colliders in the 2HDM and scan the corresponding cross sections across the updated parameter space of Type-X and Type-Y 2HDMs. The enhancement factor, defined as the ratio of the cross section for SM-like Higgs pair production in the 2HDM to the corresponding SM prediction, is analyzed across the viable regions of parameter space of Type-X and Type-Y 2HDMs. Finally, both CP-even and CP-odd Higgs pair production is investigated across the allowed parameter space of the considered models. In the Type-Y scenario, at $ \sqrt{s} = 10 $ TeV with an integrated luminosity of $ {\cal{L}} = 10000\; \text{fb}^{-1} $, both CP-even and CP-odd Higgs pair production in the $ t\bar{t}b\bar{b} $ final state, including subsequent top-quark decays into leptons and bottom quarks, can be probed with a statistical significance exceeding $ 2\sigma $ at several viable parameter points.

The outline of this paper is as follows. In Section 2, we present the 2HDM and discuss the relevant theoretical and experimental constraints. Full one-loop electroweak radiative corrections to the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ in the SM are computed in detail in Section 3. Phenomenological analyses of neutral scalar pair production via W-boson fusion at multi–TeV muon colliders are presented in Section 4. Finally, conclusions and outlook are provided in Section 5, while checks of the calculations are reported in the appendices.

In this section, we present the structure of the models under investigation along with the corresponding theoretical and experimental constraints. For further reviews of the 2HDM, including phenomenological studies, we refer the reader to Refs. [51−55]. In particular, the gauge and fermion content of the 2HDM remains identical to that of the SM. In contrast, the scalar Higgs sector is extended by introducing an additional complex scalar field with hypercharge $ Y = 1/2 $. With this extension, the most general scalar Higgs potential consistent with gauge invariance and renormalizability is expressed as follows (adopting the notation of Ref. [53]):

Following our previous work [51−54], we adopt the CP-conserving version of the 2HDM in our phenomenological analysis. Accordingly, all parameters $ m_{11}^2 $, $ m_{22}^2 $, $ m_{12}^2 $, and $ \lambda_1, \ldots, \lambda_5 $ in the potential above are taken to be real. A discrete $ Z_2 $ symmetry is imposed on the scalar potential, with only a soft-breaking term allowed, thereby forbidding tree-level flavor-changing neutral currents. Consequently, four distinct types of 2HDMs are classified according to their Yukawa coupling structures. Following Ref. [51], the general form of the Yukawa Lagrangian reads:

where v denotes the vacuum expectation value. The Cabibbo–Kobayashi–Maskawa (CKM) matrix, with elements $ V_{ij} $, is included in the above Lagrangian to describe quark mixing. The operators $ P_{L/R} = \frac{1 \mp \gamma_5}{2} $ are projection operators, and the left- and right-handed leptons are denoted by $ \ell_{L/R} $. The coupling coefficients appearing in Eq. 2 are listed in Table 2 of Ref. [56] or in Table 1.

After EWSB, the additional scalar particles in the model include a CP-even Higgs boson (H), a CP-odd Higgs boson (A), and a pair of charged Higgs bosons ($ H^\pm $). The independent parameters of the 2HDM employed in our analysis comprise the mixing parameters $ s_{\beta-\alpha} $ and $ \tan\beta $, the scalar masses $ m_H $, $ m_A $, and $ m_{H^\pm} $, as well as the soft-breaking parameter $ m_{12}^2 $. Before proceeding to the study of neutral scalar pair production at multi–TeV muon colliders, we briefly outline the updated parameter space of the model. The six independent parameters of the 2HDM are subject to both theoretical and experimental constraints. For a detailed discussion of these constraints and the corresponding results, we refer the reader to our previous works [51−53]. The parameter scan is performed over the following ranges: $ m_{h}=125.09\; {\rm{GeV}} $, $ m_{H} \in [130, 1000]\; {\rm{GeV}} $, $ s_{\beta-\alpha} \in [0.97, 1] $, $ m_{A,H^{\pm}} \in [130, 1000]\; {\rm{GeV}} $, $ \tan\beta \in [0.5, 45] $, and $ m_{12}^2 \in [0, 10^6]\; {\rm{GeV}}^2 $. Theoretical constraints are first imposed on the scanned parameter space, including perturbative unitarity, perturbativity, and vacuum stability. The resulting parameter points are then tested against electroweak precision observables (EWPOs). In this context, the electroweak oblique parameters S, T, and $ U\; $[57], which are related to the W-boson mass [66], are taken into account. The allowed parameter space is obtained by requiring S, T, and U to lie within the $ 95$% confidence level limits. All theoretical and EWPO constraints are implemented using the public code $ {\mathrm{2HDMC-1.8.0}}$ [58]. In addition, compatibility with exclusion limits from Higgs searches at LEP, the LHC, and the Tevatron is ensured using $ {\mathrm{HiggsBounds}}$-5.10.1 [59], which compares theoretical predictions with experimental upper bounds at the $ 95$% confidence level. The consistency of the SM-like Higgs boson with LHC measurements is further tested using $ {\mathrm{HiggsSignals}}$-2.6.1 [60]. Finally, constraints from flavor physics are taken into consideration, particularly those associated with B-meson observables measured by LHCb. These constraints are evaluated using the $ {\mathrm{SuperISO}}$ v4.1 package [61], and only parameter points consistent with experimental data within the $ 2\sigma $ bounds are retained.

Based on the updated parameter-space datasets provided in Refs. [51−53], we present representative scatter plots in the subsequent analysis. In Fig. 1, the parameter space of the six independent parameters is shown for the Type-X 2HDM (left) and the Type-Y 2HDM (right). The upper panels display scatter plots of $ m_A - m_H $, $ m_{H^\pm} - m_H $, and $ m_{H^\pm} - m_A $, while the lower panels show the distributions of $ m_{12}^2 $, $ m_A $, and $ \tan\beta $. The results in the upper panels indicate that the mass splittings among the scalar states are typically constrained to be below approximately $ 200 $ GeV by EWPO bounds. In the Type-X scenario, the mass difference between the charged Higgs boson and the CP-even Higgs boson can reach values up to $ \pm 600 $ GeV. Nevertheless, the dominant region in the $ m_{H^\pm}-m_A $ plane is concentrated along $ m_{H^\pm} \simeq m_A $. Similar features are observed in the Type-Y 2HDM. In this case, the allowed range of the mass difference is $ -600\; \text{GeV} \leqslant m_A - m_H \leqslant 400\; \text{GeV} $, which is narrower than the corresponding interval $ -600\; \text{GeV} \leqslant m_A - m_H \leqslant 600\; \text{GeV} $ found in the Type-X scenario. In the lower panels, we present the distributions of $ m_{12}^2 $, $ m_A $, and $ \tan\beta $ for both 2HDM types. Over the full range of CP-odd Higgs masses, the soft $ Z_2 $-breaking parameter lies within $ 3 \cdot 10^{3} \leqslant m_{12}^2 \leqslant 5 \cdot 10^{5}\; {\rm{GeV}}^2 $. In addition, the parameter $ \tan\beta $ is found to lie in the range $ 1 \leqslant \tan\beta \leqslant 20 $. We observe distinct differences in the viable parameter spaces of the Type-X and Type-Y 2HDMs for several reasons. First, these differences arise from constraints obtained by confronting the model predictions with experimental data, including high-precision measurements of the SM-like Higgs boson properties as well as exclusion limits from scalar boson searches. These measurements impose stringent restrictions on the masses of the CP-even and CP-odd scalar states. Second, additional constraints originate from flavor-physics observables, which place strong bounds on the charged Higgs boson mass and the mixing parameter $ \tan\beta $.

Figure 1.  (color online) The parameter space for the six independent parameters in the Type-X 2HDM (left) and Type-Y 2HDM (right). The upper panels show scatter plots of $ m_A - m_H $, $ m_{H^\pm} - m_H $, and $ m_{H^\pm} - m_A $, while the lower panels display scatter plots of $ m_{12}^2 $, $ m_A $, and the mixing angle $ \tan\beta $.

In this section, we present, for the first time, the full one-loop electroweak radiative corrections to the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ within the Standard Model (SM). The calculations are performed using the $ {\mathrm{GRACE-Loop}}$ system, which is described in detail in Ref. [62]. In general, the total cross section for the process $ \mu^- \mu^+ \to VV \to hh $ can be expressed as follows:

Here, $ \hat{\sigma}_{V_{\lambda_1}V_{\lambda_2} \to hh}(\tau s) $ denotes the cross section for the partonic process $ V_{\lambda_1}V_{\lambda_2} \to hh $, where $ \lambda_1 $ and $ \lambda_2 $ denote the polarization states of the vector bosons. In this work, the partonic processes are computed including one-loop radiative corrections with the $ {\mathrm{GRACE-Loop}}$ system. Moreover, the vector-boson splitting functions, which describe the probability for an initial-state muon to emit a weak vector boson, correspond to the helicity-dependent parton distribution functions (PDFs) [67−69] and are given by:

Here, $ Q^2 $ denotes the energy scale. We emphasize that the expression in Eq. 4 corresponds to the transverse polarization of the vector-boson PDF, whereas Eq. 5 describes the longitudinal polarization of the vector-boson PDF. The coupling coefficient $ g_V $ is defined as follows:

In this work, we consider several partonic subprocesses, namely $ \gamma\gamma \to hh $, $ \gamma Z \to hh $, $ ZZ \to hh $, and $ W^\pm W^\mp \to hh $. However, the first two channels are loop-induced and are therefore expected to yield smaller contributions than the latter processes. The $ ZZ $-fusion process proceeds as $ \mu^- \mu^+ \to ZZ \to \mu^- \mu^+ hh $, while the $ WW $-fusion channel is given by $ \mu^- \mu^+ \to W^\pm W^\mp \to \nu_{\mu}\bar{\nu}_{\mu} hh $. These processes lead to distinct final-state signatures. In the $ ZZ $-fusion case, di-Higgs production is accompanied by a muon pair, whereas in the $ WW $-fusion channel, the signal is associated with missing energy due to neutrinos in the final state. From kinematic considerations, the $ ZZ $-fusion contribution is significantly smaller than that of the $ WW $-fusion process. For these reasons, we focus on the calculation of the partonic process $ W^\pm W^\mp \to hh $. The relevant partonic processes are given by

In the $ {\mathrm{GRACE-Loop}}$ program, nonlinear gauge-fixing terms are implemented; see Ref. [62] for details. Consequently, the one-loop amplitude depends on the nonlinear gauge parameters as well as on the ultraviolet (UV) divergence parameter, $ C_{UV} $. This dependence cancels once all contributions from the one-loop and corresponding counterterm diagrams are included. For the process $ W^\pm W^\mp \to hh $, virtual photon exchange in the loop induces infrared (IR) divergences, leading to a dependence on the photon-mass regulator λ. This dependence is removed after including soft-photon radiation (the third term on the right-hand side of Eq. 8). However, the results still depend on the soft–hard separator $ k_c $. To obtain the complete one-loop radiative corrections, hard-photon emission must also be included, corresponding to the process $ W^\pm W^\mp \to hh\gamma $ with a hard photon in the final state (which corresponds to the last term on the right-hand side of Eq. 8). The final results are independent of all the aforementioned unphysical parameters. For validation of the calculation, we refer the reader to Appendix A, where explicit results are presented.

Using the corrected partonic cross sections, we evaluate the total cross section by convolving them with the W-boson PDFs, as defined in Eq. 3. We then investigate the impact of the full one-loop electroweak radiative corrections on the processes under consideration. Together with electroweak radiative corrections, initial-state radiation (ISR) effects are also important for the investigated processes at future colliders and are therefore included in this study. The ISR effects are evaluated as follows:

where s denotes the center-of-mass (CoM) energy, and x represents the energy fraction of a photon emitted by the initial-state muon. The radiator function, taken from Ref. [73], is given by

In the formulas, we use the relevant coefficients given by

In the subsequent phenomenological analysis, we use input parameters from the Particle Data Group [66]. We work in the $ G_{\mu} $ scheme, where the fine-structure constant is calculated as

with $ G_F=1.1663785\cdot 10^{-5} $ GeV^-2. We also apply the selection cuts $ p_T^{h}\geq 20 $ GeV and pseudorapidity $ |\eta_{h}|<2.4 $ for the following results.

In Fig. 2, the cross sections for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ are shown as functions of the center-of-mass (CoM) energy, ranging from 3 TeV to $ 30 $ TeV. The red curve corresponds to the tree-level prediction, while the blue curve represents the fully corrected result. The $ {\cal{O}}(\alpha^2) $ ISR corrections are represented by the green curve. We find that the electroweak corrections range from approximately $ 10$% to $ 5$% over this energy range, while the $ {\cal{O}}(\alpha^2) $ ISR corrections contribute at the level of about $ -5$%.

Figure 2.  (color online) The cross sections for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $, including the full electroweak radiative corrections, are shown as a function of the center-of-mass (CoM) energy.

The differential cross sections with respect to $ p_T^{h} $ and the rapidity $ \eta^{h} $ for the SM-like Higgs boson are shown in Fig. 3. In these plots, the red curve corresponds to the tree-level predictions, while the blue (green) curve represents the fully corrected (ISR-only) results. The left panels display the distributions at $ \sqrt{s} = 3 $ TeV, whereas the right panels show the corresponding results at $ \sqrt{s} = 10 $ TeV. We find that the full electroweak corrections span approximately from $ +10$% down to $ \sim -30$%, whereas the ISR corrections range from about $ -2$% to $ \sim -3$%. In the high-$ p_T $ tail, the cross sections become very small, so the corresponding corrections are not phenomenologically significant. In general, the electroweak corrections are of the order of $ 10$% for the distributions considered. A similar pattern is observed for the rapidity distributions.

Figure 3.  (color online) The differential cross sections as functions of the transverse momentum $ p_T^{h} $ and rapidity $ \eta^h $ are shown. The red line corresponds to the tree-level cross sections, while the blue line represents the fully corrected cross sections.

We have verified that the longitudinal modes $ \mu^- \mu^+ \to W^\pm_L W^\mp_L \to hh $ provide the dominant contribution in the high-energy regime of future multi–TeV colliders. This observation is consistent with the numerical results of Ref. [68], where the longitudinal channel accounts for more than $ 97$% of the total contribution. Consequently, the $ Q^2 $ dependence of the cross section is relatively weak, typically at the level of a few percent, as shown in Table B1. Similar results were also reported in Ref. [68]. Finally, the cross sections presented in Fig. 2 increase with the CoM energy, in accordance with the Goldstone boson equivalence theorem. In particular, the cross sections for $ W^\pm_L W^\mp_L \to hh $ are equivalent to those of $ G^\pm G^\mp \to hh $ in the high-energy limit, leading to a growth of the cross sections with $ \sqrt{s} $. The mild scale dependence, at the level of a few percent, further underscores the importance of the electroweak radiative corrections and the $ {\cal{O}}(\alpha^2) $ ISR effects evaluated in this work.

We evaluate neutral scalar pair production via W-boson fusion at multi–TeV muon colliders within the 2HDM framework. In this study, we restrict the calculations to tree-level cross sections for neutral scalar pair production. The corresponding amplitudes are generated using the $ {\mathrm{FeynArts/FormCalc}}$ packages [64, 65]. To validate our results, the amplitudes are computed in the general $ R_\xi $ gauge and tested for gauge invariance by varying the gauge-fixing parameter ξ. Having established gauge invariance, we proceed to analyze the phenomenological results in the following subsections. It should be noted that the $ {\cal{O}}(\alpha^2) $ ISR effects are incorporated into the tree-level cross sections for neutral scalar pair production in the following phenomenological analyses.

We first consider SM-like Higgs-pair production via W-boson fusion at multi–TeV muon colliders. For both the signal and background processes, the total cross sections for $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ are computed using Eq. 3. In our phenomenological analysis, we introduce an enhancement factor, defined as follows:

where $ {\cal{P}}_{ \rm{2HDM}} = \{ s_{\beta-\alpha}, \tan\beta, m_H, m_A, m_{H^\pm}, m_{12}^2 \} $. It should be stressed that we use tree-level cross sections for SM-like Higgs pair production in the 2HDM ($ \sigma^{ \rm{LO}}_{ \rm{2HDM} } $), whereas the corresponding cross sections in the SM include full one-loop corrections ($ \sigma^{ \rm{NLO}}_{ \rm{SM} } $). We study the enhancement factor over the viable parameter space of the 2HDM. In Fig. 4, we scan the enhancement factor $ \mu_{hh} $ over $ \tan\beta $ and $ m_H $ for the Type-X 2HDM (left panels) and the Type-Y 2HDM (right panels). The results are presented at $ \sqrt{s} = 3 $ TeV in the upper panels and at $ \sqrt{s} = 10 $ TeV in the lower panels. They indicate that the enhancement factor can reach values up to approximately $ 3 $ for several viable parameter points in the Type-X 2HDM, whereas it remains below unity in the Type-Y 2HDM. For the Type-X scenario, the enhancement factor is close to unity for $ \tan\beta \lesssim 5 $ and $ m_H \gtrsim 400 $ GeV. In contrast, in the Type-Y 2HDM, it lies within the range $ 0.91 $–$ 0.95 $. In general, this process proceeds identically in both the Type-X and Type-Y 2HDMs, as it does not involve scalar–fermion couplings at leading order. The observed differences therefore originate primarily from the distinct allowed parameter spaces, which are strongly constrained by the Yukawa interactions of the scalar sector. Furthermore, although the scalar trilinear couplings relevant to this process have identical expressions in both Type-X and Type-Y 2HDMs, the allowed parameter space in the $ (\tan\beta,\, m_H) $ plane differs significantly between the two scenarios. As a result, distinct behaviors of the enhancement factor are observed in the Type-X and Type-Y 2HDMs. These results demonstrate the potential to distinguish between these two types of 2HDMs via precision measurements of double Higgs production at future multi–TeV colliders.

Figure 4.  (color online) We evaluate the enhancement factor $ \mu_{hh} $ across the $ (\tan\beta,\, m_H) $ parameter space for the Type-X 2HDM (left panels) and the Type-Y 2HDM (right panels). Results are shown for $ \sqrt{s} = 3 $ TeV (upper panels) and $ \sqrt{s} = 10 $ TeV (lower panels).

In this subsection, we compute the cross section for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to AA \to t\bar{t} b\bar{b} $ in our analysis. As discussed below, we consider the decay channels of the CP-odd Higgs boson A, namely $ A \to t\bar{t} $ and $ A \to b\bar{b} $. The cross section for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to A(p_3) A(p_4) \to t\bar{t} b\bar{b} $ is evaluated as follows:

Here, $ \Gamma_{A\to b\bar{b}} $, $ \Gamma_{A\to t\bar{t}} $, and $ \Gamma_{A} $ denote the partial decay widths of the CP-odd Higgs boson into $ b\bar{b} $ and $ t\bar{t} $, and its total decay width, respectively. All decay widths mentioned above are computed using the publicly available code $ {\mathrm{H-COUP}}$ [63]. Note that a factor of $ 2 $ is included to account for the interchange of the decay modes $ A \to b\bar{b} $ and $ A \to t\bar{t} $, whereas an additional factor of $ 1/2 $ arises from the presence of two identical CP-odd Higgs bosons in the final state of the partonic process $ W^\pm W^\mp \to AA $. We also emphasize that off-shell CP-odd Higgs bosons (i.e., $ p_3^2 \neq m_A^2 $ and $ p_4^2 \neq m_A^2 $) are included in this procedure. The partonic process is first generated in the general $ {\cal{R}}_{\xi} $ gauge with the help of the computer algebra packages $ {\mathrm{FeynArts/FormCalc}}$ [64, 65]. The results are then checked for consistency by varying the ξ gauge parameters. The corrected amplitude for the process under consideration is convolved as in Eq. 18. We emphasize that the decay width of the CP-odd Higgs boson satisfies $ \Gamma_A / m_A \sim {\cal{O}}(5 \times 10^{-2}) $ (i.e., below the $ 5$% level). Therefore, Eq. 18 can be approximated as follows:

This approximation is used throughout our analysis.

The branching ratios for $ A \to t\bar{t} $ (left panels) and $ A \to b\bar{b} $ (right panels) are shown in the $ (m_A,\, \tan\beta) $ plane. The upper (lower) panels correspond to the Type-X (Type-Y) 2HDM. Across the parameter space, the decay $ A \to t\bar{t} $ dominates over $ A \to b\bar{b} $. We also find that the $ A \to b\bar{b} $ branching ratio in the Type-Y 2HDM is more than two orders of magnitude larger than in the Type-X 2HDM. This difference arises from the distinct Yukawa coupling structures of the two models.

Figure 5.  (color online) The branching ratios for the decays $ A \to t\bar{t} $ (left panels) and $ A \to b\bar{b} $ (right panels) are shown in the $ (m_A,\, \tan\beta) $ parameter space. The upper panels correspond to the Type-X 2HDM, while the lower panels correspond to the Type-Y 2HDM.

We now discuss event generation for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to A(p_3) A(p_4) \to t\bar{t} b\bar{b} $, scanning the parameter spaces of the Type-X and Type-Y 2HDMs. In Fig. 6, we present the expected event yields at $ \sqrt{s} = 3 $ TeV with an integrated luminosity of $ {\cal{L}} = 3000\; \text{fb}^{-1} $ and at $ \sqrt{s} = 10 $ TeV with $ {\cal{L}} = 10000\; \text{fb}^{-1} $. The results for the Type-X 2HDM are shown in the upper panels, while those for the Type-Y 2HDM are displayed in the lower panels. The events are shown in the $ (m_A,\, \tan\beta) $ parameter space. The SM background is indicated by red dotted lines and includes the processes $ \mu^- \mu^+ \to W^\pm W^\mp \to t\bar{t} b\bar{b} $ and $ t\bar{t}jj $, with light jets $ j=u,d,c,s $. The $ t\bar{t}jj $ channel accounts for the misidentification of light jets as b-jets. We compute the full set of Feynman diagrams contributing to these final states. Consequently, the SM background includes all relevant subprocesses, such as $ \mu^- \mu^+ \to W^\pm W^\mp \to t\bar{t} P $, where $ P = \gamma, Z, h, g $, etc., with subsequent decays $ P \to b\bar{b} $ and $ P \to j\bar{j} $. The background cross sections are evaluated after applying the kinematic selections described below. To reconstruct the CP-odd Higgs bosons, invariant-mass cuts are imposed on the final-state particles: $ |m_{t\bar{t}} - m_A | \leqslant 20\; \text{GeV} $ and $ | m_{jj} - m_A | \leqslant 20\; \text{GeV} $. These cuts efficiently suppress the SM background, as it is evaluated away from the dominant resonance regions associated with Z-boson mediation. In addition, the following kinematic cuts are applied to the jets: $ p_T^{j} \geq 60 $ GeV, $ |\eta^{j}| \leqslant 2.4 $, and $ R_{jj} \gt 0.5 $, for both signal and background events. The same selection cuts are applied to b-jets. In all cases, CP-odd Higgs boson pair production in the Type-Y 2HDM at a $ \sqrt{s} = 10 $ TeV muon collider represents a promising channel for future facilities. In this scenario, the SM background is expected to be relatively small and is therefore neglected for simplicity. We consider the leptonic decay mode of the top quark, $ t \to \ell \nu_\ell b $, with a total branching fraction of $ 0.332 $ [66]. The resulting final state is $ \ell \bar{\ell}\, b\bar{b}\, b\bar{b} $, accompanied by missing energy. The statistical significance for the process $ \mu^- \mu^+ \to W^\pm W^\mp \to AA \to t\bar{t} b\bar{b} \to \ell \bar{\ell}\, b\bar{b}\, b\bar{b} $, including missing energy, is estimated to be

Figure 6.  (color online) Production events at $ \sqrt{s} = 3 $ TeV with an integrated luminosity of $ {\cal{L}} = 3000\; \text{fb}^{-1} $ and at $ \sqrt{s} = 10 $ TeV with $ {\cal{L}} = 10000\; \text{fb}^{-1} $ are presented. The results for the Type-X 2HDM are shown in the upper panels, while those for the Type-Y 2HDM are displayed in the lower panels.

In Fig. 7, we scan the statistical significance across the viable parameter space as a function of the enhancement factor $ \mu_{hh} $. The signal exceeds $ 2\sigma $ at several allowed points in the Type-Y 2HDM; however, most viable points remain below the $ 5\sigma $ discovery threshold.

Figure 7.  (color online) The statistical significance of the signal process $ \mu^- \mu^+ \to W^\pm W^\mp \to AA \to t\bar{t} b\bar{b} \to \ell \bar{\ell}\, b\bar{b}\, b\bar{b} $, with missing energy included, is evaluated across the viable parameter space in relation to the enhancement factor $ \mu_{hh} $ for the Type-Y 2HDM at $ \sqrt{s} = 10 $ TeV.

We also investigate the pair production of CP-even heavy Higgs bosons at future multi–TeV muon colliders via the process $ \mu^- \mu^+ \to W^\pm W^\mp \to HH \to t\bar{t} b\bar{b} $. The production cross sections are computed using Eq. 18, with the corresponding CP-even Higgs mass and decay width. In the narrow-width limit $ \Gamma_H / m_H \ll 1 $, the production cross sections can also be approximated as follows:

This implies that, to apply the above approximation, the branching fractions for the CP-even Higgs boson decays into $ t\bar{t} $ and $ b\bar{b} $ must first be evaluated. These quantities are computed using the publicly available code $ {\mathrm{H-COUP}}$. In Fig. 8, we present the corresponding branching fractions: $ t\bar{t} $ (left panels) and $ b\bar{b} $ (right panels). The results for the Type-X 2HDM are shown in the upper panels, while those for the Type-Y 2HDM are displayed in the lower panels. Across the parameter space, the decay mode $ H \to t\bar{t} $ dominates. The branching ratio for $ H \to b\bar{b} $ ranges from $ 10^{-4} $ to $ 10^{-1} $ in the Type-Y 2HDM and from $ 10^{-5} $ to $ 10^{-2} $ in the Type-X 2HDM.

Figure 8.  (color online) The branching ratios for the decays $ H \to t\bar{t} $ (left panels) and $ H \to b\bar{b} $ (right panels) are shown in the $ (m_H,\, \tan\beta) $ plane. The Type-X 2HDM results are presented in the upper panels, while the Type-Y 2HDM results are shown in the lower panels.

Expected event yields for the considered process at $ \sqrt{s} = 3 $ TeV with an integrated luminosity of $ {\cal{L}} = 3000\; \text{fb}^{-1} $ and at $ \sqrt{s} = 10 $ TeV with $ {\cal{L}} = 10000\; \text{fb}^{-1} $ are shown in Fig. 9. The results for the Type-X 2HDM are presented in the upper panels, while those for the Type-Y 2HDM are displayed in the lower panels. In these plots, the red dotted lines represent the Standard Model background. We find that, for this channel at $ \sqrt{s} = 10 $ TeV with $ {\cal{L}} = 10000\; \text{fb}^{-1} $ in the Type-Y 2HDM, the expected yields offer a particularly promising scenario to be probed at future multi–TeV muon colliders. Accordingly, we further investigate the corresponding signal significance.

Figure 9.  (color online) Production event yields for the considered process at $ \sqrt{s} = 3 $ TeV with an integrated luminosity of $ {\cal{L}} = 3000\; \text{fb}^{-1} $ and at $ \sqrt{s} = 10 $ TeV with $ {\cal{L}} = 10000\; \text{fb}^{-1} $ are presented. The results for the Type-X 2HDM are shown in the upper panels, while those for the Type-Y 2HDM are displayed in the lower panels.

As in the previous case, we scan the viable parameter space of the Type-Y 2HDM at $ \sqrt{s} = 10 $ TeV and evaluate the statistical significance of the signal process $ \mu^- \mu^+ \to W^\pm W^\mp \to HH \to t\bar{t} b\bar{b} \to \ell \bar{\ell}\, b\bar{b}\, b\bar{b} $, including the effects of missing energy, as a function of the enhancement factor $ \mu_{hh} $. The results are shown in Fig. 10. We find that the statistical significance can exceed the $ 2\sigma $ level for several parameter points within the allowed parameter space of the Type-Y 2HDM.

Figure 10.  (color online) We evaluate the statistical significance of the signal process $ \mu^- \mu^+ \to W^\pm W^\mp \to HH \to t\bar{t} b\bar{b} \to $$ \ell \bar{\ell}\, b\bar{b}\, b\bar{b} $, with missing energy included, across the viable parameter space and examine its correlation with the enhancement factor $ \mu_{hh} $ for the Type-Y 2HDM at $ \sqrt{s} = 10 $ TeV.

Other double-scalar production channels, such as $ \mu^- \mu^+ \to W^\pm W^\mp \to AH $, $ hH $, and $ hA $, are suppressed by the softly broken $ Z_2 $ symmetry and are therefore not considered in this work.

In this work, we present the first results for the full one-loop electroweak radiative corrections to the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ in the SM. We then evaluate neutral-scalar pair production via vector-boson fusion at multi–TeV muon colliders within the 2HDM framework. In the phenomenological analysis, we scan the viable parameter space of the Type-X and Type-Y 2HDMs and evaluate the enhancement factor for SM-like Higgs-pair production relative to the SM prediction. We find that this factor can reach values up to $ 3 $ in several regions of the Type-X 2HDM parameter space, whereas it remains in the range $ 0.91 $–$ 0.95 $ throughout the entire parameter space of the Type-Y 2HDM. It is noteworthy that the enhancement factor exhibits distinct behavior in the Type-X and Type-Y 2HDMs. These results demonstrate the potential to distinguish between these two types of 2HDMs via precision measurements of double-Higgs production at future multi–TeV colliders. Finally, we investigate CP-even and CP-odd Higgs-pair production via W-boson fusion at multi–TeV muon colliders. In the Type-Y scenario, at $ \sqrt{s} = 10 $ TeV with an integrated luminosity of $ {\cal{L}} = 10000\; \text{fb}^{-1} $, both CP-even and CP-odd Higgs-pair production in the $ t\bar{t}b\bar{b} $ final state, including subsequent top-quark decays into leptons and bottom quarks, can be probed with a statistical significance exceeding $ 2\sigma $ at several viable parameter points. These results highlight the strong potential of future multi–TeV muon colliders to probe the scalar sector and provide deeper insight into the mechanism of EWSB.

In this Appendix, we present numerical consistency checks for the one-loop radiative corrections to the process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ in the SM. The $ {\mathrm{GRACE-Loop}}$ system implements nonlinear gauge-fixing (NLG) terms in the Lagrangian. As a result, the total cross section must be independent of the NLG parameters, namely α, β, $ \cdots $, and κ. In the context of one-loop renormalization, the total cross section is also free of UV divergences. Due to the presence of virtual photon exchange in the loop, the corresponding Feynman diagrams contain infrared divergences. To regularize these divergences, a fictitious photon mass λ is introduced for the virtual photon in the loop. Consequently, the one-loop cross sections depend on λ. By including soft-photon radiation, the cross sections become independent of λ (the photon mass). However, the inclusion of soft-photon radiation makes the cross section dependent on the photon energy cutoff parameter, denoted by $ k_c $. At this stage, the dependence on $ k_c $ is canceled by the contribution from hard-photon radiation. In the final result, the fully corrected cross sections are free of all the aforementioned parameters.

In Table 1, we present numerical checks of the NLG parameters α, β, $ \cdots $, and κ, as well as $ C_{UV} $ and λ. The results show that the cross sections are stable when these parameters are varied over wide ranges. We impose $ p_T^{h} \geq 20 $ GeV and pseudorapidity $ |\eta_{h}| \lt 2.4 $. The tests are performed at a CoM energy of $ 3 $ TeV at a future multi–TeV muon collider.

Table A1 presents tests of the stability of the results with respect to $ k_c $. In the first column, we vary $ k_c $ from $ 10^{-3} $ GeV to $ 10^{-5} $ GeV. The second column lists the soft-photon contributions, and the third column lists the hard-photon contributions. The last column gives their sum. All cross sections are in pb. We set NLG$ =\{0,0,0,0,0,0\} $ and $ \lambda=10^{-15} $ GeV for this test. The results show good stability with respect to variations of the photon-energy cutoff $ k_c $.

The tree-level process $ \mu^- \mu^+ \to W^\pm W^\mp \to hh $ in the SM has been extensively studied in the literature. We cross-check the tree-level cross sections computed in this work with results generated by $ {\mathrm{MadGraph5}}\_{\mathrm{aMC@NLO}}$ [68, 72]. In Table B1, the first column shows the center-of-mass (CoM) energies ranging from $ 4 $ TeV to $ 30 $ TeV. The second column presents our results, while the third column lists those from $ {\mathrm{MadGraph5}}\_{\mathrm{aMC@NLO}}$ [68, 72]. The fourth column displays the relative difference between the two results, expressed in percent. We observe good agreement, with deviations at the level of a few percent. The third row illustrates the scale dependence of the results, obtained by varying the factorization scale, taking $ Q^2 = \hat{s}/4 $ and $ Q^2 = \hat{s}/2 $. The observed differences can be attributed to this scale dependence and are also at the level of a few percent, as indicated in the table. No kinematic cuts are applied to the final-state Higgs bosons in the results presented in Table B1.