Measuring the Higgs self-coupling via Higgs-pair production at a 100 TeV p-p collider

Higgs pair production provides a unique handle for measuring the strength of the Higgs self interaction and constraining the shape of the Higgs potential. Among the proposed future facilities, a circular 100 TeV proton-proton collider would provide the most precise measurement of this crucial quantity. In this work, we perform a detailed analysis of the most promising decay channels and derive the expected sensitivity of their combination, assuming an integrated luminosity of 30 ab$^{-1}$. Depending on the assumed systematic uncertainties, we observe that the Higgs self-coupling will be measured with a precision in the range 2.9 - 5.5% at 68% confidence level.


Introduction
The steady progress of the LHC experiments keeps improving our knowledge of the Higgs properties [1,2]. The long-term prospects for the high-luminosity phase of the LHC (HL-LHC) set important precision goals [3], reaching the level of few percent for several of the Higgs couplings to gauge bosons and fermions. Beyond this, the per-mille level frontier is opened by a future generation of Higgs factories [4]. The measurement of the Higgs self-coupling, the key parameter controlling the shape of the Higgs potential, will remain however elusive for a long time. Aside from providing clues to the deep origin of electroweak (EW) symmetry breaking (EWSB), the determination of the Higgs potential has implications for a multitude of fundamental phenomena, ranging from the nature of the EW phase transition (EWPT) in the early universe [5], to the (meta) stability of the EW vacuum [6][7][8][9][10]. This measurement sets therefore a primary target among the promised guaranteed deliverables of any future collider programme. Comparative assessments of the potential of different collider options, relying on studies carried out through the years in preparation for their design studies, have recently appeared in two reports [4,11]. The ±50% precision projected for the HL-LHC [3] can be improved by a factor up to 2 at future e + e − colliders [4,12], exploiting the impact of radiative corrections induced by the Higgs self-coupling on single-H production at several energies below the onset of on-shell Higgspair (HH) production [13]. The direct measurement of HH production at √ s ≥ 1 TeV will provide stronger, and independent, measurements, reaching 10% and 9% for the ILC at √ s = 1 TeV [14] and CLIC at √ s = 3 TeV [15], respectively. These measurements will require a longer time scale, as they will be possible only at the last stage of the proposed ILC and CLIC programmes. On these timescales, comparable or even better precision could be possible via the study of HH production at a future high-energy proton-proton (pp) collider, like the 100 TeV Future Circular Collider 1 (FCC-hh [16] or the SPPC [17]).
HH production in hadronic collisions has long been considered as an ideal probe of the Higgs self-coupling [18][19][20], and much work along these lines has been done since the Higgs discovery. Some of the most recent work, in the context of future colliders, is documented in Refs. [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. The best estimates, obtained in these studies, of the sensitivity to the Higgs self-coupling at the FCC-hh have used the bbγγ decay channel, leading to an achievable precision between 5-10%, using this channel alone. A study focusing on the bbτ τ and bbbb final states [30] in the boosted regime achieved a sensitivity of 8% and 20%, respectively. The most up-to-date result, performed by the FCC-hh collaboration [16,34] quotes a precision of 5-7%, driven by the bbγγ channel.
The goal of the present study is to extend the scope of previous projections summarized in Ref. [16] and to provide a refined and comprehensive reference for the combined prospect for the Higgs self-coupling measurement at the FCC-hh. We improve on previous studies and show that further optimization of the most sensitive Higgs decay channels using multivariate techniques is possible. When interpreted in the framework of the Standard Model (SM), the combination of these measurements of HH production allows to reach a precision on the tri-linear Higgs self-coupling in the range δ κ λ = 2.9 − 5.5%, significantly improving previous estimates.
This article is organized as follows. We introduce the theoretical framework, discussing the relation between Higgs self-coupling and HH production, in Section 2, and we present in Section 3 the event generation tools used for this study. The detector modeling, event simulation and analysis frameworks are discussed in Section 4. In Section 5 we introduce the general measurement strategy and the procedure that we use for the signal extraction and to derive the expected precision on the self-coupling. The analyses of the three most sensitive decay channels bbγγ, bbτ τ and bbbb final states and their combination are presented in Section 6. Section 7 summarizes our results and our conclusions.

The theoretical framework
Perturbing the Higgs potential around its minimum, leads to the general expression: where m H is Higgs boson mass and λ 3 and λ 4 are respectively the trilinear and quartic Higgs self-couplings. In the SM the self-couplings are predicted to be λ SM where v is the vacuum expectation value (vev) of the Higgs field. The Higgs vev is known from its relation to Fermi constant, v = ( √ 2G F ) −1/2 = 246 GeV, and the discovery of the Higgs particle at the LHC [37,38] has fixed the last remaining free parameter of the SM, the Higgs mass m H [39]. Beyond the SM, corrections to λ 3 and λ 4 , as well as higher-order terms, are possible.
To this day, large departures from the SM potential are perfectly compatible with current observations [40,41]. This makes it possible, for example, to contemplate BSM models where the modified Higgs potential allows for a strong first order EW phase transition (SFOPT) in the early universe, instead of the smooth cross-over predicted in the SM (for a recent discussion of the interplay between collider observables and models with a SFOPT, see e.g. Ref. [42]). In the context of SM modifications of the Higgs properties [43] parameterized by effective-field-theories (EFTs), it is well known that changes of the Higgs potential are often correlated with changes of other couplings, such as those of the Higgs to the EW gauge bosons. In many instances, a very precise measurement of the latter can be as powerful in constraining new physics as the self-coupling measurement [44]. For example, Ref. [45] considered models for SFOPT with an extra real scalar singlet, and showed that a measurement of the HZZ coupling g HZZ with a precision of ∼ 1% can rule out most of the parameter space that could be probed by a measurement of the self-coupling with a ∼ 50% precision (see Fig. 1 of that paper). Should a deviation from the SM be observed in g HZZ , however, a large degeneracy would be present in the set of allowed parameters. For example, Fig. 1 of Ref. [45] shows that a ∼ 2% deviation in g HZZ would be compatible, in this class of models, with any value of 1 λ 3 /λ SM Another remark is in order: the relation between the Higgs self-coupling and HH production properties is unambiguous only in the SM. Beyond the SM, the HH production rate could be modified not only by a change in the Higgs self-coupling, but also by the presence of BSM interactions affecting the HH production diagrams. These could range from a modified top Yukawa coupling, to higher-order EFT operators leading to local vertices such as ggHH [46], WWHH [28] or ttHH [47,48]. The measurement of an anomalous HH production rate, therefore, could not be turned immediately into a shift of λ 3 ; rather, its interpretation should be made in the context of a complete set of measurements of both Higgs and EW observables, required to pin down and isolate the coefficients of the several operators that could contribute. In view of this, it is not possible to predict an absolute degree of precision that can be achieved on the measurement of λ 3 , since this will depend on the ultimate λ 3 value, on the specific BSM framework leading to that value, and on the ancillary measurements that will be available as additional inputs. As is customary in the literature 2 , we shall therefore focus on the context of the SM, neglecting the existence of interactions influencing the HH production, except for the presence of a pure shift in λ 3 . The precision with which λ 3 can be measured under these conditions has been for a long time the common standard by which the performance of future experiments is gauged, and we adopt here this perspective. Our results remain therefore indicative of the great potential of a hadron collider in the exploration of the Higgs potential.

The theoretical modeling of signals and backgrounds
The signal and background processes are modeled with the MadGraph5_aMC@NLO [51] and Powheg [52,53] Monte Carlo (MC) generators, using the parton distribution functions (PDF) set NNPDF3.0 [54] from the Lhapdf [55] repository. The evolution of the partonlevel events is performed with Pythia8 [56], including initial and final-state radiation (ISR, FSR), hadronization and underlying event (UE). The generated MC events are then interfaced with the Delphes [57] software to model the response of the FCC-hh detector, as described in Section 4.2. The full event generation chain is handled within the integrated FCC collaboration software (Fccsw) [58]. The event yields for the background and signal samples are normalized to the integrated luminosity of L int = 30 ab −1 .

The HH production processes
At √ s = 100 TeV, the dominant HH production modes are, in order of decreasing relative cross section, gluon fusion (ggHH), vector boson fusion (VBF HH), associated production with top pairs (ttHH) and double Higgs-strahlung (VHH). A subset of diagrams for these processes is given in Fig. 1. Single top associated production is also a possible production mode but it is neglected in this study. The cross-section calculations [59][60][61][62][63][64] for these main production mechanisms, reported also in Refs. [11,43,65], are given in Table 1  production becomes as important as vector boson fusion, and together they contribute to nearly 15% of the total HH cross section. The ggHH MC events have been generated at next-to-leading order (NLO) with the full top mass dependence using Powheg [53,66]. The VBF HH, ttH and VHH events were instead generated at leading order (LO) with MadGraph5_aMC@NLO. All the HH production mechanisms feature the interference between diagrams that depend on the self-coupling with diagrams that do not. This leads to a non-trivial total cross section dependence on λ 3 , as shown in Fig. 2(a), and has crucial implications for the self-coupling measurement strategy, as discussed in Section 5. In order to account for this non-trivial dependence of the cross section on the self-coupling, the MC samples for the signal processes have been generated for several possible values of κ λ = λ 3 /λ SM 3 within the interval κ λ ∈ [0.0,3.0]. In order to match the MC inclusive cross section prediction with the cross sections of Table 1, we correct the event normalisation by means of a constant K-factor (shown in the last column of Table 1). We note that in principle the K-factors are κ λ -dependent. However the cross sections at the highest possible accuracy are not known for values of κ λ = 1, therefore the (process dependent) K-factor is derived for κ λ = 1 and applied to correct the cross section at values of κ λ = 1 . The total cross section obtained with this procedure as a function of κ λ is shown in Fig. 2(a). The merging of the NLO parton-level configurations with the parton-shower evolution is realized in the Powheg samples with Pythia8. In Fig. 2(b) the transverse momentum of the HH system p HH T is shown as a validation of the NLO merging procedure. For the VBF HH, ttHH and VHH samples, Pythia8 simply adds the regular parton shower to the LO partonic final states.
The Higgs self-coupling can be probed via a number of different Higgs boson decay channels. Given the small cross section, at least one of the Higgs bosons is required to decay to a pair of b-quarks. Here, we consider the three most promising channels: HH → bbγγ, HH → bbτ τ and HH → bbbb. The di-Higgs system decay in the various modes is performed by the Pythia8 program and the respective branching fractions BR(HH → bbγγ) = 0.00262, BR(HH → bbτ τ ) = 0.072 and BR(HH → bbbb) = 0.33 are taken from Ref. [43], assuming m H = 125.10 GeV.

The background processes
The background processes for the channels under study can be classified in irreducible, reducible and instrumental backgrounds. Irreducible backgrounds feature the presence in the matrix element of the exact same final state as the ggHH signal process. These include for example prompt bbγγ (QCD) production, or Zbb with Z → bb(τ τ ). We define as reducible background the processes that contain the same final state particles as the signal, but also additional particles that can be used as handles for discrimination. This is the case for instance of ttH, H → γγ as a background for the HH → bbγγ channel or the tt background for the HH → bbτ τ channel. Finally, we call as instrumental the background processes that mimic the signal final state due to a mis-reconstruction of the event in the     Figure 2. (a) Cross section of the ggHH, VBF HH, ttHH, and VHH processes as a function of
detector. An instrumental background for the HH → bbγγ channel is the γ + jets process where one the jets gets accidentally reconstructed as an isolated photon. Special care has to be given to such backgrounds as they strongly depend on the details of the detector performance.
Single-Higgs production constitutes a background for all di-Higgs final states. The four main production modes, gluon fusion (ggH), vector boson fusion (VBF H), top pair associated production (ttH) and Higgs-strahlung (VH), have been simulated at LO, including up to two extra MLM-matched jets [67,68], using MadGraph5_aMC@NLO. The ggH matrix element was generated using the full top mass dependence. The rates of single-Higgs processes have been normalised to the most accurate cross-section calculations at √ s = 100 TeV [29]. The normalisation K-factor for the ggH process includes corrections up to N 3 LO, while the VBF H, ttH and VH modes include corrections up to NNLO.
Top-induced backgrounds, in particular top-pair production (tt), constitute a large background for the HH → bbτ τ final state, and to a minor degree for the HH → bbbb final state. This process was generated at LO using MadGraph5_aMC@NLO with up to two extra MLM-matched jets. The total cross section is normalised to match the NNLO prediction at √ s = 100 TeV. The Drell-Yan (Z/γ * +jets) and di-boson backgrounds are also mainly relevant for the HH → bbτ τ and HH → bbbb final states. These are generated at LO with MadGraph5_aMC@NLO by directly requiring the presence of bbτ τ (or jjbb for the bbbb channel) final state at the matrix element level. We generate the pure QCD contribution at order O(α 3 S ). The next contribution, Z/γ * +jets, corresponding to jjbb and bbτ τ was generated at order O(α 2 S α EW ). The latter includes for example the Z → bb(τ τ ) process. The final contribution, generated at order O(α S α 2 EW ), includes the pure EW processes such as ZZ and ZH. When this background is included, the single-Higgs ZH mode discussed earlier is indeed omitted. For the pure QCD contribution we simply assume a conservative K = 2 correction to the MC LO cross section. For the processes at orders O(α EW ) and O(α 2 EW ) we employ K-factors that match to the NNLO Drell-Yan and di-boson √ s = 100 TeV predictions. The last class of relevant background processes for the the HH → bbbb and the HH → bbτ τ final states are the ttZ and ttW processes. These were also generated at LO using MadGraph5_aMC@NLO and normalized the the highest accuracy NLO cross-section calculations.
The largest background contribution for the HH → bbγγ final state are QCD multijet production with one or more prompt photons in the final states, γγ +jets and γ +jets respectively. For the γγ + jets process we generated the matrix element of γγ plus two partons, where partons are generated in the 5-flavour (5F) scheme to allow for mis-reconstructed light and c-quark jets. In order to maximise the MC event efficiency in the signal region, the γγ + jets process was generated with the |m γγ − 125| < 10 GeV requirement at parton level. The γ + jets process was instead generated as γ plus three partons in the final state, again in the 5F scheme. Both these processes were generated at LO and a conservative K=2 correction factor on the LO prediction to account for higher order predictions was applied. The ttγγ process was also considered for this channel and its contribution was found to be negligible.

The experimental and analysis framework
The FCC project is described in detail in its Conceptual Design Reports [69,70]. We focus here on the 100 TeV pp collider, FCC-hh, designed to operate at instantaneous luminosities up to L = 3 × 10 35 cm −2 s −1 . For our study we adopt the reference total integrated luminosity of L int = 30 ab −1 , achieved after 20 years of operations, possibly combining the statistics of two general purpose detectors. The analysis of these data will set challenging requirements to the detector design and performance, which will reflect on the physics potential in general, and in particular on the measurement of the HH cross sections. We summarize here the main features of the current detector design, as implemented in the Delphes [57] simulation tool used for our study.

Detector requirements
A detector operating within the FCC-hh environment will have to be able to isolate the hard-scattering event from up to 1000 pile-up (PU) simultaneous collisions per bunchcrossing. Extreme detector granularity together with high spatial and timing resolution are therefore needed. In addition, to meet the high precision goal in key physics channels such as HH → bbγγ, an excellent photon energy resolution is needed. This requires a small calorimeter stochastic term 3 in an environment of large PU noise, which in turn can be achieved via a large sampling fraction and a fine transverse and longitudinal segmentation. Finally, physics processes occurring at moderate energy scales (Q = 100 GeV − 1 TeV) will be produced at larger rapidities compared to the LHC. Therefore high precision calorimetry and tracking need to be extended up to |η| < 6.
A prototype of a baseline FCC-hh detector that could fulfill the above requirements has been designed by the FCC-hh collaboration [70][71][72]. The detector has a diameter of 20 m and a length of 50 m, with dimensions comparable to the ATLAS detector. A central detector (covering a region up to |η| < 2.5) contains a silicon-based tracker, a Liquid Argon electromagnetic calorimeter (ECAL) and a Scintillating Tile Hadron calorimeter (HCAL) inside a 4 T solenoid with a free bore diameter of 10 m. The muon chambers are based on small Monitored Drift Tube technology (sMDTs). The tracking volume has a radius of 1.7 m with the outermost layer lying at 1.6 m from the interaction point (IP) in the central and the forward regions, providing the full lever arm up to |η| = 3. The ECAL has a thickness of 30 radiation lengths and provides, together with the HCAL, an overall calorimeter thickness of than 10.5 nuclear interaction lengths. The transverse segmentation of both the electromagnetic and hadronic calorimeters is ∼ 4 times finer than the present ATLAS [73] and CMS calorimeters [74]. A high longitudinal segmentation in the ECAL is needed to ensure a high sampling fraction, hence a small stochastic term and in turn the good photon energy resolution required in order to maximise the efficiency of the H → γγ reconstruction. In order to reach good performances at large rapidites (2.5 < |η| < 6), the forward parts of the detector are placed at 10 m from the interaction point along the beam axis. Two forward solenoids with an inner bore of 5 m provide the required bending power for forward tracking. The integrated forward calorimeter system (ECAL and HCAL) is fully based on LAr due to its instrinsic radiation hardness. Coverage up to |η| = 6 is feasible by placing the forward system at a distance z=16.6 m from the IP in the beam direction and at r=8 cm in the transverse direction. The FCC-hh baseline detector performance has been studied in full Geant4 [75] simulations and parameterised within the fast simulation framework Delphes [57,76].

Detector simulation and object reconstruction
The reconstruction of the MC-generated events in the FCC-hh detector is simulated with the Delphes framework. Delphes makes use of a parameterised detector response in the form of resolution functions and efficiencies. The Delphes simulation includes a track propagation system embedded in a magnetic field, electromagnetic and hadron calorimeters, and a muon identification system. Delphes produces physics objects such as tracks, calorimeter deposits and high level objects such as isolated leptons, jets, and missing energy.
Delphes also includes a particle-flow reconstruction that combines tracking and calorimeter information to form particle-flow candidates. Such particles are then used as input for jets clustering, missing energy, and isolation variables. In the following we will focus on describing the key parameters of the FCC-hh detector implementation in Delphes that are relevant for the self-coupling analysis presented here.
Unless specified otherwise, the jets are clustered by the anti-k T algorithm [77] with a parameter R=0.4. For leptons ( = e, µ) and photons (γ), the relative isolation I rel is computed by summing the p T of all particle-flow candidates in a cone around the particle of interest an dividing by the particle's p T (e, µ, γ). Isolated objects, such as photons originating from a HH → bbγγ decay, typically feature a small relative isolation. The reconstruction and identification (ID) efficiencies for leptons and photons are parameterised as function of p T and pseudo-rapidity η. Since a full-fledged event simulation and object reconstruction does not exists at this stage for the FCC-hh detector, the assumed object efficiencies result from extrapolations from the HL-LHC detectors. We simply mention here that a typical photon originating from a HH → bbγγ decay with p T ≈ 50 GeV at η ≈ 0 has a probability γ = 85% of being reconstructed. As mentioned previously, a dominant background for HH → bbγγ analysis is the γ + jets process. The probability for a jet to be misreconstructed as an isolated photon is small O(10 −3 ) in current LHC detectors, thanks to the excellent angular resolution of present calorimeters. As we noted in Section 4.1, the assumed granularity for the FCC-hh detector is a factor 2-4 better than present LHC detectors. We make however the conservative choice of assuming a j → γ fake-rate j→γ = 0.002 · e −p T /(30 GeV) , which is of the same magnitude as in the LHC detectors [78]. For leptons we neglect possible fake jets contributions since these are negligible at the momemtum scale relevant for the HH → bbτ τ final state. Delphes also provides heavy flavour tagging, in particular τ (hadronic) and b-jet identification. Both hadronic τ and b-jets modeling rely on a parameterisation of the (mis-)identification probability as a function of (p T , η). Again, since we cannot yet derive such performance from full-simulation, we assume efficiencies and mistag rates of the same order as in the HL-LHC detectors. For example, the tagging efficiency of central p T ≈ 50 GeV b-jet is b = 85% with mistag-rates respectively for light (l=u,d,s quarks) and c quarks of l→b = 1% and c→b = 5%. For central high p T hadronically decaying τ 's (τ h ), we assume instead an efficiency τ h = 80% and a mistag-rate j→τ h = 1%.
Finally, we note that the effect of pile-up is not simulated directly by overlaying minimum bias events to the hard scattering. Although Delphes allows for such possibility, including in the simulation up to 1000 pile-up interactions would result in an overly conservative object reconstruction performance for the simple reason that the current Delphes FCC-hh setup does not possess well-calibrated pile-up rejection tools that will necessarily be employed for a detector that will operate in such conditions, and so far in the future. These techniques will include the use of picosecond (ps) timing detectors as well as advanced machine learning based techniques for pile-up mitigation. For the present LHC detectors, as well as for presently approved future detectors (the ATLAS and CMS Phase II detectors) it is already the case that such techniques allow to recover the nominal detector performance in the absence of pile-up [79,80]. The level of degradation of the λ 3 measurement precision caused by the deterioration of the performance of specific physics objects (for examples the photon energy resolution, or the b or τ -tagging efficiencies) has been quantified in a previous study [34].

Signal extraction methodology
As mentioned in Section 3.1, the cross section for HH production has a non-trivial dependence on the self-coupling modifier κ λ = λ 3 /λ SM as shown in Fig. 1. In Fig. 1, (T )-diagrams appear on the left column while (S)-diagrams are shown on the right. The (S) and (T ) contributions are present in all HH-production mechanisms. Moreover, the contribution of the interference term between (S) and (T ) is highly non trivial. For the ggHH and VBF HH modes, the total cross section reaches a minimum respectively at κ λ ≈ 2.5 and κ λ ≈ 1.8, while the ttHH and VHH cross section carry little dependence on κ λ .
At first order one can write: where we define µ = σ/σ SM as the signal strength. One can measure λ 3 (or alternatively κ λ ) by measuring the total HH production cross section. It follows that: where δ κ λ and δ µ are respectively the uncertainty on the self-coupling modifier and on the signal strength. It can be noted that at first order the precision of the self-coupling measurement is determined by the slope of the cross section (or µ) at κ λ = 1 and by the uncertainty on the measurement of the total cross section. Since dµ dκ λ SM is a given parameter, in order to maximise the precision on the self-coupling, we have to maximise the precision on the cross section, or equivalently on µ. Assuming all other standard model parameters are known with better precision than the expected precision on κ λ 4 , the relative weight of the (S) and (T ) amplitudes (and their interference) is determined by the magnitude of κ λ . The magnitude of κ λ impacts not only the total HH rate, as discussed above, but also the HH production kinematic observables. Notably, the invariant mass of the HH pair m hh is highly sensitive the value of the self-coupling. This can be easily understood by noting that configurations with large m hh are mostly suppressed in the (S) amplitude (not in (T ) diagrams). Vice-versa, the phase-space region near threshold, at m hh 2m H , maximises the (S) contribution. The m hh distribution is shown for the ggHH and ttHH processes respectively in Figs. 3(a) and 3(b). For ggHH the dependence is distorted at values of κ λ ≈ 2 due to the large destructive interference between (S) and (T ). The transverse momenta of the two Higgs bosons (p T (h 1 ), and p T (h 2 )) also display a large dependence on κ λ , as shown in Figs. 4(a) and 4(b).
The general strategy for providing the best possible accuracy on the self-coupling will therefore rely on maximizing the cross-section precision by using obervables that are able to discriminate between signal and backgrounds as well as exploiting the shapes of observables that are highly sensitive to the value of κ λ . The signal over background optimisation is largely dependent on the class of background and will be addressed in the discussion specific to each channel below. However a common theme is that typically the strategy to  obtain a high S/B ratio relies heavily on the reconstruction of the mass peak of the two Higgs bosons. In addition we will make use of the m hh observable, and the Higgs particles transverse momentum (p T (h 1 ), and p T (h 2 )) differential distributions to further improve the sensitivity on κ λ .

Determination of the Higgs self-coupling
While the Higgs pair can be reconstructed in a large variety of final states, only the most promising ones are considered here: bbγγ, bbτ τ and bbbb. For each of these final states, the event kinematical properties are combined within boosted decision trees (BDTs) to form a powerful single observable that optimally discriminates between signal and backgrounds. The BDT discriminant is built using the ROOT-TMVA package [83,84]. The statistical procedure and the evaluation of the systematic uncertainties are summarized in Appendix A and B, respectively.

The bbγγ channel
Despite its small branching fraction, the HH → bbγγ channel is by far the most sensitive decay mode for measuring the self-coupling. The presence of two high p T photons in the final state, together with the possibility of reconstructing the decay products of both Higgses without ambiguities and with high resolution, provide a clean signature with a large S/B. The largest background processes are single-Higgs production and the QCD continuum γγ + jets and γ + jets. A discussion on the simulation of these processes was given in Section 3.2.

Event selection
In the bbγγ channel, events are required to contain at least two isolated photons and two btagged jets with the requirement p T (γ, b) > 30 GeV and |η(γ, b)| < 4.0. The leading photon and b-jet are further required to have p T (γ, b) > 35 GeV. The Higgs candidates 4-momenta are formed respectively from the two reconstructed b-jets and photons with the largest p T (γ, b). Since the γγ + jets process was generated with a parton-level requirement (see Section 3.2) on m γγ , we further require the events to pass the loose selection |m γγ − 125| < 7 GeV. The efficiency of the full event selection for the SM signal sample is approximately 26%. For an integrated luminosity L int = 30 ab −1 this event selection yields approximately 26k Higgs pair events, 250k single Higgs, 2.5M jjγγ and 3M γ + jets events. The trigger efficiency for the above selection is assumed to be 100% efficient.
In a future FCC-hh experiment, identification algorithms for photon and heavy-flavour will make use of the information of the invariant mass of the photon or jet candidate. Therefore we have to assume that the parameterised performance of the identification efficiency of such objects (in Delphes) already accounts for these variables. As a result, the photon and jet mass are not used as input variables in the BDT discriminant. The m γγ , m bb , and m hh observables, shown respectively in Figs. 5(a), 5(b) and 5(c) provide most of the discrimination against the background.
The QCD (γ+jets and γγ+jets) and single-Higgs background processes possess different kinematic properties, and are therefore treated in separate classes. In QCD backgrounds, the final photons and jets tend to be softer and at higher rapidity. Conversely, the photonpair candidates in single-Higgs processes often originate from a Higgs decay. As a result, while the m γγ observable is highly discriminating against QCD, it is not against single-Higgs processes. In order to maximally exploit these kinematic differences we perform a separate training for each class of backgrounds, that in turn produce two multivariate discriminants: BDT H and BDT QCD . During the training, each background within each class is weighted according to the relative cross section.
The output of the BDT discriminant is shown in the (BDT H , BDT QCD ) plane for the signal and the two background components in Figs. 6(a), 6(b) and 6(c), respectively. As expected, the signal (background) enriched region clearly corresponds to large (small) BDT H and BDT QCD values. We note that the multivariate discriminant correctly identifies the two main components (ggH and ttH) within the single-Higgs background. The ggH background, as opposed to ttH, is more "signal-like" and populates a region of high BDT H and BDT QCD .

Signal Extraction and results
The expected precision on the signal strenth µ = σ/σ SM and on the self-coupling modifier κ λ = λ 3 /λ SM 3 are obtained from a 2-dimensional fit of the (BDT H , BDT QCD ) output, following the procedure described in Appendix A. The results are shown in Figs. 7(a) and 7(b). The various lines correspond to the different systematic uncertainties assumptions described in Appendix B and summarized in Table 3. From Figure 7(b) one can extract the 68% and 95% confidence intervals for the various systematics assumptions. The expected precision for bbγγ is summarized in Table 6.1.2 for each assumption on the systematics. Depending on the assumed scenario, the Higgs self-coupling can be measured with a precision of 3.5-8.5% at 68% C.L using the bbγγ channel alone. We note that the achievable precision is largely dependent on the assumptions on the systematic uncertainties.    Table 2. Expected precision on the Higgs self coupling using the bbγγ channel at the FCC-hh with L int = 30 ab −1 .

The bbτ τ channel
The bbτ τ channel is very attractive thanks to the large branching fraction (7.3%) and the relatively clean final state. As opposed to the bbγγ channel, the HH → bbτ τ decay cannot be fully reconstructed due to the presence of τ neutrinos in the final state. We consider mainly two channels here: the fully hadronic final state bbτ h τ h , and the semi-leptonic one, bbτ h τ ( = e, µ).
As spelled out in Section 3.2, several processes act as background for the bbτ τ final state. The largest background contributions are QCD and tt. QCD is a background mainly for the bbτ h τ h decay channel. However, the absence of prompt missing energy in QCD events makes this background reducible. We have verified that it can be suppressed entirely and therefore has been safely neglected here. In order of decreasing magnitude, the largest backgrounds are Z/γ * +jets single Higgs, ttV and ttVV, where V=W,Z.

Event selection
Events are required to contain at least two b-jets with p T (b) > 30 GeV and |η(b)| < 3.0. For the bbτ h τ final state the presence is required of at least one isolated (I rel < 0.1) lepton = e, µ with p T ( ) > 25 GeV and |η( )| < 3.0 and at least one hadronically tagged τ -jet with p T (τ h ) > 45 GeV and |η(τ h )| < 3.0. For the bbτ h τ h final state, we require at least two hadronically tagged τ -jet with p T (τ h ) > 45 GeV and |η(τ h )| < 3.0. In what follows we refer to a τ -candidate as the lepton = e, µ or the τ -jet. In particular the τ 4-momentum is defined as the sum of the 4-momenta of the visible τ decay products. In order to maximally exploit the kinematic differences between the signal and the dominant tt background, we build a multivariate BDT discriminant using as an input the following kinematic properties: • The 3-vector components of the leading (τ 1 ) and subleading τ -candidate (τ 2 ): transverse momentum (p τ 1 T , p τ 2 T ), pseudo-rapidity (η τ 1 , η τ 2 ), and azimutal angle (φ τ 1 , φ τ 2 ).
• The transverse mass of each τ -candidate, computed as m T = 2p τ T p miss T − p T τ · p T miss .
The m τ τ and m T2 observables are shown in Figs. 8(a),9(a) and 8(b),9(b) for the bbτ h τ h and bbτ h τ final states respectively. These provide the largest discrimination against the tt background. The output of the BDT discriminant is shown in Figs. 8(c) and 9(c) for the bbτ h τ h and bbτ h τ final states.

Signal Extraction and results
The expected precision on the signal strength and the Higgs self-coupling are derived from a maximum likelihood fit on the BDT observable, according to the prescription described in Appendix A. The bbτ h τ h and bbτ h τ channels are considered separately with their relative set of systematic uncertainties and then combined assuming a 100% correlation on equal sources of uncertainties among the two channels. The combined expected precision on the bbτ τ channel is shown in Figs. 10(a) and 10(b). The coloured lines correspond to the different systematic uncertainties assumptions summarized in Table 3. From Fig. 10(b) one can extract the 68% and 95% confidence intervals for the various systematics assumptions. Depending on the assumed scenario, using the bbτ τ channel, the Higgs pair signal strength and Higgs self-coupling can be measured respectively with a precision of δ µ = 6% and δ κ λ = 12 − 13% at 68% C.L. Despite the large signal event rate in the bbτ τ channel, the sensitivity is limited by the overwhelming background contribution. Therefore, contrary to the bbγγ case, the bbτ τ channel is statistically dominated at the FCC-hh with L int = 30 ab −1 , and the achievable precision is only moderately dependent on the assumptions on the systematic uncertainties.  Figure 10. Expected negative log-Likelihood scan as a function the signal strenth µ = σ/σ SM (a) and trilear self-coupling modifier κ λ = λ 3 /λ SM 3 (b) in the bbτ τ channel (combination of the bbτ h τ h and bbτ h τ channels). The various lines correspond to the different systematic uncertainties assumptions summarized in Table 3.

The bbbb channel
The HH → bbbb decay mode has the largest branching fraction among all possible Higgs pair decays. Despite the presence of soft neutrinos from semi-leptonic b decays (that may impact negatively the reconstructed hadronic Higgs mass resolution), the Higgs decays into b-jets can be fully reconstructed. However due do the fully hadronic nature of this decay mode, this channel suffers from the presence of overwhelming QCD backgrounds and hence features a relatively small S/B. Moreover, a combinatorial ambiguity affects the possibility to correctly associate the four b-jets to the two parent Higgs candidates.
We consider mainly the case where the Higgs candidates are only moderately boosted, leading to four fully resolved b-jets. The boosted analysis, where the Higgs candidates are sufficiently boosted to decay into a single large radius jet [87,88], provides less sensitivity to the self-coupling measurement and was discussed in previous studies [30,34]. The main backgrounds to this final state are QCD and tt, followed by Zbb, single-Higgs production and ZZ.

Event selection
In order to fulfill our initial assumption of fully efficient online triggers, the event selection starts by requiring the presence of at least four b-jets with p T (b) > 30 GeV and |η(b)| < 4.0. The Higgs candidates are reconstructed as the pairing of b-jet pairs that minimizes the difference between the invariant masses of the two b-jet pairs. The Higgs candidate with the largest (smallest) p T is named h 1 (h 2 ).
The following variables are then used as input to a multivariate BDT discriminant to ensure an optimal discrimination versus the dominant QCD background: • The 4-vector components of the leading (h 1 ) and subleading (h 2 ) Higgs candidates: • The 4-vector components of the Higgs-pair candidate: transverse momentum (p hh T ), pseudo-rapidity (η hh ), azimutal angle (φ hh ) and invariant mass (m hh ).
The m h 1 and m hh observables are shown respectively in Figs. 11(a) and 11(b). The m h 1 distribution shows that the procedure described above correctly associates the b-jet pairs to the parent Higgs particle for the signal, and to the parent Z particle for the Zbb and ZZ backgrounds. Thanks to their resonant nature, the m h 1 and m h 2 distributions provide the largest discrimination against the QCD background. The output of the BDT discriminant is shown in Fig. 11(c) for the signal and various background contributions.

Results
The expected precision on the signal strength and the Higgs self-coupling are derived from a 1D maximum likelihood fit on the BDT discriminant, according to the prescription described in Appendix A. The combined expected precision on the bbbb channel is shown in Figs. 12(a) and 12(b). The coloured lines correspond to the different systematic uncertainties assumptions summarized in Table 3. The 68% and 95% confidence intervals on δ µ and δ κ λ for the various systematics assumptions can be extracted from Figs. 10(a) and 10(b). Depending on the assumed scenario, using the bbbb channel, the Higgs pair signal strenth and Higgs self-coupling can be measured respectively with a precision of δ µ = 9 − 10% and δ κ λ = 24 − 26% at 68% C.L. As for the bbτ τ case, due to the overwhelming QCD background this channel is statistically limited, and as such the achievable precision only moderately depends on the assumptions on the systematic uncertainties.   Table 3.

Combination procedure and results
When combining results from the various channels, the systematic uncertainties from the various sources are accounted for as follows. Lepton, jets, b-jets, and photon reconstruction and identification efficiencies are assumed to be fully correlated across the channels and analysis categories that use the same objects. All sources of systematic uncertainties are assumed to only affect the overall normalization of signal and background shapes and to not introduce a significant deformation of their shapes.
The combined expected negative log-Likelihood scan is shown in Fig. 13. The expected precision for the single channels is also shown. For completeness, we introduced in the   combination also the bbZZ (4 ) channel, which provides a sensitivity similar to the 4b channel. This decay channel was not re-optimized in this study and the result of the analysis is documented in Ref [34]. The expected combined precision on the Higgs self-coupling obtained after combining the channels bbγγ, bbτ τ , bbbb and bbZZ (4 ) can be inferred from the intersection of black curves with the horizontal 68% and 95% CL lines. The expected statistical precision, neglecting systematic uncertainties, can be read from the dashed black line in Fig. 13, and gives δ κ λ = 2.2% at 68% CL. The solid line corresponds to scenario II for systematic uncertainties, while the boundaries of the shaded area represent respectively the alternative scenarios I and III. From the shaded black curve one can infer the final precision when including systematic uncertainties. Depending on the assumptions, the expected precision for the Higgs self-coupling is δ κ λ = 2.9 − 5.5% at 68% CL.
The expected precision on the Higgs self-coupling as a function of the integrated luminosity is shown in Fig. 14, for the three scenarios of systematic uncertainty. Even with the most conservative scenario, a precision of δ κ λ = 10% can be reached with only 3 ab −1 of integrated luminosity (less than 2 ab −1 are sufficient with the central reference systematics scenario). The 10% target should therefore be achievable during the first 5 years of FCChh operations, combining the datasets of two experiments. Even including the duration of the FCC-ee phase of the project, and the transition period from FCC-ee to FCC-hh, this timescale is competitive with the time required by the proposed future linear colliders, which to achieve this precision need to complete their full programme at the highest beam energies.
As already discussed, the value of the self-coupling coupling can significantly alter both the Higgs pair production cross section and the event kinematic properties. In order to explore the sensitivity to possible BSM effects in Higgs pair production, a multivariate BDT discriminant was optimised against the backgrounds for several values of κ λ in the range 0 < κ λ < 3, in order to maximise the achievable precision for values of κ λ = 1. The BDT training has been performed only for the bbγγ channel, which dominates the overall sensitivity. In order to obtain the combined sensitivity for κ λ = 1, we re-scaled the error obtained at κ λ = 1 for the bbγγ channel only, by the ratio, evaluated at κ λ =1, of the bbγγ-channel uncertainty to the fully combined uncertainty. The obtained precision as a function of κ λ is shown in Fig. 15 5 . When ignoring systematic uncertainties, it can be seen that the overall precision follows the behaviour of the HH production cross section. The best precision (δ κ λ ≈ 0.02) is reached at κ λ = 0, while the maximum uncertainty at κ λ = 2.4 corresponds to the minimum of the total HH cross section (δ κ λ ≈ 0.10). It can also be noticed that when switching on the systematic uncertainties the precision at small κ λ degrades compared to the SM case. This reflects the fact that the HH kinematics at κ λ ≈ 0 are similar to the single-Higgs background. Large correlated uncertainties (such as b-tagging and photon identification efficiencies) between double and single-Higgs production thus have a larger impact at κ λ ≈ 0 as opposed to the SM case.

Discussion
The combined precision improves considerably if compared to the expected precision in the dominant channel bbγγ. Given the large amount of expected data at the FCC, for some processes, such as the QCD contribution to the bbbb channel, the likelihood fit is able to constrain the uncertainty on the overall normalization to smaller levels than the assumed pre-fit uncertainties affecting the given process. This is especially evident in background rich regions, i.e at low values of the BDT score. In these regions, the signal contribution is small, and the fit to the overall normalization is entirely driven by the background rate, which is therefore precisely determined by the data. The resulting effect is a sizeable reduction of the total contribution of systematic uncertainties on the expected precision compared to the sensitivity obtained with simple error propagation using a weighted average procedure.
More concretely the most relevant source of systematic uncertainty in the bbγγ channel is the uncertainty on b-tagging, which can be as large as 2% for each b-jet. The available statistics in this channel are not sufficient to be able to constrain this dominant source of uncertainty. Nevertheless, the much larger expected statistics in the background-dominated phase space region of the bbbb channel (at small BDT values) are responsible for reducing the value of the b-tag uncertainty down to just 5% of its original value. Indeed a (fully correlated) 8% variation on the 4 b-jets background normalization would be much larger than any possible statistical fluctuation, given the amount of data expected by the end of the FCC-hh data taking program. Since systematic uncertainties are correlated across all channels when performing the combined fit, this constraint is carried from the bbbb channel to all other channels, and in particular to the bbγγ channel that dominates the final sensitivity, which in turn translates in a significant overall reduction of uncertainty.

Conclusions and perspectives
The precise measurement of the Higgs self-coupling must be a top priority of future highenergy collider experiments. Previous studies on the potential of a 100 TeV pp collider have discussed the sensitivity of various decay channels, often based on simple rectangular cut-based analyses 6 . In the present study the measurement strategy has been optimized in the bbγγ, bbτ τ and bbbb channels using machine learning techniques. For the first time, a precise set of assumptions of possible sources of systematic uncertainties has been defined and used to derive the achievable precision. Consistently with our previous findings, the bbγγ channel drives the final sensitivity, with an expected precision of δ κ λ = 3.5 − 8.5%. The bbτ τ and bbbb channels provide instead a less precise single channel measurement, respectively of δ κ λ = 12% and δ κ λ = 25%. Contrary to naive expectations, the contribution of the least sensitive channel (bbbb) in the combined expected precision is far from negligible. Thanks to the large available statistics in background enriched regions, these channels significantly help in constraining and reducing the overall impact of systematic uncertainties (that are correlated among all channels) such as theoretical uncertainties on the production cross section, the luminosity, and heavy-flavour jet identification efficiency.
The final combined sensitivity across all considered channels leads to an expected precision at the FCC-hh on the Higgs self-coupling δ κ λ = 2.9−5.5% with an integrated luminosity of L int = 30 ab −1 . The 10% threshold can be achieved with ∼ 2 − 3 ab −1 , depending on the systematics scenario, corresponding to ∼ 3 − 5 years of early running at the start-up luminosity of 5 × 10 34 cm −2 s −1 .
This work shifts the perspective on the Higgs self-coupling ultimate precision at FCC from being statistics dominated to systematics dominated. This is a crucial new development: on one side it gives us confidence that the design parameters of FCC-hh are well tailored to reach, unique among all proposed future collider facilities, the few-percent level of statistical precision. On the other hand, it calls for a more thorough assessment of all systematic uncertainties. For example, we should validate the estimates presented in this work through full simulations of more realistic detector designs in the presence of pile-up, and explore all other possible handles to further reduce them. The huge FCC statistics will provide multiple control samples, well beyond what discussed in our paper, that could be used for these purposes, and in particular to pin down the background rates with limited reliance on theoretical calculations. At this level of precision, however, the theoretical uncertainties on the HH signal will play an important role in the extraction of the self-coupling from the measured production rate. As indicated in Appendix B, we assumed in our study an uncertainty on the HH cross sections ranging from 0.5 to 1.5%. This would need the theoretical predictions to improve relative to today's knowledge, requiring the extension of the perturbative order by at least one order beyond today's known NLO with full top mass dependence [62], and possibly beyond the N 3 LO in the m top → ∞ limit, recently achieved [89,90]. This will be very challenging, and it is impossible today to estimate the asymptotic reach in theoretical precision. Nevertheless, the innovative technical progress we have witnessed in the recent years encourages us to assume that the necessary improvements are possible within the several decades that separate us from the first FCC-hh run.
In conclusion, this study strengthens the evidence that a 100 TeV pp collider, with integrated luminosity above 3 ab −1 , can measure the Higgs self-coupling more precisely than any other proposed project, on a competitive time scale.

A Statistical procedure
The statistical methodology used in this paper relies on the strategy adopted by the AT-LAS and CMS Collaborations, and described in Ref. [91]. A detailed description of the procedures used in this paper are described in more detail in Refs. [38,92]. The Combine software package [93] has been used as statistical and fitting tool to produce the final results.
Combine is based on the standard LHC data modeling and handling toolkits RooFit [94] and RooStats [95] and it is developed and maintained by the CMS collaboration.
The parameter of interest (POI) tested in these results are either the trilinear coupling modifier κ λ = λ 3 /λ SM 3 or the double Higgs signal strength µ = σ/σ SM , defined as the ratio between the (expected) measured double Higgs yield and its SM expectation.
In the model, the POI α = k λ or α = µ is estimated with its corresponding confidence intervals using a profile likelihood ratio test statistic q(α) [91,96], in which experimental or theoretical uncertainties are incorporated via nuisance parameters (NP). Given a of POI α that depends on the set of NP θ, q(α) is defined as: (A.1) An individual NP represents a single source of systematic uncertainty. Its effect is therefore considered fully correlated between all of the the final states included in the fit that share a dependency on it, as will be discussed later in this section.
The quantitiesα andˆ θ denote the unconditional maximum likelihood estimates of the parameter value, whileˆ θ α denotes the conditional maximum likelihood estimate for fixed values of the POI α. The likelihood functions in the numerator and denominator of Eq. (A.1) are constructed using products of probability density functions (PDFs) of signal and background for the various discriminating variables used in the input analyses, as well as constraint terms for the NPs. The PDFs are built from the BDT distributions described in Section 6. It should be noted that while the signal shape depends on the value on κ λ , that dependence is relatively mild. Given the expected precision of O(10%) on the measurement of k λ at the FCC-hh, the effect of its variation on the signal lineshape is minimal when performing the measurement at a given value of κ λ and can be safely neglected. The effects on acceptance and selection efficiencies of varying k λ or µ are instead considered in the fit. The expected precision on κ λ and µ is assessed by performing a likelihood fit on a a pseudo-data set that has been constructed assuming µ = 1 and k λ = 1, using the asymptotic approximation as described in [96,97].  Table 3. Summary of the various sources of systematics. Upper part: contributions to the uncertainty in the measurement of the cross sections of the processes listed in the last column. Lower row: theoretical uncertainty of the total HH rate assumed for the extraction of µ and κ λ from the measured cross section.

B Systematic uncertainties
Systematic uncertainties can play a major role on the expected sensitivity of the selfcoupling measurement. Several assumptions have been made on the possible evolution of theoretical and experimental sources of uncertainties in order to present a realistic estimate of the physics potential of FCC-hh for the channels considered here. In particular, for each uncertainty source, we defined three possible scenarios. An intermediate scenario that we use as a reference point (II), and an optimistic (I) and a conservative (III) scenario. We note that the intermediate assumptions are almost equivalent to those made for HL-LHC projections [3,4]. A detailed list of the systematic uncertainties considered is presented in Table 3 for all the channels, together with the processes affected by each uncertainty. The numbers in the Table refer to the individual contributions to the overall yield uncertainty. In particular, we consider uncertainties on: • background normalisation, which we assume to be dominated by the uncertainty on the experimental measurement of the tt production and is varied between 0.5% and 1.5%.
• luminosity. We assume that the integrated luminosity will be known at FCC-hh at least as well as at the LHC. For this reason, we assume a conservative estimate of 2% and an optimistic (intermediate) estimate of 0.5% (1%).
• experimental uncertainties on objects reconstruction and identification efficiencies: b-jets: for each b-jet, we assume a 0.5%, 1.0%, and 2.0% uncertainty for the optimistic, intermediate and conservative scenarios, respectively.
τ -jets: for each jet originated from the hadronic decay of a τ -jet we assume an uncertainty of 1.0%, 2.0% and 5.0% for the optimistic, intermediate and conservative scenarios, respectively.
leptons: we assume the same uncertainty on the lepton identification and reconstruction efficiency for electrons and muons: a 0.5%, 1.0%, and 1.5% uncertainty for the optimistic, intermediate and conservative scenarios, respectively.
photons: we assume that photons performances will be comparable to the electron ones. For this reason we assign a systematic uncertainty to photon reconstruction of 0.5%, 1.0%, and 2.0% for the optimistic, intermediate and conservative scenarios, respectively.
Furthermore, when interpreting the results in terms of µ and κ λ , we include the additional contribution from the assumed theoretical uncertainty on the HH cross section, as specified in the bottom row of Table 3. We assume that several backgrounds will be measured with high statistical accuracy from "side bands". This is the case for example for the QCD and non-single-Higgs backgrounds that dominate the bbbb and bbγγ channels background contributions. In these cases, while there is no uncertainty associated to the normalization of the backgrounds, the statistical uncertainty due to the possible fluctuation of the number of events in the sidebands is considered in the fit.
When performing the fit for the combination across different channels, systematic uncertainties of same physical origin are considered fully correlated across processes and final states. Otherwise they are considered as completely uncorrelated.