Triple Higgs boson production at a 100 TeV proton-proton collider

We consider triple Higgs boson production at a future 100 TeV proton-proton collider. We perform a survey of viable final states and compare and contrast triple production to Higgs boson pair production. Focussing on the $hhh \rightarrow (b\bar{b}) (b\bar{b}) (\gamma \gamma)$ final state, we construct a baseline analysis for the Standard Model scenario and simple deformations, demonstrating that the process merits investigation in the high-luminosity phase of the future collider as a new probe of the self-coupling sector of the Higgs boson.

We consider triple Higgs boson production at a future 100 TeV proton-proton collider. We perform a survey of viable final states and compare and contrast triple production to Higgs boson pair production. Focussing on the hhh → (bb)(bb)(γγ) final state, we construct a baseline analysis for the Standard Model scenario and simple deformations, demonstrating that the process merits investigation in the high-luminosity phase of the future collider as a new probe of the self-coupling sector of the Higgs boson.

I. MULTI-HIGGS BOSON PRODUCTION AT HADRON COLLIDERS
With the exception of very few interactions, most of the terms that comprise the Standard Model (SM) Lagrangian have been measured or constrained, their strengths found to be suggestively close to the expected ones. An important category of interactions not directly observed are those of the the Higgs boson with itself. The so-called 'self-couplings' and their energy dependence are crucial in determining the stability of the vacuum. Current observations suggest that our Universe may be sitting at a metastable false vacuum [1][2][3][4][5][6][7][8] and measurements of these couplings will illuminate this fact further.
At colliders, these terms, i.e. those proportional to h n , h being the Higgs boson scalar field, can be directly probed through the simultaneous production of (n − 1) Higgs bosons. Unfortunately, the production rates for processes with n ≥ 3, i.e. more than one Higgs boson, are small, mainly due to the relatively large invariant mass of the final state system. In particular, at the Large Hadron Collider (LHC) with 14 TeV proton-proton centre-ofmass energy, gluon-fusion Higgs boson pair production is expected to have a cross section of ∼40 fb [9][10][11][12][13][14][15][16][17], whereas triple production is expected to have a rather dwarfish rate, with a cross section of O(0.1 fb) [15]. Hence, even though there is optimism that Higgs boson pair production will provide important information and constraints through LHC measurements , any direct measurement of SM-like triple Higgs boson production will be essentially impossible at the LHC, even at the end of the high-luminosity phase (HL-LHC) [76,77]. However, with a significant increase in the collision energy, a Future Circular hadron-hadron Collider (FCC-hh), colliding protons at 100 TeV, stands a good chance at observing and constraining the self-coupling of the Higgs bosons through Higgs boson pair production [64,71,[78][79][80], the cross section rising to ∼1.6 pb [81]. Additionally, at the FCC-hh one may also get the chance to observe three onshell Higgs bosons being produced, since the total cross section rises to ∼5 fb [15]. The evaluation of this possibility is the main object of the present article.
Concretely, the part of the Higgs boson potential which includes the self-interactions, may be written as: where v 246 GeV is the vacuum expectation value (vev), m h 125 GeV is the measured Higgs boson mass and c 3 and d 4 parametrize possible deviations from the standard model expectation (i.e. the SM is recovered for c 3 = d 4 = 0). Figure 1 shows some of the Feynman diagrams contributing to triple Higgs boson production. It is clear that the production cross section depends on both c 3 and d 4 parameters. This should be contrasted to double Higgs boson production, which does not depend on d 4 . In Ref. [76] the dependence of the triple Higgs boson cross section on the parameters c 3 and d 4 was investigated at 14 TeV and 200 TeV proton-proton colliders for a Higgs boson mass m h = 120 GeV. We produce an equivalent result for proton-proton collisions at 100 TeV, for m h = 125 GeV, shown in Fig. 2. The conclusions are similar to those drawn in [76]: the cross section dependence on d 4 is mild, the deviations due to d 4 = ±1 being at most ±20% for c 3 = 1. Hence modifications of the d 4 coefficient itself will be very challenging to probe. This is also demonstrated in the contour plot of Fig. 3(a), which shows the cross section normalised to the SM value, on the c 3 − d 4 parameter space. On this plane, one can observe that the dependence along d 4 is much weaker than that along c 3 .
In terms of constraining c 3 , triple Higgs boson production cannot be superior to double Higgs boson production due to its small production cross section. On the other hand, triple production would be the best process to constrain d 4 , although, as we will demonstrate, even the FCC-hh with 30 ab −1 of integrated luminosity can only provide O(1) constraints on d 4 , because its dependency of the cross section is very modest. However, observing the triple Higgs boson production process is an interesting task in its own right, and as will be seen, indeed challenging at the FCC-hh. The goal of this article is to provide a first baseline study of Standard Model-like triple Higgs boson production via gluon fusion (ggF), at a future 100 TeV proton-proton collider. Furthermore, we investigate triple Higgs production in two scenarios where it is affected by new physics: (i) in the SM augmented by a single higher-dimensional operator in an effective field theory approach and (ii) the generic case on the (c 3 − d 4 )-plane.
The article is organised as follows: in Section I A we investigate an explicit scenario that contains a single higher-dimensional operator. In Section II we list, for future reference, the final states that could be interesting in the study of Higgs boson triple production. The Monte Carlo event generation, simulation of b-jet and photon tagging are described in Section III. Differential distributions at parton level for triple Higgs boson production at 100 TeV, compared to those of Higgs boson pair production and the analysis of the channel (bb)(bb)(γγ) is described in Section IV. We use this analysis to provide constraints in two scenarios. Finally, we provide discussion and conclusions in Section V.
A. The self-coupling in D = 6 EFT In the framework of the dimension-6 operator extension to the Standard Model (D = 6 EFT), one can compare the sensitivity of multi-Higgs production to variations of the operator Wilson coefficients [50]. Here we consider, as an illustrative example, a simplified mode with the assumption that the effect of all coefficients apart from a single one, originating from an operator of the form O 6 ∼ |H| 6 , where H is the Higgs doublet scalar before electroweak symmetry breaking: where µ 2 and λ are the conventional parameters employed in the SM potential for the Higgs doublet H.
The changes in the quartic and the triple Higgs couplings, defined in Eq. 1, are related via [50]: * Due to the relation appearing in Eq. 3, the cross section for triple Higgs boson production is a quartic polynomial in c 6 , i.e. it contains terms up to c 4 6 . Such terms come from squared matrix elements of diagrams containing two triple Higgs couplings, such as the one shown in Fig. 1(d).
In Fig. 3(b) we show the variation of the inclusive leading-order cross sections for ggF hh and hhh with respect to the SM (c 6 = 0). The fit as a function of c 6 for * Note that, in general, c 3 and d 4 would be multiplied by v 2 /Λ 2 in D = 6 EFT. We have set Λ = v for simplicity here. the two cases, at 100 TeV, is: The line d 4 = 6c 3 is also shown as a dissection on the c 3 − d 4 plane in Fig. 3(a).

II. TRIPLE HIGGS PRODUCTION FINAL STATES
We list the dominant Higgs boson triple production final states, i.e. those that yield N events > 10 with 30 ab −1 of integrated luminosity at a proton collider at 100 TeV centre-of-mass energy, in Table I.  and their branching ratios (BR). The subscript " x " denotes the number of leptons x in the final state, originating from the di-bosons. The cross section used for pp → hh at 100 TeV is σNLO = σLO × 2.0 = 5.78 fb, where a K-factor K = 2.0 has been applied to obtain an estimate of the NLO cross section. The number of events has been rounded to the nearest integer. If we apply further requirements to the final states listed in Table I: • to possess greater than 100 events at 30 ab −1 of integrated luminosity, • and all gauge bosons fully decay to leptons, then we are left with the following interesting final states: In particular, the expected combined number of events in the multi-b-jet and multi-τ final states is ∼45000 over the lifetime of the FCC-hh, and will most likely provide valuable information on the triple Higgs boson process. In the present study we focus on the rare but clean final state (bb)(bb)(γγ).

A. Detector simulation
In the hadron-level analysis that follows, we consider all particles within a pseudorapidity of |η| < 5 and p T > 400 MeV. We reconstruct jets using the anti-k t algorithm available in the FastJet package [82,83], with a radius parameter of R = 0.4. We only consider jets with p T > 40 GeV within |η| < 3.0 in our analysis. We consider photons within |η| < 3.5 and p T > 40 GeV and 100% reconstruction efficiency. The jet-to-photon misidentification probability is taken to be P j→γ = 10 −3 , flat over all momenta above the p T cut and over all pseudorapidities. † We also consider the mis-tagging of two light jets to bottom-quark-initiated jets with a flat probability of 1% for each mis-tag, corresponding to a flat b-jet identification rate of 80% and demand that they lie within |η| < 3.0. We demand all photons to be isolated, an isolated photon having i p T,i less than 15% of its transverse momentum in a cone of ∆R = 0.2 around it. Finally, no detector-smearing effects have been considered.

B. Event generation
Events for the hhh signal samples have been generated via the loop-induced module of the MadGraph 5/aMC@NLO package [84][85][86][87][88]. The SM loop model present in MadGraph 5/aMC@NLO was modified to allow for deformations of the Higgs boson triple and quartic selfcouplings away from the SM values. All tree-level and next-to-leading order (i.e. matched via the MC@NLO method [89]) background processes have been generated using MadGraph 5/aMC@NLO, apart from the di-Higgs plus jets (hh + jets) background, which was simulated using HERWIG++ in conjunction with the OpenLoops † Note that the HL-LHC expectation has the approximate form P j→γ = 0.0093 × e −0.036p T j /GeV [78]. For a p T ∼ 40 GeV, this gives approximately P j→γ ∼ 2 × 10 −3 . Thus, the value employed here is expected to be a reasonable approximation to future detector performance.  1)). On the right-hand panel one can see the variation with respect to the SM in a theory where the SM is extended with a O6 ∼ |H| 6 operator as in Eq. 2, for both Higgs boson pair production (hh) and Higgs boson triple production (hhh). For both calculations, the NNPDF23 nlo as 0119 parton density function set was used.
matrix-element generator [32,90]. The default parton density functions were used in each case: for the signal and tree-level backgrounds (including hh+jets) the NNPDF23 nlo as 0119 set was used, whereas for the NLO samples the NNPDF23 nlo as 0118 qed set was employed [91].
Due to the large cross sections and high-multiplicity final states present at a 100 TeV collider, we only generate the tree-level processes to include true photons and true b-quarks at parton level. This implies that light extra jets for these processes will be generated by the parton shower, for which we employ the HERWIG++ generalpurpose event generator [92][93][94][95]. ‡ Inevitably this introduces an uncertainty to the results presented herein, rendering any observables related to these light jets leadinglog accurate. § We do not expect this, however, to alter the main conclusions of this first, baseline, study. Furthermore, generation-level cuts that anticipate the analysis cuts at hadron level are imposed on the b quarks and the photons. In the case of decaying resonances (i.e. h and Z bosons) no cuts are imposed. The phase-space cuts applied on the samples bbbb, bbbbγ, bbbbγγ, bbγγ are shown in Table II. ‡ Simulation of hadronization and the underlying event were also included. [96]. No simulation of pile-up events was considered. § The hh+jets process is the only exception, with the first jet being leading-order accurate [32].  At this point one should stress that even though NLO event generation matched to the parton shower has been largely automated, NLO calculations for the highmultiplicity final states, particularly with many coloured particles and complicated phase space cuts, remain challenging at present. We hence apply a conservatively large flat K-factor of K = 2.0 to all the processes calculated at tree level, as well as the hhh and hh+jets loop-induced processes. This is a crucial point that should be addressed in future studies at higher-energy hadron colliders, as such final states will become increasingly common.
The analysis of the signal and backgrounds generated for the final state (bb)(bb)(γγ) is presented in section IV B.

A. Differential distributions
We investigate the shape of the differential distributions in Higgs triple production in the Standard Model.
Here we keep the Higgs bosons stable and include parton shower effects. We compare the shape of the hhh distributions to those coming from the more familiar case of Higgs boson pair production (hh) at 100 TeV. Figure 4(a) shows the transverse momentum of any single Higgs boson either in hh or hhh production, p T,h . Evidently, the transverse momentum of a Higgs boson in hhh is softer than that of hh, peaking at ∼ 100 GeV instead of ∼ 150 GeV.
In Fig. 4(b) we show the the spectrum of the transverse momentum of the Higgs boson "system", p T,h n , i.e. the triplet of Higgs bosons in hhh, and the two Higgs bosons in hh. One can observe that the p T,h n is harder in hhh than that of the pair in hh.
We examine the distance between two Higgs bosons, ∆R(h, h), in hh and hhh production in Fig. 4(c). In the case of triple production the distance is calculated between any two Higgs bosons. The Higgs bosons in hh are found to be more back-to-back than those in hhh, as expected.
Finally, in Fig. 4(d) we show the the invariant mass of all Higgs bosons in hh or hhh production, M h n . The invariant mass distribution in hhh peaks just above M h 3 ∼ 600 GeV, whereas that in Higgs pair production, just above M h 2 ∼ 400 GeV.

B. hhh → (bb)(bb)(γγ)
The hhh → (bb)(bb)(γγ) process is expected to be relatively clean and simple to reconstruct. ¶ The excellent resolution of the di-photon invariant mass, that has contributed to the Higgs boson discovery at the LHC's Run 1, can be exploited to facilitate background rejection.
The present analysis follows a simple path, using the R = 0.4 anti-k t jets as described in Section III. Note, however, that an analysis utilising the jet substructure of boosted Higgses to a bottom-anti-bottom pairs, e.g. as in [98], could assist in signal-background separation. We defer this task to future work.
We ask for four b-jets, or light jets mis-identified as b-jets, within |η| < 3.0, possessing transverse momenta p T,b {1,2,3,4} > {80, 50, 40, 40} GeV, where the subscripts 1, 2, 3, 4 denote the first, second, third and fourth hardest b-jets respectively. We ask for two photons, or mis-identified jets as photons, within |η| < 3.0 and ¶ Note that this final state has been considered in [97], in the context of the two-Higgs doublet model hH → hhh final state. Here we consider the SM case. p T,γ {1,2} > {70, 40} GeV. Due to the fact that, for the majority of b-jets we cannot identify whether they originated from a b-quark or an anti-b-quark, there exists a 3-fold combinatorial ambiguity in combining b-jets into the two Higgs boson candidates. As a simple choice, we take the highest-p T b-jet and pair it with the closest bjet in ∆R = ∆η 2 + ∆φ 2 , and pair the other two remaining b-jets together. We thus construct the paired b-jet invariant mass, respectively, m close,1 bb and m close,2 bb , for which we demand m close,1 bb ∈ [100, 160] GeV and m close,2 bb ∈ [90, 170] GeV. The rather large mass windows are chosen to maintain high signal efficiency given the small initial cross section. Moreover, we construct the distance between the highest-p T b-jet and the corresponding paired one, and impose ∆R close,1 bb ∈ [0.2, 1.6]. * * For the photon pair, we simply construct the invariant mass and impose a strong window on the measured Higgs boson mass m γγ ∈ [124, 126] GeV. † † We also restrict the distance between the two photons to ∆R γγ ∈ [0.2, 4.0]. We collect these selection cuts in Table III  We show a summary of the processes considered in the analysis in Table IV. The most significant backgrounds in our set-up turn out to be the SM bbbbγγ and those coming from Higgs boson pair production in association with extra jets. Specifically, the latter emulates the signal well, as the di-photon mass window is expected to have similar efficiency to the signal. Moreover, as we have pointed out at the beginning of the section, the Higgs bosons in hh are harder on average than than those in We have verified explicitly that an alternative method based on minimization of the squared sum of (m bb − m h ) from each combination yields results that differ by O(1%) compared to the simpler ∆R method. * * The distance between the other paired b-jets was not found to have significant discriminatory power. † † This cut implies that the di-photon resolution should be better than ∼ 1 GeV at the FCC-hh. The current resolution at the LHC is 1-2 GeV, [99,100] and thus it is not unreasonable to expect an improvement at the detectors of the future collider. hhh, thus passing transverse momentum cuts easily. This background could be tackled in future studies via h → bb tagging using jet substructure techniques that exploit the decay versus the g → bb branching that produces the additional bb pair in hh+jets. ‡ ‡ ‡ ‡ Note that the additional two b-jets in hh+jets and hZ have been generated by gluon splitting into bb, performed by the shower Monte Carlo.

C. Sensitivity in D = 6 EFT
Despite the rather large backgrounds, a signal-tobackground ratio of O(1) can be obtained for the SM case. To summarise the results of the analysis, we present in the first two columns of Table V, respectively, the number of expected hhh events and the total expected number of events, for the SM, as well as for the two simple deformations obtained by including the D = 6 operator O 6 , with coefficient values c 6 = ±1. The third column of Table V indicates that, if one assumes that the SM is the underlying theory, then c 6 = ±1 can be excluded at 95% C.L. or better, using hhh → (bb)(bb)(γγ) at the 'high-luminosity' phase of the FCC-hh.
Furthermore, we show in Fig. 5 the expected exclusion region on the c 6 coefficient, as well as the expected  The processes considered in the analysis of the (bb)(bb)(γγ) final state. The parton-level cross section, including the cuts given in the main text is given (if any), the analysis efficiency and the expected number of events at 30 ab −1 are given. A flat K-factor of K = 2.0 has been applied to all tree-level processes (including hh+jets) as an estimate of the expected increase in cross section from LO to NLO. The hZZ, hhZ and hZ processes have been produced at NLO and hence no K-factor is applied. Even though the hhZ process depends on c6, we only consider the SM case, as it was found to be negligible after cuts.
number of events after cuts, at 30 ab −1 . The theoretical uncertainty on the expected number of events for the hh and the hh+jets processes was taken to be 40% and uncorrelated between the two. The analysis efficiencies for hhh and hh+jets were individually fitted using points in the region c 6 ∈ [−3.0, 4.0]. § § We assume that there is negligible uncertainty on the 'other' backgrounds, which are taken to consist of the bbbbγγ and bbbbγ+jets processes. 3.5. Note that the analysis can be optimised for each value of c 6 to obtain a higher significance, but in light of the many sources of uncertainties we do not pursue this here. Such optimisation could substantially alter the shape of the hhh and hh+jets curves in Fig. 5.

D. Sensitivity on the (c3 − d4)-plane
Higgs boson triple production can be used to place constraints on the (c 3 − d 4 )-plane. This can subsequently be used to impose constraints on arbitrary relations between the triple and quartic coefficients in explicit models. We approximate the hhh signal efficiency over the whole plane by calculating its average value for cuts along the plane. The standard deviation on the efficiency obtained this way was found to be ∼ 20% along this direction in the given interval. Considering the magnitude of the uncertainties on the signal and background predictions, we consider this to be adequate at present. For the hh+jets background we use the efficiency fit calculated for the D = 6 EFT case. We show the projected constraints on the (c 3 − d 4 )-plane an integrated luminosity of 30 ab −1 in Fig. 6. As a sanity check, we draw the

V. DISCUSSION AND CONCLUSIONS
Evidently, discovering Standard Model-like triple Higgs boson production will be a challenging task. Our analysis of the hhh → (bb)(bb)(γγ) channel has demonstrated that the process merits serious investigation at a future collider running at 100 TeV proton-proton centreof-mass energy. It is important at this point to emphasise the defining points and caveats that lead this phenomenological analysis to this conclusion: • The detector of an FCC-hh needs to have excellent photon identification and resolution, so that a di-photon invariant mass window of width 2 GeV around the Higgs boson mass can imposed. As we already mentioned, the current resolution at the LHC is 1-2 GeV, [99,100]. Moreover, the projections for photon identification efficiency at the high-luminosity LHC are at O(80%) [101]. It is not unreasonable to expect an improvement in both of these parameters at the FCC-hh, to a resolution of 1 GeV or photon identification of 90%.
• Tagging of b-jets should be extremely good, at least in the range of 70-80%, with excellent light jet rejection of O(1%) over a wide range of transverse momenta and pseudorapidities. Reducing the tagging probability from 80% to 70% would reduce the final number of events in 'true' 4-b-jet final states by about 40%. We note that the expected performance of the b-tagging algorithms for the LHC Run 2 is already at this ballpark [102].
• Any analysis of triple Higgs production that includes bb pairs will also benefit from a very good forward coverage, allowing identification of b-jets up to pseudo-rapidities of |η| ∼ 3.0. Good forward coverage for photons to |η| ∼ 3.5 would also benefit the analysis. For example, the fraction of signal events with two b-jets falling in |η b | ∈ [2.5, 3.0] is ∼ 15% and the fraction of events with two photons falling in |η γ | ∈ [2.5, 3.5] is ∼ 5%. These two are approximately uncorrelated, and thus an LHC-like coverage of |η b | < 2.5, |η γ | < 2.5 would cause a ∼ 20% reduction in signal efficiency compared to the analysis presented in this article.
• Predictions of the triple Higgs boson production cross section, as for the case of double production, posses large theoretical uncertainties at present, due to the unknown higher-order corrections. The best available calculation includes only exact real emission diagrams in combination with 'low-energy theorem' results [15]. A full next-to-leading order calculation will reduce this and allow one to use the process to extract constraints on various models of new physics.
• Crucially, the Monte Carlo event generation of multiple coloured partons (4-6) at next-to-leading or-der, with complicated phase-space cuts, matched to the parton shower, is essential. Technical improvements in this direction, along with increase in computing power, will allow us to perform predictions with reduced theoretical uncertainties, as well as perform analyses of more hhh final states, such as those mentioned in Section II (as well as other processes that involve multiple Higgs bosons).
• Due to the aforementioned theoretical and technical limitations, as well as the unknown characteristics of the future collider, we have not attempted to fully quantify the theoretical uncertainties permeating our results. We expect that future improvements in all of these aspects would allow one to obtain a more reliable quantitative result, including a reasonable expectation of uncertainty.
We note here that our event selection is optimised for the assumed detector performance, and if some of these assumptions are changed, the event selection should also be changed to optimise the signal acceptance and background rejection. Moreover in the scenario that the FCChh performance is substantially worse than what we have assumed, other channels could come into play, such as hhh → (bb)(bb)(τ + τ − ) or hhh → (bb)(τ + τ − )(τ + τ − ).
In conclusion, the study of triple Higgs production should be an important aspect of any future collider pro-gramme. It could provide complementary information on the nature of the Higgs boson and its role in electroweak symmetry breaking, as well as extensions of the Higgs boson sector beyond the standard model. This first baseline study resurrects this process and prompts further investigation into how it can be put into use.