Heavy Higgs Bosons in 2HDM at a Muon Collider

We study the discovery potential of the non-Standard Model (SM) heavy Higgs bosons in the Two-Higgs-Doublet Models (2HDMs) at a multi-TeV muon collider and explore the discrimination power among different types of 2HDMs. We find that the pair production of the non-SM Higgs bosons via the universal gauge interactions is the dominant mechanism once above the kinematic threshold. Single Higgs boson production associated with a pair of heavy fermions could be important in the parameter region with enhanced Yukawa couplings. For both signal final states, $\mu^+\mu^-$ annihilation channels dominate over the vector boson fusion (VBF) processes, except at high center of mass energies where the VBF processes receive large logarithmic enhancement with the increase of energies. Single Higgs boson $s$-channel production in $\mu^+\mu^-$-annihilation via the radiative return can also be important for the Type-L 2HDM in the very large $\tan\beta$ region, extending the kinematic reach of the heavy Higgs boson mass to the collider energy. Considering both the production and decay of non-SM Higgs bosons, signals can be identified over the Standard Model backgrounds. With the pair production channels via annihilation, 95\% C.L. exclusion reaches in the Higgs mass up to the production mass threshold of $\sqrt{s}/2$ are possible when channels with different final states are combined. Including single production modes can extended the reach further. Different types of 2HDMs can be distinguishable for moderate and large values of $\tan\beta$.

Indirect limits on the heavy Higgs masses can also be obtained via the precision measurements of the SM-like Higgs couplings at future Higgs factories, when loop corrections are included [6][7][8][9].
Recently, there have been renewed interests for muon colliders operating at high energies in the range of multi-TeV [10][11][12]. This would offer great physics opportunity to open unprecedented new energy threshold for new physics, and provide precision measurements in a clean environment in leptonic collisions [13][14][15][16][17][18][19][20][21]. Recent studies indeed demonstrated the impressive physics potentials exploring the EW sector, including precision Higgs boson coupling measurements [19], the electroweak dark matter detection [20], and discovery of other BSM heavy particles [21].
In this article, we explore the discovery potential for the non-SM heavy Higgs bosons at a high-energy muon collider in the framework of 2HDMs [22]. We adopt the commonly studied four categories according to the assignments of a discrete Z 2 symmetry, which dictates the pattern of the Yukawa couplings. We identify the relevant parameters of the Higgs masses and couplings and predict the decay branching fractions. We take a conservative approach in the alignment limit for the mixing parameter so that there are no large corrections to the SM Higgs physics.
We consider the benchmark energies for the muon colliders [10] in the range of √ s = 3−30 TeV, with the integrated luminosity scaled as We study both the heavy Higgs boson pair production as well as single production associated with two heavy fermions. Both µ + µ − annihilation channels and Vector Boson Fusion (VBF) processes are considered, which are characteristically different. We also analyze the radiative return s-channel production of a heavy Higgs boson in µ + µ − annihilation, given the possible enhancement of the muon Yukawa couplings in certain models. We further consider the Higgs bosons to decay to heavy fermions and study their signatures and the SM backgrounds.
We design appropriate cuts to select the signals according to the different final states and kinematics, while effectively suppress the backgrounds. Combining together the production channels and the decay patterns, we also show how the four different types of 2HDMs can be distinguished.
The rest of the paper is organized as follows. In Sec. II, we briefly introduce the 2HDMs, identify the relevant parameters of the Higgs masses and couplings and predict the decay branching fractions. In Sec. III, we present the results for the heavy Higgs boson pair production and discuss the signal observability above the SM backgrounds. We also discuss how different types of 2HDMs can be distinguished by studying the production and decays of heavy Higgs bosons. In Sec. IV, we present the results for single heavy Higgs boson production in association with a pair of fermions. In Sec. V, we present the radiative return production of a non-SM Higgs boson and compare that to the other production mechanisms.
In Sec. VI, we summarize and draw our conclusion.
In this table, the interactions are classified into suppressed (∝ c β−α ) and un-suppressed (∝ s β−α ) categories. This is motivated by the scaling properties of the coupling strengths in the limit of cos(β − α) = 0, the so-called alignment limit. Under this limit, the couplings of the SM-like Higgs boson h with two gauge bosons restore the couplings in the SM, while the couplings of pair of SM gauge bosons to the non-SM Higgs bosons vanish. Given the Higgs coupling measurements at the LHC [6,24], the 2HDM parameter spaces have already been constrained to be near the alignment limit [7,8,25,26] . In our analyses below, we will assume the alignment limit of cos(β − α) = 0 when h will have the same tree-level couplings to gauge bosons and fermions as in the SM, and the other heavy Higgs bosons exhibit the universal gauge interactions. Even in the large tan β region of the Type-I 2HDM when the largest deviation of cos(β − α) from 0 is possible (at tree-level): cos(β − α) 0.3, the production cross sections will not be suppressed too much comparing to the results under the exact alignment limit as we presented in this paper. Furthermore, additional production channels could appear for non-zero cos(β − α), for example, µ + µ − → hA, which could be used to extend the reach of BSM Higgs bosons beyond the pair production threshold.
Note that for φ 1 φ 2 V and φ 1 φ 2 V 1 V 2 couplings, interactions involving two non-SM Higgs bosons are un-suppressed under the alignment limit, while interactions involving one SM-like Higgs h and one non-SM Higgs are always suppressed. This leads to the pair production of non-SM Higgs bosons via the full gauge interaction strength at a muon collider, while the single non-SM Higgs production in association with a SM gauge boson is suppressed by the small mixing.
The Lagrangian of Yukawa couplings is whereΦ = iσ 2 Φ * and Φ u,d,e are either Φ 1 or Φ 2 , depending on the Z 2 charge assignments.
Expanding Eq. (4) in terms of the physical Higgs fields, the interactions of Higgs bosons with fermions can be expressed as where u = (u, c, t), d = (d, s, b), V ij is the CKM matrix, and P L/R ≡ (1 ∓ γ 5 )/2 are the projection operators for the left-/right-handed fermions. In this expression, factors ξ are the coupling strengths normalized to the corresponding SM value, which are shown below at Eq. (6) under the alignment limit of cos(β − α) = 0.
There are four different types of Yukawa coupling structure depending on how the two Higgs doublets are coupled to the leptons and quarks. The Z 2 charge assignments and the corresponding Yukawa interactions are listed in Table II. The tree-level expressions of various ξs can be found in Ref. [22]. Under the alignment limit of cos(β − α) = 0, the Higgs is the sine (cosine) of the weak mixing angle θ W , and c 2W ≡ cos(2θ W ), couplings of H and A to the SM fermions normalized to the SM values are Type-II: ξ Huu = ξ Auu = cot β, −ξ Hdd = ξ Add = tan β, −ξ H = ξ A = tan β; Type-L: ξ Huu = ξ Auu = cot β, ξ Hdd = −ξ Add = cot β, −ξ H = ξ A = tan β; Type-F: ξ Huu = ξ Auu = cot β, −ξ Hdd = ξ Add = tan β, ξ H = −ξ A = cot β. (6) While the top-quark Yukawa coupling is always enhanced at small tan β and suppressed at large tan β, the Yukawa couplings to the bottom quark and tau lepton can be either suppressed or enhanced in difference regions of tan β, depending on the types of 2HDMs.
Theoretical considerations, such as the requirements of vacuum stability, perturbativity and unitarity impose additional constraints on the 2HDM model parameters, which can be translated to the masses and their splittings of the non-SM Higgs bosons. Detailed analyses can be found in Refs. [27][28][29][30][31]. The constraints are the weakest for λv 2 = m 2 H −m 2 12 /(s β c β ) = 0 under the degenerate mass assumptions, when all values of tan β are allowed. The allowed range of tan β shrinks when the value of λv 2 deviates from zero. For example, for λv 2 = (300 GeV) 2 , the allowed range of tan β is between 0.4 and 2.5 under the alignment and degenerate mass limits [6]. Furthermore, ranges of the scalar mass splittings also depend on the values of λv 2 . For our numerical analyses of non-SM Higgs production below, we set λv 2 = 0. Non-zero values of λv 2 will not change the cross sections significantly.
There have been extensive searches for non-SM Higgs bosons at the LHC [1,24]. For the neutral Higgs bosons, A/H → τ + τ − sets the strongest direct search limits. For Type-II 2HDM with enhanced bb and τ τ Yukawa couplings at large tan β, the lower limit on m A/H is 1 TeV (1.5 TeV) for tan β = 10 (50) [1,24]. At the small tan β region, both Type-I and Type-II get strong constraints from A/H → γγ, A/H → tt and four tops non-resonance search channels. For cos(β − α) ∼ 0, m H/A < 4m t with tan β < 1 is currently excluded, and the exclusion reaches 1 TeV for tan β < 0.2.
Electroweak precision measurements also provide indirect constraints on the mass splitting between H/A and H ± . Given the current 95% C.L. range of the oblique parameters [7,8], the charged Higgs boson mass is constrained to be close to one of the neutral Higgs masses: m H ± ∼ m H or m H ± ∼ m A . Precision measurements at the Z-pole and the Higgs factory could further limit the range of the mass splittings, as studied in Refs. [7,8].
In addition, since the mass differences between the BSM Higgses are controlled by the quartic couplings in the Higgs potential, theoretical considerations (vacuum stability, unitarity and perturbativity) result in ∆m 200 (100) GeV for m Φ = 1 (2) TeV [2]. In our analyses below, we take the degenerate mass assumption of m Φ ≡ m H = m A = m H ± . Numerically, the pair production cross of H + H − or Higgs produced in associated with a pair of heavy fermions via annihilation will not change since it only depends on the corresponding BSM Higgs mass. Other production process could have different production cross sections, due to the change of the masses of either the intermediate or the final state Higgses. The numerical difference, however, is not expected to be large given the viable ranges of the mass differences mentioned above. It is also worth mentioning that once the mass difference is larger than m W,Z , additional exotic decay modes open, for example, A → HZ, H ± → HW ± , which could be used to enhance the reach [2].
Flavor constraints, such as those from B 0 d −B 0 d mixing, b → sγ, B and D meson/baryon decays also bound the 2HDM parameter space, in particular, on m H ± and regions with Yukawa coupling enhancement. Most notably, the measured value for BR(b → sγ) imposes a lower limit on the charged Higgs mass to be larger than ∼ 600 GeV in the Type-II 2HDM [32,33].

B. Higgs boson decays
Under the alignment limit with degenerate heavy Higgs boson masses, only the decays of the heavy Higgs bosons to a pair of SM fermions are relevant. Since the couplings to the fermions are proportional to their masses, the leading decay channels are to the heavy fermions H/A → tt, bb, and τ + τ − , H ± → tb, and τ ν τ .
For illustration, we take m Φ ≡ m H = m A = m H ± = 2 TeV, under the alignment limit cos(β − α) = 0. We calculate the decay branching fractions in all the four types of 2HDMs, using 2HDMC [34], and the results are presented in Fig. 1  H ± decay, the suppression of tb decay caused by the enhancement of other decay channels is less obvious given that H ± → tb also enhances at large tan β in both Type-II and Type-F.
H ± → τ ν τ , however, could be dominant in the large tan β region for Type-L.
One noteworthy point is, while the decay branching fraction of the dominant decay mode for the different types of 2HDMs degenerates at small tan β, it is quite distinct at large tan β, which would allow the discrimination between different types of 2HDMs by examining the decays of heavy Higgs bosons.

SON FUSION
A high energy muon collider would have the capability to open a new energy threshold at the energy frontier. While the µ + µ − annihilation will be most efficient in exploiting the available c.m. energy for heavy particle production, it has been argued that the VBF processes will become increasingly more important at higher energies and offer a variety of   production channels due to the initial state spectrum.

A. Production cross sections
Once crossing the pair production threshold, the heavy Higgs bosons can be produced in pair via the µ + µ − annihilation The Feynman diagrams of the leading contributions are shown in Fig. 2. In the alignment limit of cos(β − α) = 0, the production is fully governed by the EW gauge interactions, TeV (solid curves), 2 TeV (dashed curves) and 5 TeV (dotted curves). Red and green curves are for H + H − and HA productions. We see the threshold behavior for a scalar pair production in a P-wave as σ ∼ β 3 , with β = 1 − 4m 2 H /s. Well above the threshold, the cross sections asymptotically approach σ ∼ α 2 /s, which is insensitive to the heavy Higgs mass. The excess of the H + H − production cross section over that of HA is attributed to the γ * -mediated process. The vertical axis on the right shows the corresponding events for a 10 ab −1 integrated luminosity, yielding a sizeable event sample.
At high c.m. energies, the VBF processes become increasingly important. For a fusion process of the initial state vector boson partons V i and V j , we write the production cross section of an exclusive final state F and the unspecified remnants X in terms of the parton luminosity dL ij /dτ and the corresponding partonic sub-process cross sectionσ where τ =ŝ/s with √ŝ being the parton-level c.m. energy. The production threshold is In this expression, f i (ξ, Q 2 ) stands for the electroweak parton distribution function (EW PDF) of particle V i radiated off the initial muon beam carrying an energy fraction x at a factorization scale Q. In our study we adopt the leading-order EW PDFs [37,38]. Recently the EW PDFs have been calculated with the DGLAP evolution at the double-log accuracy [11]. The numerical difference from the leading-order result is typical less than 10%. In Fig. 4, we present the vector-boson luminosity dL ij /dτ for µ + µ − collisions at 14 TeV. The QED γγ luminosity in the left panel is the largest at the lowŝ from the enhancement of log(Q 2 /m 2 µ ) versus log(Q 2 /m 2 V ) for massive vector bosons, while the difference becomes much smaller at higher energies. Such logarithmic enhancements are absent for the longitudinal massive gauge bosons, which leads to the suppressed partonic luminosity, as seen for W L W L luminosity by the dashed black curve in the right panel. The smallness of the ZZ luminosity is related to the accidentally small vector-like coupling of the neutral current, proportional to (1/2 − 2 sin 2 θ W ) for the un-polarized muon beam.

Heavy Higgs boson pair production via VBF is via
where We see the expected logarithmic enhancement over the energy log 2 (s/m 2 µ ) (or log 2 (s/m 2 V )) for initial state photons (weak bosons). Unlike the production by the µ + µ − annihilation, the cross section for the VBF processes are very sensitive to the heavy Higgs masses. The decrease of the cross sections with large m Φ is primarily from the suppression of EW PDF threshold ∼ 1/M 2 F . Again, the dominance of H + H − production cross section over the neutral ones is due to the exclusive contribution from the γγ fusion. HA production cross section is much smaller comparing to other processes since the first diagram of 4-point interaction in Fig. 5 without a heavy propagator suppression gives the leading contribution, which is absent for the HA production.
In general, the µ + µ − annihilation process yields more Higgs pairs than the VBF process, except for the H + H − production, when VBF takes over from One of the advantages for adopting the EW PDF approach for the calculations is the effective separation of the individual contributions from the fusion sub-processes. We illustrate this by presenting the contributing channels in Fig. 6 for H + H − (left panel) and HH (right panel) production. We see that the γγ (red) fusion sub-process contributes dominantly to the H + H − production through electromagnetic interaction, which is explained by the abundant availability of γγ in Fig. 4. Among the W W fusions of different polarizations, W T W T (purple) is much more copious than W L W L (blue) because the scaling behavior of log(Q 2 /m 2 W ) for the transversely polarized gauge bosons is absent for the longitudinally polarized gauge bosons. The small contribution from the ZZ fusion sub-process is related to its smallness in the partonic luminosity. For the HH production, all sub-process initiated with γ fusion are absent, thus the production is mainly initiated by the W T W T fusion and the total cross section is much smaller. Compared to H + H − , the ZZ fusion in HH process is more prominent, which can be explained by the fact that the ZH + H − coupling involved in H + H − process is smaller than the ZHA coupling involved in the t-channel contribution to HH process by a factor of c 2W (see Table I).
We note that, although a VBF channel has two accompanying leptons associated with the initial state gauge bosons, they are mostly unobservable due to their forward-backward collinear nature. As such, the µ + µ − annihilation and VBF both lead to the same observable Higgs pair final states. However, the invariant mass distributions of the Higgs pair system present a qualitatively different feature for these two processes, which may serve as an effective discriminator to separate these two processes. The invariant mass of the Higgs pair is approximately equal to the collider c.m. energy m φ 1 φ 2 ≈ √ s for the direct annihilation process, while peaked near the threshold m φ 1 φ 2 ≈ m φ 1 + m φ 2 for the VBF process, as shown    7 for H + H − production process. The long tail in the low invariant mass region for the annihilation process is due to the ISR effects taken into account in our calculations. For the rest of the analyses, we assume that the µ + µ − annihilation process and those from VBF are readily separable by kinematics.

B. Signals and backgrounds
Given the pair-production of the heavy Higgs bosons and focusing on leading decay channels as shown in the last section, the signal will be four heavy third-generation fermions.
Thus, the leading irreducible backgrounds are tttt, ttbb and bbbb for third generation quark final states, dominated by the EW-QCD mixed processes of quark pair production, followed by a gluon radiation and splitting to another pair of quarks, typically µ + µ − → ttg * , bbg * followed by g * → bb. For the sake of illustration, we consider the signal for a heavy Higgs mass m Φ = 2 TeV and the corresponding backgrounds. We note that the signal kinematics for the heavy Higgs pair production is quite different from the background processes, which would help for our signal reconstruction and background suppression. First, the decay products of the heavy Higgs bosons possess a Jacobian peak in the transverse momentum p T ≈ m Φ /2 and they are more central in the angular distribution, while the fermions in the σ √ s ttbb tttt bbbb  background tend to be softer and much more forward-backward. We thus impose the basic acceptance cuts where the choice for the angular cut is to simulate the detector coverage [39].
Depending on the specific production channels and decays, the signals will be characterized by the mass reconstruction. We therefore further require for HA channel : m tt , m bb > 0.9M H/A , θ tt , θ bb < 150 • .
to reconstruct the resonance masses, where θ is the opening angle between the two fermions in the µ + µ − lab frame. The angular cut above is to suppress the background processes with the fermion pair back-to-back production, while the Higgs decay products are more collinear due to the potentially boosted Higgs bosons. With those selection cuts above, the four heavy fermion backgrounds are highly suppressed as shown in Table IV, resulting in the background rates about 10 −3 fb or below, while the signal efficiencies shown in Table V are kept at 60% − 80%, depending on the collider energy. Here we only calculated the signal efficiency for the µ + µ − -annihilation process since the production through VBF is much smaller at √ s = 14 TeV (see Fig. 3). For the bbbb final states, the slightly lower signal efficiencies as seen in Table V    There are other kinematic features that could help further purify the signal samples from the backgrounds. For instance, at high energies, the Higgs bosons are produced back-to-back θ φ 1 φ 2 ≈ π, so that two fermion pairs may form a large angle. In contrast, the background is primarily from a back-to-back fermion pair production, followed by a collinear gluon splitting to the second fermion pair. Consequently, three quarks tend to cluster close by, going against another single energetic quark.
It is important to note that there are two classes of kinematic topologies for the signal of Higgs pair production, namely, the µ + µ − annihilation at high invariant mass and the VBF near the Higgs pair threshold, as seen in Fig. 7. When needed, the invariant mass variable of the Higgs pair system can serve as a discriminator for the production mechanisms. The kinematics of the decay products, however, would look largely the same since it is governed by the heavy Higgs mass.

C. Distinguishing 2HDMs
The pair production rates for all four types of 2HDMs are the same, since those only involve gauge coupling structures, or tri-Higgs couplings, which are independent of the Yukawa coupling structures. The decay branching fractions into fermion pairs, however, are different, which are determined by their Yukawa couplings characterized by tree-level parameter tan β, as shown in Fig. 1. We will focus on the leading decay channels for the non-SM Higgs bosons as in Eq. (7). In • For low values of tan β < 5, the four models cannot be distinguished since the dominating decay channels are the same: H/A → tt, H ± → tb.
• For large values of tan β > 10, the decay modes of H/A → τ + τ − , H ± → τ ν become substantial in Type-L, which can be used to separate Type-L from the other three types of 2HDMs.
• For tan β > 5, the enhancement of the bottom Yukawa coupling with tan β in Type-II/F leads to the growing and even the leading of H/A → bb decay branching fraction, which can be used to separate Type-II/F from the Type-I 2HDM. Adding other channels will increase the reach, in particular, beyond the pair production mass threshold of √ s/2.
A representative Feynman diagram of the dominant contributions is shown in Fig. 9. The calculation is performed with tree-level diagrams. However, we include the large higherorder effects for the running of the Yukawa couplings (Y u,d,e in Eq. (4)) to the corresponding scale µ = m Φ by solving the renormalization Group Equations (RGEs) [40]. All the input parameters listed in Sec. II as well as the quark/lepton masses for the RGEs are given at µ = m Z [41]. For tan β = 1 at m Z , the running Yukawa couplings at m Z , 1 TeV and 2 TeV are listed in Table VII. Effectively, compared with results using parameters at a fixed scale m Z , we will have a suppression about 10% and 14% − 30% for top and bottom quark processes, respectively.
To simulate the detector acceptance in our partonic-level calculations, we first apply the simple cuts on the final state fermions p f T > 50 GeV and 10 • < θ f < 170 • .
A veto cut of 0.8m Φ < m f f < 1.2m Φ is applied to the associated fermions to remove contributions from resonant Higgs decays. The signal cross sections are shown in the left m Z 1 TeV 2 TeV m Z 1 TeV 2 TeV m Z 1 TeV 2 TeV  for m H ± =1 TeV. ttH production cross section is about factor of 3 smaller. bbH and τ + τ − H cross sections are further reduced by a factor of (m b,τ /m t ) 2 . However, this hierarchical structure could be altered by the choice of tan β in different types, which will be discussed later in Sec. IV C. We see the advantage of the accessible events below the pair production threshold of √ s < 2m Φ . The production rates well above the threshold are smaller than that of the pair production via the µ + µ − annihilation in Sec. III by a factor of few due to the 3-body kinematics. The mass dependence of the single production cross sections is also stronger comparing to that of the pair production processes via annihilation, as shown in the solid curves for m Φ = 1 TeV and dashed curves for m Φ = 2 TeV in Fig. 10, comparing to that in the left panel of Fig. 3.
The fermion associated single heavy Higgs can also be produced via the VBF processes.
In addition to the charge-neutral states in Eq. (14), the fusions of W ± γ/Z give rise to rich charged final states. The complete set of the dominant ones are Some representative Feynman diagrams of the dominant contributions are shown in Fig. 11.
There also exist the pure ν τ associated final states such as ν τντ H/A and ν τντ H ± via VBF processes. However, given the absence of neutrino Yukawa couplings, their cross sections are about two orders of magnitude smaller than the corresponding production with τ . The simple acceptance cuts as in Eq. (15) are again applied, which help to regularize the singularities of the outgoing fermions in the forward regions. The cross sections for the quark associated production as a function of the c.m. energy √ s for tan β = 1 are shown in the right panel of Fig. 10. While the similar hierarchical structure of production cross section is apparent, the production cross sections also manifest the rising trend with the increasing of the c.m. energy typical to the VBF processes. Compared to the charge-neutral final states, the charged final states have comparable cross sections. This is due to the fact that the partonic luminosities of W γ and γγ are about the same at higher energies, which dominate the production of charged and charge-neutral final states, respectively. The VBF production cross sections also exhibit sensitive mass dependence rough as 1/m 2 Φ , as shown by the solid (m Φ = 1 TeV) and dashed (m Φ = 2 TeV) lines. The relative size of various VBF processes, however, could vary as m Φ varies. For m Φ = 1 TeV, tbH ± is larger than ttH ± , while the order is flipped for m Φ = 2 TeV. Similarly behavior occurs for bbH ± versus tbH/A. There are several factors contributing to this, for example, the difference in the dominant V V contributions in difference processes, parton luminosity differences between the γ, V L and V T , and the contributions of top and bottom Yukawa couplings that enter differently in different processes, as well as the chiral suppression of the bottom quark mass. The cross sections of the leading production channels are summarized in Table VIII. tbH ± has the largest production cross section for the annihilation, while both tbH ± and ttH ± contribute for the VBF processes.
For the neutral Higgs production, ttH/A via annihilation is important for lower c.m. energies, while tbH/A, bbH/A and ttH/A via VBF could be important for higher c.m. energies and low scalar masses. One of the advantages for adopting the EW PDF approach is to appreciate the underlying contributions from the individual sub-processes. In Fig. 12, we plot the contributions from individual VBF processes to the charge-neutral final state process tbH ± , ttH and bbH, and charged final state process ttH ± with the benchmark heavy Higgs mass m Φ = 2 TeV and tan β = 1. By comparing the top two plots, we find that the production of ttH ± through W γ fusion is slightly larger than tbH ± through γγ fusion. That is because the contributions from diagram (a) in Fig. 13 is twice as large for W γ fusion in ttH ± production with H/A in the internal line, comparing to γγ fusion in tbH ± production with H ± in the internal line.
γγ fusion contribution in ttH is much smaller since contribution from diagram (a) is absent.

B. Signals and backgrounds
Following the discussions of the signal construction and background suppression as in Sec. III B, we present the scheme to identify the Higgs boson signal from the associated production with a pair of heavy fermions, again illustrated for m Φ = 2 TeV and √ s = 14 TeV. We first note that the four-fermion background processes are of the same origin as in Sec. III B. As for the signal, there is only one heavy Higgs boson in the events. As an illustration, we focus on the two leading production channels tbH ± and ttH. We first impose the basic acceptance cuts of p T and θ cuts in Eq. (15) on all the final state fermions, and then propose the appropriate cuts to suppress the irreducible backgrounds for the dominate decay modes.

tbH ± → tbtb
We first consider the signal production via the µ + µ − annihilation. For the signal reconstruction, we find the pair of tb which gives the closest invariant mass to the hypothetical m H ± . The other pair of tb unlikely to be from the Higgs decay is denoted by In the upper left panel of Fig. 15, we present the distribution for the scatter events Comparing it with the left panel, we find that they exhibit very distinct kinematic features.
To further separate the signal and background, in the bottom two panels of Fig. 15, we depict their distributions in the plane of ∆R(tt )-∆R(bb ). Again, we see the back-to-back feature of the two tb pairs for the signal, namely, most signal events concentrate at θ tt = π σ (fb)   and θ bb = π, because t and b originating from the H ± decay tends to fly in the same direction due to the high energy boost. However background events own two concentrations in the lower right panel, the one with bb from the gluon splitting, and the one with tb from W decay. Based on these distributions, we thus propose the following additional cuts to suppress the ttbb background for the tbH ± channel for the benchmark Higgs mass of 2 TeV: where t 1 and b 1 are the top and bottom quarks with leading p T . We note that the leading p T fermions may not be from the heavy Higgs decay in this production mechanism. With the above cuts, we can suppressed the ttbb background down to the level of 10 −3 fb, while retaining the signal at the level of 10 −2 fb. The cuts are more efficient at high c.m. energy since they're proposed based on the specific features of the event distribution manifested at high enough energy. At √ s = 14 TeV, we can achieve a signal rate up to 72%. However, at low c.m. energies, the distributions shown in Fig. 15 would be less distinguishable. The optimal cut selection should be properly tuned according to the hypothetical Higgs mass and the c.m. energy. The background cross sections after the cuts are given in Table IX and the corresponding signal rates are shown in Table X.
For the VBF processes, the system is close to the production threshold and thus the final state quarks are less boosted. The angular cuts should be adjusted accordingly. For the top and bottom quarks that are paired from the heavy Higgs decay, we require For the other pair of top and bottom quarks, we require Then for the two tops and two bottoms, we require: We again achieve efficient signal-background separation, as presented in Table X and Ta-ble IX for the signal and background, respectively.

ttH → ttbb
We again first consider the annihilation production. One advantage of this decay mode is that the heavy Higgs reconstruction is free of the combinatorial problem. Based on a similar analysis as in the previous section, in additional to the basic acceptance in Eq. (15), we implement the following cuts to suppress the background where t 2 and b 2 are top and bottom quarks sub-leading in p T .
For the VBF process, the cuts are accordingly adjusted as The achieved signal and background separation, as presented in Table X and Table IX for the signal and background, respectively.

ttH → tttt
Tops sorted by p T from high to low are labeled by t 1 , t 2 , t 3 , t 4 . The top pair chosen to reconstruct heavy Higgs are denoted by tt while the other two tops are denoted by t t .
We reconstruct heavy Higgs based on invariant mass close to the hypothetical m Φ in the analyses. For the µ + µ − annihilation production, the cuts we propose to suppress the tttt background are With those cuts, we find that the overall event efficiencies are slightly lower than the ttbb final states due to the combinatorics, ranging from 62% compared to 72% at √ s = 14 TeV and 60% compared to 71% at √ s = 30 TeV.
For the VBF production, the cuts are slightly different, The achieved signal and background separation, as presented in Table X and Table IX for the signal and background, respectively.
In summary, we have demonstrated that it is quite conceivable to archive very desirable signal identification over the SM backgrounds for the associated heavy Higgs production for both direct µ + µ − annihilation and the VBF channels.

C. Distinguishing 2HDMs
The production of heavy Higgs bosons in association with fermions is achieved through Yukawa couplings, thus the cross sections are sensitive to tan β, which behave differently for different types of 2HDMs. In Fig. 16, we demonstrate the tan β dependence of quarkassociated production channels in the top two panels and τ associated production channels in the bottom two panels. The left panels are for the production through µ + µ − annihilation and the right panels are the production through VBF. As a benchmark point, we choose degenerate heavy Higgs mass m Φ = 2 TeV (Φ = H, A, H ± ) and the collider c.m. energy √ s = 14 TeV.
The tan β dependence for the production through µ + µ − annihilation is directly related to the Yukawa couplings in Eq. (6). For the quark associated production, they're uniformly proportional to 1/ tan 2 β in Type-I/L. In Type-II/F, the tt and bb-associated production is proportional to (Y t / tan β) 2 and (Y d tan β) 2 , respectively. The charged Higgs production For the production through VBF, diverse production diagrams in Fig. 11 make the tan β dependence more complicated, but the overall behaviours can be explained by the large contributions from diagram (a) in Fig. 13, unless it is severely suppressed by the small Yukawa coupling, or the absence of the V φ 1 φ 2 coupling.
For the τ associated production, the cross sections scale as 1/ tan 2 β in Type-I/F, and tan 2 β in Type-II/L. They are only large enough for observation in Type-II/L at large tan β region.
In Table XI we summarized the leading signal channels of the Higgs associated production with fermions in four types of 2HDMs in different regimes of tan β. Several observations can be made: • In the small tan β < 5 region, all six production channels have sizable production cross sections. However, it is hard to distinguish different types of 2HDMs since they all lead to the same final states.
• In the large tan β > 10 region, all the production channels for the Type-I are suppressed, while Type-II/F have sizable production in tbH ± , bbH ± , bbH/A and tbH/A channels. Type-II and Type-F can be further separated by studying the sub-dominant decay channels of H ± → τ ν τ and H/A → τ + τ − in the Type-II. Same final states of Type-II can also be obtained via τ ν τ H ± , τ + τ − H ± , τ + τ − H/A and τ ν τ H/A production. enough to compensate, resulting in a rather low signal production rate. A rich set of final states, however, are available given the various competing decay modes of H ± and H/A.

V. RADIATIVE RETURN
While the cross sections for heavy Higgs pair production are un-suppressed under the alignment limit, the cross section has a threshold cutoff at m Φ ∼ √ s/2. The resonant production for a single heavy Higgs boson may further extend the coverage to about m Φ ∼ √ s, as long as the coupling strength to µ + µ − is big enough. The drawback for the resonant production is that the collider energy would have to be tuned close to the mass of the heavy Higgs, which is less feasible at future muon colliders. A promising mechanism is to take advantage of the initial state radiation (ISR), so that the colliding energy is reduced to a lower value for a resonant production, thus dubbed the "radiative return" [42], as shown in Fig. 17. This mechanism can be characterized by the process where γ can be a mono-photon observed in the detector, or unobserved along the beam as the collinear radiation. We first calculate the cross section of the mono-photon process for m H = 1, 5, 15 TeV at tan β = 1. 10 • < θ γ < 170 • is imposed for the photon detection acceptance.
For a single photon production, its energy is mono-chromatic E γ = (s − m 2 H )/2 √ s. The results are given in the left panel of Fig. 18 by the dashed curves.
As a comparison, we calculate the µ + µ − → H process with ISR spectrum where x is the energy fraction carried by the muon after the ISR. The partonic cross section To compare with process in Eq. (25), we calculate the cross section to the first order of α by convoluting the ISR spectrum to one muon beam, The results are given in the left panel of Fig. 18 by the solid curves. The cross section increases as the heavy Higgs mass approaches the collider c.m. energy, closer to the s-channel resonant production.
The right panel of Fig. 18 shows the tan β dependence of the cross section for √ s = 14 TeV and m H = 12 TeV. While the cross section at tan β = 1 is much smaller than the other production channels we considered earlier, the cross section scales like tan 2 β in Type-II/L, which could be sizable at large tan β. It could even be the dominant production for heavy Higgs in the large tan β region of Type-L, when pair production is kinematically forbidden and quark associated productions are suppressed.

VI. SUMMARY
High energy muon colliders offer new opportunities for the direct production of heavy particles. In this paper, we studied the discovery potential of the heavy Higgs bosons in Two-Higgs-Doublet Models (2HDMs) at a high-energy muon collider. We take the alignment limit so that the interactions between the Higgs bosons and the SM gauge bosons are of the universal gauge interactions. We explored the pair production of non-SM Higgs bosons, and single non-SM Higgs production in association with a pair of heavy fermions from both µ + µ − -annihilation and the VBF mechanism, as well as radiative return production of a single non-SM Higgs boson directly coupled to µ + µ − . We considered the heavy Higgs boson decays to heavy fermions such as the t-quark, b-quark, and a τ lepton. With appropriate cuts on the invariant mass, transverse momenta, and angular separation between heavy fermions, the dominant SM backgrounds can be effectively suppressed to a negligible level.
We found that the pair production of the heavy Higgs bosons is the dominant mechanism for √ s > 2m Φ , while the single non-SM Higgs production associated with a pair of heavy fermions could be important for heavier masses, and in regions of tan β with Yukawa coupling enhancement. We also compared the annihilation production versus the VBF production, and found that VBF processes could be dominating at large center of mass energy and low scalar masses. Radiative return for the single Higgs boson production, in particular, could be important in the large tan β region of Type-L, extending the mass coverage to the kinematic limit √ s ∼ m Φ . With the pair production channels via annihilation, 95% C.L. exclusion reaches in the Higgs mass up to the production mass threshold of √ s/2 are possible when channels with different final states are combined. Including single production modes can extend the reach further. We reiterate that the discovery coverage at a muon collider would be quite complementary to that at future hadron colliders, where the signal-to-background ratio is low and the signal identification depends heavily on the final states from the heavy Higgs decays and on the sophisticated multiple variable analyses [2][3][4][5].
We also assessed the discrimination power of a muon collider on different types of 2HDMs.
With the combination of both the production mechanisms and decay patterns, we found that while it is challenging to distinguish different types of 2HDMs at the low tan β region, the intermediate and large tan β values offer great discrimination power to separate Type-I and Type-L from Type-II/F. To further identify either Type-II or Type-F, we need to study the subdominant channels with τ final states, which could be sizable in the signal rate in Type-II.
Our analyses were performed on the four general types of 2HDMs, in which a Z 2 symmetry is imposed such that one fermion species only couples to one Higgs doublet to avoid treelevel FCNC. There exist scenarios in the literature [43][44][45][46][47] in which flavor alignment or other mechanisms are adopted to suppress the dangerous FCNC effects. As a result, some of the