Multi Higgs production via photon fusion at future multi-TeV muon colliders

Multi-TeV muon colliders promise an unprecedented potential for exploring the particle-physics energy frontier, and, at the same time, can probe with fantastic accuracy the precise structure of the Standard Model, and in particular of the Higgs boson sector. Here we consider the possibility to measure the loop-suppressed single-, double-, triple-Higgs boson production mediated by the collinear-photon scattering in the channels $\mu^+\mu^-\to\mu^+\mu^-\gamma^\ast \gamma^\ast \to \mu^+\mu^- [H,HH,HHH]$. We study total rates and kinematical distributions in the Standard Model, and compare them with the corresponding vector-boson-fusion processes $V^\ast V^\ast\to H,HH,HHH$ at muon collisions center-of-mass energies in the range between 1.5 and 100 TeV. Possible strategies for enhancing the $\gamma\gamma\to H,HH,HHH$ signal over the dominant vector-boson-fusion production are proposed. The sensitivity of total rates to possible anomalies in the Higgs-boson self-couplings is also discussed.

Multi-TeV muon colliders could then provide a possible new path to expand the energy frontier in accelerator physics [23][24][25][26], potentially allowing the direct production of new heavy states, even beyond the reach of alternative future high-energy colliders possibly following LHC [1,[27][28][29].Many possible phenomenological implications for beyond-SM searches have recently been considered for multi-TeV muon colliders in .At the same time, multi-TeV muon collisions can offer a brand new laboratory to test with high precision the SM predictions for processes involving highmomentum transfer [82][83][84].This holds true especially for the electroweak (EW) sector of the SM whose phenomenology is moderately affected by QCD backgrounds at lepton colliders 1 .In particular, Higgs boson physics can be tested in new energy regimes [27,83,[87][88][89][90] with the enhancement of extremely challenging production channels like the triple Higgs production which is sensitive to the quartic Higgs self-coupling [91] .
Benefitting from the multi-TeV regime, here we want to analyse the multiple Higgs-boson production initiated by collinear (almost-on-shell) photons radiated by the initial muon beams, and mediated by loops of heavy particles In particular, we consider the production of one, two and three Higgs bosons in the channels 1 PDFs effects in muon beams have been considered in [85,86].The spectrum of quasi real photons emitted by an energetic charged lepton beam can be described by the leading-order effective photon approximation (EPA) where, for an energy fraction x of the charged lepton of initial energy E, the Weizsäcker-Williams spectrum is given by [92,93] with the splitting functions P γ/ℓ (x) = (1 + (1 − x)2 )/x for ℓ → γ, and P ℓ/ℓ (x) = (1 + x 2 )/(1 − x) for ℓ → ℓ.
We first analyze the total rates and kinematical features of the partonic processes γγ → H, HH, HHH.We then convolute them according to the EPA to obtain total rates and kinematical features of the corresponding processes µ + µ − → µ + µ − (H, HH, HHH) , mediated by collinear photons.In EPA, the final muons will be forward and completely lost in the detector.
Our benchmark scenarios for muon collision c.m. energies √ s and corresponding integrated luminosities L , according to a luminosity scaling L ∼ 10 ab −1 ( √ s/10 TeV) 2 [1], will be √ s ≃ 3, 6, 10, 14, 30, 100 TeV, (7) and respectively 2 .In order to put the quite moderate photon-induced multi-Higgs production in context, we will compare the γγ-induced µ + µ − → µ + µ − (H, HH, HHH) rates and distributions with the dominant ones for multi-Higgs production at tree-level [27,28,83,87,91] We stress that our numerical predictions for the processes in Eqs. ( 9)- (10) include complete tree-level matrix elements, i.e. the contribution of both vector-boson fusion (VBF) and s-channel production channels.However, since the VBF contribution is by far the dominant one at the c.m. energies in Eq. ( 7), we will dub the contributions in Eq. ( 9) and Eq.(10) as WW and ZZ fusion channels, respectively.The final neutrinos replacing muons in the WW -fusion channel will then mimic the photon induced process but with different kinematical features of the Higgs system.On the other hand, in the ZZ-fusion process, the final muons will be produced at an average p µ T ∼ M Z , which will differentiate the ZZfusion signal from the γγ-induced one.Indeed, we will show that in both the WW and ZZ fusion channels the transversemomentum distribution for the H, HH, HHH systems has a maximum at p T values of the order of the W, Z masses which is typical of vector-boson fusion induced processes (see, e.g., [94], and [95] for single Higgs production in VBF).This will give a handle to separate the latter from the γγinduced signal that in the EPA in Eq. ( 6) is characterized by We will also compare the photon-induced loop production rates with other γγ multi-Higgs cross sections, where the Higgs bosons are radiated at tree-level by the intermediate amplitudes γγ → W ( * ) W ( * ) and γγ → t ( * ) t( * ) , in the processes respectively, with forward (untagged) final muons.The two extra final W 's and top quarks in Eqs. ( 13)-( 14) will of course distinguish these channels from the γγ loop production.We will see that, at particularly large √ s values, the production mechanism in Eq. ( 13) becomes dominant over the one-loop γγ fusion.On the other hand, the contributions in Eq. ( 14) are always smaller than the loop-induced ones in the setups under consideration.
The plan of the paper is the following.In Section 2, we describe the computational method used to obtain the photon-photon one-loop amplitudes and cross sections for multiple Higgs production.In Section 3, we study the partonic γγ → (nH) cross sections and kinematical distributions versus √ s γγ .The convoluted µ + µ − → µ + µ − γ * γ * → µ + µ − (nH) cross sections (kinematical distributions) are considered in Section 4 (5), and compared with the corresponding VBF processes.In Section 6, we discuss a possible strategy to enhance the photon induced signal over the VBF production.We then study the sensitivity of the µ + µ − → µ + µ − γ * γ * → µ + µ − HH cross section to a trilinear anomalous Higgs self-coupling, and of the µ + µ − → µ + µ − γ * γ * → µ + µ − HHH cross section to a trilinear or quartic anomalous Higgs self-coupling in Section 7. Our conclusions are presented in Section 8.

Cross section computational method
In this section we describe the procedure we followed to obtain the amplitudes and corresponding cross sections for the different processes described below.
For the phase-space integration, we used the Monte Carlo integrator already employed in [104,105] for diboson production in hadronic collisions.In the Monte Carlo program, the proton PDFs have been replaced by a photon PDF describing the collinear µ → µγ splitting as in Eq. (6), where m ℓ is the muon mass and E is the nominal energy of the muon beam.
The same procedure has been followed to compute the cross sections for the tree-level processes γγ → (nH)W + W − and γγ → (nH)t t (n = 1, 2, 3), relevant for the channels in Eqs. ( 13)- (14), and for the subsequent convolution with the photon PDF in the muon beams.
The tree-level results for µ + µ − → (H, HH, HHH)ν µ νµ via WW fusion, and µ + µ − → (H, HH, HHH)µ + µ − via ZZ fusion have been obtained by using the Whizard Monte Carlo event generator [106,107].In order to prevent singularities arising from s-channel diagrams, we include the gauge-boson widths in a gauge-invariant way by means of the complexmass scheme [108][109][110].We checked that, for these two classes of processes, the effect of s-channel diagrams is essentially dominated by the Z resonance region and it is numerically relevant only for the 3 TeV setup in Eq. ( 7), in particular for HHH production, while it becomes completely negligible when increasing the collider energy.

Partonic γγ cross sections
In this section, we discuss the total rates and kinematical distributions for the partonic channels γγ → H, HH, HHH, in the γγ c.m. system, before the convolution with the photon energy spectrum in the muon beams.
In the present approximation, the single Higgs channel γγ → H is characterized by a production rate resonating at the Higgs mass, σ ∼ δ ( √ s − m H ), with the Higgs boson at rest in the γγ c.m. system.Once convoluted with the collinear photon spectrum in muon collisions, the final Higgs will still have zero transverse momentum, and a rapidity distribution that will just reflect the one of the initial γγ system.
As for the partonic γγ → HH, HHH productions, we present in Figure 4 the corresponding total cross sections versus the γγ c.m. energy √ s γγ , in the range 2m H < ∼ √ s γγ < ∼ 5 TeV.The two-Higgs cross section reaches its maximum value σ max (γγ → HH) ∼ 0.28 fb for √ s γγ ∼ 470 GeV, (15) while the three-Higgs cross section has a maximum The Higgs transverse-momentum, rapidity and angular distributions for the γγ → HH process in the γγ c.m. system are shown in Figure 5, with the two Higgs bosons having the same energy E H = √ s γγ /2, and a back-to-back configuration.
Figure 6 shows the γγ → HHH case (including the nontrivial Higgs energy distribution), by ordering the three final Higgs-bosons labels according to their transverse momentum.

HH, HHH) cross sections
In this section, we present the µ + µ − total cross sections for the γγ-induced one-loop (H, HH, HHH) production after the convolution of the partonic cross section discussed in Section 3 with the photon spectrum described by the EPA in Eq. ( 6).
In Table 1, the inclusive µ + µ − cross sections (in ab) are presented versus √ s, in the range 1.5 TeV < ∼ √ s < ∼ 100  8), would correspond to about 16.000 events.At √ s ∼ 3 TeV with 1 ab −1 , one has σ (γγ→H) ∼ 0.37 fb, corresponding to about 400 events.These rates are about 3 and 2 orders of magnitude less than, respectively, the WW → H and ZZ → H ones for all √ s values.As for double Higgs boson production, the γγ → HH cross section varies in the range σ (γγ→HH) ∼ [0.4 − 8] ab.At √ s ∼ 14 TeV, σ (γγ→HH) ∼ 3 ab, which, with 20 ab −1 , would correspond to about 60 events.At √ s ∼ 3 TeV with 1 ab −1 , one has σ (γγ→HH) ∼ 0.9 ab, corresponding to about 1 event.Similarly to single Higgs production, these HH rates are, about 3 and 2 orders of magnitude less than, respectively, the WW → HH and ZZ → HH ones for all √ s values.Finally, for triple Higgs boson production, the γγ → HHH cross section varies in the range σ (γγ→HHH) ∼ [0.07 − 2.3] × 10 −2 ab, and requires the maximum c.m. energy to become observables.Indeed, at √ s ∼ 14 TeV, σ (γγ→HHH) ∼ 8.4 × 10 −3 ab, which, with 20 ab −1 , would correspond to no event produced, just as in the case of √ s ∼ 3 TeV (with 1 ab −1 ), where one has σ (γγ→HHH) ∼ 2 × 10 −3 ab.On the other hand, at √ s ∼ 30 TeV, σ (γγ→HHH) ∼ 1.3 ×10 −2 ab, which, with 90 ab −1 , would correspond to about one event.Hence, only going to the very upper √ s range considered here, one could have an observable event number.I particular, at √ s ∼ 100 TeV with 103 ab −1 , σ (γγ→HHH) ∼ 2.3 × 10 −2 ab, one would give about 23 events.Note that, at low √ s values, the HHH rates are less depleted than in the case of H and HH production with respect to VBF production, being only about 2 and 1 orders of magnitude less than the WW → HHH and ZZ → HHH ones, respectively.
The behavior of µ + µ − total rates versus √ s is more clearly shown in Figure 7, where for comparison also VBF cross sections are plotted.In the left plot the µ + µ − → Hµ + µ − production (blue solid line) is compared with the same process when replacing µ ± → e ± (blue dashed line), and with WW -fusion (red solid line), ZZ-fusion (green solid line) single-Higgs production, for √ s < ∼ 400 GeV.The difference of about a factor three in the rates for the electron and muon initiated processes for γγ fusion stems from the double logarithm arising from Eq. 6 applied twice to initial beams.The dashed red and green lines represent the integrated crosssections for single Higgs productions in pure WW and ZZ fusion, respectively, in µ + e − scattering: in practice, since both electrons and muons in our simulations are treated as massless particles (with the only exception of the mass parameter in the photon spectrum in Eq. ( 6)), these lines correspond to the predictions for H production in WW and ZZ fusion when only the t-channel production mechanism is considered.The difference between red (green) solid and dashed lines is large at the opening of the Z resonance ( √ s ∼ M H + M Z ) and decreases for larger values of √ s.In the right plot of Figure 7, µ + µ − collision rates for the H (solid lines), HH (dashed lines), HHH (dot-dashed lines) one-loop γγ production (blue lines) are shown, and compared with the corresponding rates for tree-level WW -fusion (red lines), and ZZ fusion (green lines) production, for 1 TeV < ∼ √ s < ∼ 100 TeV.In this partonic center of mass range, the µ + µ − → nHν ν over µ + µ − → nHµ + µ − cross-section ratio is essentially flat, ranging from about 10 for single-Higgs production to approximately 6 for the triple-Higgs case: this is mainly a consequence of the differences in the HVV , HHVV , and V f f couplings for V = W, Z with the leading contribution coming from the latter interaction as documented in Ref. [95] for single-Higgs production.

Higgs boson distributions in γγ-induced HH and HHH production
In this section, we go through the main kinematical features of the Higgs bosons in γγ-induced µ + µ − → µ + µ − (HH, HHH) production.We first present the inclusive Higgs kinematical distributions.Then we include the effect of a reduced acceptance in the Higgs kinematics due to a possible shielding of the detector in the forward regions, where the beam-induced background (BIB) might become prohibitive [111][112][113][114][115][116][117][118].To this end, we will assume a two-body decay of the produced Higgs (giving the dominant contribution to the Higgs decay width, through the channels H → b b, c c, gg, ττ), and put a minimum transverse-momentum cut and a maximum rapidity cut on the Higgs decay products (applying no decay branching ratio) 3 .In particular, we simulate the effect of introducing "shielding nozzles" in the detector, by assuming a detector coverage for all kind of particles only for transverse momentum and rapidity p T > 20 GeV, |y| < 3.
As for the single Higgs production, in the present approximation, the Higgs kinematics simply reflects the rapidity distribution of the initial γγ system as provided by the collinear photon spectrum described by the EPA in Eq. ( 6), which gives vanishing Higgs transverse momentum.
The inclusive Higgs transverse-momentum and rapidity distributions for the two Higgs production in µ + µ − → µ + µ − HH mediated by γ * γ * → HH are presented in Figure 8, for different values of √ s.The same distributions for the three Higgs production in µ + µ − → µ + µ − HHH mediated by γ * γ * → HHH are shown in Figure 9, where H 1 labels the highest-p T Higgs, and H 2 the second highest-p T Higgs.
In Figure 10, we present instead the photon-induced µ + µ − → µ + µ − HH Higgs transverse-momentum and rapidity distributions after applying the selection cuts p b T > 20 GeV, |y b | < 3 on the H → b b decay products (blue lines), at √ s = 3, 10, 30 TeV.For comparison, the same distributions for the treelevel WW channel (red lines) and ZZ channel (green lines) are also shown.For the process µ + µ − → µ + µ − HH , calculated with EPA, both Higgs bosons have the same transverse momentum and therefore only the blue solid line is shown.
6 Enhancing the γγ → n H signal over the WW, ZZ → n H background The previous discussion clearly shows that one-loop γγ → H, HH, HHH production in µ + µ − collisions is quite suppressed with respect to the WW, ZZ → H, HH, HHH production.Independently of that, a γγ → HHH signal can be considered statistically out of reach even for the most optimistic luminosities scenarios.On the other hand, the higher statistics channels γγ → H, HH, being suppressed by orders of magnitudes with respect to WW, ZZ → H, HH production, will be in general overwhelmed by the latter, making the detection of loop-induced events quite hard in general.
In this section we propose a strategy that could dramatically enhance a γγ → H, HH, HHH signal over the VBF ones.
It is well know that VBF single Higgs production is characterized by a typical Higgs transverse momentum of the order of the M W,Z [Eq.(11)].This feature is shared by VBF multiple Higgs production, and is connected to the mechanism of vector-boson radiation from the initial beams.
In Figure 12, we show the transverse-momentum distribution of the HH (solid lines) and HHH (dashed lines) systems in VBF-mediated HH and HHH production, respectively, in µ + µ − collisions, at various √ s.The transversemomentum distribution for the HH and HHH systems indeed keeps the single Higgs VBF production main feature, with a maximum at p T (HH, HHH) ∼ M W,Z .This should be confronted with the γγ one-loop process where, in the approximation of Eq. ( 6), one has p T (H, HH, HHH) ∼ 0.
In order to enhance the γγ signal over the VBF one, one can then impose, depending on the momentum resolution ∆ p T (n H) of the detector, a cut on large p T (n H) events (that is events characterized by p T (n H) > ∼ ∆ p T (n H)), which leaves unsuppressed the photon-induced p T (n H) ∼ 0 events.In this way, the actual background from VBF events is reduced to the sub-dominant fraction of events characterized by a small p T (n H), hopefully as small as much less than O(M W ).
Figure 13 shows the effects on µ + µ − total rates of a cutflow selecting low n H transverse-momentum events.The n H transverse-momenta are reconstructed by summing up the vectorial momenta of the H → b b decay products, assuming the finite acceptance p b T > 20 GeV, |y b | < 3 (not including any decay BR) .Then, the inclusive cross sections for the n H processes (solid lines) are compared with the production rates that pass the H → b b basic acceptance cuts (dashed lines), and with the rates where a further constrain on the maximum n H transverse-momentum (as reconstructed from the Higgs decay products) is applied (dotdashed lines).Two cases are analyzed for the maximum n H transverse-momentum, which are p T (n H) < ∼ 10 GeV (left plot), and 30 GeV (right plot).
Figure 13 clearly shows how a cut on large transverse momentum of the Higgs system can upset the hierarchy of the rates for the γγ and VBF processes, hence allowing to control the VBF background to multiple Higgs γγ production.In particular, in the p T (n H) < ∼ 10 GeV case (left plots), one can see that after the acceptance cuts the γγ single Higgs production gets comparable to the ZZ channel at all √ s, while both the double and triple Higgs γγ productions get the upper-hand over the ZZ channels.The WW channel is maintained moderately dominant over γγ for both H and HH productions, and gets comparable to γγ in HHH production.The VBF background suppression gets of course less dramatic by relaxing the p T (n H) to 30 GeV (Figure 13, right plots).
The presence of a veto on the total transverse momentum of the nH system clearly affects also the differential distributions of the background studied in Figures 10 and 11.On the one hand, the main effect of such cut is to induce a suppression of the µ + µ − → HH(H)ν ν and µ + µ − → HH(H)µ + µ − processes similar to the one observed at the integrated crosssection level.On the other hand, the veto on p T (nH) also introduces some shape distortions for the background processes, in particular for the Higgs-boson p T distributions where the impact of this cut becomes more severe for increasing Higgs p T values.

Cross section sensitivities to anomalous Higgs self-couplings
Before concluding, we want to briefly discuss the sensitivity of the γγ → HH, HHH cross sections to possible anomalies in the Higgs boson self-interactions.In particular, we will assume a deviation in either the trilinear Higgs coupling λ 3 (1 + δ 3 )HHH or the quartic Higgs coupling λ 4 (1 + δ 4 )HHHH, where λ 3,4 are the self-coupling SM values, and δ 3,4 ̸ = 0 would reveal possible deviations BSM in the Higgs potential.We leave a more sophisticated analysis performed in terms of effective-field theories for future work.In Figure 14, left plot, we detail the ratio of the γγ → HH total cross section with and without anomalous Higgs self- couplings versus δ 3 at different muon collision √ s, while, in the right plot, we show the same quantity after applying the acceptance cuts on Higgs decays.At the level of inclusive quantities, there is a good sensitivity to anomalous trilinear coupling of the double Higgs production.In particular, a δ 3 ≃ ±1 variation increases the cross sections by about 15-20%.After imposing the acceptance cuts, the sensitivity is in general slightly reduced In Figure 15 we show the δ 3 dependence for δ 4 = 0 (upper plots), and δ 4 dependence for δ 3 = 0 (lower plots) for the ratio of the cross sections for triple Higgs production with and without anomalous couplings via photon scattering at different muon √ s.The results corresponding to the setups with (right plots) and without (left plots) acceptance cuts are shown.Note that for the reference integrated luminosities in Eq. ( 8), the tiny γγ → HHH production cross section will not allow a precision cross-section measurement.However, as in the case of double Higgs production, there is a sizeable dependence on anomalous Higgs self-couplings, in particular for coupling variations of the order of δ 3(4) ∼ ±1, at δ 4(3) ≃ 0. Furthermore, γγ → HHH cross section has a minimum around its SM value for δ 4 ∼ 0 and variable δ 3 , while the δ 4 dependence is monotonic in the δ 4 range considered in Figure 15 (lower plots).We checked that the latter monotonic δ 4 dependence is the outcome of a nontrivial interplay between the top-loop and the W -loop amplitudes providing the final sensitivity to anomalous δ 4 .
Fig. 13 Effects of the cut-flow selecting low n H transverse-momentum events.The µ + µ − inclusive cross section for the n H processes (solid lines) are compared with the production rates that pass the Higgs decays basic acceptance (dashed lines), and with the rates where a further constrain on the maximum n H transverse-momentum (as reconstructed from Higgs decay products) is applied (dot-dashed lines), which is 10 GeV in the left plot, and 30 GeV in the right plot.Fig. 14 Ratio of the integrated cross sections with and without anomalous Higgs self-couplings for the process γγ → HH in the inclusive setup (left plot) and after imposing acceptance cuts on the Higgs decay products (right plot) as a function of δ 3 .

Conclusions
Single, double and triple Higgs-boson production mediated by photon fusion via loops of heavy particles (top quarks and W bosons) have been studied at multi-TeV lepton colliders.For integrated collision luminosities scaling as L ∼ 10 ab −1 ( √ s/10 TeV) 2 , at √ s ∼ 14(3) TeV one has inclusive cross sections corresponding to 16.000 (400) single Higgs events and 60 (1) double Higgs events, while triple Higgs production is quite suppressed, and needs √ s > ∼ 30 TeV for getting more than 1 HHH event.We have studied kinematical distributions at different √ s, and analyzed the effects of a finite acceptance on the Higgs decay products in order to control BIB effects via dead cones along the beams.The comparison with the dominant VBF mechanism for producing single, double and triple Higgs final states has been carried out, and possible strategies for pinpointing γγ-induced events on the basis of the total transverse momentum of the Higgs system have been suggested.Finally, the inclusive cross sections sensitivity to triple and quartic Higgs selfcouplings have been discussed and found sizeable.
Multi-TeV muon colliders offer a unique possibility to probe such processes.Apart from single Higgs production, which has been discussed at future e + e − colliders [119], we are not aware of any similar study at multi-TeV future colliders [120] allowing to compare the corresponding potential for multi-Higgs production induced by quasi-real photon collisions.
Our study proceeded under the strict assumption of collinear initial photons described by means of the EPA, which is valid in first approximation.A more refined analysis should relax this assumption, and consider in γγ one-loop multi-Higgs production also interference and coherence effects with Z-boson mediated processes described in terms of electroweak PDFs [94].This would imply including at the amplitude level both VBF ZZ fusion channels (possibly including higherorder EW corrections) and one-loop γZ-induced channels.However, the effects of the latter two processes should affect mainly non-vanishing p T (nH) final states, and the description in terms of collinear EW PDFs might need to be improved to better describe the transverse momenta of the initial-state vector bosons.We then leave an extensive discussion of this issue to further dedicated studies.Fig. 15 Ratio of the integrated cross sections with and without anomalous Higgs self-couplings for the process γγ → HHH in the inclusive setup (left plots) and after imposing acceptance cuts on the Higgs decay products (right plots) as a function of δ 3 (upper plots) and δ 4 (lower plots).

Fig. 1
Fig. 1 Feynman diagrams contributing to the process γγ → H, involving top-quark and W -boson loops.

Fig. 2
Fig. 2 Feynman diagrams contributing to the process γγ → HH, involving top-quark and W -boson loops.

Fig.
Fig. Transverse-momentum distribution of the HH (solid lines) and HHH (dashed lines) systems in vector-boson-fusion mediated double and triple Higgs production, respectively, in muon collisions, at different √ s.