1 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 \(\pm 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 Higgs-pair (HH) production [13]. The direct measurement of HH production at \(\sqrt{s}\ge 1\) TeV will provide stronger, and independent, measurements, reaching 10% and 9% for the ILC at \(\sqrt{s}=1\) TeV [14] and CLIC at \(\sqrt{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.Footnote 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,37,38,39,40,41]. The best estimates, obtained in these studies, of the sensitivity to the Higgs self-coupling at the FCC-hh have used the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) decay channel, leading to an achievable precision between 5 and 10%, using this channel alone. A study focusing on the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) final states [31] 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, 37] quotes a precision of 5–7%, driven by the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) 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 multi-variate 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 \(\,\delta _{\,\kappa _{\lambda }} =3.4{-}7.8\%\), 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 Sect. 2, and we present in Sect. 3 the event generation tools used for this study. The detector modeling, event simulation and analysis frameworks are discussed in Sect. 4. In Sect. 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 \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) final states and their combination are presented in Sect. 6. Section 7 summarizes our results and our conclusions.

2 The theoretical framework

Perturbing the Higgs potential around its minimum, leads to the general expression:

$$\begin{aligned} {\mathcal {L}}_{h} = \frac{1}{2}m_\mathrm {H} ^2 H^2 + \,\lambda _3 H^3 + \,\lambda _4 H^4, \end{aligned}$$
(2.1)

where \(m_\mathrm {H} \) is Higgs boson mass and \(\lambda _3\) and \(\lambda _4\) are respectively the trilinear and quartic Higgs self-couplings. In the SM the self-couplings are predicted to be \(\lambda _3^{SM}={m_\mathrm {H} ^2}/{2v}\), \(\lambda _4^{SM} = {m_\mathrm {H} ^2}/{8v^2}\), where v is the vacuum expectation value (vev) of the Higgs field. The Higgs vev is known from its relation to Fermi constant, \(v=(\sqrt{2}G_F)^{-1/2}=246\) GeV, and the discovery of the Higgs particle at the LHC [42, 43] has fixed the last remaining free parameter of the SM, the Higgs mass \(m_\mathrm {H}\) [44]. Beyond the SM, corrections to \(\,\lambda _3 \) and \(\,\lambda _4 \), as well as higher-order terms, are possible.

To this day, large departures from the SM potential are perfectly compatible with current observations [45, 46]. 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. [47]). In the context of SM modifications of the Higgs properties [48] 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 [49]. For example, Ref. [50] 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 \(\sim 1\%\) can rule out most of the parameter space that could be probed by a measurement of the self-coupling with a \(\sim 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. [50] shows that a \(\sim 2\%\) deviation in \(g_{HZZ}\) would be compatible, in this class of models, with any value of \(1\lesssim \,\lambda _3/\,\lambda _3 ^{SM}\lesssim 2\). A precise direct measurement of \(\,\lambda _3 \) is therefore necessary, independently of what other observables could possibly probe, and is an indispensable component of the Higgs measurement programme.

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 [51], WWHH [29] or \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) [52, 53]. The measurement of an anomalous HH production rate, therefore, could not be turned immediately into a shift of \(\,\lambda _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 \(\,\lambda _3 \), since this will depend on the ultimate \(\,\lambda _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,Footnote 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 \(\,\lambda _3 \). The precision with which \(\,\lambda _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.

Fig. 1
figure 1

Diagrams contributing to Higgs pair production: a gluon fusion, b vector-boson fusion, c double Higgs-strahlung and d double Higgs bremsstrahlung off top quarks. The trilinear Higgs self-coupling is marked in red

3 The theoretical modeling of signals and backgrounds

The signal and background processes are modeled with the \(\textsc {MadGraph5}\)_aMC@NLO [56] and \(\textsc {Powheg}\) [57, 58] Monte Carlo (MC) generators, using the parton distribution functions (PDF) set \(\textsc {NNPDF3.0}\) [59] from the \(\textsc {Lhapdf}\) [60] repository. The evolution of the parton-level events is performed with \(\textsc {Pythia8}\) [61], including initial and final-state radiation (ISR, FSR), hadronization and underlying event (UE). The generated MC events are then interfaced with the \(\textsc {Delphes}\) [62] software to model the response of the FCC-hh detector, as described in Sect. 4.2. The full event generation chain is handled within the integrated FCC collaboration software (\(\textsc {Fccsw}\)) [63]. The event yields for the background and signal samples are normalized to the integrated luminosity of  \(\mathcal {L}_{int}=30\text { ab}^{-1}\).

3.1 The HH production processes

At  \(\sqrt{s}=\text {100 TeV}\), the dominant HH production modes are, in order of decreasing cross section, gluon fusion (\(\mathrm {g g H H}\)), vector boson fusion (\(\mathrm {VBF\,HH}\)), associated production with top pairs (\(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\)) and double Higgs-strahlung (\(\mathrm {V H H}\)). 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 [64,65,66,67,68,69,70,71,72,73] for these main production mechanisms, reported also in Refs. [11, 48, 74], are given in Table 1. We note that the relative rate of the sub-dominant modes (\(\mathrm {VBF\,HH}\), \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) and \(\mathrm {V H H}\)) increases significantly from  \(\sqrt{s}=\text {14 TeV}\) to  \(\sqrt{s}=\text {100 TeV}\). In particular, associated top pair production becomes as important as vector boson fusion, and together they contribute to nearly 15% of the total HH cross section.

Table 1 Signal cross sections (\(\sigma \), in fb) for HH production, including the QCD corrections recommended by the LHC Higgs Cross Section Working Group [48, 74]. For each process, scale variations have been symmetrized and added in quadrature to PDF+ \(\alpha _\text {S}\) uncertainties. For the \(\mathrm {g g H H}\) process, we added in quadrature also the dominant uncertainty induced by the finite \(m_{top}\) corrections. The cross sections of \(\mathrm {W}^{-} \mathrm {H H}\), \(\mathrm {W}^{+} \mathrm {H H}\) and \(\mathrm {Z H H}\) processes have been summed together in a single \(\mathrm {V H H}\) line and their uncertainties have been summed in quadrature

The \(\mathrm {g g H H}\) MC events have been generated at next-to-leading order (NLO) with the full top mass dependence using \(\textsc {Powheg}\) [58, 75]. The \(\mathrm {VBF\,HH}\), \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) and \(\mathrm {V H H}\) events were instead generated at leading order (LO) with \(\textsc {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  \(\lambda _3\), as shown in Fig. 2a, and has crucial implications for the self-coupling measurement strategy, as discussed in Sect. 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  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) within the interval \(\,\kappa _{\lambda } \in \) [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  \(\kappa _{\lambda }\)-dependent. In this work, the (process dependent) K-factor is derived for \(\,\kappa _{\lambda } =1\) and applied to correct the cross section at values of \(\,\kappa _{\lambda } \ne 1\). This is justified by the explicit calculation of the \(\hbox {N}^3\)LO corrections at \(\,\kappa _{\lambda } \ne 1\) for the VBF production channel [69], and by the study of the \(\,\kappa _{\lambda } \) dependence of the NNLO/NLO ratio for \(\mathrm {g g H H}\) in Ref. [76]. In the latter case, the shape variation of kinematical distributions for \(\,\kappa _{\lambda } \ne 1\) from NLO to NNLO is small compared to the overall size of the difference between \(\,\kappa _{\lambda } \ne 1\) and \(\,\kappa _{\lambda } = 1\). The total cross section obtained with this procedure as a function of  \(\kappa _{\lambda }\) is shown in Fig. 2a. The merging of the NLO parton-level configurations with the parton-shower evolution is realized in the \(\textsc {Powheg}\) samples with \(\textsc {Pythia8}\). In Fig. 2b the transverse momentum of the HH system \(p_{\text {T}}^{\mathrm {HH}}\) is shown as a validation of the NLO merging procedure. For the \(\mathrm {VBF\,HH}\), \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) and \(\mathrm {V H H}\) samples, \(\textsc {Pythia8}\) simply adds the regular parton shower to the LO partonic final states.

Fig. 2
figure 2

a Cross section of the \(\mathrm {g g H H}\), \(\mathrm {VBF\,HH}\), \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\), and \(\mathrm {V H H}\) processes as a function of  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\). b Transverse momentum spectrum of the HH system in \(\mathrm {g g H H}\) NLO events after parton-shower merging for \(\,\kappa _{\lambda } =0\), \(\,\kappa _{\lambda } =1\), \(\,\kappa _{\lambda } =2\) and \(\,\kappa _{\lambda } =3\)

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: \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\). The di-Higgs system decay in the various modes is performed by the \(\textsc {Pythia8}\) program and the respective branching fractions \(\,\text {BR}({\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma }) =0.00262\), \(\,\text {BR}({\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau }) =0.072\) and \(\,\text {BR}({\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}}) =0.33\) are taken from Ref. [48], assuming \(m_\mathrm {H}=125.10\,\,\text {Ge}\text {V}\).

3.2 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 \(\mathrm {g g H H}\) signal process. These include for example prompt \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) (QCD) production, or \(\mathrm {Z b}\bar{\mathrm {b}}\) with \(\mathrm {Z}\rightarrow \mathrm {b}\bar{\mathrm {b}} (\,\tau \tau )\). 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 \(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\), \(\mathrm {H} \rightarrow \gamma \gamma \) as a background for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel or the \(\mathrm {t}\bar{\mathrm {t}}\) background for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) 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 detector. An instrumental background for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel is the \(\gamma +\mathrm {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 (\(\mathrm {ggH}\)), vector boson fusion (\(\mathrm {VBF\,H}\)), top pair associated production (\(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\)) and Higgs-strahlung (\(\mathrm {VH}\)), have been simulated at LO, including up to two extra MLM-matched jets [77, 78], using \(\textsc {MadGraph5}\)_aMC@NLO. The \(\mathrm {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  \(\sqrt{s}=\text {100 TeV}\) [30]. The normalisation K-factor for the \(\mathrm {ggH}\) process includes corrections up to \(\mathrm {N}^3 \mathrm {LO} \), while the \(\mathrm {VBF\,H}\), \(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\) and \(\mathrm {VH}\) modes include corrections up to \(\mathrm {NNLO} \).

Top-induced backgrounds, in particular top-pair production (\(\mathrm {t}\bar{\mathrm {t}}\)), constitute a large background for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) final state, and to a minor degree for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) final state. This process was generated at LO using \(\textsc {MadGraph5}\)_aMC@NLO with up to two extra MLM-matched jets. The total cross section is normalised to match the \(\mathrm {NNLO} \) prediction at  \(\sqrt{s}=\text {100 TeV}\). The Drell–Yan (\(\text {Z}/\gamma ^*\text {+jets}\)) and di-boson backgrounds are also mainly relevant for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) final states. These are generated at LO with \(\textsc {MadGraph5}\)_aMC@NLO by directly requiring the presence of \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) (or \(\mathrm {jj}\mathrm {b}\bar{\mathrm {b}}\) for the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel) final state at the matrix element level. The pure QCD contribution, \(\mathrm {jj}\mathrm {b}\bar{\mathrm {b}}\)  has also been generated with \(\textsc {MadGraph5}\)_aMC@NLO at order \({\mathcal {O}}(\,\alpha _\text {S} ^3)\). The next contribution, \(\text {Z}/\gamma ^*\text {+jets}\), corresponding to \(\mathrm {jj}\mathrm {b}\bar{\mathrm {b}}\) and \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) was generated at order \({\mathcal {O}}(\,\alpha _\text {S} ^2 \,\alpha _\text {EW})\). The latter includes for example the \(\mathrm {Z}\rightarrow \mathrm {b}\bar{\mathrm {b}} (\,\tau \tau )\) process. The final contribution, generated at order \({\mathcal {O}}(\,\alpha _\text {S} \,\alpha _\text {EW} ^2)\), 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 \({\mathcal {O}}(\,\alpha _\text {EW})\) and \({\mathcal {O}}(\,\alpha _\text {EW} ^2)\) we employ K-factors that match to the \(\mathrm {NNLO} \) Drell-Yan and di-boson  \(\sqrt{s}=\text {100 TeV}\) predictions. The last class of relevant background processes for the the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) and the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) final states are the \(\mathrm {ttZ}\) and \(\mathrm {ttW}\) processes. These were also generated at LO using \(\textsc {MadGraph5}\)_aMC@NLO and normalized to the highest accuracy NLO cross-section calculations.

The largest background contribution for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) final state are QCD multijet production with one or more prompt photons in the final states, \(\gamma \gamma +\mathrm {jets}\) and \(\gamma +\mathrm {jets}\) respectively. For the \(\gamma \gamma +\mathrm {jets}\) process we generated the matrix element of \(\gamma \gamma \) plus two partons with \(\textsc {MadGraph5}\)_aMC@NLO, 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 \(\gamma \gamma +\mathrm {jets}\) process was generated with the \(|m_{\gamma \gamma }- 125| < 10\) GeV requirement at parton level. The \(\gamma +\mathrm {jets}\) process was instead generated as \(\gamma \) plus three partons in the final state, again in the 5F scheme. Both these processes were generated at LO and a conservative \(\hbox {K}=2\) correction factor on the LO prediction to account for higher order predictions was applied. The \(\mathrm {t}\bar{\mathrm {t}}\gamma \gamma \) process was also considered for this channel and its contribution was found to be negligible.

4 The experimental and analysis framework

The FCC project is described in detail in its Conceptual Design Reports [79, 80]. We focus here on the 100 TeV pp collider, FCC-hh, designed to operate at instantaneous luminosities up to \({\mathcal {L}}=3\times 10^{35}~\text{ cm}^{-2}~\text{ s}^{-1}\). For our study we adopt the reference total integrated luminosity of  \(\mathcal {L}_{int}=30\text { ab}^{-1}\), to be 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 \(\textsc {Delphes}\) [62] simulation tool used for our study.

4.1 Detector requirements

A detector operating within the FCC-hh environment will have to isolate the hard-scattering event from up to 1000 pile-up (PU) simultaneous collisions per bunch-crossing. 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 \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), an excellent photon energy resolution is needed. This requires a small calorimeter stochastic termFootnote 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\,\,\text {Ge}\text {V}{-}1\,\,\text {Te}\text {V}\)) will be produced at larger rapidities compared to the LHC. Therefore high precision calorimetry and tracking need to be extended up to \(|\eta | <6\).

A prototype of a baseline FCC-hh detector that could fulfill the above requirements has been designed for the FCC CDR [80,81,82]. 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 \(|\eta | < 2.5\)) contains a silicon-based tracker, a Liquid Argon (LAr) 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 \(|\eta | = 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 \(\sim 4\) times finer than the present ATLAS [83] and CMS calorimeters [84]. 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 \(\mathrm {H} \rightarrow \gamma \gamma \) reconstruction. In order to reach good performances at large rapidites (\(2.5< |\eta | < 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 \(|\eta | = 6\) is feasible by placing the forward system at a distance \(\hbox {z}=16.6\) m from the IP in the beam direction and at \(\hbox {r}=8\) cm in the transverse direction. The FCC-hh baseline detector performance has been studied in full \(\textsc {Geant4}\) [85] simulations and parameterised within the fast simulation framework \(\textsc {Delphes}\)  [62, 86].

4.2 Detector simulation and object reconstruction

The reconstruction of the MC-generated events in the FCC-hh detector is simulated with the \(\textsc {Delphes}\) framework. \(\textsc {Delphes}\) makes use of a parameterised detector response in the form of resolution functions and efficiencies. The \(\textsc {Delphes}\) simulation includes a track propagation system embedded in a magnetic field, electromagnetic and hadron calorimeters, and a muon identification system. \(\textsc {Delphes}\) produces physics objects such as tracks, calorimeter deposits and high level objects such as isolated leptons, jets, and missing energy. \(\textsc {Delphes}\) also includes a particle-flow reconstruction that combines tracking and calorimeter information to form particle-flow candidates, i.e. charged hadrons, neutral hadrons and photons. Such particles are then used as input for jet clustering, missing energy, and isolation variables. In the following we will focus on describing the key parameters of the FCC-hh detector implementation in \(\textsc {Delphes}\) that are relevant for the self-coupling analysis presented here.

Jets are clustered by the  \(\text {anti-}k_{\text {T}} \) algorithm [87] with a parameter \(\hbox {R}=0.4\). For leptons (\(\ell =e,\mu \)) and photons (\(\gamma \)), the relative isolation  \(I_\text {rel}\) is computed by summing the \(p_{\text {T}}\) of all particle-flow candidates in a cone around the particle of interest an dividing by the particle’s \(p_{\text {T}} (e,\mu ,\gamma )\). Isolated objects, such as photons originating from a \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) decay, typically feature a small value of  \(I_\text {rel}\). The reconstruction and identification (ID) efficiencies for leptons and photons are parameterised as function of \(p_{\text {T}}\) and pseudo-rapidity \(\eta \).

Table 2 Performance of physics objects for the various scenarios. Objects efficiencies and mistag rates are given for a representative \(p_{\text {T}} \approx 50\,\text {Ge}\text {V}\). For b and \(\tau \)-tagging (and their respective mistag rates) numbers for two different working points are given (medium and tight)

We note that the effect of pile-up is not simulated directly by overlaying minimum bias events to the hard scattering. Although \(\textsc {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 \(\textsc {Delphes}\) FCC-hh setup does not possess the well-calibrated pile-up rejection tools that will necessarily be employed for a detector operating 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 [88, 89]. The level of degradation of the  \(\lambda _3\) measurement precision caused by the deterioration of the performance of specific physics objects (for example the photon energy resolution and reconstruction efficiency or the b-tagging efficiency) has been quantified in previous studies [37]. The impact of degrading the photon energy resolution due to pile-up contamination was studied in full simulation with up 1000 pile-up interactions in Ref. [81]. The degraded resolution was then propagated in \(\textsc {Delphes}\) and the effect on the  \(\lambda _3\) precision was found to be approximately 1% (or 20% in relative terms). We stress however that this level of degradation should be considered as a worse case scenario given that a simple sliding window algorithm was used, and timing information was not exploited. In Ref. [37] we have also studied the impact of degrading the photon reconstruction or, equivalently, of increasing the jet-to-photon probability, which also showed an effect of 1% on the  \(\lambda _3\) precision. Since a full-fledged event simulation and object reconstruction does not exist at this stage for the FCC-hh detector, the assumed object efficiencies result from extrapolations from the LHC detectors. In order to account for a possible degradation of the detector performance in the presence of pile-up, we define 3 baseline scenarios:

  • Scenario I: Optimistic – target detector performance, similar to Run 2 LHC conditions.

  • Scenario II: Realistic – intermediate detector performance.

  • Scenario III: Conservative – pessimistic detector performance, assuming extrapolated HL-LHC performance using present-day algorithms.

Table 3 Summary of the sources of systematic uncertainties in the 3 scenarios. The last column indicates the processes that are affected by the corresponding source of uncertainty. For each given object (b-jet, \(\tau \)-jet, \(\gamma \), lepton), the quoted uncertainty on reconstruction and identification efficiency is applied as many times as the object appears in the final-state

The assumptions on the performance of various physics objects for each baseline scenario, are summarized in Table 2. As mentioned previously, a dominant background for \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) analysis is the \(\gamma +\mathrm {jets}\) process. The probability for a hard scattering jet to be mis-reconstructed as an isolated photon is small, \({\mathcal {O}}(10^{-3})\), in current LHC detectors, thanks to the excellent angular resolution of present calorimeters. As we noted in Sect. 4.1, the assumed granularity for the FCC-hh detector is a factors 2–4 better than present LHC detectors. We make however the conservative choice of assuming a \(\text {j}\rightarrow \gamma \) fake-rate \(\,\epsilon _{j \rightarrow \gamma } = 0.0007 \cdot e^{-p_{\text {T}} [GeV]/187}\), which is of the same magnitude as in the HL-LHC detectors [90]. In addition, we account for the probability for a pile-up jet to be reconstructed as photon by further multiplying the fake rate by a factor 2. This factor has been derived by simulating 1000 PU collision with the CMS Phase-II detector in Delphes, applying a pile-up ID mistag rate of 10 % (from Ref. [89]) and applying the fake-rate probability for calibrated pile-up jets given in Figure 5b in Ref. [90]. This procedure is used for Scenario I. For Scenario II and III we multiply the above fake-rate by factors of 2 and 4 respectively. For leptons we neglect possible fake jets contributions since these are negligible at the momemtum scale relevant for the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) final state. \(\textsc {Delphes}\) also provides heavy flavour tagging, in particular \(\tau \) (hadronic) and b-jet identification. Both hadronic \(\tau \)’s and b-jets are reconstructed using the total visible 4-momentum of the jet. The tagging efficiencies rely on a parameterisation of the (mis-)identification probability as a function of (\(p_{\text {T}}\), \(\eta \)). Again, since we cannot yet derive such performance from full-simulation, we assume efficiencies and mistag rates of the same order as for the (HL-)LHC detectors. For Scenario I, the efficiencies for \(\tau \) and b-jets are modelled after the CMS performance given in Refs. [91, 92]. For scenarios II and III a degradation of the efficiencies for a constant mistag-rate probability is assumed. It should be noted that this is a conservative assumption, since present estimates for heavy-flavour tagging in LHC Phase II conditions project a similar performance as in present conditions, due to superior tracking and high-precision timing detectors [93, 94]. For \(\tau \) and b-jets we also consider two definitions, a “Medium” (M) and a “Tight” (T) working point, in order to operate at an optimal signal-to-background rejection in each decay channel. As shown in Table 2, we also consider the impact of photon and b-jet energy resolutions. Both are relevant since all di-Higgs processes are resonant and the reconstructed Higgs mass directly affects the final sensitivity. The di-photon resolution for scenarios I, II and III was directly determined from full-simulation respectively with 0, 200, and 1000 pile-up interactions (see Ref. [81]). For the di-jet invariant mass, for scenario I we assume the invariant mass resolution as obtained with a multi-variate regression technique in CMS Run 2 (in Ref. [95]), whereas for scenario II and III respectively we assume a factor 1.5 and 2 degradation compared to scenario I. The \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel is a special case, since two di-jet invariant masses are reconstructed (see Sect. 6.3). In that case, as a conservative assumption, we assume that only the Higgs candidate with the largest \(p_{\text {T}}\) is affected by the \(m_{\mathrm {bb}}\) resolution assumed for the scenarii I–III, while the sub-leading Higgs candidate mass has the resolution of Scenario III.

4.3 Systematic uncertainties

Systematic uncertainties can play a major role on the expected sensitivity of the self-coupling 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, following the general principle introduced in Sect. 4.2. 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:

  • Theoretical cross-section affecting the single-Higgs and ZZ backgrounds. Due to their moderate yields we assume these backgrounds to be estimated from Monte Carlo at the FCC-hh. We also assume these two processes to be well known and well reproduced by Monte Carlo simulations at the FCC-hh, with an overall uncertainty varied between 0.5 and 1.5% depending on the scenario. Furthermore, we include a similar theoretical uncertainty on the HH cross section, affecting the interpretation of the HH rate measurement in terms of \(\mu \) and  \(\kappa _{\lambda }\).

  • 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%), reflecting future opportunities to extract the luminosity from hard processes like Z production. As for the theoretical uncertainties, the luminosity affects both the signal and the single-Higgs and ZZ backgrounds.

  • 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. Since we expect this to be one of the dominant uncertainties, it is applied during event simulation accounting for the \(p_{\text {T}}\) dependence of the b-jet efficiency uncertainty (taken from Ref. [96]). This procedure allows to take into account the effect of the uncertainty on the shape of the resulting BDT distribution.

    • \(\tau \)-jets: For each jet originated from the hadronic decay of a \(\tau \)-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 the photon related uncertainties will be comparable to electrons. 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.

The uncertainty on the energy resolution of photons has been studied in Ref. [37]. Degrading the photon resolution by a relative 100% has an effect of order 1% on the self-coupling precision. According to Ref. [97], the uncertainty on the resolution, parametrised as an additional constant term, amounts to 0.4% in the barrel. Assuming a photon energy \(\hbox {E}=60\) GeV, a nominal constant term of 0.8% and a stochastic term of 10%, this results in a relative difference of less than 5% (in relative termsFootnote 4) on the photon energy resolution. Such a degradation has an effect, at first order, of less than 0.1% of the self-coupling measurement (assuming that a degradation of more than 100% has an effect of 1% on the self-coupling precision). A similar argument applies to the effect of the jet energy resolution, which is measured by ATLAS with a relative precision of 7% according to Ref. [98]. We therefore neglect such source of systematic uncertainties. As far as the scale uncertainty goes, for both photons and jets, the above references show that they amount respectively to 0.1% and 1%. We have then explicitly verified that the effect of shifting the scale of both photon and jets by 1% results in a relative change in the significance of order 0.5%, translating at worst into a relative change in the self-coupling precisionFootnote 5 of 1%. We conclude therefore that also the photon and jet (and necessarily hadronic \(\tau \)’s) energy scale have a negligible impact on the self-coupling precision. Moreover, we stress that the effect of a degradation of the detector performance, especially in terms of photons and b-jet energy resolution, is probed by studying the various scenarios given in Table 2. The absolute performance degradation considered in Scenarios II and III largely overcomes the corresponding systematic uncertainties on the object resolution that were neglected.

We assume that several backgrounds will be measured with high statistical accuracy from “side bands” or “control regions”. This is the case for example for the \(\mathrm {t}\bar{\mathrm {t}}\), QCD, and non-single-Higgs backgrounds (with the exception of ZZ, that we assume to be predicted by the Monte Carlo) that dominate the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) and \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channels background contributions. In these cases, while there is no uncertainty associated to the normalisation of the backgrounds, the statistical uncertainty due to the possible fluctuations of the number of events in the side bands is considered in the fit.

Fig. 3
figure 3

Higgs pair invariant-mass distribution in \(\mathrm {g g H H}\) (a) and \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) (b) events for \(\,\kappa _{\lambda } =0\), \(\,\kappa _{\lambda } =1\), \(\,\kappa _{\lambda } =2\) and \(\,\kappa _{\lambda } =3\)

Fig. 4
figure 4

Transverse momentum spectra of the leading (a) and sub-leading (b) Higgs boson in \(\mathrm {g g H H}\) events for \(\,\kappa _{\lambda } =0\), \(\,\kappa _{\lambda } =1\), \(\,\kappa _{\lambda } =2\) and \(\,\kappa _{\lambda } =3\)

When performing the fit for the combination across different channels, systematic uncertainties with the same physical origin are considered fully correlated across processes and final states. Otherwise they are considered as completely uncorrelated, with the notable exception of the b-tagging efficiency (and mistag rate) uncertainty, which is correlated across channels but not across processes. We therefore consider separate shape b-tagging and mistag uncertainties for each process (single H, ZZ, and HH), correlated across the various channels. For example, the b-tagging uncertainty affecting the single-Higgs production is correlated across all channels, but is uncorrelated from the the b-tagging uncertainty affecting double Higgs, ZZ, and so on. This reflects possible differences in the properties of b-tagged jets created in different processes.

5 Signal extraction methodology

As mentioned in Sect. 3.1, the cross section for HH production has a non-trivial dependence on the self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) due to the presence, at LO, of diagrams that contain the trilinear interaction vertex ( \((S)\)) as well as diagrams that do not ( \((T)\)), 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 \(\mathrm {g g H H}\) and \(\mathrm {VBF\,HH}\) modes, the total cross section reaches a minimum respectively at \(\,\kappa _{\lambda } \approx \) 2.5 and \(\,\kappa _{\lambda } \approx \) 1.8, while the slopes of the \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) and \(\mathrm {V H H}\) cross sections carry little dependence on  \(\kappa _{\lambda }\) (see Fig. 2a).

At first order one can write:

$$\begin{aligned} \mu (\,\kappa _{\lambda }) = 1 + (\,\kappa _{\lambda }-1)\,\frac{d\mu }{d\,\kappa _{\lambda }}\Bigr |_{\text {SM}}, \end{aligned}$$
(5.1)

where we define  \(\mu =\sigma / \sigma _{\mathrm {SM}}\) as the signal strength. One can measure  \(\lambda _3\) (or alternatively  \(\kappa _{\lambda }\)) by measuring the total HH production cross section. It follows that:

$$\begin{aligned} \,\delta _{\,\kappa _{\lambda }} = \frac{\,\delta _{\mu }}{\,\frac{d\mu }{d\,\kappa _{\lambda }}\Bigr |_{\text {SM}}}, \end{aligned}$$
(5.2)

where  \(\delta _{\,\kappa _{\lambda }} \) and  \(\delta _{\mu } \) 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 \(\mu \)) at \(\,\kappa _{\lambda } =1\) and by the uncertainty on the measurement of the total cross section. Since  \(\frac{d\mu }{d\,\kappa _{\lambda }}\Bigr |_{\text {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 \(\mu \). Assuming all other standard model parameters are known with better precision than the expected precision on  \(\kappa _{\lambda }\),Footnote 6 the relative weight of the  \((S)\) and  \((T)\) amplitudes (and their interference) is determined by the magnitude of  \(\kappa _{\lambda }\).

The magnitude of  \(\kappa _{\lambda }\) 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_\mathrm {hh}\) is highly sensitive to the value of the self-coupling. This can be easily understood by noting that configurations with large \(m_\mathrm {hh}\) are mostly suppressed in the  \((S)\) amplitude (not in  \((T)\) diagrams). Vice-versa, the phase-space region near threshold, at \(m_\mathrm {hh}\) \(\gtrsim \) 2\(m_\mathrm {H}\), maximises the  \((S)\) contribution. The \(m_\mathrm {hh}\) distribution is shown for the \(\mathrm {g g H H}\) and \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) processes respectively in Fig. 3a, b. For \(\mathrm {g g H H}\) the dependence is distorted at values of  \(\kappa _{\lambda }\) \(\approx \) 2 due to the large destructive interference between  \((S)\) and  \((T)\). The transverse momenta of the two Higgs bosons (\(p_{\text {T}} (h_1) \), and \(p_{\text {T}} (h_2) \)) also display a large dependence on  \(\kappa _{\lambda }\), as shown in Fig. 4a, 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  \(\kappa _{\lambda }\). 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  \(\mathrm {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_\mathrm {hh}\) observable, and the Higgs particles transverse momentum (\(p_{\text {T}} (h_1) \), and \(p_{\text {T}} (h_2) \)) differential distributions to further improve the sensitivity on  \(\kappa _{\lambda }\).

6 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: \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\). 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 [101, 102]. The statistical procedure and the evaluation of the systematic uncertainties are summarized in Appendix A and 4.3, respectively.

For a similar analysis in the case of HL-LHC, see Ref. [103] and the studies by the ATLAS and CMS collaborations [104, 105], contributed to Ref. [3].

6.1 The \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel

Despite its small branching fraction, the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel is by far the most sensitive decay mode for measuring the Higgs self-coupling. The presence of two high \(p_{\text {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  \(\mathrm {S/B}\). The largest background processes are single-Higgs production and the QCD continuum \(\gamma \gamma +\mathrm {jets}\) and \(\gamma +\mathrm {jets}\). A discussion of the simulation of these processes was given in Sect. 3.2.

6.1.1 Event selection

In the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel, events are required to contain at least two isolated photons and two b-tagged jets with the requirement \(p_{\text {T}} (\gamma ,\mathrm {b})\) > 30 GeV and \(|\eta (\gamma ,\mathrm {b})\)| < 4.0. The leading photon and b-jet are further required to have \(p_{\text {T}} (\gamma ,\mathrm {b})\) > 35 GeV. The Higgs candidates 4-momenta are formed respectively from the two reconstructed b-jets and photons with the largest \(p_{\text {T}} (\gamma ,\mathrm {b})\). The b-jets are identified with the “Medium” working point criterion, defined in Table 2. Since the \(\gamma \gamma +\mathrm {jets}\) process was generated with a parton-level requirement (see Sect. 3.2) on \(m_{\gamma \gamma }\), we further require the events to pass the loose selection \(|m_{\gamma \gamma }- 125| < 7\) GeV. The efficiency of the full event selection for the SM signal sample is approximately 26%. For an integrated luminosity  \(\mathcal {L}_{int}=30\text { ab}^{-1}\)   this event selection yields approximately 10k Higgs pair events, 125k single Higgs, 2.6M \(\mathrm {j j} \gamma \gamma \) and 7M \(\gamma +\mathrm {jets}\) events for Scenario I. The trigger efficiency for the above selection is assumed to be 100% efficient.

Fig. 5
figure 5

Invariant mass spectra of the \(\mathrm {H} \rightarrow \gamma \gamma \)  (a), \(\mathrm {H}\rightarrow \mathrm {b}\bar{\mathrm {b}}\) (b), HH (c) candidates after applying the event pre-selection. The SM Higgs pair process is normalized to 50 times the expected yield with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\)

Fig. 6
figure 6

Spectrum of SM signal (a), the QCD (b) and single Higgs (c) backgrounds in the ( \(\mathrm {BDT}_\mathrm {H}\)\(\mathrm {BDT}_{\mathrm {QCD}}\)) plane

Fig. 7
figure 7

Expected negative log-Likelihood scan as a function the signal strenth  \(\mu =\sigma / \sigma _{\mathrm {SM}}\) (a) and trilinear self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) (b) in the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel. The various lines correspond to the different systematic uncertainty assumptions summarized in Table 3. The black dashed line shows the likelihood profile when only the statistical uncertainty is included under scenario I

In order to maximally exploit the kinematic differences between signal and background, a boosted decision tree (BDT) is trained using most of the available kinematic information in the event:

  • The 3-vector components of the leading ( \(\gamma _{1}\)) and subleading photon ( \(\gamma _{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\gamma _{1}}}\), \(p_{\text {T}}^{\mathrm {\,\gamma _{2}}}\)), pseudo-rapidity (\(\eta _{\,\gamma _{1}}\), \(\eta _{\,\gamma _{2}}\)), and azimuthal angle (\(\phi _{\,\gamma _{1}}\), \(\phi _{\,\gamma _{2}}\)).

  • The 3-vector components of the leading ( \(\mathrm {b}_{1}\)) and subleading b-jet ( \(\mathrm {b}_{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{2}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {b}_{1}}\), \(\eta _{\,\mathrm {b}_{2}}\)), and azimuthal angle (\(\phi _{\,\mathrm {b}_{1}}\), \(\phi _{\,\mathrm {b}_{2}}\)).

  • The 3-vector components of the leading ( \(\mathrm {j}_{1}\)) and subleading additional reconstructed jets in the event ( \(\mathrm {j}_{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{2}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {j}_{1}}\), \(\eta _{\,\mathrm {j}_{2}}\)), and azimuthal angle (\(\phi _{\,\mathrm {j}_{1}}\), \(\phi _{\,\mathrm {j}_{2}}\)). If no additional jets are found, dummy values are given to these variables.

  • The 4-vector components of the \(\mathrm {H} \rightarrow \gamma \gamma \) candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\gamma \gamma }}\)), pseudo-rapidity (\(\eta _{\,\gamma \gamma }\)), azimuthal angle (\(\phi _{\,\gamma \gamma }\)) and invariant mass (\(m_{\gamma \gamma }\)).

  • The 4-vector components of the \(\mathrm {H}\rightarrow \mathrm {b}\bar{\mathrm {b}}\) candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {bb}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {bb}}\)), azimuthal angle (\(\phi _{\,\mathrm {bb}}\)) and invariant mass (\(m_{\mathrm {bb}}\)).

  • The 4-vector components of the Higgs pair candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {hh}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {hh}}\)), azimuthal angle (\(\phi _{\,\mathrm {hh}}\)) and invariant mass (\(m_\mathrm {hh}\)).

  • The number of reconstructed b-jets \(N_b\), the number of light jets \(N_l\) and the total number of jets \(N_j = N_b + N_l\).

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 \(\textsc {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_{\gamma \gamma }\), \(m_{\mathrm {bb}}\), and \(m_\mathrm {hh}\) observables, shown respectively in Fig. 5a–c provide most of the discrimination against the background.

The QCD (\(\gamma +\mathrm {jets}\) and \(\gamma \gamma +\mathrm {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 photon-pair candidates in single-Higgs processes often originate from a Higgs decay. As a result, while the \(m_{\gamma \gamma }\) 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, producing two multivariate discriminants:  \(\mathrm {BDT}_\mathrm {H}\) and  \(\mathrm {BDT}_{\mathrm {QCD}}\). During the training, each background within each class is weighted according to the relative cross section.

The discrimination against single-Higgs backgrounds ( \(\mathrm {BDT}_\mathrm {H}\)) is largely driven by the \(m_{\mathrm {bb}}\) variable, followed by \(m_\mathrm {hh}\) and \(p_{\text {T} \,\mathrm {bb}}\). Additional separation power, in particular against the \(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\) background, is provided by the number of reconstructed jets \(N_j\) as well as \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{1}}}\) and \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{2}}}\). On the other hand, the discrimination against QCD is driven by the \(m_{\gamma \gamma }\) and \(m_{\mathrm {bb}}\) variables, followed by \(p_{\text {T} \,\gamma \gamma }\), \(p_{\text {T} \,\mathrm {bb}}\), \(p_{\text {T}}^{\mathrm {\,\gamma _{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{1}}}\), \(m_\mathrm {hh}\), \(p_{\text {T}}^{\mathrm {\,\gamma _{2}}}\) and \(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{2}}}\).

The output of the BDT discriminant is shown in the  \((\,\mathrm {BDT}_\mathrm {H},\,\,\mathrm {BDT}_{\mathrm {QCD}})\) plane for the signal and the two background components in Fig. 6a–c, respectively. As expected, the signal (background) enriched region clearly corresponds to large (small)  \(\mathrm {BDT}_\mathrm {H}\) and  \(\mathrm {BDT}_{\mathrm {QCD}}\) values. We note that the multivariate discriminant correctly identifies the two main components (\(\mathrm {ggH}\) and \(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\)) within the single-Higgs background. The \(\mathrm {ggH}\) background, as opposed to \(\mathrm {t}\bar{\mathrm {t}}\mathrm {H}\), is more “signal-like” and populates a region of high  \(\mathrm {BDT}_\mathrm {H}\) and  \(\mathrm {BDT}_{\mathrm {QCD}}\).

6.1.2 Signal extraction and results

The expected precision on the signal strength  \(\mu =\sigma / \sigma _{\mathrm {SM}}\) and on the self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) are obtained from a 2-dimensional fit of the  \((\,\mathrm {BDT}_\mathrm {H},\,\,\mathrm {BDT}_{\mathrm {QCD}})\) output, following the procedure described in Appendix A. The results are shown in Fig. 7a, b. The various lines correspond to the different scenarios described in Sects. 4.2 and 4.3. From Fig. 7b one can extract the symmetrized 68% and 95% confidence intervals for the various scenarios. The expected precision for \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) is summarized in Table 4 for each of these assumptions. Depending on the assumed scenario, the Higgs self-coupling can be measured with a precision of 3.8–10% at 68% C.L using the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel alone. We note that the achievable precision is largely dependent on the assumptions on the detector configuration and the systematic uncertainties. Such result needs to be compared against the precision of  \(\delta _{\,\kappa _{\lambda }} \) \(=3.4{-}7.4\%\), obtained using statistical uncertainties alone. It is clear that the performance of the detector. In particular degrading the photon and jet energy resolution can have a substantial impact on the achievable precision.

6.2 The \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel

The \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel is very attractive thanks to the large branching fraction (7.3%) and the relatively clean final state. As opposed to the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel, the \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) decay cannot be fully reconstructed due to the presence of \(\tau \) neutrinos in the final state. We consider mainly two channels here: the fully hadronic final state \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\), and the semi-leptonic one, \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \)(\(\ell =e,\mu \)).

Table 4 Expected precision at 68% CL on the di-Higgs production signal strength and Higgs self coupling using the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel at the FCC-hh with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\). The symmetrized value \(\delta =(\delta ^++\delta ^-)/2\) is given in %
Table 5 Expected precision at 68% CL on the di-Higgs production cross-section and Higgs self coupling using the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel at the FCC-hh with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\). The symmetrized value \(\delta =(\delta ^++\delta ^-)/2\) is given in %

As spelled out in Sect. 3.2, several processes act as background for the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) final state. The largest background contributions are QCD and \(\mathrm {t}\bar{\mathrm {t}}\). QCD is a background mainly for the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _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. Moreover, analyses using CMS data [106] show that QCD is overall a subdominant background at the LHC, and negligible in the signal region. As a result recent CMS Phase II projections neglect this background altogether [105]. In order of decreasing magnitude, the largest backgrounds are \(\text {Z}/\gamma ^*\text {+jets}\)   single Higgs, ttV and ttVV, where \(\hbox {V}=\hbox {W},\hbox {Z}\).

6.2.1 Event selection

Events are required to contain at least two b-jets with \(p_{\text {T}} (\mathrm {b})\) > 30 GeV and \(|\eta (\mathrm {b})\)| < 3.0. We require at least, and not exactly, two bjets in order not to suppress the \(\mathrm {t}\bar{\mathrm {t}}\mathrm {HH}\) signal contribution. For the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) final state the presence is required of exactly one isolated (\(\,I_\text {rel} ~<~0.1\)) lepton \(\ell =e,\mu \) with \(p_{\text {T}} (\ell )\) > 25 GeV and \(|\eta (\ell )\)| < 3.0 and at least one hadronically tagged \(\tau \)-jet with \(p_{\text {T}} (\tau _h)\) > 45 GeV and \(|\eta (\tau _h)\)| < 3.0. For the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) final state, we require at least two hadronically tagged \(\tau \)-jet with \(p_{\text {T}} (\tau _h)\) > 45 GeV and \(|\eta (\tau _h)\)| < 3.0. The hadronic \(\tau \) is identified according to the “Tight” criterion defined in Table 2. This ensures a highly efficient rejection of the QCD background (at the cost of a smaller \(\tau _h\) efficiency) which strengthens the solidity of our assumption of neglecting the QCD background altogether. In what follows we refer to a \(\tau \)-candidate as the lepton \(\ell =e,\mu \) or the \(\tau \)-jet. In particular the \(\tau \) 4-momentum is defined as the sum of the 4-momenta of the visible \(\tau \) decay products. In order to maximally exploit the kinematic differences between the signal and the dominant \(\mathrm {t}\bar{\mathrm {t}}\) background, we build a multivariate BDT discriminant using as an input the following kinematic properties:

  • The 3-vector components of the leading ( \(\tau _{1}\)) and subleading \(\tau \)-candidate ( \(\tau _{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\tau _{1}}}\), \(p_{\text {T}}^{\mathrm {\,\tau _{2}}}\)), pseudo-rapidity (\(\eta _{\,\tau _{1}}\), \(\eta _{\,\tau _{2}}\)), and azimuthal angle (\(\phi _{\,\tau _{1}}\), \(\phi _{\,\tau _{2}}\)).

  • The 3-vector components of the leading ( \(\mathrm {b}_{1}\)) and subleading b-jet ( \(\mathrm {b}_{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {b}_{2}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {b}_{1}}\), \(\eta _{\,\mathrm {b}_{2}}\)), and azimuthal angle (\(\phi _{\,\mathrm {b}_{1}}\), \(\phi _{\,\mathrm {b}_{2}}\)).

  • The 3-vector components of the leading ( \(\mathrm {j}_{1}\)) and subleading additional reconstructed jets in the event ( \(\mathrm {j}_{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{2}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {j}_{1}}\), \(\eta _{\,\mathrm {j}_{2}}\)), and azimuthal angle (\(\phi _{\,\mathrm {j}_{1}}\), \(\phi _{\,\mathrm {j}_{2}}\)). If no additional jets are found, dummy values are given to these variables.

  • The 4-vector components of the \(\mathrm {H} \rightarrow \tau \tau \) candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\tau \tau }}\)), pseudo-rapidity (\(\eta _{\,\tau \tau }\)), azimuthal angle (\(\phi _{\,\tau \tau }\)) and invariant mass (\(m_{\tau \tau } \)).

  • The 4-vector components of the \(\mathrm {H}\rightarrow \mathrm {b}\bar{\mathrm {b}}\) candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {bb}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {bb}}\)), azimuthal angle (\(\phi _{\,\mathrm {bb}}\)) and invariant mass (\(m_{\mathrm {bb}}\)).

  • The 4-vector components of the Higgs pair candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {hh}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {hh}}\)), azimuthal angle (\(\phi _{\,\mathrm {hh}}\)) and invariant mass (\(m_\mathrm {hh}\)).

  • The transverse missing energy \(p_{\text {T}} ^{\text {miss}} \).

  • The transverse mass of each \(\tau \)-candidate, computed as \(\,m_\text {T}= \sqrt{2p_{\text {T}}^{\mathrm {\tau }} p_{\text {T}} ^{\text {miss}}- \vec {p_{\text {T}}}^{\tau } \cdot \vec {p_{\text {T}}}^{\text {miss}}}\).

  • The event s-transverse mass  \(m_\text {T2}\) as defined in Refs.  [107, 108].

  • The number of reconstructed b-jets \(N_b\), the number of light jets \(N_l\) and the total number of jets \(N_j = N_b + N_l\).

The \(m_{\tau \tau } \) and  \(m_\text {T2}\) observables are shown in Figs. 8a, 9a and 8b, 9b for the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) and \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) final states respectively. These provide, together with the \(m_{\mathrm {bb}}\) observable, the largest discrimination against the \(\mathrm {t}\bar{\mathrm {t}}\) background. Additional discrimination against the \(\mathrm {t}\bar{\mathrm {t}}\) background is provided by \(p_{\text {T}}^{\mathrm {\,\tau \tau }}\) and \(p_{\text {T}}^{\mathrm {\,\mathrm {bb}}}\) and the \(p_{\text {T}} ^{\text {miss}} \)  followed by \(m_\mathrm {hh}\), \(\eta _{\,\tau \tau }\), \(\eta _{\,\mathrm {bb}}\) and the \(p_{\text {T}}\) and \(\eta \) of  \(\tau _{2}\),  \(\tau _{1}\),  \(\mathrm {b}_{2}\),  \(\mathrm {b}_{1}\),  \(\mathrm {j}_{2}\) and  \(\mathrm {j}_{1}\), and finally \(N_j\). The output of the BDT discriminant is shown in Figs. 8c and 9c for the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) and \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) final states.

Fig. 8
figure 8

Distributions in the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) final state of the invariant mass of the  \(\tau \tau \) pair (left),  \(m_\text {T2}\) (center), and output of the BDT multi-variate discriminant (right)

Fig. 9
figure 9

Distributions in the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) final state of the invariant mass of the  \(\tau \tau \) pair (left),  \(m_\text {T2}\) (center), and output of the BDT multi-variate discriminant (right)

6.2.2 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 \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) and \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) 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 \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel is shown in Fig. 10a, b. The various lines correspond to the different scenarios described in Sects. 4.2 and 4.3. From Fig. 10b one can extract the 68% and 95% confidence intervals for the various systematics assumptions. Depending on the assumed scenario, using the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel, the Higgs pair signal strength and Higgs self-coupling can be measured respectively with a precision of \(\,\delta _{\mu } =5.8{-}7.9\%\) and \(\,\delta _{\,\kappa _{\lambda }} =9.8{-}13.8\%\) at 68% C.L.Footnote 7 Despite the large signal event rate in the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel, the sensitivity is limited by the large background. Therefore, the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel is statistically dominated at the FCC-hh with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\), and the achievable precision is only moderately dependent on the assumptions on the systematic uncertainties. We note that Ref. [31] quotes a precision of \(\,\delta _{\,\kappa _{\lambda }} =13\%\) using the resolved semi-leptonic final state, which is consistent with the result presented here of \(\,\delta _{\,\kappa _{\lambda }} =14{-}18\%\). The same study shows that further precision can be obtained using the boosted topology, which has not been considered here.

Fig. 10
figure 10

Expected negative log-likelihood scan as a function the signal strength  \(\mu =\sigma / \sigma _{\mathrm {SM}}\) (a) and trilinear self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) (b) in the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) channel (combination of the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) and \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) channels). The various lines correspond to the different systematic uncertainties assumptions summarized in Table 3. The black dashed line shows the likelihood profile when only the statistical uncertainty is included under scenario I

6.3 The \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel

The \(\mathrm {HH} \rightarrow {\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) 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 very large QCD backgrounds and hence features a relatively small  \(\mathrm {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 [109, 110], provides less sensitivity to the self-coupling measurement and was discussed in previous studies [31, 37]. The main backgrounds to this final state are QCD and \(\mathrm {t}\bar{\mathrm {t}}\), followed by \(\mathrm {Z b}\bar{\mathrm {b}}\), single-Higgs production and ZZ.

6.3.1 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_{\text {T}} (\mathrm {b})\) > 30 GeV and \(|\eta (\mathrm {b})\)| < 4.0. The b-jets are identified with the “Medium” working point defined in Table 2. 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_{\text {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 3-vector components of the four leading b-jets in the event (\(b_1\), \(b_2\), \(b_3\), \(b_4\)): transverse momenta (\(p_{\text {T}}^{\mathrm {b_i}}\)), pseudo-rapidities (\(\eta _{b_i}\)), and azimuthal angles (\(\phi _{b_i}\), \(\text {i}=1..4\)).

  • The 3-vector components of the leading ( \(\mathrm {j}_{1}\)) and subleading additional reconstructed jets in the event ( \(\mathrm {j}_{2}\)): transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{2}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {j}_{1}}\), \(\eta _{\,\mathrm {j}_{2}}\)), and azimuthal angle (\(\phi _{\,\mathrm {j}_{1}}\), \(\phi _{\,\mathrm {j}_{2}}\)). If no additional jets are found, dummy values are given to these variables.

  • The 4-vector components of the leading (\(h_1\)) and subleading (\(h_2\)) Higgs candidates: transverse momentum (\(p_{\text {T}}^{\mathrm {h_1}}\), \(p_{\text {T}}^{\mathrm {h_2}}\)), pseudo-rapidity (\(\eta _{h_1}\), \(\eta _{h_2}\)), azimuthal angle (\(\phi _{h_1}\), \(\phi _{h_2}\)) and invariant mass (\(m_{h_1}\), \(m_{h_2}\)).

  • The 4-vector components of the Higgs-pair candidate: transverse momentum (\(p_{\text {T}}^{\mathrm {\,\mathrm {hh}}}\)), pseudo-rapidity (\(\eta _{\,\mathrm {hh}}\)), azimuthal angle (\(\phi _{\,\mathrm {hh}}\)) and invariant mass (\(m_\mathrm {hh}\)).

The \(m_{h_1}\) and \(m_\mathrm {hh}\) observables are shown respectively in Fig. 11a, 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 \(\mathrm {Z b}\bar{\mathrm {b}}\) 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. A substantial discrimination against the QCD background is provided, in decreasing order of importance, by \(p_{\text {T}}^{\mathrm {b_2}}\), \(p_{\text {T}}^{\mathrm {b_4}}\), \(p_{\text {T}}^{\mathrm {b_3}}\) and \(p_{\text {T}}^{\mathrm {b_1}}\), followed by \(m_\mathrm {hh}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {hh}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{1}}}\), \(p_{\text {T}}^{\mathrm {\,\mathrm {j}_{2}}}\), and the \(\eta _{b_i}\). The output of the BDT discriminant is shown in Fig. 11c for the signal and various background contributions.

Fig. 11
figure 11

Distributions in the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) final state of the highest \(p_{\text {T}}\) reconstructed Higgs invariant mass (left), the Higgs pair invariant mass (center), and the output of the BDT multi-variate discriminant (right)

Fig. 12
figure 12

Expected Negative log-Likelihood scan as a function the signal strength  \(\mu =\sigma / \sigma _{\mathrm {SM}}\) (a) and trilear self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) (b) in the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel. The various lines correspond to the different systematic uncertainties assumptions summarized in Table 3. The black dashed line shows the likelihood profile when onlythe statistical uncertainty is included under scenario I

Fig. 13
figure 13

Expected negative log-Likelihood scan as a function of the trilinear self-coupling modifier  \(\,\kappa _{\lambda } =\,\lambda _3/ \,\lambda _3 ^{\mathrm {SM}}\) in all channels, and their combination. The solid line corresponds to the scenario II for systematic uncertainties. The band boundaries represent respectively scenarios I and III. The dashed line represents the sensitivity obtained including statistical uncertainties only, under the assumptions of scenario I

Table 6 Expected precision at 68% CL on the di-Higgs production cross-section and Higgs self coupling using the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel at the FCC-hh with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\). The symmetrized value \(\delta =(\delta ^++\delta ^-)/2\) is given in %

6.3.2 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 \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel is shown in Fig. 12a, b. The coloured lines correspond to the different systematic uncertainties assumptions summarized in Table 3. The 68% and 95% confidence intervals on  \(\delta _{\mu } \) and  \(\delta _{\,\kappa _{\lambda }} \) for the various systematics assumptions can be extracted from Fig. 12a, b and the results are summarized in Table 6. Depending on the assumed scenario, using the \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channel, the Higgs pair signal strenth and Higgs self-coupling can be measured respectively with a precision of \(\,\delta _{\mu } =8{-}18\%\) and \(\,\delta _{\,\kappa _{\lambda }} =18{-}32\%\) at 68% C.L. As for the \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) case, due to the huge QCD background this channel is statistically limited. We note that, despite the large statistical uncertainty, the systematic uncertainties play a large role in this channel. This is the result of the significant contamination in the signal region from single-Higgs background events. Since this background is assumed to be estimated from Monte Carlo, its uncertainty, even though at the percent level, is reflected in a larger signal systematics.

6.4 Combined precision

When combining results from the various channels, the systematic uncertainties from the various sources affecting those processes that we assume to be estimated from Monte Carlo simulations (HH, single Higgs, and ZZ) are accounted for as follows.

Fig. 14
figure 14

Expected precision on the Higgs self-coupling as a function of the integrated luminosity

Lepton (e/\(\mu \), \(\tau \)) uncertainties are correlated across all process and across the \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _h\) and \({\mathrm {b}\bar{\mathrm {b}}} \tau _h \tau _\ell \) channels. In the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel, the photon uncertainty for the single and double Higgs processes are correlated. The luminosity uncertainty is correlated across all these processes and all channels. The same applies, for each process independently, to the overall normalisation uncertainty. The uncertainties on lepton ID, luminosity and normalisation are assumed to affect only the overall normalization of signal and background shapes and to not introduce a significant deformation of their shapes. For the b-jets ID, we take into account both the shape and normalisation uncertainty due to the b-jets systematic uncertainties. This is achieved by shifting the efficiency for each jet by the (\(p_{\text {T}}\)-dependent) values reported in Table 2, and re-computing the resulting BDT distributions. As discussed in Sect. 4.3, this uncertainty is correlated across different channels when it affects the same process. Overall normalisation uncertainties are cancelled out when a background is estimated from a control region. Moreover, we expect all control regions to be well-populated at the FCC-hh. For these reasons, we assume the systematic uncertainties affecting those backgrounds that are most likely to be estimated from well populated control regions or side bands, such as QCD, Zbb, photon(s) + jets, and \(\mathrm {t}\bar{\mathrm {t}}\), to be negligible and we do not include them in the fit procedure.

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 \({\mathrm {b}\bar{\mathrm {b}}ZZ~(4}\ell \mathrm {)}\) 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.  [37]. The expected combined precision on the Higgs self-coupling obtained after combining the channels \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \), \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) and \({\mathrm {b}\bar{\mathrm {b}}ZZ~(4}\ell \mathrm {)}\) can be inferred from the intersection of black curves with the horizontal 68% and 95% CL lines. The expected statistical precision for Scenario I, neglecting systematic uncertainties, can be read from the dashed black line in Fig. 13, and gives \(\,\delta _{\,\kappa _{\lambda }} = 3.0\)% at 68% CL. The solid line corresponds to scenario II, 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 \(\,\delta _{\,\kappa _{\lambda }} = 3.4{-}7.8\)% at 68% CL. The signal strength and self-coupling precision for the combination are summarized in Table 7.

Table 7 Combined expected precision at 68% CL on the di-Higgs production cross- and Higgs self coupling using all channels at the FCC-hh with  \(\mathcal {L}_{int}=30\text { ab}^{-1}\). The symmetrized value \(\delta =(\delta ^++\delta ^-)/2\) is given in %

The expected precision on the Higgs self-coupling as a function of the integrated luminosity is shown in Fig. 14, for the three scenarios, with and without systematic uncertainties. With the most aggressive scenario I, a precision of \(\,\delta _{\,\kappa _{\lambda }} =10\%\) can be reached with only 3 \(\hbox {ab}^{-1}\) of integrated luminosity, whereas approximately 20 \(\hbox {ab}^{-1}\) are required for the most conservative scenario III. Therefore, assuming scenario I, the 10% target should therefore be achievable during the first 5 years of FCC-hh 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.

Fig. 15
figure 15

Expected precision on the Higgs self-coupling as a function of the  \(\kappa _{\lambda }\) value for each scenario. To improve readability, the positions of scenario I (III) bands are slightly offset in the negative (positive) direction along the  \(\kappa _{\lambda }\) axis

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  \(\kappa _{\lambda }\) in the range \(0<\,\kappa _{\lambda } <3\), in order to maximise the achievable precision for values of \(\,\kappa _{\lambda } \ne 1\). The BDT training has been performed only for the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel, which dominates the overall sensitivity, whereas for the other channels we conservatively employ the BDT trained at \(\,\kappa _{\lambda } =1\). The obtained precision as a function of \(\,\kappa _{\lambda } \) is shown in Fig. 15.Footnote 8

It can be seen that the overall precision follows the behaviour of the HH production cross section as function of  \(\kappa _{\lambda }\) given in Fig. 2a. The best precision,  \(\delta _{\,\kappa _{\lambda }} \) \(\approx \) 2%, is reached at \(\,\kappa _{\lambda } =0\) where the value  \(\frac{d\mu }{d\,\kappa _{\lambda }} \) is large. Conversely, the maximum uncertainty  \(\delta _{\,\kappa _{\lambda }} \) \(\approx \) 60% is obtained at \(\,\kappa _{\lambda } \approx 2.5\), and corresponds to the minimum of the total HH cross section, where  \(\frac{d\mu }{d\,\kappa _{\lambda }} \) \(\rightarrow 0\) . As can be expected, the likelihood function presents a broad second minimumFootnote 9 in correspondence of the minimum of the HH cross section at \(\,\kappa _{\lambda } =2.5\). The presence of this minimum is the reason behind the asymmetric behaviour of the uncertainties for the points near \(\,\kappa _{\lambda } =2.5\). If the measurement is performed close enough to \(\,\kappa _{\lambda } =2.5\) the likelihood falls in the second minimum before reaching the 68% C.L. threshold, thus enlarging the measurement uncertainty in one direction. It should be noted that, while the HH cross section is roughly symmetric around \(\,\kappa _{\lambda } =2.5\), we do not expected the uncertainties to be symmetric as well, as the kinematic behaviour of the HH system are quite different between \(\,\kappa _{\lambda } <2.5\) and \(\,\kappa _{\lambda } >2.5\). It can also be noticed that when switching on the systematic uncertainties the precision at small  \(\kappa _{\lambda }\) degrades compared to the SM case. This reflects the fact that the HH kinematics at \(\,\kappa _{\lambda } \approx 0\) are similar to the single-Higgs background.

7 Conclusions and perspectives

The precise measurement of the Higgs self-coupling must be a top priority of future high-energy 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.Footnote 10 In the present study the measurement strategy has been optimized in the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \), \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channels using machine learning techniques. For the first time, a precise set of assumptions of detector performances and possible sources of systematic uncertainties has been defined and used to derive the achievable precision. Consistently with our previous findings, the \({\mathrm {b}\bar{\mathrm {b}}} \gamma \gamma \) channel drives the final sensitivity, with an expected precision of \(\,\delta _{\,\kappa _{\lambda }} = 3.8{-}10.0\%\) depending on the detector and systematic assumptions. The \({\mathrm {b}\bar{\mathrm {b}}} \,\tau \tau \) and \({\mathrm {b}\bar{\mathrm {b}}\mathrm {b}\bar{\mathrm {b}}}\) channels provide instead a less precise single channel measurement, respectively of \(\,\delta _{\,\kappa _{\lambda }} = 10{-}14\%\) and \(\,\delta _{\,\kappa _{\lambda }} = 22{-}32\%\).

The final combined sensitivity across all considered channels leads to an expected precision at the FCC-hh on the Higgs self-coupling \(\,\delta _{\,\kappa _{\lambda }} = 3.4{-}7.8\%\) with an integrated luminosity of  \(\mathcal {L}_{int}=30\text { ab}^{-1}\). By considering the most agressive detector and systematics assumptions, the 10% threshold can be achieved with \(\sim 3\) \(\hbox {ab}^{-1}\), corresponding to \(\sim 3{-}5\) years of early running at the start-up luminosity of \(5\times 10^{34}\) \(\text{ cm}^{-2}~\text{ 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 Sect. 4.3, 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 [68], and possibly beyond the \(\hbox {N}^3\)LO in the \(m_{top}\rightarrow \infty \) limit, recently achieved [111, 112]. 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 \(\hbox {ab}^{-1}\), can measure the Higgs self-coupling more precisely than any other proposed project, on a competitive time scale.