## Higgs boson production at µ µ + + colliders

## Yu Hamada, a,b Ryuichiro Kitano, c,d Ryutaro Matsudo, c Shohei Okawa, c Ryoto Takai, c,d Hiromasa Takaura e and Lukas Treuer c,d

- a Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
- b Research and Education Center for Natural Sciences, Keio University, 4-1-1 Hiyoshi, Yokohama, Kanagawa 223-8521, Japan
- c KEK Theory Center, Tsukuba 305-0801, Japan
- d The Graduate University for Advanced Studies (SOKENDAI), Tsukuba 305-0801, Japan
- e Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

E-mail: yu.hamada@desy.de , ryuichiro.kitano@kek.jp , matsudo@post.kek.jp , shohei.okawa@kek.jp , rtakai@post.kek.jp , hiromasa.takaura@yukawa.kyoto-u.ac.jp , ltreuer@post.kek.jp

Abstract: We study Higgs boson production at µ + µ + colliders at high energy. Since both initial-state particles are positively charged, there is no W boson fusion at the leading order, as it requires a W W + -pair. However, we find that the cross section of the higher-order, γ -and Z -mediated W boson fusion process is large at high center-of-mass energies √ s , growing as (log s ) 3 . This is in contrast to the log s behavior of the leading-order W boson fusion. Thus, even though it is a higher-order process, the rate of Higgs boson production for 10 TeV energies at µ + µ + colliders with polarized beams can be as high as about half of the one at µ + µ -colliders, assuming the same integrated luminosity. To calculate the cross section of this process accurately, we carefully treat the collinear emission of the photon in the intermediate state. The thereby obtained large cross section furthermore shows the significance of Higgs production with an extra W boson in the final state also at µ + µ -and e + e -colliders.

## Contents

| 1   | Introduction                                                    | Introduction                                                    |   1 |
|-----|-----------------------------------------------------------------|-----------------------------------------------------------------|-----|
| 2   | γ - and Z -mediated W boson fusion process at µ + µ + colliders | γ - and Z -mediated W boson fusion process at µ + µ + colliders |   3 |
|     | 2.1                                                             | Triple logarithms in the Higgs boson production                 |   3 |
|     | 2.2                                                             | Contributions from the Z -mediated process                      |   6 |
|     | 2.3                                                             | Numerical calculation with MadGraph                             |   7 |
| 3   | Calculation with controlled theoretical uncertainties           | Calculation with controlled theoretical uncertainties           |  10 |
|     | 3.1                                                             | Low- p T cross sections from EPA                                |  11 |
|     | 3.2                                                             | High- p T cross sections from MadGraph                          |  13 |
|     | 3.3                                                             | Summing up the low- and high- p T cross sections                |  13 |
| 4   | Event shape of the W boson fusion process                       | Event shape of the W boson fusion process                       |  15 |
| 5   | Summary                                                         | Summary                                                         |  17 |
| A   | Equivalent photon approximation (EPA)                           | Equivalent photon approximation (EPA)                           |  19 |
| B   | Error in the massless lepton approximation                      | Error in the massless lepton approximation                      |  22 |
| C   | Semi-automatic MadGraph implementation of p (cut) T calculation | Semi-automatic MadGraph implementation of p (cut) T calculation |  23 |

## 1 Introduction

Higgs boson factories are the best-motivated future colliders, enabling us to explore the nature of electroweak symmetry breaking in depth. Among others, muon colliders at 10 TeV energies have been discussed as one of the most attractive possibilities, since the size of the accelerator facilities can be as compact as O (10) km in circumference, and they can potentially probe physics up to scales as high as O (100) TeV [1-3]. By comparison, proton colliders with similar reach require a circumference of at least O (100) km.

Recently, muon colliders based on ultra-slow muon technology [4] have been proposed [5]. This is because ultra-slow muons from laser-ionized muonium, a µ + e -bound state, can be used to create µ + beams with excellent emittances, i.e., very collimated beams that facilitate high luminosities. As such, the corresponding proposal is to build on this technology, which was developed for the muon g -2/EDM experiment at J-PARC [6], leading to an estimated luminosity seemingly good enough for Higgs factories based on µ + e -and µ + µ + colliders.

At such colliders, even though there are no s -channel annihilation processes that are present in µ + µ -or e + e -colliders, Higgs boson production can still occur via vector boson

fusion processes, which dominate over s -channel annihilations at energies beyond O (1) TeV. This is due to a logarithmic enhancement of the cross section for vector boson fusion as a function of center-of-mass energy. At µ + e -colliders, while W boson fusion is possible, the center-of-mass energy is limited by the beam energy of the electron. On the other hand, for µ + µ + colliders, one can consider, e.g., O (10) TeV beam energies, but there is no W boson fusion process at leading order since both of the antimuons can only emit positively charged W bosons. Although Z boson fusion is possible at leading order since it is independent of the muon charge, its cross section is about an order of magnitude smaller than W boson fusion at µ + µ -colliders due to an unfortunate suppression of the coupling between the Z boson and leptons. This motivates us to look beyond the leading order.

In this paper, we thus investigate W boson fusion at µ + µ + colliders at higher order in perturbation theory. The pertinent process, which we show in Fig. 1, begins with emission of a photon or Z boson from one of the antimuons, followed by the photon or Z boson splitting into a W W + -pair. The W -from this splitting then collides with the W + emitted from the other antimuon to produce a Higgs boson. For large center-of-mass energies √ s , we find that this sub-leading-order process involves a (log s ) 3 factor, in contrast to the single log s appearing in the leading-order W boson fusion at µ + µ -and e + e -colliders [2, 7]. Due to this large enhancement, the cross section for Higgs boson production at µ + µ + colliders becomes comparable to that at µ + µ -colliders for O (10) TeV beam energies.

To compute the part of the process involving photon emission, we must carefully address infrared divergences, which are physically cut off by the muon mass. However, directly using numerical codes such as the event generator MadGraph [8] leads to instabilities in the numerical phase-space integration, since the muon mass is either set to zero, or significantly smaller than the center-of-mass energy. Therefore, we also discuss how to reliably and accurately compute fixed-order cross sections for such infrared-divergent processes.

We furthermore discuss the size of subsequent higher-order processes to determine if the fixed-order computation is meaningful. We find that at 10 TeV energies, the processes with further emissions of extra gauge bosons are much smaller than the process of our interest, such that fixed-order calculations are still valid.

Lastly, it is important for actual experiments whether the final-state W + has a transverse momentum large enough to be visible. To assess this, we generate sample events using MadGraph and find that for √ s = 2 TeV or 10 TeV, the vast majority of jets from the W + decays can be detected. This means that it is possible to analyze the W boson fusion process with a W h + final state separately from Z boson fusion to perform precision measurements of the Higgs boson couplings.

At the same time, the large size of the γ -and Z -mediated process means that an analogous process gives a large contribution to the Higgs production cross section also at µ + µ -and e + e -colliders. One should thus include this higher-order process in the simulation and analysis of the Higgs boson production for coupling measurements, as it may contaminate the events of the pure W boson fusion process.

This paper is organized as follows. In Sec. 2, we estimate the cross section of the γ -and Z -mediated Higgs boson production process at µ + µ + colliders analytically at the leading-logarithm level, as well as numerically, and discuss its high-energy behavior. We

then explain a more reliable method to obtain the cross section in Sec. 3, where we confirm the large enhancement of the cross section. Next, since the γ - and Z -mediated process gives an extra W boson in the final state, we simulate collider events and discuss the possibility of identifying the process in Sec. 4. Finally, we summarize our findings in Sec. 5.

In Appendix A, we review the derivation of the formula of the equivalent photon approximation, and discuss the uncertainties of the approximation. We also discuss the error associated with the treatment of neglecting the muon mass in the numerical calculation in Appendix B. Lastly, we detail a semi-automatic numerical method to implement the calculation explained in Sec. 3 in Appendix C.

## 2 γ - and Z -mediated W boson fusion process at µ + µ + colliders

In this section, we discuss how W boson fusion processes are possible at µ + µ + colliders, and why they are important at high energies. Therein, we derive an analytic formula for the cross section at the leading-logarithm level. We also calculate the cross section numerically using the event generator MadGraph [8] with a setting for parton distribution functions of the photon in one of the muons. Note that both methods have uncertainties related to their scale settings, and the results should therefore be viewed as estimations. We will thus discuss a method to obtain the cross section with controlled uncertainties in the next section. Lastly, in this section, we compare the approximate results with those obtained by the more reliable and accurate computation method, and show that the γ -and Z -mediated W boson fusion process is indeed important.

## 2.1

## Triple logarithms in the Higgs boson production

At high-energy lepton colliders, a large amount of Higgs bosons can be produced through vector boson fusion. When two antimuons collide with center-of-mass energy √ s ≳ 1 TeV, Z boson fusion is enhanced by a factor of log s , similarly to Z and W boson fusion processes at µ + µ -and e + e -colliders.

By contrast, there is no coupling between µ + and W -in the Standard Model at the first order in the electromagnetic coupling e or weak coupling g . This makes W boson fusion difficult, since it requires a W W + -pair. In particular, while the W boson fusion process may take place at higher orders in perturbation theory, these are seemingly suppressed by additional powers of the small couplings. Therefore, we na¨ ıvely do not expect a satisfyingly large contribution to Higgs production from W boson fusion at µ + µ + colliders. However, this expectation turns out to be incorrect, as we will see below.

We show two representative Feynman diagrams for Higgs boson production via such a higher-order W boson fusion in Fig. 1. Here, one antimuon provides a W + , and the other a W -through pair production by an intermediate photon or Z boson. When beam energies are much larger than the electroweak scale given by the vacuum expectation value of the Higgs field, v ≃ 246 GeV, these two diagrams are identical up to terms of the order of m /s 2 Z , where m Z is the mass of the Z boson. Note that one needs to, however, account for the different couplings of the photon and Z boson, as well as the infrared scales of their respective emission.

Figure 1 . We show representative Feynman diagrams of single Higgs boson production through quasi-real photon (left) or Z boson (right) emission, splitting, and subsequent W boson fusion.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-8f78d8c2-886a-4337-89a6-e1cee9e5b8ca.jpg)

First, we focus on the contribution of the photon, ignoring the Z boson. Quantum electrodynamics enables us to factorize the contribution of the photon into a parton distribution function with respect to the antimuon beam, and a partonic cross section with an on-shell photon in the initial state. The corresponding parton distribution function is given by

$$f _ { \gamma / \mu } ( x, \mu _ { f } ^ { 2 } ) = \frac { \alpha } { 2 \pi } \frac { 1 + ( 1 - x ) ^ { 2 } } { x } \log \frac { \mu _ { f } ^ { 2 } } { m _ { \mu } ^ { 2 } }$$

at the leading-logarithm level [9, 10], where α = e / 2 4 π is the fine structure constant, x is the longitudinal momentum fraction carried by the photon, µ f is the factorization scale, and m µ the muon mass. The large logarithm originates from an integral over the transverse momentum p T of the differential parton distribution function

$$\frac { \mathrm d f _ { \gamma / \mu } } { \mathrm d p _ { \mathrm T } ^ { 2 } } = \frac { \alpha } { 2 \pi } \frac { 1 + ( 1 - x ) ^ { 2 } } { x } \frac { 1 } { p _ { \mathrm T } ^ { 2 } } \,,$$

which comprises a p -4 T from the propagator of the photon and a p 2 T from the splitting amplitude. Here, we adopt the muon mass m µ as infrared scale of the emission, which physically cuts off the associated collinear divergence.

Similarly, one can also derive the parton distribution function of the W boson. For its longitudinal component, this yields [11, 12]

$$f _ { W _ { L } ^ { + } / \mu _ { R } ^ { + } } ( x ) = \frac { \alpha } { 2 \pi \sin ^ { 2 } \theta _ { W } } \frac { 1 - x } { x } \quad \text{and} \quad f _ { W _ { L } ^ { + } / \mu _ { L } ^ { + }$$

where θ W is the weak mixing angle, also known as Weinberg angle. The Nambu-Goldstone bosons eaten by the longitudinally polarized W bosons have a typical scale v , which modifies the dependence of the differential parton distribution function on the transverse momentum from 1 /p 2 T to ∼ v /p 2 4 T [13]. The integral over p 2 T thus does not lead to a large logarithm in the parton distribution function of Eq. (2.3).

To estimate the full cross section of the process as seen on the left in Fig. 1 using these parton distribution functions, we also need the cross section of the subprocess γW + →

W h + , which is given by

$$\tau \alpha _ { y } \quad & \sigma _ { \gamma W } ( s ) = \frac { \pi \alpha ^ { 2 } } { m _ { W } ^ { 2 } \sin ^ { 2 } \theta _ { W } } + \mathcal { O } \left ( \frac { 1 } { s } \right ) \quad & \quad \quad \quad$$

in the high energy limit s ≫ s min ≡ ( m h + m W ) 2 [14]. Here, m W and m h are the masses of W and Higgs bosons, respectively. We obtain this cross section by averaging over the polarizations of the initial photon, and summing over the ones of the final-state W + , while taking the initial-state W + to be longitudinally polarized. This is because the highenergy cross section is dominated by the contribution of order s 0 , which comes from the longitudinally polarized W boson in the initial state, since its polarization vector is almost proportional to its momentum. By contrast, the contribution from a transversely polarized W boson falls off as 1 /s at high energies. Therefore, it does not contribute in the highenergy limit, and we may assume a longitudinally polarized W boson in the initial state.

Finally, the photon contribution to the process µ + µ + → µ + ν µ W h + is thus approximately given by

$$\text{analy given by} \\ \sigma _ { \gamma } ( s ) & \simeq 2 \quad \int _ { \substack { x y s > s _ { \min } \\ x y s > s _ { \min } } } \, \text{d} x d y \, f _ { \gamma / \mu ^ { + } } ( x, \mu _ { f }$$

at the leading-logarithm level, where P µ + denotes the polarization of the antimuon beams. The factor of two in the first line accounts for photon or W + emission by either antimuon. For this estimate, we take the factorization scale to be the center-of-mass energy at the parton level, µ 2 f = xys .

We thus find three large logarithms in the final formula. Compared with the cross section of the leading-order Z boson fusion process [15, 16],

$$\sigma _ { Z \text{BF} } ( s ) \sim \frac { \alpha ^ { 3 } ( 1 - 4 \sin ^ { 2 } \theta _ { W } + 8 \sin ^ { 4 } \theta _ { W } ) ^ { 2 } } { 6 4 m _ { Z } ^ { 2 } \sin ^ { 6 } \theta _ { W } \cos ^ { 6 } \theta _ { W } } \log \frac { s } { m _ { h } ^ { 2 } } \,,$$

we observe a rapid growth of the cross section as a function of the collider energy. Furthermore, the cross section of the Z boson fusion process has an unfortunate suppression factor in the numerator. Because of this, at √ s ∼ O (10) TeV, the sub-leading, γ -mediated W boson fusion cross section becomes much larger than for the leading-order Z boson fusion process. Even when compared to the cross section of the leading-order W boson fusion process at µ + µ -or e + e -colliders [15, 16],

$$\sigma _ { W B F } ( s ) \sim \frac { \alpha ^ { 3 } } { 1 6 m _ { W } ^ { 2 } \sin ^ { 6 } \theta _ { W } } \log \frac { s } { m _ { h } ^ { 2 } } \,,$$

the γ - (and Z -) mediated process becomes important at high energy as two extra logarithmic factors can compensate the O ( α/π ) suppression. For reference, the logarithmic factors of Eq. (2.5) and Eq. (2.7) at 10 TeV evaluate to approximately 1100 and 10, respectively.

This may resemble a breakdown of perturbation theory of fixed-order computations. We will therefore return to the discussion of further higher-order processes later in this section.

## 2.2 Contributions from the Z -mediated process

At high energies, the contribution of the Z boson becomes important. However, since the γ - and Z -mediated diagrams interfere, we cannot simply add the Z -mediated cross section, estimated using a parton distribution function of the Z boson, to the γ -mediated one. Instead, a matrix form of parton distribution functions needs to be used [13, 17, 18]. Here, we therefore discuss how to estimate the additional contributions at the leading-logarithm order.

If one could treat the Z boson as a massless particle just like photon, the amplitudes in Fig. 1 should have the same form, i.e.,

$$\mathcal { M } = \mathcal { M } _ { \gamma } + \mathcal { M } _ { Z } = g _ { \gamma } \mathcal { M } _ { 0 } + g _ { Z } \mathcal { M } _ { 0 }$$

for each chirality of the antimuons, where g γ and g Z are the respective couplings of the photon and Z boson to the antimuons. The ratios of γ and Z couplings to antimuons and W bosons are given by

$$\frac { g _ { Z } ^ { \mu _ { R } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } = \frac { \frac { 1 } { 2 } - \sin ^ { 2 } \theta _ { W } } { \sin \theta _ { W } \cos \theta _ { W } } \simeq 0. 6 6 8 \,, \quad \frac { g _ { Z } ^ { \mu _ { L } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } = \frac { - \sin ^ { 2 } \theta _ { W } } { \sin \theta _ { W } \cos \theta _ { W } } \simeq - 0. 5 3 5 \,,$$

and

$$\frac { g _ { Z } ^ { W } } { g _ { \gamma } ^ { W } } = \frac { 1 } { \tan \theta _ { W } } \simeq 1. 8 7$$

for sin 2 θ W ≃ 0 222, where we have taken the value used by MadGraph for consistency with . later computations. Since the full matrix element can be expressed in terms of a shared amplitude M 0 when taking the masslessZ limit, we can approximate the squared matrix element as the sum of three squared amplitudes as

$$| \mathcal { M } | ^ { 2 } = | g _ { \gamma } | ^ { 2 } | \mathcal { M } _ { 0 } | ^ { 2 } + | g _ { Z } | ^ { 2 } | \mathcal { M } _ { 0 } | ^ { 2 } + 2 \text{Re} ( g _ { \gamma } g _ { Z } ) | \mathcal { M } _ { 0 } | ^ { 2 } + \mathcal { O } ( m _ { Z } ^ { 2 } / s ) \,. \quad \ \ ( 2. 1 1 )$$

This result allows us to estimate the cross section by using the three parton distribution functions of the photon, the Z boson, and their mixing for the corresponding terms. These functions satisfy

$$\frac { \mathrm d f _ { Z _ { \mathrm T } / \mu _ { R } ^ { + } } } { \mathrm d p _ { \mathrm T } ^ { 2 } } = \left ( \frac { g _ { Z } ^ { \mu _ { R } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \right ) ^ { 2 } \frac { \mathrm d f _ { \gamma / \mu } } { \mathrm d p _ { \mathrm T } ^ { 2 } }, \quad \frac { \mathrm d f _ { Z _ { \mathrm T } / \mu _ { L } ^ { + } } } { \mathrm d p _ { \mathrm T } ^ { 2 } } = \left ( \frac { g _ { Z } ^ { \mu _ { L } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \right ) ^ { 2 } \frac { \mathrm d f _ { \gamma / \mu } } { \mathrm d p _ { \mathrm T } ^ { 2 } } \,,$$

and

$$\frac { \mathrm d f _ { [ \gamma Z _ { \mathrm T } ] / \mu _ { R } ^ { + } } } { \mathrm d p _ { \mathrm T } ^ { 2 } } = \frac { g _ { Z } ^ { \mu _ { R } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { \mathrm d f _ { \gamma / \mu } } { \mathrm d p _ { \mathrm T } ^ { 2 } } \,, \quad \frac { \mathrm d f _ { [ \gamma Z _ { \mathrm T } ] / \mu _ { L } ^ { + } } } { \mathrm d p _ { \mathrm T } ^ { 2 } } = \frac { g _ { Z } ^ { \mu _ { L } ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { \mathrm d f _ { \gamma / \mu } } { \mathrm d p _ { \mathrm T } ^ { 2 } } \,.$$

Here, we have neglected the longitudinal components of the Z boson since its parton distribution function does not have logarithmic enhancement, analogously to the longitudinal W boson. (See the first expression in Eq. (2.3).) Since the Z boson contribution starts from p 2 T ∼ m 2 Z , we can obtain the respective correction by integrating these differential equations with f Z /µ T = f [ γZ T ] /µ = 0 at p 2 T = m 2 Z as the boundary condition.

At the leading-logarithm order, we thus obtain

$$\sigma ( s ) = \sigma _ { \gamma } ( s ) + \sigma _ { Z _ { \mathrm T } } ( s ) + \sigma _ { [ \gamma Z _ { \mathrm T } ] } ( s ) \,,$$

where and

$$\sigma _ { Z _ { \Gamma } } ( s ) & = \sigma _ { \gamma } ( s ) \, g ( s ) \left [ \frac { 1 + P _ { \mu ^ { + } } } { 2 } \left ( \frac { g _ { Z } ^ { \mu ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { g _ { Z } ^ { W } } { g _ { \gamma } ^ { W } } \right ) ^ { 2 } + \frac { 1 - P _ { \mu ^ { + } } } { 2 } \left ( \frac { g _ { Z } ^ { \mu ^ { + } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { g _ { Z } ^ { W } } { g _ { \gamma } ^ { W } } \right ) ^ { 2 } \right ] \\ & \simeq \sigma _ { \gamma } ( s ) \, g ( s ) \, \left [ 1. 2 8 + 0. 2 8 1 \, P _ { \mu ^ { + } } \right ] \,,$$

$$\sigma _ { [ \gamma Z _ { T } ] } ( s ) = 2 \, \sigma _ { \gamma } ( s ) \, g ( s ) \left [ \frac { 1 + P _ { \mu ^ { + } } } { 2 } \frac { g _ { Z } ^ { \mu ^ { + } _ { R } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { g _ { Z } ^ { W } } { g _ { \gamma } ^ { W } } + \frac { 1 - P _ { \mu ^ { + } } } { 2 } \frac { g _ { Z } ^ { \mu ^ { + } _ { L } } } { g _ { \gamma } ^ { \mu ^ { + } } } \frac { g _ { Z } ^ { W } } { g _ { \gamma } ^ { W } } \right ]$$

$$\simeq \sigma _ { \gamma } ( s ) \, g ( s ) \left [ 0. 2 5 0 + 2. 2 5 \, P _ { \mu ^ { + } } \right ] \,.$$

Here, we define the ratio of logarithmic factors originating from the p 2 T integrals as

$$g ( s ) = \frac { \log ( s / m _ { Z } ^ { 2 } ) - \frac { 2 } { 3 } \log ( s / s _ { \min } ) } { \log ( s / m _ { \mu } ^ { 2 } ) - \frac { 2 } { 3 } \log ( s / s _ { \min } ) } \,,$$

which, e.g., gives overall factors of about 0.189, 0.238 and 0.268 for √ s = 2 , 10 and 30 TeV, respectively.

Thus, by adding up the contributions involving the Z boson, we finally find

$$\sigma ( s ) \, \simeq \, \sigma _ { \gamma } ( s ) \, \cdot \, \left [ 1 + \left ( 1. 5 3 + 2. 5 3 \, P _ { \mu ^ { + } } \right ) \, g ( s ) \right ],$$

which is significantly enhanced for positively polarized antimuons due to positive interference effects. Together with the enhancement in Eq. (2.5), positive polarization of the µ + beams can enhance the cross section by more than a factor of two. It is important to note that a polarized muon beam may be available, e.g., at µ + µ + colliders with ultra-slow muon technology. In the following, we present the results in both the unpolarized case with P µ + = 0, and the polarized case with P µ + = +0 8 as suggested in [5]. .

## 2.3 Numerical calculation with MadGraph

The large enhancement by (log s ) 3 motivates us to calculate the cross section more accurately and reliably, which we discuss in the next section. However, we first numerically estimate the cross sections using the event generator MadGraph, and compare them to the results of the full calculation and the cross sections of Higgs boson production via the leading W and Z boson fusion processes at µ + µ -and e + e -colliders.

We note here in particular that the γ -and Z -mediated Higgs production process also exists at µ + µ -and e + e -colliders with the same cross section as for µ + µ + colliders. Therefore, its large enhancement means that there are large corrections to the leadingorder cross section of W boson fusion at high energies.

When estimating the cross section of the photon-mediated process in MadGraph, it is possible to use a parton distribution function for the photon in one of the antimuons. We use one of these parton distribution function options, the Improved Weizs¨cker-Williams a (IWW) setting [19], for the computations presented here. Note that we do not use the parton distribution function approximation for W + from the other antimuon, and instead treat γµ + → ν µ W h + as the hard process to be convoluted with the IWW parton distribution function of the photon in the antimuon. We then take the contribution from the Z -mediated process into account by subsequently multiplying the obtained result by the factor in Eq. (2.18).

In calculations with parton distribution functions, it is necessary to set the scale for the hard process, also referred to as factorization scale µ f . Among several options implemented in MadGraph, we take two extremal choices: the default setting -1 , given by the transverse mass after k T clustering to a 2 → 2 system; and the setting 4 , given by the partonic centerof-mass energy as µ f = √ xs [20], which is a commonly used choice for the factorization scale. Note that these are settings for the dynamical scale choice in MadGraph, which is then used for the factorization scale. Here, dynamical means that the scale is determined on an event-by-event basis, e.g., via the respective momentum fraction x in µ f = √ xs of the event. One may also set a fixed factorization scale, but we do not use this here.

In Fig. 2, we show the cross section of the γ - and Z -mediated Higgs production process, given by µ + µ + → µ + ν µ W h + . For the estimate using MadGraph, the upper band represents the results with P µ + = +0 8, while the lower one corresponds to . P µ + = 0. We obtain the upper and lower ends of each band with different choices of factorization scales for the hard process with an initial-state photon, as mentioned above. We then include the contributions from Z boson exchanges by multiplying these results by the factor in Eq. (2.18) to obtain the bands shown in the figure. The solid curves near the bands are the results from the accurate calculation, which we explain in the next section. We see that the calculation with the IWW parton distribution function multiplied by the Z -boson factor gives reasonable estimates for the cross section.

The leading-logarithm formula of Eq. (2.5) with the Z -boson factor in Eq. (2.18) is depicted as lines without data points (solid and dashed lines for P µ + = +0 8 and . P µ + = 0, respectively), and seems to overestimate the cross sections by a factor of two to three even at √ s ∼ 10 TeV. This discrepancy stems largely from the sub-leading logarithmic terms ignored in Eq. (2.5). These originate partly from the lows behavior of the partonic cross section σ γW ( s ) near the threshold, which effectively shifts s min to a larger value; and partly from the terms dropped after convolution of the partonic cross section with the parton distribution functions. However, we have confirmed that including the sub-leading terms brings down the cross section, and gives numerically consistent results compared with the ones obtained from MadGraph, even for √ s ∼ O (1) TeV. Therefore, while the subleading logarithms give noticeable corrections at O (10) TeV energies, we can thus confirm the parametric (log s ) 3 behavior.

We furthermore overlaid the cross section of the leading-order W boson fusion process at µ + µ -colliders, where we do not assume beam polarizations. We see that at energies of 10 TeV or above, the cross section of the higher-order process with polarized beams

Cross Section vs. Center-of-Mass Energy for Single Higgs Production via W Boson Fusion at + +  Colliders Comparison of Different Processes, Calculation Methods, and Polarizations P +

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-42609baa-5af7-4dc9-9206-dc20e02dd52e.jpg)

Figure 2 . We show the cross section for single Higgs boson production at µ + µ + colliders via W boson fusion as a function of center-of-mass energy. We compare the results for this process obtained using different calculation methods and polarizations P µ + for both beams, as well as different processes for single Higgs production. First, we show the W boson fusion at µ + µ -colliders, via µ + µ -→ ν µ ν µ h , with unpolarized beams, i.e., P µ ± = 0 (orange, solid with diamonds). Second, the leading-logarithm approximation of Eq. (2.5) (including the Z -boson factor of Eq. (2.18)) for W boson fusion at µ + µ + colliders, via µ + µ + → µ + ν µ W h + , with both beams +0.8 polarized, i.e., P µ + = +0 8 (dark blue, solid), . and unpolarized (dark green, dashed). Third, we show the full result for W boson fusion at µ + µ + colliders with polarized (blue, solid with dots), and unpolarized (green, dashed with dots) beams. Fourth, the bands for the Improved Weizs¨cker-Williams (IWW) a calculation in MadGraph for W boson fusion at µ + µ + colliders, multiplied by the Z -boson factor of Eq. (2.18), with polarized (purple, solid with upwards-pointing triangles) and unpolarized (yellowgreen, dashed with downwards-pointing triangles) beams. The bounds of these bands are obtained by setting the dynamical scale in MadGraph to the standard setting -1 , which yields the lower bounds, and to the partonic center-of-mass energy denoted by the setting 4 , giving the upper bounds. Lastly, we also show the Z boson fusion process in µ + µ ± collisions, via µ + µ ± → µ + µ ± h , with polarized (light blue, dash-dotted with squares) and unpolarized (light green, dotted with squares) beams.

becomes about half as large as the leading one at µ + µ -colliders. See also Tab. 1, where we list the cross sections for representative collider energies. Therefore, the disadvantage due to missing W boson fusion at the leading order is remedied at the next-to-leading order at ∼ 10 TeV energies. Importantly, this large contribution to Higgs production is also present for µ + µ -and e + e -colliders. This will thus add a new significant process at any such lepton collider.

Lastly, we also show in Fig. 2 the cross section of Z boson fusion at µ + µ + colliders for P µ + = +0 8 and . P µ + = 0. We see that the cross section of the sub-leading W boson fusion process becomes larger than for the leading Z boson fusion beyond a few TeV.

The observed large correction to W boson fusion at the next-to-leading order may indicate that we are about to lose perturbativity. We thus investigate higher-order pro-

Table 1 . We show the total cross section in the unit of fb for single Higgs production via W boson fusion in µ + µ -collisions ( µ + µ -→ ν µ ν µ h ) with unpolarized beams, W boson fusion in µ + µ + collisions ( µ + µ + → µ + ν µ W h + ) with both beams +0.8 polarized or unpolarized, and via Z boson fusion in µ + µ + collisions ( µ + µ + → µ + µ + h ) with both beams +0.8 polarized or unpolarized for different center-of-mass energies.

| Center-of-Mass Energy [TeV]                     | 1      | 2      | 3      | 10     |   30 |
|-------------------------------------------------|--------|--------|--------|--------|------|
| σ ( µ + µ - → ννh ) [fb], P µ ± = 0             | 211    | 385    | 498    | 842    | 1165 |
| σ ( µ + µ + → µ + νW + h ) [fb], P µ + = +0 . 8 | 15 . 6 | 61 . 0 | 109    | 371    |  799 |
| σ ( µ + µ + → µ + νW + h ) [fb], P µ + = 0      | 7 . 05 | 27 . 7 | 49 . 9 | 172    |  374 |
| σ ( µ + µ + → µ + µ + h ) [fb], P µ + = +0 . 8  | 29 . 0 | 53 . 8 | 70 . 4 | 121    |  168 |
| σ ( µ + µ + → µ + µ + h ) [fb], P µ + = 0       | 20 . 9 | 39 . 0 | 50 . 9 | 87 . 6 |  122 |

cesses with additional powers of log s , which are, e.g., µ + µ + → µ + µ + W W h µ µ + -, + + → µ + ν µ W + γ/Z h and µ + µ + → ν µ ν µ W W h + + . In particular, the first process has a parametric (log s ) 5 dependence, while the others have (log s ) 4 , and we thus expect it to have the largest contribution at high energies among them. Using the method we describe in Sec. 3, we calculate the cross section of this process to be σ ≃ 14 fb at √ s = 20 TeV for P µ + = +0 8, which is notably smaller than the cross section . σ ≃ 600 fb of the lower-order process µ + µ + → µ + ν µ W h + . Therefore, at least for Higgs boson production at O (10) TeV energies, the higher-order processes are not yet in the non-perturbative regime. Nevertheless, in the analysis of actual events, one should include the higher-order processes since they may contribute to the signals defined by lower-order ones.

## 3 Calculation with controlled theoretical uncertainties

We saw in the previous section that a large enhancement of the cross section occurs at high energies. We furthermore examined its leading-logarithm approximation and a numerical method implementing parton distribution functions; however, it should be possible to calculate the cross section accurately without uncertainties, as we are considering a fixed-order, tree-level process. One could try to calculate the full fixed-order process in MadGraph, but as previously mentioned, in that case one encounters numerical instabilities in the integration of the phase space for small transverse momenta p T of the antimuon emitting the photon-this is why we previously opted for using a parton distribution function instead.

Here, we therefore discuss a method to calculate the cross section accurately and reliably, and with controlled theoretical uncertainties. The method is simply to split the integral over the p T of the final-state antimuon into two parts with p T &lt; p (cut) T and p T ≥ p (cut) T , respectively, for some value of p (cut) T . We calculate the first part semi-analytically using the Weizs¨cker-Williams method, also referred to as Equivalent Photon Approximation a (EPA) [9, 10, 21-24], and the second part numerically using MadGraph, since the lower

bound on p T mitigates the numerical instabilities. This means we calculate

$$\sigma = \sigma \Big | _ { p _ { T } < p _ { T } ^ { ( \text{cut} ) } } + \sigma \Big | _ { p _ { T } \geq p _ { T } ^ { ( \text{cut} ) } }$$

$$\simeq \sigma _ { \text{EPA} } \Big | _ { p _$$

where in the second line, we approximate the first term by ignoring higher-order terms in p (cut) T , and the second term, by ignoring the muon mass. While splitting the phase-space integral itself is completely independent of p (cut) T no matter its value, the premise of EPA and MadGraph being reliable in their respective regions does depend on the p (cut) T at which we split the phase space. As such, the accuracy and applicability of the different calculation methods relies on choosing a suitable value for p (cut) T . Therefore, when calculating the split total cross section using this method, we need to check whether the sum of the two cross sections of Eq. (3.2) exhibits a plateau-i.e., a constant region- as a function of p (cut) T . This is the aforementioned insensitivity of Eq. (3.1) with respect to the value of p (cut) T . If we find such a plateau, we know on one hand that the combined calculation scheme works, and on the other hand in which range of p (cut) T it is applicable. As we will see in the case of Higgs production at µ + µ + colliders, this range is physically well-motivated, which makes the application to other processes straightforward.

We note that splitting the phase space and calculating its parts with a parton distribution function in one part and the full matrix element in another part is not an entirely new approach. It has for instance been used in Ref. [25] under the name of matrix-element matching to show the validity of the parton distribution functions of transversely polarized weak bosons.

## 3.1 Lowp T cross sections from EPA

In our EPA calculation, we use the distribution function of the photon in the lepton ℓ , given by

$$f _ { \gamma / \mu } ^ { p _ { T } ^ { ( \text$$

where we put an upper bound on the transverse momentum via p (cut) T , and m ℓ is the mass of the lepton ℓ . See Appendix A for the derivation. The first expression is obtained by rewriting Eq. (20) in Ref. [19], which is originally derived for kinematical regions with small virtuality q 2 of the photon. We rewrite it using the relation of the photon virtuality and the transverse momentum, given by

$$p _ { \text{T} } ^ { 2 } = - q ^ { 2 } ( 1$$

## Cross Section  vs. p (cut) T Values for Representative Center-of-Mass Energies s Data from EPA (lowp T ) Calculation for + + + W + h

Figure 3 . We show the EPA (lowp T ) cross section for µ + µ + → µ + ν µ W h + as a function of p (cut) T for representative center-of-mass energies (top to bottom): 30 TeV (blue, solid), 10 TeV (orange, dashed), 3 TeV (green, dash-dotted), and 1 TeV (red, dotted). The vertical dashed lines indicate the muon mass m µ and the Z boson mass m Z . We calculate the error bars from the cross section value and the relative error of Eq. (3.6), combined with the supplied error of the numerical convolution, and the error supplied by MadGraph for the calculation of the partonic cross section after propagating it through the convolution.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-5c4a889c-b623-49fb-a7f9-456733bb2632.jpg)

as we also discuss in Appendix A. The second line of Eq. (3.3) is a further-approximated formula obtained when p (cut) T is much larger than the lepton mass m ℓ . This simplified formula can be found in Ref. [24], again as a function of the photon virtuality. The first expression of the distribution function consistently takes into account the finite lepton mass, up to corrections of O ( m /s 2 ℓ ). Note that although we use the first expression of Eq. (3.3) for our numerical calculations, the second expression should actually be enough for our purposes by taking p (cut) T ≫ m ℓ . An important point is that we have a non-logarithmic term, unlike the leading-logarithm formula in Eq. (2.1), for which there is no uncertainty associated with the scale of the hard process as long as the lepton mass is negligible-i.e., the second line of Eq. (3.3).

Using this parton distribution function, the total cross section for the lowp T region is given by

$$\sigma _ { \text{EPA} } \Big | _ { p _ { \text{$$

where σ γµ ( xs ) is the cross section of the γµ + → ν µ W h + process at the center-of-mass energy √ xs . The lower end of the x integration is given by s min /s , where s min is the smallest invariant mass of the final-state particles in the subprocess, i.e., s min = ( m h + m W ) 2 , since the subprocess' center-of-mass energy needs to be large enough for it to take place, with √ xs ≥ √ s min . The cross section obtained from this formula is the leading-order value of

the systematic expansion in terms of p (cut) T . The relative uncertainty of this method is of the order of

$$\begin{array} {$$

Note that EPA here is not the leading-logarithm approximation. The uncertainty above can be much smaller than [ log ( p (cut) T /m µ )] -1 , and the formula is valid for small p (cut) T as long as p (cut) T ≫ m µ . We explain the derivation of the formula and the associated uncertainties in Appendix A.

The partonic cross section σ γµ can be calculated using MadGraph since there are no intermediate light particles in this subprocess, and it can thus be evaluated numerically. Therefore, we scan over center-of-mass energies to obtain the cross section as a function of energy, and then convolute it with the parton distribution function, as described above. For a small enough p (cut) T compared to s min (and large enough compared to m µ ), the cross section can thus be calculated with small theoretical uncertainties. However, since we do not take into account the Z -boson-exchange diagrams using this method, we expect it to be reliable for p (cut) T ≪ m Z .

We show in Fig. 3 the lowp T cross section as a function of p (cut) T for representative center-of-mass energies √ s = 1, 3, 10, and 30 TeV. The error bars in the figure are calculated primarily using the relative uncertainty of Eq. (3.6), which is invisible for small p (cut) T .

## 3.2 Highp T cross sections from MadGraph

The calculation for the highp T part can be performed numerically using MadGraph by specifying a lower bound on p T , given by p (cut) T . Thusly, we can avoid the collinear divergences (or numerical instabilities) in the phase-space integral, which are associated with the muon mass being either set to zero, or significantly smaller than the center-of-mass energy.

In the default set-up, MadGraph calculates the cross section with the muon mass set to zero. This approximation results in relative uncertainties of the order of

$$\frac { m _ { \mu } ^ { 2 } } { \left ( p _ { \mathrm T } ^ { ( \text{cut} ) } \right ) ^ { 2 } } \, \left [ \log \frac { s } { \left ( p _ { \mathrm T } ^ { ( \text{cut} ) } \right ) ^ { 2 } } \, \log \frac { s } { s _ { \min } } - \frac { 2 } { 3 } \left ( \log \frac { s } { s _ { \min } } \right ) ^ { 2 } \right ] ^ { - 1 }.$$

See Appendix B for the derivation and discussion of the error associated with this massless approximation. Therefore, by taking p (cut) T ≫ m µ , we both avoid numerical instabilities, and have negligibly small uncertainties associated with the finite muon mass.

We show the obtained cross section as a function of p (cut) T in Fig. 4 for the center-of-mass energies √ s = 1, 3, 10, and 30 TeV, analogously to the lowp T cross section.

## 3.3 Summing up the low- and highp T cross sections

We show in Fig. 5 the sum of the low- and highp T cross sections, calculated using the EPA scheme and MadGraph, respectively. As expected, we clearly observe the plateau of the

Cross Section  vs. p (cut) T Values for Representative Center-of-Mass Energies s Data from MadGraph (highp T ) Calculation for + + + W + h

Figure 4 . We show the MadGraph (highp T ) cross section for µ + µ + → µ + ν µ W h + as a function of p (cut) T for representative center-of-mass energies (top to bottom): 30 TeV (blue, solid), 10 TeV (orange, dashed), 3 TeV (green, dash-dotted), and 1 TeV (red, dotted). The vertical lines indicate the muon mass m µ and the Z boson mass m Z . We calculate the error bars from the cross section value and the relative error of Eq. (3.7), combined with the error supplied by MadGraph for the calculation of the cross section.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-b1d96dbc-e90b-4060-94e4-378459952083.jpg)

summed cross sections as a function of p (cut) T . Based on the discussions of the uncertainties in both calculations, we find an optimal choice for p (cut) T to be around the geometric mean of the muon and Z boson masses, p (cut) T = √ m m µ Z ≃ 3 10 GeV, where we indeed find the . plateau. Nevertheless, for too small or too large p (cut) T , the summed cross sections deviate from their plateau values.

On the one hand, for small p (cut) T (i.e., ≲ m µ ), the cross section calculated by MadGraph is too large due to numerical instabilities, as previously mentioned. However, the error bars in the small p (cut) T region represent the uncertainties from the massless approximation, whereas the stark increase of the cross section calculated by MadGraph seems to be largely independent of this approximation. Indeed, this behavior persists even when including the finite muon mass in MadGraph. This is because the muon mass, while non-zero, is significantly smaller than the center-of-mass energies considered, leading again to numerical instabilities. Thus, even when including a finite muon mass, one cannot obtain a reliable result from MadGraph for, e.g., p (cut) T = 0.

On the other hand, for large p (cut) T (i.e., ≳ m Z ), the contribution from Z -bosonexchange diagrams, as well as the higher-order terms in the p (cut) T expansion of EPA become important for the lowp T cross section. Neglecting these thus leads to a deficit in the summed cross section. Note that there are also diagrams without collinear divergences, which are taken into account in the MadGraph (highp T ) cross section. Increasing p (cut) T

Cross Section  vs. p (cut) T Values for Representative Center-of-Mass Energies s Sum of Data from EPA (lowp T ) and MadGraph (highp T ) Calculations for + + + W + h

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-00c7b9d9-f855-4e83-8330-3b0765bc58d9.jpg)

Figure 5 . We show the sum of the EPA (lowp T ) and MadGraph (highp T ) cross sections for µ + µ + → µ + ν µ W h + as a function of p (cut) T for representative center-of-mass energies (top to bottom): 30 TeV (blue, solid), 10 TeV (orange, dashed), 3 TeV (green, dash-dotted), and 1 TeV (red, dotted). The vertical lines indicate the muon mass m µ and the Z boson mass m Z . The error bars are determined by the root-square-sum of the EPA and MadGraph errors, which are primarily obtained from Eqs. (3.6) and (3.7), respectively. The violet rectangle marks the approximate region of the plateau, i.e., the p (cut) T -independent region of the summed cross section, which we take as 0 2 GeV . ≲ p (cut) T ≲ 20 GeV.

too far thus also neglects part of their contribution.

Note that the results discussed here for unpolarized beams apply in the same way to the case with polarized beams.

By reading off the cross sections in the plateau region, we can now reliably calculate the fixed-order cross section of the µ + µ + → µ + ν µ W h + process as shown in Fig. 6. The upper and lower curves are the cross sections with beam polarizations P µ + = +0 8 and . P µ + = 0, respectively. We also show error bands associated with the choice of p (cut) T = √ m m µ Z for the range √ m m / µ Z 4 &lt; p (cut) T &lt; 4 √ m m µ Z . However, they are almost invisible. We thereby estimate the uncertainty associated with the choice of p (cut) T to be less than one percent.

For convenience, we present in Appendix C a semi-automatic method to implement the p (cut) T scheme in MadGraph, in particular the lowp T calculation.

## 4 Event shape of the W boson fusion process

As seen so far, we obtained reliable results for the higher-order Higgs production by dividing the cross section into two regions, corresponding to low and high p T of the final-state antimuon, which we calculate using EPA and MadGraph, respectively. Compared to the leading-order W boson fusion, this process contains an additional W boson in the final state. Thus, to distinguish event signals from background, the produced Higgs and W bosons should both be detected.

## Cross Section vs. Center-of-Mass Energy for + + + W + h Full Result with Different Polarizations P +  for p (cut) T = m mZ 3.10 GeV  with p (cut) T Error Bands

Figure 6 . We show the full, summed cross section as a function of center-of-mass energy for p (cut) T = √ m m µ Z ≃ 3 10 GeV, the geometric mean of the muon mass . m µ and the Z boson mass m Z . We obtain the error bands by varying p (cut) T in the range √ m m / µ Z 4 &lt; p (cut) T &lt; 4 √ m m µ Z . The upper curve (blue, dashed) and its error band (red, diagonal lines) correspond to a polarization of P µ + = +0 8 for both beams, and the lower curve (green, dash-dotted) and its band (magenta, . solid) to P µ + = 0, i.e., unpolarized beams.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-f625a6fa-0378-41b3-a502-9f7bd4310c58.jpg)

However, because the µ + beams induce undesired backgrounds such as µ + -decay products and incoherent e + e -pairs (see Ref. [26] for a review), we need to put shielding nozzles around the interaction point to reduce these beam-induced backgrounds. This makes placement of detectors quite non-trivial, leading to constraints of the detectors' coverage angle. Therefore, to estimate event efficiency, it is necessary to study the scattering angles of final-state particles and to compare them with the coverage angle.

In the following, we present event histograms of the distribution of the scattering angles. As one of the main decay channels, we consider the produced h and W + decaying into bb ¯ , and light quarks ( ud ¯ ), respectively. We generate the events using MadGraph with p T ≥ 0 1 GeV, which corresponds to the high-. p T region of the previous sections. Here, we neglect the lowp T region because it does not have a significant number of events for this choice of lowerp T cut. We show in Figs. 7 and 8 the histograms for the intermediate h and W + bosons and their decay products b , ¯ , b u , and ¯ for d √ s = 2 TeV and 10 TeV, respectively. The histograms are at the parton level without smearing of the energies or momenta of the final-state particles. Since a µ + µ + collider is symmetric in the two initial-state particles, we concentrate on the angles 0 ◦ ≤ θ ≤ 90 . ◦

From these plots, we observe that most decay products have sufficiently large scattering angles, or transverse momenta, and are hence visible. We therefore now estimate how many events can be caught by detectors. Although there is not yet a concrete detector design for µ + µ + colliders, we assume the coverage of the hadron calorimeter to be around 10 ◦ [27, 28]. Thus, we find that around 83% of events can be caught by the detector, and reconstructing

Figure 7 . We show event histograms for the scattering angles of final-state particles in the W boson fusion process with √ s = 2 TeV, where the calculation is performed by MadGraph with a minimump T cut 0 1 GeV. The number of the generated events is 5 . × 10 . 4

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-b2b7abc4-b751-4538-b96e-e55e447410d6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-fc14b7b1-f98b-4fd5-90c2-d4c18046fcf0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-205c2430-4d00-499f-969f-f99f8c60e7eb.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-adc48854-a979-4c34-8b09-d7607148690c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-b6967fa1-a2e4-42be-ba7d-98bbcc09ebe5.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-c0a7f924-76ff-4b29-9c30-1249f2193a1c.jpg)

the W and Higgs bosons is possible for both √ s = 2 TeV and 10 TeV. We expect this reconstruction to help to distinguish these signal events from background events. Further studies, such as detector simulation, will be performed elsewhere.

## 5 Summary

We studied Higgs boson production at µ + µ + colliders, where W boson fusion is not possible at the leading order. This is because a µ + cannot directly emit the W -necessary to create

Figure 8 . We show event histograms for scattering angles of final-state particles in the W boson fusion process with √ s = 10 TeV, where the calculation is performed by MadGraph with a minimump T cut 0 1 GeV. The number of the generated events is 5 . × 10 . 4

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-7ba04cee-aa93-4950-808c-5cfdf3d6b024.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-a4708f94-3a88-4fad-8e4a-aeed247bb098.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-097df63f-8c88-4176-855d-877d918d400d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-814d6a16-1eef-432f-aace-21bea56e1f1c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-38ee1aac-4a3b-4f3c-94c0-41e94ec9aaea.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.01068v1-figure-a6b78f5c-f913-40d6-87b8-6fb0381d7a0e.jpg)

the W W + -pair fusing to a Higgs boson. Nevertheless, we show that at high-energy µ + µ + colliders, a γ -and Z -mediated W boson fusion is possible, with its cross section growing as (log s ) 3 with the center-of-mass energy √ s . This means that at high energies, such colliders can produce almost as many Higgs bosons as µ + µ -or e + e -colliders, since their production cross section only increases as log s .

We calculated the cross section by carefully treating the divergence associated with collinear emission of the intermediate photon from the antimuon. Therein, we split the

integration region of the transverse momentum p T of the antimuon into the sum of a lowand a highp T region. We then calculate the lowp T part using the equivalent photon approximation (EPA), and the highp T part directly using MadGraph. EPA gives the leading-order value in a systematic expansion in the maximal value of p T , and Z boson contributions to the lowp T region become important for p T ≳ m Z . Therefore, we can control the systematic uncertainties associated with this approximation by choosing a small enough value p (cut) T at which we split the low- and highp T regions. On the other hand, to avoid the collinear divergences or numerical instabilities from MadGraph, we must have p (cut) T ≳ m µ . Thereby, we obtained a reliable result for the cross section by verifying its independence on the p T value at which we split the regions. This is well satisfied, e.g., for the choice p (cut) T = √ m m µ Z .

The cross section we obtained in this manner is indeed enhanced at high energy; note, however, that the leading-logarithm formula that only includes (log s ) 3 terms seems to overestimate it. Nevertheless, even compared with the leading-order W boson fusion process at µ + µ -colliders, the cross section of the γ - and Z -mediated process with polarized beams can be as large as half of the one for the leading-order process at 10 TeV energies. Therefore, while beam polarization capabilites are important, a µ + µ + collider is as good a Higgs boson factory as a µ + µ -collider.

## Acknowledgements

We would like to thank Koji Nakamura, Sayuka Kita and Mitsuhiro Yoshida for useful discussions. This work was in part supported by JSPS KAKENHI Grant Numbers JP22K21350 (R.K., R.M. and S.O.), JP21H01086 (R.K.), JP19H00689 (R.K.), JP24KJ1157 (R.T.), JP19K14711 (H.T.) and JP23K13110 (H.T.). H.T. is the Yukawa Research Fellow supported by Yukawa Memorial Foundation. This work is also supported by the Deutsche Forschungsgemeinschaft under Germany's Excellence Strategy - EXC 2121 Quantum Universe - 390833306.

## A Equivalent photon approximation (EPA)

Repeating a similar calculation to Ref. [19], we derive the EPA formula applicable to the case where a p T cut is applied. This section includes a review of Ref. [19]. Let us consider the photon-mediated process ℓ p ( ) + f ( k ) → ℓ p ( ′ ) + X , where ℓ , f , and X denote a lepton, massless parton ( k 2 = 0), and a generic final-state system, respectively. The four-momenta of the leptons and the virtual photon are defined as

$$p = E ( 1, 0, 0, \beta ), \ \ p ^ { \prime } = E ^ { \prime } ( 1, 0, \beta ^ { \prime } \sin \theta, \beta ^ { \prime } \cos \theta ) \ \text{ and } \ q = p - p ^ { \prime } \text{ \quad \ \ (A.1)}$$

with

$$\beta = \sqrt { 1 - \frac { m _ { \ell } ^ { 2 } } { E ^ { 2 } } } \quad \text{and} \quad \beta ^ { \prime } = \sqrt { 1 - \frac { m _ { \ell } ^ { 2 } } { E ^ { \prime 2 } } } \,,$$

and m ℓ the lepton mass.

The cross section is given by

$$\mathrm d \sigma _ { \text{EPA} } = \frac { 1 } { 8 k \cdot p } \frac { e ^ { 2 } W ^ { \mu \nu } T _ { \mu \nu } } { q ^ { 4 } } \frac { \mathrm d ^ { 3 } p ^ { \prime } } { ( 2 \pi ) ^ { 3 } 2 E ^ { \prime } } = \frac { 1 } { 8 k \cdot p } \frac { \alpha } { 4 \pi } \frac { W ^ { \mu \nu } T _ { \mu \nu } } { q ^ { 4 } } \, \mathrm d q ^ { 2 } \mathrm d x \,,$$

where α is the fine-structure constant, and the two independent kinematic variables, photon virtuality and longitudinal momentum fraction, are given by

$$q ^ { 2 } = 2 m _ { \ell } ^ { 2 } - 2 E E ^ { \prime } ( 1 - \beta \beta ^ { \prime } \cos \theta ) \ \text{ and } \ x = 1 - \frac { E ^ { \prime } ( 1 + \beta ^ { \prime } \cos \theta ) } { E ( 1 + \beta ) } \,, \quad \quad ( \text{A.4} )$$

respectively. Furthermore, x satisfies x = 1 -k p · ′ k p · = k q · k p · . The leptonic tensor T µν and the hadronic tensor W µν of the electromagnetic current are given by

$$T _ { \mu \nu } = 4 \left ( \frac { 1 } { 2 } q ^ { 2 } g _ { \mu \nu } + p _ { \mu } p ^ { \prime } _ { \nu } + p _ { \nu } p ^ { \prime } _ { \mu } \right )$$

and

$$W ^ { \mu \nu } & = W _ { 1 } ( q ^ { 2 }, k \cdot q ) \left ( - g ^ { \mu \nu } + \frac { q ^ { \mu } q ^ { \nu } } { q ^ { 2 } } \right ) - \frac { q ^ { 2 } W _ { 2 } ( q ^ { 2 }, k \cdot q ) } { ( k \cdot q ) ^ { 2 } } \left ( k ^ { \mu } - \frac { k \cdot q } { q ^ { 2 } } q ^ { \mu } \right ) \left ( k ^ { \nu } - \frac { k \cdot q } { q ^ { 2 } } q ^ { \nu } \right ).$$

Contracting them then yields

$$W ^ { \mu \nu } T _ { \mu \nu } & = - 4 \left [ 2 m _ { \ell } ^ { 2 } W _ { 1 } ( q ^ { 2 }, k \cdot q ) + q ^ { 2 } \left ( W _ { 1 } ( q ^ { 2 }, k \cdot q ) + \frac { 2 ( 1 - x ) } { x ^ { 2 } } W _ { 2 } ( q ^ { 2 }, k \cdot q ) \right ) \right ] \\ & \simeq - 4 W _ { 1 } ( 0, k \cdot q ) \left [ 2 m _ { \ell } ^ { 2 } + q ^ { 2 } \cdot \frac { 1 + ( 1 - x ) ^ { 2 } } { x ^ { 2 } } \right ],$$

where we use

$$W _ { 1 } ( q ^ { 2 }, k \cdot q ) = W _ { 1 } ( 0, k \cdot q ) + \mathcal { O } ( q ^ { 2 } ) \ \text{ and } \ W _ { 2 } ( q ^ { 2 }, k \cdot q ) = W _ { 1 } ( 0, k \cdot q ) + \mathcal { O } ( q ^ { 2 } ) \,, \quad ( \text{A.8} )$$

obtained by requiring that W µν be analytic in q 2 as q 2 → 0.

Therefore, the cross section can be approximated by

$$d \sigma _ { \text{EPA} } = f _ { \gamma / \ell } ( x ) \sigma _ { \gamma f } ( q, k ) d x \,,$$

with [19]

$$f _ { \gamma / \ell } ( x ) = \frac { \alpha } { 2 \pi } \left [ 2 m _ { \ell } ^ { 2 } x \left ( \frac { 1 } { q _ { \max } ^ { 2 } } - \frac { 1 } { q _ { \min } ^ { 2 } } \right ) + \frac { 1 + ( 1 - x ) ^ { 2 } } { x } \log \frac { q _ { \min } ^ { 2 } } { q _ { \max } ^ { 2 } } \right ] \quad \text{(A.10)}$$

and

$$\sigma _ { \gamma f } ( q, k ) = \frac { W _ { 1 } ( 0, k \cdot q ) } { 4 k \cdot q } \,.$$

We obtain these expressions by integrating Eq. (A.3) with Eq. (A.7) over q 2 from q 2 min to q 2 max .

Ref. [19] further presents the formula for the case where a θ cut is applied, with 0 ≤ θ &lt; θ c with θ c ≪ 1. In this case, the range of q 2 reads [19]

$$q _ { \max } ^ { 2 } = - \frac { m _ { \ell } ^ { 2 } x ^ { 2 } } { 1 - x } \quad \text{and} \quad q _ { \min } ^ { 2 } = - \frac { m _ { \ell } ^ { 2 } x ^ { 2 } } { 1 - x } - E ^ { 2 } ( 1 - x ) \theta _ { c } ^ { 2 } \,.$$

Substituting these results into Eq. (A.10), one obtains f γ/ℓ ( x ) in the θ -cut scheme.

Instead of a θ -cut, we give f γ/ℓ ( x ) for the case where a p T cut is applied, with 0 ≤ p T &lt; p (cut) T . In particular, we find a simple relation between p 2 T and q 2 :

$$p _ { \text{T} } ^ { 2 } = - q ^ { 2 } ( 1 - x ) - m _ { \ell } ^ { 2 } x ^ { 2 } \,,$$

as seen from p T = E β ′ ′ sin θ . This is especially useful, since it does not depend on E , in contrast to the case with θ . Therefore, the range of q 2 reads

$$q _ { \max } ^ { 2 } = - \frac { m _ { \ell } ^ { 2 } x ^ { 2 } } { 1 - x } \quad \text{and} \quad q _ { \min } ^ { 2 } = - \frac { 1 } { 1 - x } \left ( \left ( p _ { \text{T} } ^ { ( \text{cut} ) } \right ) ^ { 2 } + m _ { \ell } ^ { 2 } x ^ { 2 } \right ).$$

Substituting these results into Eq. (A.10), we obtain f γ/ℓ ( x ) in the p T cut scheme, as shown in Eq. (3.3).

Next, we discuss the uncertainty of the EPA formula. The EPA calculation utilizes the leading term of the expansion of W q ,k i ( 2 · k ) in q 2 :

$$W _ { i } ( q ^ { 2 }, k \cdot q ) = \sum _ { n \geq 0 } A _ { i } ^ { ( n ) } ( k \cdot q ) \left ( \frac { q ^ { 2 } } { ( k \cdot q ) } \right ) ^ { n }$$

for i = 1, 2. Therefore, we study the impact of higher-order terms by again following the corresponding calculation in Ref. [19]. Due to

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \\ W ^ { \mu \nu } T _ { \mu \nu } = & - 4 W _ { 1 } ( 0, k \cdot q ) \left [ 2 m _ { \ell } ^ { 2 } + q ^ { 2 } \frac { 1 + ( 1 - x ) ^ { 2 } } { x ^ { 2 } } \right ] - 8 m _ { \ell } ^ { 2 } \left ( A _ { 1 } ^ { ( 1 ) } \frac { q ^ { 2 } } { k \cdot q } + A _ { 1 } ^ { ( 2 ) } \frac { q ^ { 4 } } { ( k \cdot q ) ^ { 2 } } \right ) \\ & - 4 q ^ { 2 } \left ( A _ { 1 } ^ { ( 1 ) } \frac { q ^ { 2 } } { k \cdot q } + \frac { 2 ( 1 - x ) } { x ^ { 2 } } A _ { 2 } ^ { ( 1 ) } \frac { q ^ { 2 } } { k \cdot q } \right ) + \mathcal { O } ( q ^ { 6 } ) \,, \\. & \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots$$

the error caused by the higher-order terms is given by

## ∆d σ

$$& \Delta d \sigma _ { \text{EPA} } \\ & \simeq - \frac { \alpha } { 2 \pi } \frac { 1 } { 4 k \cdot p } \left \{ \frac { 2 m _ { \ell } ^ { 2 } A _ { 1 } ^ { ( 1 ) } } { k \cdot q } \frac { 1 } { q ^ { 2 } } + \frac { 2 m _ { \ell } ^ { 2 } A _ { 1 } ^ { ( 2 ) } } { ( k \cdot q ) ^ { 2 } } + \frac { A _ { 1 } ^ { ( 1 ) } } { k \cdot q } + \frac { 2 ( 1 - x ) } { x ^ { 2 } } \frac { A _ { 2 } ^ { ( 1 ) } } { k \cdot q } \right \} d q ^ { 2 } d x \\ & \simeq - \frac { \alpha } { 2 \pi } \left \{ \frac { m _ { \ell } ^ { 2 } } { s } \frac { A _ { 1 } ^ { ( 1 ) } } { k \cdot q } \log \frac { q _ { \max } ^ { 2 } } { q _ { \min } ^ { 2 } } + \left [ \frac { 2 m _ { \ell } ^ { 2 } A _ { 1 } ^ { ( 2 ) } } { s } \frac { 1 } { k \cdot q } + \frac { 1 } { 2 } \frac { A _ { 1 } ^ { ( 1 ) } } { k \cdot q } + \frac { A _ { 2 } ^ { ( 1 ) } } { k \cdot q } \frac { 1 - x } { x ^ { 2 } } \right ] \frac { q _ { \max } ^ { 2 } - q _ { \min } ^ { 2 } } { s } \right \} d x \,.$$

In the second equality, we performed the q 2 -integral, then used the relation

$$k \cdot q = x k \cdot p = \frac { x s } { 2 } + \mathcal { O } \left ( \frac { m _ { \ell } ^ { 2 } } { s } \right ),$$

and finally combined 1 / k ( · q ) with A ( n ) i to take

$$\frac { A _ { i } ^ { ( n ) } } { k \cdot q } \sim \mathcal { O } ( 1 ) \times \sigma _ { \gamma f }.$$

This identification is analogous to what we did for the n = 0 case. Noting that

$$q _ { \max } ^ { 2 } - q _ { \min } ^ { 2 } = \begin{cases} & E ^ { 2 } ( 1 - x ) \theta _ { c } ^ { 2 } & \text{for a $\theta$ cut}, \\ & \frac { 1 } { 1 - x } \left ( p _ { T } ^ { ( \text{cut} ) } \right ) ^ { 2 } & \text{for a $p_{T}$ cut}, \end{cases}$$

and also

$$\int _ { x _ { \min } } \frac { d x } { x ^ { 2 } } \sim \frac { 1 } { x _ { \min } } \,,$$

where x min is given by s min /s , we see that the dominant error stems from the last term inside the square brackets of Eq. (A.17), and we thus find:

$$\frac { \Delta \sigma _ { \text{EPA} } \Big | _ { p _ { \text{T} } < p _ { \text{T} } ^ { ( \text{cut} ) } } } { \sigma _ { \text{EPA} } \Big | _ { p _ { \text{T} } < p _ { \text{T} } ^ { ( \text{cut} ) } } } \sim \begin{cases} \mathcal { O } ( E ^ { 2 } \theta _ { c } ^ { 2 } / s _ { \min } ) & \text{for a } \theta \text{ cut}, \\ \mathcal { O } ( ( p _ { \text{T} } ^ { ( \text{cut} ) } ) ^ { 2 } / s _ { \min } ) & \text{for a } p _ { \text{T} } \text{ cut}. \end{cases}$$

## B Error in the massless lepton approximation

We now discuss the error in the calculation with MadGraph for the highp T cross section in Eq. (3.2) when the muon mass is neglected. For this, we split σ MG ∣ ∣ p T ≥ p (cut) T as

$$\sigma _ { \text{MG} } \Big | _ { p _ { \text{T} } \geq p _ { \text{T} } ^ { ( \text{cut} ) } } = \sigma _ { \text{MG} } \Big | _ { p _ { \text{T} } ^ { * } \geq p _ { \text{T} } \geq p _ { \text{T} } ^ { ( \text{cut} ) } } + \sigma _ { \text{MG} } \Big | _ { p _ { \text{T} } > p _ { \text{T} } ^ { * } }$$

by using an additional cut p ∗ T . The dominant error caused by the massless approximation should stem from the first cross section of the sum, where p T is smaller. We then set the new cut p ∗ T such that the EPA formula is also valid for the cross section with p ∗ T ≥ p T ≥ p (cut) T . (The plateaus we observe in Fig. 5 indicate that such a p ∗ T exists.) Thus, we can estimate the error by

$$\Delta \sigma _ { \text{MG} } \Big | _ { p _ { T } \geq p _ { T } ^ { ( \text{cut)} } } & \simeq \sigma _ { \text{MG} } \Big | _ { p _ { T } ^ { * } \geq p _ { T } \geq p _ { T } ^ { ( \text{cut)} }, m _ { \mu } ^ { 2 } = 0 } \\ & \simeq \sigma _ { \text{EPA} } \Big | _ { p _ { T } ^ { * } \geq p _ { T } \geq p _ { T } ^ { ( \text{cut)} } } - \sigma _ { \text{EPA} } \Big | _ { p _ { T } ^ { * } \geq p _ { T } \geq p _ { T } ^ { ( \text{cut)} }, m _ { \mu } ^ { 2 } = 0 } \\ & = \int d x \left [ f _ { \gamma / \mu } ^ { p _ { T } ^ { ( \text{cut)} }, p _ { T } ^ { * } } ( x ) - f _ { \gamma / \mu } ^ { p _ { T } ^ { ( \text{cut)} }, p _ { T } ^ { * } } ( x ) \Big | _ { m _ { \mu } ^ { 2 } = 0 } \right ] \sigma _ { \gamma \mu } ( x s ) \,, \\ & \text{(cut)} \ \.$$

where the function f p (cut) T , p ∗ T γ/µ ( x ) denotes the EPA parton distribution function in Eq. (A.10) with

$$q _ { \min } ^ { 2 } = - \frac { 1 } { 1 - x } \left ( p _ { \text{T} } ^ { * \, 2 } + m _ { \mu } ^ { 2 } x ^ { 2 } \right ) \quad \text{and} \quad q _ { \max } ^ { 2 } = - \frac { 1 } { 1 - x } \left ( \left ( p _ { \text{T} } ^ { ( \text{cut} ) } \right ) ^ { 2 } + m _ { \mu } ^ { 2 } x ^ { 2 } \right ), \quad \text{(B.)}$$

and the function f p (cut) T , p ∗ T γ/µ ( x ) ∣ ∣ ∣ m 2 µ =0 denotes the analogous expression with m µ → 0.

Assuming m 2 µ ≪ ( p (cut) T ) 2 ≪ p ∗ 2 T , the difference of these functions is

$$& \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \julia \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \lq \\ & \quad f _ { \gamma / \mu } ^ { p _ { T } ^ { ( \text{cut)} }, p _ { T } ^ { * } } ( x ) - f _ { \gamma / \mu } ^ { p _ { T } ^ { ( \text{cut)} }, p _ { T } ^ { * } } ( x ) \Big | _ { m _ { \mu } ^ { 2 } \rho _ { \mu } ^ { 2 } } \left ( \frac { 1 - x } { m _ { \mu } ^ { 2 } x ^ { 2 } + p _ { T } ^ { * } } - \frac { 1 - x } { m _ { \mu } ^ { 2 } x ^ { 2 } + ( p _ { T } ^ { ( \text{cut)} } ) ^ { 2 } } \right ) + \frac { 1 + ( 1 - x ) ^ { 2 } } { x } \log \frac { 1 + m _ { \mu } ^ { 2 } x ^ { 2 } / p _ { T } ^ { * } } { 1 + m _ { \mu } ^ { 2 } x ^ { 2 } / \Big ( p _ { T } ^ { ( \text{cut)} } ) ^ { 2 } } \right ) ^ { 2 } \right ] \\ & \simeq - \frac { m _ { \mu } ^ { 2 } } { \Big ( p _ { T } ^ { ( \text{cut)} ) ^ { 2 } } \frac { \alpha } { 2 \pi } x ( 2 - x ) ^ { 2 }, } \\ & \quad \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \ rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \rho \\ \end{matrix}$$

and therefore the error size is controlled by m / p 2 µ ( (cut) T ) 2 . Furthermore, we observe that Eq. (B.4) is not singular at x = 0, and loses the logarithmic factor present in f p (cut) T γ/µ ( x ). This indicates that Eq. (B.2) has two fewer logarithmic factors compared with σ ∣ ∣ p T ≥ p (cut) T . Noting also that the infrared cutoff for σ ∣ ∣ p T ≥ p (cut) T is p (cut) T , rather than m µ , we hence estimate the relative uncertainty by

$$\frac { \Delta \sigma _ { \text{MG} } \Big | _ { p _ { \text{T} } \geq p _ { \text{T} } ^ { ( \text{cut} ) } } } { \sigma \Big | _ { p _ { \text{T} } \geq p _ { \text{T} } ^ { ( \text{cut} ) } } } \sim \frac { m _ { \mu } ^ { 2 } } { \left ( p _ { \text{T} } ^ { ( \text{cut} ) } \right ) ^ { 2 } } \left [ \log \frac { s } { \left ( p _ { \text{T} } ^ { ( \text{cut} ) } \right ) ^ { 2 } } \log \frac { s } { s _ { \min } } - \frac { 2 } { 3 } \left ( \log \frac { s } { s _ { \min } } \right ) ^ { 2 } \right ] ^ { - 1 }.$$

## C Semi-automatic MadGraph implementation of p (cut) T calculation

As previously mentioned, we also implemented a semi-automatic calculation of the p (cut) T scheme for convenience. Here, semi-automatic means that both the low- and highp T cross sections can be calculated in MadGraph directly; however, the two parts still need to be summed by the user. In particular, we implement the calculation of the lowp T cross section using EPA, since the highp T part can already be computed directly. In this appendix, we therefore discuss how the calculation is implemented, how to use it, and its accuracy.

First, we briefly mention how MadGraph calculates EPA, or IWW, via the iww setting in the run card. With this setting, MadGraph uses the parton distribution function formula of Eq. (A.10) with q 2 max = -m x / 2 ℓ 2 (1 -x ), as given in Eq. (A.12); we also note here that the lepton mass in q 2 max is hard-coded, and specified through the collider type set in the run card of the respective run. However, q 2 min is not given by a fixed expression, but is instead provided to the parton distribution function as external input, together with the momentum fraction x . In particular, q 2 min takes on the role of (squared) factorization scale, which is set either by a fixed scale for all events or a dynamical scale determined on an event-by-event basis, depending on the user specification in the run card. This is because the factorization scale is the largest scale included in the parton distribution function description, whose square is precisely the maximum absolute value of the photon virtuality, ∣ ∣ q 2 min ∣ ∣ . Thus, based on the dynamical scale choice specified by the user, q 2 min is calculated by MadGraph for each event, and then passed on to the IWW parton distribution function. For instance, the setting 4 sets the dynamical (factorization) scale to the partonic center-of-mass energy, which is a common choice, and thus ∣ ∣ q 2 min ∣ ∣ = xs . Note that this is calculated using the

four-momenta of the incoming partons of the respective event; and if applicable, x may be a product of momentum fractions of the incoming partons.

Since MadGraph supports user-defined dynamical scales, we implement the parton distribution function of Eq. (3.3) through such a custom dynamical scale, which can then be imported to MadGraph and used for the lowp T , EPA calculation. The implementation is partly based on the original MadGraph code, and works in the following way: First, we import the collider information, in particular the beam energies and beam types. For this implementation, the beams are restricted to electrons, muons, and their respective anti-particles. Second, we import the kinematical information of the respective parton event and define the quantities necessary for the calculation. Therein, the value of p (cut) T is hard-coded, and thus needs to be changed by hand depending on which value is required for the process of interest. Third, we calculate and set the pertinent variables from the kinematical data. In particular, we set the lepton mass based on the information from the run card, and then calculate the momentum fraction x using

$$x = \frac { \hat { s } } { s } + \mathcal { O } \left ( \frac { m _ { \text{beam} } ^ { 2 } } { s } \right ) = \frac { \left ( p _ { \text{parton} } 1 + p _ { \text{parton} } 2 \right ) ^ { 2 } } { \left ( E _ { \text{beam} } 1 + E _ { \text{beam} } 2 \right ) ^ { 2 } } + \mathcal { O } \left ( \frac { m _ { \text{beam} } ^ { 2 } } { s } \right ), \quad \text{(C.1)}$$

where √ s is the center-of-mass energy of the collider, √ ˆ that of s the partonic process, m beam1 2 , are the masses of the beam particles set in the run card, p parton 1 2 , are the fourmomenta of the incoming partons, and E beam1 2 , are the respective beam energies. Here, we used that ˆ is determined by the incoming partons' four-momenta and is equal to the s product xs , up to corrections of O ( m 2 beam1 2 , /s ). Lastly, we calculate the dynamical scale via Eq. (A.14) as

$$( \text{dynamical scale} ) = \sqrt { | q _ { \min } ^ { 2 } | } = \sqrt { \frac { 1 } { 1 - x } \left ( \left ( p _ { \text{T} } ^ { ( \text{cut} ) } \right ) ^ { 2 } + m _ { \ell } ^ { 2 } x ^ { 2 } \right ) } \,, \quad \text{(C.2)}$$

and return the result, which will then be used in the iww parton distribution function.

To use this implementation in practice, the following steps are necessary. First, save the code as .f (Fortran) text file on your machine, and copy the file path. Second, adjust the p (cut) T value in GeV according your process of interest and save the changes; the default setting is the geometric mean of the muon and Z boson mass, √ m m µ Z ≃ 3 10 GeV, since . this was the pertinent value for our considerations. The corresponding variable in the file is denoted as pt\_cut\_iww and marked by a box around it for ease of use. Third, run MadGraph as usual, set the photon-emitting beam as ± 3 for e ∓ or ± 4 for µ ∓ , and the parton distribution function for the photon to iww . Fourth, set the dynamical\_scale\_choice parameter to 0 , which means that a custom dynamical scale is to be used. Note, in particular, that no fixed factorization scale should be set. Lastly, set the parameter custom\_fcts to the file path of the .f file containing the code to calculate √∣ ∣ q 2 min ∣ ∣ . After this, the run can be continued in the usual manner.

Note that to obtain the correct lowp T cross section for the case of identical beam particles, the obtained result needs to be multiplied by a factor of two, as is also the case for a manual convolution.

In terms of accuracy, the results thereby obtained by the MadGraph implementation of the lowp T cross section are consistently about 5% smaller on average, when compared with the manual calculation for our Higgs production process at center-of-mass energies up to 30 TeV. The deviation is a bit larger for the lower center-of-mass energies, and a bit smaller for the higher ones, with a respective O (1)% change in both directions. Nevertheless, we observe an approximate plateau for similar p (cut) T values as in the manual calculation, irrespective of beam polarizations and energies. Since the size of the deviation of the cross section is largely independent of beam polarizations and center-of-mass energies, we assume that it is a systematic error. Numerically, it can be mitigated by adjusting the scalefact parameter in the run card to the value ∼ 1 3, . which modifies the dynamical scale as (new dynamical scale) = scalefact × (old dynamical scale) . The results obtained from this modification are within about ± 0 4% of the manual calculation, across center-of-. mass energies and polarizations.

Finally, note that due to the calculation of x as ˆ s/s , and since the dynamical scale is calculated based on Eq. (A.14) specifically for the photon, we recommend using partonic processes where a parton distribution function is only used for the photon.

## References

- [1] D. Buttazzo, R. Franceschini and A. Wulzer, Two Paths Towards Precision at a Very High Energy Lepton Collider , JHEP 05 (2021) 219 [ 2012.11555 ].
- [2] H. Al Ali et al., The muon Smasher's guide , Rept. Prog. Phys. 85 (2022) 084201 [ 2103.14043 ].
- [3] K. M. Black et al., Muon Collider Forum report , JINST 19 (2024) T02015 [ 2209.01318 ].
- [4] K. Nagamine, Y. Miyake, K. Shimomura, P. Birrer, J. P. Marangos, M. Iwasaki et al., Ultraslow Positive-Muon Generation by Laser Ionization of Thermal Muonium from Hot Tungsten at Primary Proton Beam , Phys. Rev. Lett. 74 (1995) 4811.
- [5] Y. Hamada, R. Kitano, R. Matsudo, H. Takaura and M. Yoshida, µ TRISTAN PTEP , 2022 (2022) 053B02 [ 2201.06664 ].
- [6] M. Abe et al., A New Approach for Measuring the Muon Anomalous Magnetic Moment and Electric Dipole Moment , PTEP 2019 (2019) 053C02 [ 1901.03047 ].
- [7] ILC collaboration, The International Linear Collider Technical Design Report - Volume 2: Physics , 1306.6352 .
- [8] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations , JHEP 07 (2014) 079 [ 1405.0301 ].
- [9] C. F. v. Weizs¨cker, a Ausstrahlung bei st¨ßen sehr schneller elektronen o , Zeitschrift f¨r Physik u 88 (1934) 612.
- [10] E. J. Williams, Correlation of certain collision problems with radiation theory , Kong. Dan. Vid. Sel. Mat. Fys. Med. 13N4 (1935) 1.
- [11] S. Dawson, The Effective W Approximation , Nucl. Phys. B 249 (1985) 42.

- [12] G. L. Kane, W. W. Repko and W. B. Rolnick, The Effective W ± , Z 0 Approximation for High-Energy Collisions , Phys. Lett. B 148 (1984) 367.
- [13] J. Chen, T. Han and B. Tweedie, Electroweak Splitting Functions and High Energy Showering , JHEP 11 (2017) 093 [ 1611.00788 ].
- [14] K. Hagiwara, I. Watanabe and P. Zerwas, Higgs boson production in e γ collisions , Physics Letters B 278 (1992) 187.
- [15] R. N. Cahn and S. Dawson, Production of Very Massive Higgs Bosons , Phys. Lett. B 136 (1984) 196.
- [16] R. N. Cahn, Production of Heavy Higgs Bosons: Comparisons of Exact and Approximate Results , Nucl. Phys. B 255 (1985) 341.
- [17] P. Ciafaloni and D. Comelli, Electroweak evolution equations , JHEP 11 (2005) 022 [ hep-ph/0505047 ].
- [18] C. W. Bauer, N. Ferland and B. R. Webber, Standard Model Parton Distributions at Very High Energies , JHEP 08 (2017) 036 [ 1703.08562 ].
- [19] S. Frixione, M. L. Mangano, P. Nason and G. Ridolfi, Improving the Weizsacker-Williams approximation in electron - proton collisions , Phys. Lett. B 319 (1993) 339 [ hep-ph/9310350 ].
- [20] V. Hirschi and O. Mattelaer, Automated event generation for loop-induced processes , JHEP 10 (2015) 146 [ 1507.00020 ].
- [21] R. B. Curtis, Meson production by electrons , Phys. Rev. 104 (1956) 211.
- [22] R. H. Dalitz and D. R. Yennie, Pion production in electron-proton collisions , Phys. Rev. 105 (1957) 1598.
- [23] S. J. Brodsky, T. Kinoshita and H. Terazawa, Dominant colliding beam cross-sections at high-energies , Phys. Rev. Lett. 25 (1970) 972.
- [24] N. Arteaga-Romero, A. Jaccarini, P. Kessler and J. Parisi, Photon-photon collisions, a new area of experimental investigation in high-energy physics , Phys. Rev. D 3 (1971) 1569.
- [25] R. Ruiz, A. Costantini, F. Maltoni and O. Mattelaer, The Effective Vector Boson Approximation in high-energy muon collisions , JHEP 06 (2022) 114 [ 2111.02442 ].
- [26] M. Casarsa, D. Lucchesi and L. Sestini, Experimentation at a muon collider , 2311.03280 .
- [27] N. V. Mokhov, Y. I. Alexahin, V. V. Kashikhin, S. I. Striganov and A. V. Zlobin, Muon Collider Interaction Region and Machine-Detector Interface Design , Conf. Proc. C 110328 (2011) 82 [ 1202.3979 ].
- [28] N. V. Mokhov and S. I. Striganov, Detector Background at Muon Colliders , Phys. Procedia 37 (2012) 2015 [ 1204.6721 ].
Yu Hamada Ryuichiro Kitano Ryutaro Matsudo Shohei Okawa Ryoto Takai Hiromasa Takaura Lukas Treuer
Higgs boson production at $μ^+ μ^+$ colliders
We study Higgs boson production at $\mu^+ \mu^+$ colliders at high energy.
Since both initial-state particles are positively charged, there is no $W$
boson fusion at the leading order, as it requires a $W^+ W^-$ pair. However, we
find that the cross section of the higher-order, $\gamma$- and $Z$-mediated $W$
boson fusion process is large at high center-of-mass energies $\sqrt s$,
growing as $(\log s)^3$. This is in contrast to the $\log s$ behavior of the
leading-order $W$ boson fusion. Thus, even though it is a higher-order process,
the rate of Higgs boson production for 10 TeV energies at $\mu^+ \mu^+$
colliders with polarized beams can be as high as about half of the one at
$\mu^+ \mu^-$ colliders, assuming the same integrated luminosity. To calculate
the cross section of this process accurately, we carefully treat the collinear
emission of the photon in the intermediate state. The thereby obtained large
cross section furthermore shows the significance of Higgs production with an
extra $W$ boson in the final state also at $\mu^+ \mu^-$ and $e^+ e^-$
colliders.