Boosting Higgs pair production in the [Formula: see text] final state with multivariate techniques.

The measurement of Higgs pair production will be a cornerstone of the LHC program in the coming years. Double Higgs production provides a crucial window upon the mechanism of electroweak symmetry breaking and has a unique sensitivity to the Higgs trilinear coupling. We study the feasibility of a measurement of Higgs pair production in the [Formula: see text] final state at the LHC. Our analysis is based on a combination of traditional cut-based methods with state-of-the-art multivariate techniques. We account for all relevant backgrounds, including the contributions from light and charm jet mis-identification, which are ultimately comparable in size to the irreducible 4b QCD background. We demonstrate the robustness of our analysis strategy in a high pileup environment. For an integrated luminosity of [Formula: see text] ab[Formula: see text], a signal significance of [Formula: see text] is obtained, indicating that the [Formula: see text] final state alone could allow for the observation of double Higgs production at the High Luminosity LHC.

Introduction

The measurement of double Higgs production will be one of the central physics goals of the LHC program in its recently started high-energy phase, as well as for its future high-luminosity upgrade (HL-LHC) which aims to accumulate a total integrated luminosity of 3 ab [1, 2]. Higgs pair production [3] is directly sensitive to the Higgs trilinear coupling and provides crucial information on the electroweak symmetry breaking mechanism. It also probes the underlying strength of the Higgs interactions at high energies, and it can be used to test the composite nature of the Higgs boson [4, 5]. While Standard Model (SM) cross sections are small, many Beyond the SM (BSM) scenarios predict enhanced rates for double Higgs production; therefore searches have already been performed by ATLAS and CMS with Run I data [6-10] and will continue at Run II. The study of Higgs pair production will also be relevant to any future high-energy collider, either at a 100 TeV circular machine [11-14] or at a linear or circular electron–positron collider [15].

Analogously to single Higgs production [16], in the SM the dominant mechanism for the production of a pair of Higgs bosons at the LHC is gluon fusion (see [3, 17] and references therein). For a center-of-mass energy of  TeV, the next-to-next-to-leading order (NNLO) total cross section is approximately 40 fb [18], which is increased by a further few percent once next-to-next-to-leading logarithmic (NNLL) corrections are accounted for [19]. Feasibility studies in the case of a SM-like Higgs boson in the gluon-fusion channel at the LHC have been performed for different final states, including [20-22], [23-26], [25, 27] and [21, 23, 25, 28, 29]. While these studies differ in their quantitative conclusions, the consistent picture emerges that the ultimate precision in the determination of the Higgs trilinear coupling requires the full integrated luminosity of the HL-LHC, ab and should rely on the combination of different final states. The interplay between kinematic distributions for the extraction of from the measured cross sections and the role of the associated theoretical uncertainties have been intensely scrutinised recently [17, 30–37].

In addition to the gluon-fusion channel, Higgs pairs can also be produced in the vector-boson fusion channel hhjj [5, 26, 38, 39], the associated production modes hhW and hhZ [3, 40, 41] (also known as Higgs-Strahlung), and also in association with top quark pairs [42]. All these channels are challenging due to the small production rates: at 14 TeV, the inclusive total cross sections are 2.0 fb for VBF hhjj [43], 0.5 fb for W(Z)hh [3] and 1.0 for [42].

While the SM production rates for Higgs pairs are small, they are substantially enhanced in a variety of BSM scenarios. Feasibility studies of Higgs pair production in New Physics models have been performed in a number of different frameworks, including Effective Field Theories (EFTs) with higher-dimensional operators and anomalous Higgs couplings [14, 44–50], resonant production in models such as extra dimensions [51-54], and Supersymmetry and Two Higgs Doublet models (2HDMs) [55-61]. Since BSM dynamics modify the kinematic distributions of the Higgs decay products, for instance boosting the di-Higgs system, different analysis strategies might be required for BSM Higgs pair searches as compared to SM measurements.

Searches for the production of Higgs pairs have already been performed with 8 TeV Run I data by ATLAS in the [7] and [8] final states, and by CMS in the same [9] and [10] final states. In addition, ATLAS has presented [6] a combination of its di-Higgs searches in the , and bbbb final states. Many other exotic searches involve Higgs pairs in the final state, such as the recent search for heavy Higgs bosons H [62].

In the context of SM production, the main advantage of the final state is the enhancement of the signal yield from the large branching fraction of Higgs bosons into pairs, [16]. However, a measurement in this channel needs to deal with an overwhelming QCD multi-jet background. Recent studies of Higgs pair production in this final state [28, 29] estimate that, for an integrated luminosity of  ab, a signal significance of around can be obtained. In these analysis, irreducible backgrounds such as 4b and are included, however, the reducible components, in particular bbjj and jjjj, are neglected. These can contribute to the signal yield when light and charm jets are mis-identified as b-jets. Indeed, due to both selection effects and b-quark radiation in the parton shower, the contribution of the 2b2j process is as significant as the irreducible 4b component.

In this work, we revisit the feasibility of SM Higgs pair production by gluon fusion in the final state at the LHC. Our strategy is based upon a combination of traditional cut-based methods and multivariate analysis (MVA). We account for all relevant backgrounds, including the contribution from mis-identified light and charm jets. We also assess the robustness of our analysis strategy in an environment with high pileup (PU). Our results indicate that the final state alone should allow for the observation of double Higgs production at the HL-LHC.

The structure of this paper proceeds as follows. In Sect. 2 we present the modeling of the signal and background processes with Monte Carlo event generators. In Sect. 3 we introduce our analysis strategy, in particular the classification of individual events into different categories according to their topology. Results of the cut-based analysis are then presented in Sect. 4. In Sect. 5 we illustrate the enhancement of signal significance using multivariate techniques, and we assess the robustness of our results against the effects of PU. In Sect. 6 we conclude and outline future studies to estimate the accuracy in the determination of the trilinear coupling and to provide constraints in BSM scenarios.

Modeling of signal and background processes

In this section we discuss the Monte Carlo generation of the signal and background process samples used in this analysis. We shall also discuss the modeling of detector resolution effects.

Higgs pair production in gluon fusion

Higgs pair production is simulated at leading order (LO) using MadGraph5_aMC@NLO [63]. We use a tailored model [34] for gluon-fusion Higgs boson pair production which includes mass effects from the exact form factors for the top-quark triangle and box loops [64]. Equivalent results can be obtained using the recently available functionalities for the calculation of loop-induced processes [65] in MadGraph5_aMC@NLO. The calculation is performed in the scheme, accounting for b-quark mass effects. The renormalisation and factorisation scales are taken to be , withthe scalar sum of the transverse masses of all final-state particles. For the input parton distribution functions (PDFs) we adopt the NNPDF 3.0 LO set [66] with , interfaced via LHAPDF6 [67]. The Higgs boson couplings and branching ratios are set to their SM values, and its mass is taken to be  GeV [68-70]. In the SM, the Higgs trilinear coupling is given by , with  GeV the Higgs vacuum expectation value.

In Fig. 1 we show representative Feynman diagrams for LO Higgs pair production in gluon fusion. The non-trivial interplay between the heavy quark box and the triangle loop diagrams can lead to either constructive or destructive interference and complicates the extraction of the trilinear coupling from the measurement of the Higgs pair production cross section. Higher-order corrections [17, 18] are dominated by gluon radiation from either the initial-state gluons or from the heavy quark loops.

Fig. 1

Representative Feynman diagrams for Higgs pair production in gluon fusion at leading order. Only the fermion triangle loop diagram (right) is directly sensitive to the Higgs trilinear coupling . In the SM, the fermion loops are dominated by the contribution from the top quark

The total inclusive cross section for this processes is known up to NNLO [18]. Resummed NNLO+NNLL calculations for Higgs pair production are also available [19], leading to a moderate enhancement of the order of few percent as compared to the fixed-order NNLO calculation. To achieve the correct higher-order value of the integrated cross section, we rescale our LO signal sample to match the NNLO+NNLL inclusive calculation. This corresponds to a K-factor , as indicated in Table 1.

Table 1

Details of the signal and background Monte Carlo samples used in this work. Also provided are the inclusive K-factors which are applied to reproduce the known higher-order results

Process	Generator	N_evt	σ_LO (pb)	K-factor
pp→hh→4b	MadGraph5_aMC@NLO	1M	6.2×10^-3	2.4 (NNLO+NNLL [18, 19])
pp→bb¯bb¯	SHERPA	3M	1.1×10^3	1.6 (NLO [63])
pp→bb¯jj	SHERPA	3M	2.7×10^5	1.3 (NLO [63])
pp→jjjj	SHERPA	3M	9.7×10^6	0.6 (NLO [77])
pp→tt¯→bb¯jjjj	SHERPA	3M	2.5×10^3	1.4 (NNLO+NNLL [78])

Parton-level signal events are then showered with the Pythia8 Monte Carlo [71, 72], version v8.201. We use the default settings for the modeling of the underlying event (UE), multiple parton interactions (MPI), and PU, by means of the Monash 2013 tune [73], based on the NNPDF2.3LO PDF set [74, 75].

Backgrounds

Background samples are generated at leading order with SHERPA [76] v2.1.1. As in the case of the signal generation, the NNPDF 3.0 LO set with strong coupling is used for all samples, and we use as factorisation and renormalisation scales . We account for all relevant background processes that can mimic the signal process. This includes QCD 4b multi-jet production, as well as QCD 2b2j and 4j production, and top-quark pair production. The latter is restricted to the fully hadronic final state, since leptonic decays of top quarks can be removed by requiring a lepton veto. Single Higgs production processes such as and (see Appendix A) along with electroweak backgrounds e.g. , are much smaller than the QCD backgrounds [28, 29] and are therefore not included in the present analysis.

The LO cross sections for the background samples have been rescaled so that the integrated distributions reproduce known higher-order QCD results. For the 4j sample, we rescale the LO cross section using the BLACKHAT [77] calculation, resulting in an NLO/LO K-factor of 0.6. For the 4b and 2b2j samples NLO/LO K-factors of 1.6 and 1.3, respectively, have been determined using MadGraph5_aMC@NLO [63]. Finally, the LO cross section for production has been rescaled to match the NNLO+NNLL calculation of Ref. [78], leading to a K-factor of 1.4. The K-factors that we use to rescale the signal and background samples are summarised in Table 1.

At the generation level, the following loose selection cuts are applied to background events. Each final-state particle in the hard process must have  GeV, and be located in the central rapidity region with . At the matrix-element level all final-state particles must also be separated by a minimum . We have checked that these generator-level cuts are loose enough to have no influence over the analysis cuts. From Table 1 we see that the and QCD 4b cross sections are of the same order of magnitude. However, the former can be efficiently reduced by using top quark reconstruction criteria. The bbjj cross section is more than two orders of magnitude larger than the 4b result, but it will be suppressed by the light and charm jet mis-identification rates, required to contribute to the 4b final state.

As a cross-check of the SHERPA background cross sections reported in Table 1, we have produced leading-order multi-jet samples using MadGraph5_aMC@NLO, benchmarked with the results for the same processes reported in Ref. [63]. Using common settings, we find agreement, within scale uncertainties, between the MadGraph5_aMC@NLO and SHERPA calculations of the multi-jet backgrounds.

Modeling of detector resolution

While it is beyond the scope of this work to perform a full detector simulation, it is important to include an estimate of detector effects in the analysis, particularly for the finite energy and angular resolutions which directly degrade the reconstruction of important kinematic variables, such as the invariant mass of the Higgs candidates. Here we simulate the finite energy resolution of the ATLAS and CMS hadronic calorimeters by applying a Gaussian smearing of the transverse momentum with mean zero and standard deviation for all final-state particles before jet clustering, that is,with a univariate Gaussian random number, different for each of the particles in the event. We take as a baseline value for the transverse-momentum smearing a factor of .

To account for the finite angular resolution of the calorimeter, the plane is divided into regions of , and each final-state particle which falls in each of these cells is set to the same and values of the center of the corresponding cell. Finally, the energy of each final-state particle is recalculated from the smeared , and values to ensure that the resulting four-momentum is that of a light-like particle, since we neglect all jet constituent masses in this analysis.

Our modeling of detector simulation has been tuned to lead to a mass resolution of the reconstructed Higgs candidates consistent with the hadronic mass resolutions of the ATLAS and CMS detectors [79-81], as discussed in Sect. 3.5.

Analysis strategy

In this section we describe our analysis strategy. First of all we discuss the settings for jet clustering and the strategy for jet b-tagging. Following this we discuss the categorisation of events into different topologies, and how the different topologies may be prioritised. We motivate our choice of analysis cuts by comparing signal and background distributions for representative kinematic variables. Finally, we describe the simulation of PU and validate the PU-subtraction strategy.

Jet reconstruction

After the parton shower, final-state particles are clustered using the jet reconstruction algorithms of FastJet [82, 83], v3.1.0. Here we use the following jet definitions:For the boosted and intermediate categories, which involve the use of large-R jets, we use jet substructure variables [90, 91] to improve the significance of the discrimination between signal and background events in the MVA. In particular we consider the following substructure variables:

Small-R jets. These are jets reconstructed with the anti- clustering algorithm [84] with radius. These small-R jets are required to have transverse momentum  GeV and pseudo-rapidity , within the central acceptance of ATLAS and CMS, and therefore within the region where b-tagging is possible.

Large-R jets. These jets are also constructed with the anti- clustering algorithm, now using a radius. Large-R jets are required to have  GeV and lie in a pseudo-rapidity region of . The more restrictive range in pseudo-rapidity as compared to the small-R jets is motivated by mimicking the experimental requirements in ATLAS and CMS related to the track-jet based calibration [85, 86]. In addition to the basic and acceptance requirements, large-R jets should also satisfy the BDRS mass-drop tagger (MDT) [87] conditions, where the FastJet default parameters of and are used. Before applying the BDRS tagger, the large-R jet constituents are reclustered with the Cambridge/Aachen (C/A) algorithm [88, 89] with . In the case of the analysis including PU, a trimming algorithm [106] is applied to all large-R jets to mitigate the effects of PU, especially on the jet mass. For further details, see Sect. 3.5.

Small-R subjets. All final-state particles are clustered using the anti- algorithm, but this time with a smaller radius parameter, namely . The resulting anti- (AKT03) jets are then ghost-associated to each large-R jets in order to define its subjets [7]. These AKT03 subjets are required to satisfy  GeV and , and they will be the main input for b-tagging in the boosted category.

The -splitting scale [87, 92]. This variable is obtained by reclustering the constituents of a jet with the algorithm [93], which usually clusters last the harder constituents, and then taking the distance measure between the two subjets at the final stage of the recombination procedure, with and the transverse momenta of the two subjets merged in the final step of the clustering, and the corresponding angular separation.

The ratio of 2-to-1 subjettiness [94, 95]. The N-subjettiness variables are defined by clustering the constituents of a jet with the exclusive algorithm [96] and requiring that N subjets are found, where is the of the constituent particle k and the distance from subjet i to constituent k. In this work we use as input to the MVA the ratio of 2-subjettiness to 1-subjettiness, namely which provides good discrimination between QCD jets and jets arising from the decay of a heavy resonance.

The ratios of energy correlation functions (ECFs) [97] and [98]. The ratio of energy correlation functions is defined as while is instead defined as a double ratio of ECFs, that is, The energy correlation functions are defined in [97] with the motivation that -point correlators are sensitive to N-prong substructure. The free parameter is set to a value of , as recommended by Refs. [97, 98].

Tagging of b-jets

In this analysis we adopt a b-tagging strategy along the lines of current ATLAS performance [91, 99], though differences with respect to the corresponding CMS settings [100, 101] do not modify qualitatively our results. For each jet definition described above, a different b-tagging strategy is adopted:For the b-tagging probability , along with the b-mistag probability of light () and charm () jets, we use the values , and .

Small-R jets. If a small-R jet has at least one b-quark among their constituents, it will be tagged as a b-jet with probability . In order to be considered in the b-tagging algorithm, b-quarks inside the small-R jet should satisfy  GeV [99]. The probability of tagging a jet is not modified if more than one b-quark is found among the jet constituents. If no b-quarks are found among the constituents of this jet, it can be still be tagged as a b-jet with a mistag rate of , unless a charm quark is present instead, and in this case the mistag rate is . Only jets that contain at least one (light or charm) constituent with  GeV can induce a fake b-tag. We attempt to b-tag only the four (two) hardest small-R jets in the resolved (intermediate) category. Attempting to b-tag all of the small-R jets that satisfy the acceptance cuts worsens the overall performance as the rate of fake b-tags increases substantially.

Large-R jets. Large-R jets are b-tagged by ghost-associating anti- (AKT03) subjets to the original large-R jets [7, 91, 102, 103]. A large-R jet is considered b-tagged if both the leading and the subleading AKT03 subjets, where the ordering is done in the subjet , are both individually b-tagged, with the same criteria as the small-R jets. Therefore, a large-R jet where the two leading subjets have at least one b-quark will be tagged with probability . As in the case of small-R jets, we only attempt to b-tag the two leading subjets, else one finds a degradation of the signal significance. The treatment of the b-jet mis-identification from light and charm jets is the same as for the small-R jets.

Event categorisation

The present analysis follows a strategy similar to the scale-invariant resonance tagging of Ref. [51]. Rather than restricting ourselves to a specific event topology, we aim to consistently combine the information from the three possible topologies: boosted, intermediate and resolved, with the optimal cuts for each category being determined separately. This approach is robust under variations of the underlying production model of Higgs pairs, for instance in the case of BSM dynamics, which can substantially increase the degree of boost in the final state.

The three categories are defined as follows:

Boosted category. An event which contains at least two large-R jets, with the two leading jets being b-tagged. Each of these two b-tagged, large-R jets are therefore candidates to contain the decay products of a Higgs boson.

Intermediate category. An event with exactly one b-tagged, large-R jet, which is assigned to be the leading Higgs candidate. In addition, we require at least two b-tagged, small-R jets, which must be separated with respect to the large-R jet by an angular distance of . The subleading Higgs boson candidate is reconstructed by selecting the two b-tagged small-R jets that minimise the difference between the invariant mass of the large-R jet with that of the dijet obtained from the sum of the two small-R jets.

Resolved category. An event with at least four b-tagged small-R jets. The two Higgs candidates are reconstructed out of the leading four small-R jets in the event by considering all possible combinations of forming two pairs of jets and then choosing the configuration that minimises the relative difference of dijet masses.

Fig. 2

Comparison of the distributions of the leading (left) and subleading (right) large-R jets in the boosted category, for signal and background events. Distributions have been normalised to unity. The total background is the sum of all components listed in Table 1

Once a Higgs boson candidate has been identified, its invariant mass is required to lie within a fixed window of width around the nominal Higgs boson mass of  GeV. Specifically we require the conditionwhere is the invariant mass of each of the two reconstructed Higgs candidates. This cut is substantially looser than the corresponding cut used in the typical ATLAS and CMS analyses [79, 80]. The motivation for such a loose cut is that further improvements of the signal significance will be obtained using an MVA. Only events where the two Higgs candidates satisfy Eq. (8) are classified as signal events.

These three categories are not exclusive: a given event can be assigned to more than one category, for example, satisfying the requirements of both the intermediate and the resolved categories at the same time. The exception is the boosted and intermediate categories, which have conflicting jet selection requirements.

This is achieved as follows. First of all we perform an inclusive analysis, and optimise the signal significance in each of the three categories separately, including the MVA. We find that the category with highest significance is the boosted one, followed by the intermediate and the resolved topologies, the latter two with similar significance. Therefore, when ascertaining in which category an event is to be exclusively placed: if the event satisfies the boosted requirements, it is assigned to this category, else we check if it suits the intermediate requirements. If the event also fails the intermediate category requirements, we then check if it passes the resolved selection criteria. The resulting exclusive event samples are then separately processed through the MVA, allowing for a consistent combination of the significance of the three event categories.

Motivation for basic kinematic cuts

We now motivate the kinematic cuts applied to the different categories, comparing representative kinematic distributions between signal and background events. First of all, we present results without PU, and then discuss the impact of PU on the description of the kinematic distributions. In the following, all distributions are normalised to their total integral.

In Fig. 2 we show the distributions of the leading and subleading large-R jets in the boosted category. We observe that the background distribution falls off more rapidly as a function of than the di-Higgs signal. On the other hand, the cut in cannot be too strong to avoid a substantial degradation of signal selection efficiency, specially taking into account the subleading large-R jet. This comparison justifies the cut of  GeV for the large-R jets that we impose in the boosted category.

Another selection requirement for the boosted category is that the two leading AKT03 subjets of the large-R jet should satisfy  GeV. To motivate this cut, in Fig. 3 we show the distribution in of the leading and subleading AKT03 subjets in the subleading large-R jet in events corresponding to the boosted category. It is clear from the comparison that the subjet spectrum is relatively harder in the signal with respect to the background. On the other hand, considering the subleading AKT03 subjet, this cut in cannot be too harsh to maintain a high signal selection efficiency. Therefore, as for the previous distribution, the chosen cut value is a compromise between suppressing backgrounds but keeping a large fraction of signal events is crucial.

Fig. 3

Same as Fig. 2 for the leading (left) and subleading (right) AKT03 subjets in the subleading Higgs candidate large-R jet

Turning to the resolved category, an important aspect to account for in the selection cuts is the fact that the distribution of the four leading small-R jets of the event can be relatively soft, especially for the subleading jets. As noted in [29], this is due to the fact that the boost from the Higgs decay is moderate; therefore the selection cuts for the small-R jets cannot be too large. In Fig. 4 we show the distribution in of the four leading small-R jets in signal and background events: we observe that both distributions peak at  GeV, with the signal distribution falling off less steeply at large . The feasibility of triggering on four small-R jets with a relatively soft distribution is one of the experimental challenges for exploiting the resolved category in this final state, and hence the requirement that  GeV for the small-R jets. In Fig. 4 we also show the rapidity distribution of the small-R jets in the resolved category. As expected, the production is mostly central, and more so in the case of signal events, since backgrounds are dominated by QCD t-channel exchange; therefore the selection criteria on the jet rapidity are very efficient.

Fig. 4

Same as Fig. 2, now for the and rapidity distributions of the small-R jets corresponding to the resolved selection

One of the most discriminating selection cuts is the requirement that the invariant mass of the Higgs candidate (di)jets must lie within a window around the nominal Higgs value, Eq. (8). In Fig. 5 we show the invariant mass of the leading reconstructed Higgs candidates, before the Higgs mass window selection is applied, for the resolved and boosted categories. While the signal distribution naturally peaks at the nominal Higgs mass, the background distributions show no particular structure. The width of the Higgs mass peak is driven both from QCD effects, such as initial-state radiation (ISR) and out-of-cone radiation, as well as from the four-momentum smearing applied to final-state particles as part of our minimal detector simulation.

Fig. 5

Same as Fig. 2 for the invariant mass distribution of the leading Higgs candidates in the resolved (left) and boosted (right) selections

The invariant mass of the di-Higgs system is another important kinematic distribution for this process. The di-Higgs invariant mass is a direct measure of the boost of the system, which in BSM scenarios can be substantially enhanced, for instance due to specific EFT operators [14]. One important advantage of the final state for di-Higgs production is that it significantly increases the reach in as compared to other channels with smaller branching ratios, such as or . In Fig. 6 we show the invariant mass distribution of the reconstructed Higgs pairs, comparing the resolved and the boosted categories.

Fig. 6

Same as Fig. 2 for the invariant mass distribution of the di-Higgs system , in the resolved (left) and boosted (right) categories

In the resolved case, we see that the distribution in is rather harder for the signal as compared to the background, and therefore one expects that cutting in would help signal discrimination. For the boosted category the overall trend of the distribution is different because of the selection criteria, and the distribution now peaks at higher values of the invariant mass. In this case, signal and background distributions are not significantly differentiated. Note that at parton level the distribution for signal events has a kinematic cut-off at  GeV, which is smeared due to parton shower and detector resolution effects.

In Fig. 7 we show the transverse momentum of the di-Higgs system, , for the resolved and boosted categories. Once more we see that the background has a steeper fall-off in than the signal, in both categories, therefore this variable should provide additional discrimination power, motivating its inclusion as one of the inputs for the MVA. In our LO simulation the distribution is generated by the parton shower, an improved theoretical description would require merging higher-multiplicity matrix elements [35] or matching to the NLO calculation [17],

Fig. 7

Same as Fig. 2 for the transverse momentum distribution of the di-Higgs system

We shall now investigate the discrimination power provided by jet substructure quantities. In Fig. 8 we show the distributions of representative substructure variables for the boosted category: the splitting scale , Eq. (3), the ECF ratio , Eq. (6), and the 2–to–1 subjettiness ratio , Eq. (5), all for the leading Higgs candidates, and also for the subleading Higgs candidates.

Fig. 8

Distribution of representative substructure variables in the boosted category at the end of the cut-based analysis, to be used as input to the MVA. From top to bottom and from left to right we show the splitting scale , the energy correlation ratio and the subjettiness ratio for the leading Higgs. In the case of the distributions for the subleading Higgs are also given

From Fig. 8 we observe how for these substructure variables the shapes of the signal and background distributions reflect the inherent differences in the internal structure of QCD jets and jets originating from Higgs decays. Signal and background distributions peak in rather different regions. For example, the splitting scale peaks around 80 GeV (40 GeV) for signal (background) events, while the distribution of the ECF ratio is concentrated at small values for signal and is much broader for background events. From Fig. 8 we also see the distributions of the subjettiness ratio are reasonably similar for both the leading and the subleading jets.

Impact of pileup

Now we turn to discuss how the description of kinematic distributions for signal and background processes are modified in the presence of pileup. To study the impact of PU, Minimum Bias events have been generated with Pythia8, and then superimposed to the signal and background samples described in Sect. 2. We have explored two scenarios, one with a number of PU vertices per bunch crossing of , and another with . In the following we adopt as our baseline, and denote this scenario by PU80. We have verified that the combined signal significance is similar if is adopted instead.

In order to subtract PU in hadronic collisions, a number of techniques are available [87, 102, 104–114].1 In this work, PU is subtracted with the SoftKiller (SK) method [111], as implemented in FastJet, whose performance has been shown to improve the commonly used area-based subtraction [104]. The idea underlying SoftKiller consists of eliminating particles below a given cut-off in their transverse momentum, , whose value is dynamically determined so that the event-wide transverse-momentum flow density vanishes, where is defined asand where the median is computed over all the regions i with area and transverse momentum in which the plane is partitioned.

From its definition in terms of the median, it follows that the value of will be dynamically raised until half of the regions have . The size and number of these regions is a free parameter of the algorithm—here we will use square regions with length . We restrict ourselves to the central rapidity region, , for the estimation of the flow density . The SoftKiller subtraction is then applied to particles at the end of the parton shower, before jet clustering.

In addition, jet trimming [106], as implemented in FastJet, is applied to large-R jets. The trimming parameters are chosen such that the constituents of a given jet are reclustered into subjets with . Subjets with transverse momentum less than 5 % of the total transverse momentum of the large-R jet are then removed. The use of trimming in addition to PU removal with SoftKiller is necessary to correct the jet mass in the boosted category, which is particularly susceptible to soft, wide-angle contaminations. No trimming is applied to the small-R jets and to the case without PU.

Fig. 10

The transverse momentum of the leading Higgs candidate (upper plots) and of the invariant mass of the di-Higgs system (lower plots) in the resolved (left) and boosted (right) categories. We compare the results without PU with those with PU80 and SK+Trim subtraction, as explained in the text

In Fig. 9 we show the invariant mass distributions of the Higgs candidates for signal events in the resolved and boosted categories. In the resolved category, we compare the results without PU with those with PU80, with and without SK subtraction. If PU is not subtracted, there is a large shift in the Higgs mass peak, by more than 30 GeV. Once SK subtraction is performed, we recover a distribution much closer to the no PU case, with only a small shift of a few GeV and a broadening of the mass distribution. In the boosted case, the comparison is performed between no PU, PU with only SK subtraction, and PU with both SK and trimming. We find that the mass distribution for jets to which no trimming is applied peaks at around 160 GeV, even after PU subtraction with SoftKiller. When trimming is applied in addition to SoftKiller, the distribution peaks close to the nominal Higgs mass, as in the case of the resolved category.

Fig. 9

The invariant mass distributions of Higgs candidates in signal events in the resolved (left) and boosted (right) categories. In the resolved category, we compare the results without PU with those with PU80 with and without SK subtraction. In the boosted case, the comparison is performed between no PU, PU with only SK subtraction, and PU with both SK and trimming

In Fig. 10 we compare the transverse momentum of the leading Higgs candidate, and the invariant mass of the di-Higgs system , in both the boosted and the resolved categories, between the no PU and the PU+SK+Trim cases. In the case of the distribution, the differences between the selection criteria for the resolved and boosted categories is reflected in the rightward shift of the latter. After subtraction, the effects of PU are small in the two categories. A similar behaviour is observed in the di-Higgs invariant mass distribution.

We can also assess the impact of PU on the substructure variables that will be used as input to the MVA in the boosted and intermediate categories. In Fig. 11 we show the 2-to-1 subjettiness ratio , Eq. (5), and the ratio of energy correlation functions , Eq. (6), for the leading Higgs candidate. We observe that the shapes of both substructure variables are reasonably robust in an environment including significant PU. Therefore we can consider the PU subtraction strategy as validated for the purposes of this study, although further optimisation should still be possible, both in terms of the SoftKiller and of the trimming input settings.

Fig. 11

Same as Fig. 10 for the substructure variables (left) and (right) for the leading Higgs candidates in the boosted category

It is also interesting to quantify how the relative differences between signal over background distributions are modified by the inclusion of PU. Considering the boosted category initially, in Fig. 12 we compare various kinematic distributions for signal and background events, with and without PU for the leading Higgs candidate: the transverse-momentum distribution , the of the leading AKT03 subjet, the 2-to-1 subjettiness ratio , and the splitting scale . We verify that the relevant qualitative differences between signal and background distributions are maintained in the presence of PU. This is especially noticeable for the substructure variables, which exhibit a similar discriminatory power both with and without PU.

Fig. 12

Comparison of kinematic distributions for the leading Higgs candidate, in the boosted category, for signal and background events in the case of PU subtraction with SK+Trim: its transverse momentum , the of its leading AKT03 subjet, and the substructure variables and

We can also perform a similar comparison for the resolved category. In Fig. 13 we compare the kinematic distributions for signal and background events, with and without PU, for the invariant mass and the transverse momentum of the leading Higgs candidate. Again, the PU-subtracted background distributions appear reasonably close to their counterparts without PU, and thus the distinctive features between signal and background are maintained after PU subtraction.

Fig. 13

Same as Fig. 12 for the resolved category

It is illustrative to determine the mass resolution obtained for the reconstructed Higgs candidates in the various cases considered in the present study. In Table 2 we indicate the shift of the fitted invariant mass peak as compared to the nominal Higgs mass, , and the corresponding width of the distribution, , obtained from fitting a Gaussian to the mass distributions of leading and subleading Higgs candidates in the resolved and boosted categories. We show results for three cases: without PU, with PU80 but without subtraction (only for the resolved category), and the same with SK+Trim subtraction.

Table 2

Resolution of the invariant mass distribution of reconstructed Higgs candidates in the resolved and boosted categories. We show three cases: no PU, with PU80 without subtraction (only for resolved), and the same with SK+Trim subtraction. We indicate the shift of the fitted invariant mass peak for the Higgs candidates as compared to the nominal Higgs mass , as well as the fitted Gaussian width

	mhreco-m_h (GeV)	σ_m_h (GeV)
Resolved category
No PU
Leading h	−3.8	8.5±0.2
Subleading h	−5.8	9.1±0.3
PU80
Leading h	+33	8.8±1.5
Subleading h	+31	11.7±3.3
PU80+SK
Leading h	+3.9	10.7±0.3
Subleading h	+2.1	10.5±0.3
Boosted category
No PU
Leading h	+2.0	8.2±0.5
Subleading h	+1.0	8.8±0.5
PU80+SK+Trim
Leading h	−2.2	8.7±0.7
Subleading h	−4.9	9.0±0.8

In both categories, we find a mass resolution of around 9 GeV in the case without PU. In the case of PU with SK+Trim subtraction, in the resolved category the mass resolution worsens only slightly to around 11 GeV, while in the boosted category we find the same resolution as in the no PU case. We also note that after SK+Trim subtraction, the peak of the invariant mass distributions of Higgs candidates coincides with the nominal values of within a few GeV for the two categories.

Pre-MVA loose cut-based analysis

In this section we present the results of the pre-MVA loose cut-based analysis described in the previous section, and provide cut flows for the different analysis steps. We study how the signal significance is affected if only the 4b component of the QCD multi-jet background is taken into account. This section presents the results in an environment without pileup; the following one contains those obtained including significant PU.

Cut flow and signal significance

Here we compare the cross sections for signal and background events at various stages of the analysis. We consider all relevant backgrounds (see Sect. 2), and discuss how results are modified in the case where only the 4b background is considered. In Table 3 the different steps of the cut flow in the present analysis are summarised, separated into the boosted, intermediate, and resolved topologies. The different analysis steps proceed as follows:Signal and background events satisfying all the analysis cuts up to the C2 level are then used as input for the MVA training, to be described next in Sect. 5.

Table 3

Definition of the cuts imposed successively for the three selections

	Boosted	Intermediate	Resolved
C1a	NjetsR10≥2	NjetsR04≥2, NjetsR10=1	NjetsR04≥4
		+p_T cuts and rapidity cuts
C1b	+N_MDT≥2	+NjetsR10=1 with MDT	+Higgs reconstruction
		+Higgs reconstruction
C1c	+m_h window cut
C2	+b-tagging

C1a: check that we have at least two large-R jets (in the boosted case), one large-R jet and at least 2 small-R jets (in the intermediate case) and at least four small-R jets (in the resolved case). In addition, require that these jets satisfy the corresponding thresholds;  GeV for large-R jets and  GeV for small-R jets, as well as the associated rapidity acceptance constraints.

C1b: the two leading large-R jets must be mass-drop tagged in the boosted category. In the intermediate category, the large-R jet must also be mass-drop tagged.

C1c: after the two Higgs candidates have been reconstructed, their invariant masses are required to lie within a window around , in particular between 85 and 165 GeV, Eq. (8).

C2: the b-tagging conditions are imposed (see Sect. 3.2), and the event is categorised exclusively into one of the three topologies, according to the hierarchy determined in Sect. 3.3.

In Table 4 we collect the values for the signal and background cross sections at the different analysis steps. Results are divided into the resolved, intermediate and boosted categories, and they are inclusive up to the C2 level, where exclusivity is imposed. In Table 4 we also provide the signal over background ratio, S / B, and the signal significance, , corresponding to an integrated luminosity of ab. These are computed either taking into account all the background components or the 4b QCD background only. We find that after b-tagging, the 2b2j component is of the same order of magnitude as the 4b component in all categories. This implies that the signal significance at the end of the cut-based analysis is degraded due to the contribution of light and charm jets being mis-identified as b-jets.

Table 4

The cross sections for the signal and the background processes at different steps of the analysis (see Table 3), for the resolved (upper), intermediate (middle) and boosted (lower table) categories, for the analysis without PU. For each step, the signal over background ratio S / B, and the signal significance for ab are also provided, considering either the total background, or only the 4b component

	hh4b	Total bkg	Cross section [fb]				S / B		S/B
			4b	2b2j	4j	tt¯	Tot	4b	Tot	4b
HL-LHC, resolved category, no PU
C1a	9	2.2×10^8	6.9×10^4	1.5×10^7	2.0×10^8	2.1×10^5	4.0×10^-8	1.3×10^-4	0.03	1.9
C1b	9	2.2×10^8	6.9×10^4	1.5×10^7	2.0×10^8	2.1×10^5	4.0×10^-8	1.3×10^-4	0.03	1.9
C1c	2.6	4.4×10^7	1.6×10^4	3.2×10^6	4.1×10^7	8.8×10^4	6.1×10^-8	1.6×10^-4	0.02	1.1
C2	0.5	4.9×10^3	1.7×10^3	2.9×10^3	2.1×10^2	47	1.1×10^-4	2.9×10^-4	0.4	0.6
HL-LHC, intermediate category, no PU
C1a	2.8	8.4×10^7	2.1×10^4	5.3×10^6	7.9×10^7	3.3×10^4	3.4×10^-8	1.3×10^-4	0.02	1.1
C1b	2.6	5.8×10^7	1.4×10^4	3.6×10^6	5.5×10^7	3.0×10^4	4.5×10^-8	1.9×10^-4	0.02	1.2
C1c	0.5	3.5×10^6	8.7×10^2	2.1×10^5	4.3×10^7	8.8×10^3	1.6×10^-7	6.1×10^-4	0.02	1.0
C2	0.09	1.8×10^2	56	96	22	3.1	5.3×10^-4	1.6×10^-3	0.4	0.6
HL-LHC, boosted category, no PU
C1a	3.9	4.6×10^7	1.1×10^4	2.9×10^6	4.3×10^7	2.4×10^4	8.2×10^-8	3.4×10^-4	0.03	2.0
C1b	2.7	3.7×10^7	7.5×10^3	2.1×10^6	3.5×10^7	2.2×10^4	7.4×10^-8	3.7×10^-4	0.03	1.7
C1c	1.0	3.9×10^6	8.0×10^2	2.3×10^5	3.7×10^6	7.1×10^3	2.6×10^-7	1.3×10^-3	0.03	2.0
C2	0.16	2.5×10^2	53	1.9×10^2	13	1.6	5.7×10^-4	2.7×10^-3	0.5	1.1

In the boosted category, at the end of the loose cut-based analysis, we find that around 500 events are expected at the HL-LHC, with a large number, , of background events. This leads to a pre-MVA signal significance of and a signal over background ratio of . From Table 4 it is also possible to compute the corresponding pre-MVA expectations for the LHC Run II with  fb: one expects in the boosted category around 50 signal events, with signal significance dropping down to . Such signal significances could have been enhanced by applying tighter selection requirements, but our analysis cuts have been left deliberately loose so that such optimisation may be performed by the MVA.

The resolved category benefits from higher signal yields, but this enhancement is compensated for by the corresponding increase in the QCD multi-jet background. In both resolved and intermediate categories the signal significance is , similar to that of the boosted category. A further drawback of the resolved case is that S / B is substantially reduced as compared to the boosted and intermediate cases.

Combining the results from the boosted, intermediate and resolved categories, we obtain an overall pre-MVA significance for the observation of the Higgs pair production in the final state at the HL-LHC of .

The role of light and charm jet mis-identification

One of the main differences in the present study as compared to previous work is the inclusion of both irreducible and reducible background components, which allows us to quantify the impact of light and charm jet mis-identification. Two recent studies that have also studied the feasibility of SM Higgs pair production in the final state are from the UCL group [28] and from the Durham group [29]. The UCL study is based on requiring at least four b-tagged anti- jets in central acceptance with  GeV, which are then used to construct dijets (Higgs candidates) with  GeV,  GeV and between the two components of the dijet. In addition to the basic selection cuts, the constraints from additional kinematic variables are included by means of a Boosted Decision Tree (BDT) discriminant. The backgrounds included are the 4b and 2b2c QCD multijets, as well as , Zh, and . For the HL-LHC, a signal significance of is obtained.

The Durham group study [29] requires events to have two C/A jets with  GeV, and in addition two b-tagged subjets inside each large-R jet with  GeV each. To improve the separation between signal and background, both the BDRS method and the Shower Deconstruction (SD) [116, 117] technique are used. The backgrounds considered are QCD 4b as well as , hZ and hW. At the HL-LHC, their best result is obtained by requiring two SD-tagged large-R jets, which leads to . Using the BDRS tagger results in slightly poorer performance.

From our results in Table 4, we observe that the signal significance for the boosted, intermediate, and resolved categories is increased to 1.1, 0.6 and 0.6, respectively, when only the QCD 4b background is included. Combining the signal significance in the three categories, we obtain , twice as large as the result found when all background components are included. Note the importance of the combination of the three exclusive event topologies, as opposed the exploitation of a single specific category. Taking into account the loose selection cuts, we see that our pre-MVA results including only the 4b background are consistent with those reported in previous studies.

From Table 4 we can compare the interplay between the reducible and irreducible components of the QCD backgrounds. In all cases, the 4b and 2b2j components have comparable magnitudes within the uncertainties from missing higher-order corrections. On the other hand, the 4j component is always substantially smaller. So while the 4j component can be safely neglected, the inclusion of the 2b2j component is essential to assess the feasibility of measuring Higgs pairs in this final state robustly, especially in the boosted category. This has the important consequence that a promising avenue to improve the prospects of this measurement would be to reduce, as much as possible, the light and charm jet mis-identification rate.

In Fig. 14 we show a comparison of the shapes of the 4b and 2b2j components of the QCD background for the transverse momentum of the leading Higgs candidate and for invariant mass of the reconstructed di-Higgs system in the resolved and boosted categories. The two components possess a rather similar shape for the two distributions, albeit with some differences. In the boosted category, the 4b component exhibits a less steep fall-off of the distribution at large , while in the resolved case the 2b2j component has a slightly harder distribution of the invariant mass . We also observe that the 2b2j distributions are affected by somewhat larger Monte Carlo fluctuations as compared to 4b, despite the large size of the initial sample.

Fig. 14

Upper plots Comparison of the shapes of the 4b and 2b2j components of the QCD background for the of the leading Higgs candidate in the resolved (left plot) and boosted (right plot) categories. Lower plots Same comparison for the invariant mass of the reconstructed di-Higgs system

In the resolved category, the cross section before b-tagging is two orders of magnitude larger in the 2b2j sample as compared to the 4b sample. After b-tagging, a naive assessment would suggest a suppression of the 2b2j cross section by a factor , as compared to the 4b component, since a total of four b-tags are required to classify the event as a Higgs candidate. In this case the ratio of 2b2j over 4b would be around , and therefore negligible. While we have checked that this expectation is borne out at the parton level, we find that when parton shower effects are accounted for the situation is different, due both to radiation of pairs and from selection effects. Due to these, the number of b quarks in the final state is increased substantially in the 2b2j component as compared to the parton level, while at the same time the number of events in the 4b sample with 4 b-jets passing selection cuts is reduced.

We can make these statements more quantitative in the following way. To first approximation, neglecting the contribution from charm mis-identification, the overall efficiency of the b-tagging requirements in the resolved category will be given by the following expression:with being the fraction of events satisfying all the selection requirements, where j jets out of the leading four jets of the event contain b quarks (with  GeV). Similar expressions can be derived for the boosted and intermediate categories.

The naive expectation is that all events in the 4b sample have and for , while the events in the 2b2j sample should have and zero otherwise. This leads to the ratio of overall b-tagging selection efficienciesHowever, after the parton shower, the above estimate is no longer accurate. First of all, we will have a non-negligible fraction with also in the 2b2j sample, due to b-quark pair radiation during the shower. Second, not all events in the 4b sample will lead to four small-R b-jets, due to a combination of selection cuts and parton shower effects.

In Table 5 we collect the values of for the signal and the three QCD background samples. We find that rather than the estimate Eq. (11), the correct ratio of b-tagging selection efficiencies is insteadThis suppression factor is of the same order as the ratio of 4b to 2b2j cross sections in the resolved category before b-tagging. This explains why the 2b2j contribution cannot be neglected as compared to the irreducible 4b component of the QCD background. A similar calculation from the numbers in Table 5 shows that, on the other hand, the 4j component of the background can be neglected.

Table 5

The relative fractions of events for the resolved selection for which out of the four leading small-R jets of the event, j jets contain at least one b-quark with  GeV. This information is provided for the di-Higgs signal events and for the three QCD background samples. The last column indicates the overall selection efficiency as defined in Eq. (10)

	n0(b-jet) (%)	n1(b-jet) (%)	n2(b-jet) (%)	n3(b-jet) (%)	n4(b-jet) (%)	EFF_b-tag (%)
hh→4b	0.1	3	25	53	20	8.5
QCD 4b	1	8	27	44	20	8.4
QCD 2b2j	9	42	49	1	0.1	0.04
QCD 4j	96	3.5	0.5	0.01	3×10^-4	2×10^-4

Multivariate analysis

At the end of the loose cut-based analysis, by combining the three event topologies, we obtain a signal significance of with all backgrounds (only QCD 4b) considered. This section describes how this signal significance can be enhanced when the cut-based analysis is complemented by multivariate techniques. These are by now a mature tool in high-energy physics data analysis, opening new avenues to improve the performance of many measurements and searches at high-energy colliders. In particular, the classification of events into signal and background processes by means of MVAs is commonly used in LHC applications [28, 46, 80, 118–120].

In this section, first we present the specific MVA that we use, based on feed-forward multi-layer neural networks. Then we introduce the input variables that are used in the MVA, including the jet substructure variables, and then present the signal significance obtained by applying the MVA. Then we assess the robustness of the MVA strategy in the case of significant contamination from pileup.

Deep artificial neural networks

The specific type of MVA that we use to disentangle signal and background events is a multi-layer feed-forward artificial neural network (ANN), known as a perceptron.2 This family of ANNs are also known as deep neural networks, due to their multi-layered architecture. The MVA inputs are a set of kinematic variables describing the signal and background events which satisfy the requirements of the cut-based analysis. The output of the trained ANNs also allows for the identification, in a fully automated way, of the most relevant variables in the discrimination between signal and background.

In this work, the ANN that we use has the following architecture.where represents the number of input variables for the MVA, which is different in the resolved, intermediate, and boosted categories. All neural-network layers use a sigmoid activation function, allowing for a probabilistic interpretation of the ANN output. In Fig. 15 we show an illustrative example of an ANN used in this work, corresponding to the case of the boosted category (thus , as we explain below).

Fig. 15

Schematic of the Artificial Neural Network (ANN) used for the analysis of the boosted category, with input variables and thus the same number of neurons in the first layer. The color code in the neuron connections (the weights) is a heat map obtained at the end of the Genetic Algorithms training, with red indicating larger values and black indicating smaller values

The training of the ANN for the signal/background classification task proceeds as follows. Given a set of kinematic variables associated with the event i and a set of neural-network weight parameters , we interpret the neural-network output (the activation state of the neuron in the last layer) as the probability that the event i originates from the signal process,where represents the true classification of the event i, i.e., for signal and for background events. With this interpretation, our general classification probability including background events is given byconsequently we can define an error function to be minimised during the ANN training. In this case, the error function is the cross-entropy function, defined aswhere is the number of Monte Carlo events that are used for the ANN training. The ANN is trained both on the signal and background MC events, so it is important to ensure that the input MC sample is large enough to avoid contamination from MC statistical fluctuations.

The training of the neural networks therefore consists of the minimisation of the cross-entropy error, Eq. (16), which in this work is achieved using a Genetic Algorithm (GA). GAs [125-128] are non-deterministic minimisation strategies suitable for the solution of complex optimisation problems, for instance when a very large number of quasi-equivalent minima are present. GAs are inspired on natural selection processes that emulate biological evolution. In our case, the GA training is performed for a very large number of generations, , to avoid the risk of under-training. We have verified that if a much larger number of generations are used, the results are unchanged.

In addition, in order to avoid the possibility of over-fitting, we have used a cross-validation stopping criterion, in particular the same one as that used in the NNPDF3.0 analysis [66]. This cross-validation proceeds by dividing the input MC dataset into two disjoint sets, using one for training the ANN and the other for validation: the optimal stopping point is then given by the minimum of the error function Eq. (16) to the validation sub-sample. This indicates the point where the ANN begins to train upon statistical fluctuations in the input MC samples, rather than learning the underlying (smooth) physical distributions.

Input kinematic variables

In this work we use different sets of input variables for the three categories. In the case of large-R jets, we exploit the available information on jet substructure. For the three categories, boosted, intermediate and resolved, the following common variables are used as input to the MVA:In addition, in the boosted category we use the transverse momenta of the leading, and and subleading, and , Higgs candidate AKT03 subjets. In the resolved category instead, the corresponding variables are the transverse momenta of the four leading b-tagged small-R jets in the event. In the intermediate category, we use the transverse momenta of the subjets from the large-R jet and and the transverse momenta of the two leading b-tagged small-R jets. Therefore, we have 13 variables which are common to the three categories.

The transverse momenta of the leading and subleading Higgs, and .

The transverse momentum of the reconstructed Higgs pair, .

The invariant masses of the leading and subleading Higgs candidates, and .

The invariant mass of the reconstructed Higgs pair, .

The separation in the – plane between the two Higgs candidates, .

The separation in between the two Higgs candidates, .

The separation in between the two Higgs candidates, .

In the boosted and intermediate categories, we also include the jet substructure variables introduced in Sect. 3 for the large-R jets: the splitting scales , the ratio of 2-to-1 subjettiness , and the ratios of energy correlation functions and . This leads to a total of and 21 variables for the resolved, intermediate and boosted categories, respectively.

Given that the MVA is able to identify the most discriminatory variables in an automated way, and to suppress those which have little effect, it is advantageous to include a wide array of input variables. This is one of the main advantages of ANNs in this context: their inherent redundancy means that adding additional information, even if carries very little weight, should not degrade the classification power of the MVA.

MVA results

We now present the results of the MVA, first without PU, and then later including the effects of PU. First of all, in Fig. 16 we show the distribution of the ANN output at the end of the GA minimisation, separately for the boosted, intermediate and resolved categories. All distributions are normalised so that their integral adds up to one. The separation between signal and background is achieved by introducing a cut, , on the ANN output, so that MC events with are classified as signal events, and those with as background events. Therefore, the more differentiated the distribution of the ANN output is for signal and background events, the more efficient the MVA discrimination will be.

Fig. 16

The distributions, at the end of the GA training, for the signal and background MC events in the three categories: boosted (upper plot), intermediate (lower left plot) and resolved (lower right plot), as a function of the ANN output

From Fig. 16 we see that in the boosted category the MVA can produce a clear discrimination between signal and background, with the two distributions forming peaks at their respective optimal limits. This indicates that introducing a suitable cut in the ANN output will substantially reduce the background, while keeping a reasonable signal efficiency. The performance of the MVA discrimination is similar, although slightly worse, in the intermediate and resolved categories.

The results for the signal selection efficiency and the background rejection rate as a function of the cut in the ANN output define the so-called Receiver-Operating Characteristic (ROC) curve, shown in Fig. 17. It is clear that we can achieve high signal efficiency by using a small value of , but such a choice would be affected by poor background rejection. Conversely, using a higher value of the cut will increase background rejection at the cost of dropping signal efficiency. As could already be inferred from the distribution of neural-networks output in Fig. 16, we find that our MVA is reasonably efficient in discriminating signal over background. The performance is best in the case of the boosted category, and then slightly worse in the resolved and intermediate categories, consistent with the distributions of the ANN outputs in Fig. 16.

Fig. 17

Left ROC curve for the background rejection rate as a function of the signal selection efficiency, as the cut in the ANN output is varied. Right Number of signal (dashed) and background (solid) events expected at the HL-LHC as a function of the

It is useful to estimate, for each value of the cut in the ANN output , how many signal and background events are expected at the HL-LHC with  ab. This comparison is shown in Fig. 17. We observe that in the boosted category, for a value we end up with around 300 signal events and background events. Similar results are obtained in the intermediate and resolved categories: in the former we find 130 () signal (background) events for (0.60), and in the latter 630 () signal (background) events for . Therefore, the MVA achieves a substantial background suppression with only a moderate reduction of signal efficiency.

A useful property of MVAs such as the one used in our analysis is that they can provide direct physical insight about which of the input variables contribute to the separation between signal and background. In the case of ANNs, this can be quantified by computing the sum of the absolute values of all the weights connected to a given input neuron i, that is,with the value of the weight connecting the kth neutron of the second layer with the ith neuron of the first (input) layer, and the number of neurons in the second layer. Those input variables with a larger value of will be those that play a more significant role in enhancing the signal discrimination using the MVA. We note, however, that the estimate provided by Eq. (17) is necessarily qualitative.

In Fig. 18 we show the distribution of the total associated weight, Eq. (17) for each of the input variables of the three categories, using the notation for the kinematic variables as in Sect. 5.2. In the resolved category, the variables that carry a higher discrimination power are the transverse momentum of the two reconstructed Higgs candidates and their invariant masses. In the case of the boosted category, the invariant mass distribution of the Higgs candidates is also the most discriminatory variable, followed by the subjet distributions and substructure variables such as and .

Fig. 18

Distribution of the total associated weight, Eq. (17) for each of the input variables of the resolved (upper left), intermediate (upper right) and boosted (lower plot) categories

The results for the signal significance and the signal over background ratio S / B as a function of for the three categories are given in Fig. 19. The values for correspond to those at the end of the loose cut-based analysis. We observe how in the three categories there is a marked improvement in signal significance as compared to the pre-MVA results. We also observe a substantial enhancement in S / B, arising from the background suppression achieved by the MVA, reaching values of 1, 6 and 3.5 % in the resolved, intermediate and boosted categories. This improvement in S / B is crucial to ensure the feasibility of this measurement, since it allows systematic uncertainties in the background determination to be at most of a similar size.

Fig. 19

The values of the signal significance, , and of the signal over background ratio, S / B, for the boosted, intermediate and resolved categories as a function of the cut in the ANN output. The results are those at the end of the cut-based analysis

The optimal value of the cut in the ANN output, , can be determined from the maximisation of , ensuring that the number of signal events expected at the HL-LHC does not become too low. In addition, we require that the number of MC events used to define the signal category (events with ) is sufficiently large in order to avoid the biases and statistical fluctuations associated to a small training sample. In Table 6 we quote, for the optimal value of in each category, the number of signal and background events expected at the HL-LHC, as well as and S / B. For completeness, we also include the corresponding pre-MVA results.

Table 6

Post-MVA results, for the optimal value of the ANN discriminant in the three categories, compared with the corresponding pre-MVA results (). We quote the number of signal and background events expected for ab, the signal significance and the signal over background ratio S / B. The pre-MVA results correspond to row C2 in Table 4

HL-LHC, no PU
Category	N_ev signal	N_ev back	S/B	S / B
Boosted
y_cut=0	440	7.6×10^5	0.5	6×10^-4
y_cut=0.90	290	1.2×10^4	2.7	0.03
Intermediate
y_cut=0	280	5.3×10^5	0.4	5×10^-4
y_cut=0.85	130	3.1×10^3	2.3	0.04
Resolved
y_cut=0	1500	1.5×10^7	0.4	1×10^-4
y_cut=0.60	630	1.1×10^5	1.9	0.01

From Table 6 we see that following the application of the MVA, the signal significance in the boosted category increases from 0.5 to 2.7, with S / B increasing from 0.06 to . For the intermediate and resolved categories, increases from 0.4 to 2.3 and 1.9, respectively, with the signal over background ratio raising from 0.05 and to 4 and 1 %. Combining the three categories, taking into account all background components, we obtain the overall signal significance:The signal significance for ab is thus well above the threshold for the observation of Higgs pair production. However, given that the HL-LHC will be a high-PU environment, which will affect the description of the various kinematic distributions used as input to the MVA, it is essential to quantify the robustness of these results in a realistic environment including the effects of significant PU.

It should be emphasised that MVAs such as the ANNs used in this work can always be understood as a combined set of correlated cuts. Once the ANNs have been trained, it is possible to compare kinematical distributions after and before the ANN cut to verify its impact. This information would allow one in principle to perform a cut-based analysis, without the need of using ANNs, and finding similar results.

To illustrate this point, in Fig. 20 we show the distribution of the leading AKT04 small-R jets and the invariant mass of reconstructed Higgs candidates in the resolved category, comparing the pre-MVA results () with the post-MVA results () for signal and background events. The distributions are not normalised, to better visualise the effect of the MVA cut. Unsurprisingly, the ANN cut effectively selects events which lead to similar kinematical distributions between signal and background events. In the case of the small-R jets distribution, the ANN cuts favours the high- region, while for the invariant mass distribution only the region around the Higgs mass peak is selected for background events.

Fig. 20

The distribution of the leading AKT04 small-R jets (upper plots) and the invariant mass of reconstructed Higgs candidates (lower plots) in the resolved category, comparing the pre-MVA results () with the post-MVA results () for signal (left) and background (right plot) events. In this case the distributions are not normalised, to better visualise the effects of the MVA cut

A particularly challenging aspect of our analysis is the modeling of the 2b2j and 4j background, especially for the latter, which require extremely large MC samples. In the analysis reported here, out of the original 3M 4j generated events, only around 100 survive the analysis cuts, and thus these low statistics have associated a potentially large uncertainty in the calculation of the post-MVA 4j cross section. On the other hand, since the 4j cross sections are always quite smaller than the sum of the 4b of the 2b2j components, these low statistics should not modify qualitatively our conclusions above. To verify explicitly this expectation, and obtain a more robust estimate of the background cross section from mis-identified jets, we have increased by a factor 10 the size of the 2b2j and 4j background samples, up to a total of 30M each. Processing these events though our analysis, including retraining the MVA, we find , consistent with Eq. (18), indicating that the low statistics of the 4j background is not a limiting factor.

Impact of PU in the MVA

In this section we study how the MVA results are modified when the analysis is performed including significant PU. The loose cut-based analysis and the subsequent MVA optimisation have been performed using the same settings as in the case without PU. In Table 7 we provide the pre-MVA cut flow in the case of PU80, the corresponding version without PU being Table 4. The interplay between the signal cross sections and the various background components is qualitatively unchanged as compared to the no PU case.

Table 7

Same as Table 4, now for the case of PU80+SK+Trim

	hh4b	Total bkg	Cross section [fb]				S / B		S/B
			4b	2b2j	4j	tt¯	Tot	4b	Tot	4b
HL-LHC, resolved category, PU+SK with n_PU=80
C1a	11	4.4×10^8	1.5×10^5	3.0×10^7	4.1×10^8	2.6·10^5	2.4×10^-8	7.2×10^-5	0.03	1.5
C1b	11	4.4×10^8	1.5×10^5	3.0×10^7	4.1×10^8	2.6·10^5	2.4×10^-8	7.2×10^-5	0.03	1.5
C1c	3	1.1×10^8	4.2×10^4	7.7×10^6	9.9×10^7	1.1·10^5	2.8×10^-8	7.4×10^-5	0.02	0.8
C2	0.6	9.0×10^3	3.5×10^3	5.1×10^3	3.1×10^2	50	6.5×10^-5	1.7×10^-4	0.4	0.5
HL-LHC, intermediate category, PU+SK+Trim with n_PU=80
C1b	2.7	8.1×10^7	2.1×10^4	5.2×10^6	7.6×10^7	3.0×10^4	3.4×10^-8	1.3×10^-4	0.02	1.0
C1c	2.6	6.2×10^7	1.5×10^4	3.9×10^6	5.8×10^7	2.8×10^4	4.1×10^-8	1.7×10^-4	0.02	1.1
C1d	0.5	2.8×10^6	7.9×10^2	1.9×10^5	2.7×10^6	6.5×10^3	1.8×10^-7	6.2×10^-4	0.02	1.0
C2	0.09	2.6×10^2	47	1.8×10^2	30	2.2	3.4×10^-4	1.8×10^-3	0.3	0.7
HL-LHC, boosted category, PU+SK+Trim with n_PU=80
C1a	3.5	4.1×10^7	1.0×10^4	2.7×10^6	3.8×10^7	2.0×10^4	8.6×10^-8	3.4×10^-4	0.03	1.9
C1b	2.5	3.2×10^7	6.8×10^3	1.9×10^6	3.0×10^7	1.9×10^4	7.8×10^-8	3.6×10^-4	0.02	1.6
C1c	0.8	2.2×10^6	5.4×10^2	1.4×10^5	2.0×10^6	4.8×10^3	3.8×10^-7	1.6×10^-3	0.03	2.0
C2	0.14	1.5×10^2	40	86	22	1.8	9.0×10^-4	3.5×10^-3	0.6	1.2

In Table 8 we compare the results for the PU80+SK+Trim case between the pre-MVA loose cut-based analysis and the post-MVA results for the optimal values of the ANN output cut . As in Table 6, we also quote the number of signal and total background events expected for ab and the values of and S / B. We observe that the pre-MVA signal significance is close to the results of the simulations without PU for the three categories. We now find values for of 0.4, 0.3 and 0.6, in the resolved, intermediate and boosted categories, respectively, to be compared with the corresponding values without PU, namely 0.4, 0.4 and 0.5. The number of selected signal events in each category at the end of the cut-based analysis is only mildly affected by PU. The slight pre-MVA improvement in for the boosted case arises from a reduction in the number of background events that are classified in this category as compared to the case without PU.

Table 8

Same as Table 6, now for the case of PU80+SK+Trim

HL-LHC, PU80+SK+Trim
Category	N_ev signal	N_ev back	S/B	S / B
Boosted
y_cut=0	410	4.5×10^5	0.6	10^-3
y_cut=0.8	290	3.7×10^4	1.5	0.01
Intermediate
y_cut=0	260	7.7×10^5	0.3	3×10^-4
y_cut=0.75	140	5.6×10^3	1.9	0.03
Resolved
y_cut=0	1800	2.7×10^7	0.4	7×10^-5
y_cut=0.60	640	1.0×10^5	2.0	0.01

Once the MVA is applied, the signal significance in the resolved, intermediate and boosted categories increases to 2.0, 1.9 and 1.5 respectively, to be compared with the corresponding values without PU, namely 1.9, 2.3 and 2.7. Therefore, the post-MVA effect of PU on is a moderate degradation of the boosted and intermediate categories, especially for the former, while the resolved category is largely unchanged.3 We also observe that, due to the MVA, the signal over background ratio is increased from 0.007, 0.03 and 0.1 % up to 1, 3 and 1 % in the resolved, intermediate and boosted categories, respectively. This indicates that, while this measurement is still highly challenging, requiring a careful extraction of the QCD background from the data, it should be within reach.

In Fig. 21 we show the number of signal and background events that are expected for  ab as a function of , together with the corresponding ROC curve. The slight degradation of the boosted category in the case of PU can be seen by comparing with the corresponding results without PU in Fig. 17. In Fig. 22 we show the signal significance, , and the signal over background ratio, S / B, accounting now for the effects of PU. The corresponding results in the case without PU were shown in Fig. 19. As can be seen, the MVA-driven enhancement remains robust in the presence of PU, with only moderately degraded. Therefore, the qualitative conclusions drawn in the case without PU also hold when the analysis is performed in a high-PU environment. Since no specific effort has been made to optimise PU subtraction, for instance by tuning the values of the patch length a in SoftKiller or the threshold during jet trimming, we believe that there should be still room for further improvement.

Fig. 21

Same as Fig. 17 for the PU80+SK+Trim case

Fig. 22

Same as Fig. 19 for the PU80+SK+Trim case

It is useful to quantify which of the MVA input variables carry the highest discrimination power in the case of PU, by means of Eq. (17), and compare this with the corresponding results without PU shown in Fig. 18. We have verified that the relative weight of the different input variables to the MVA is mostly unchanged in the case of PU. In the resolved category, the highest total associated weight is carried by the Higgs candidates and invariant mass, as well as by the of the individual small-R jets. For the boosted category, the highest weight is carried by the Higgs invariant mass, followed by the Higgs ,  , the of the AKT03 subjets and the substructure variables, with a similar weighting among them.

In Table 9 we provide the post-MVA number of signal and background events expected for  ab. For the backgrounds, we quote both the total number, , and the QCD 4b component only, . We quote results for the no PU and PU80+SK+Trim cases. We also quote in each case the corresponding values for the signal significance and the signal over background ratio. Note that the MVA is always trained to the inclusive background sample, though differences in the kinematic distributions of the 4b and 2b2j processes are moderate; see Fig. 14. From Table 9 one observes that all categories exhibit a marked improvement from eliminating the contamination from light and charm jet mis-identification. For instance, in the intermediate category, increases from 2.3 to 3.3 (1.9 to 2.9) in the no PU (PU80) case, with similar improvements in the resolved and boosted categories.

Table 9

Post-MVA number of signal and background events with  ab. For the backgrounds, both the total number, , and the 4b component only, , are shown. Also provided are the values of the signal significance and the signal over background ratio, both separated in categories and for their combination. We quote the results without PU and for PU80+SK+Trim

Category	Signal	Background		S/B_tot	S/B_4b	S/B_tot	S/B_4b
	N_ev	Nevtot	Nev4b
Boosted
No PU	290	1.2×10^4	8.0×10^3	2.7	3.2	0.03	0.04
PU80+SK+Trim	290	3.7×10^4	1.2×10^4	1.5	2.7	0.01	0.02
Intermediate
No PU	130	3.1×10^3	1.5×10^3	2.3	3.3	0.04	0.08
PU80+SK+Trim	140	5.6×10^3	2.4×10^3	1.9	2.9	0.03	0.06
Resolved
No PU	630	1.1×10^5	5.8×10^4	1.9	2.7	0.01	0.01
PU80+SK	640	1.0×10^5	7.0×10^4	2.0	2.6	0.01	0.01
Combined
No PU				4.0	5.3
PU80+SK+Trim				3.1	4.7

In Table 9 we also provide the results for obtained by combining the three categories. Taking into account all background components, we obtain for the case of an overall signal significance ofindicating that a measurement of Higgs pair production in the final state at the HL-LHC should be above the threshold for observation, even when realistic PU conditions are accounted for. A similar signal significance is obtained in the case of . Under the assumption that the only relevant background would be the irreducible QCD 4b component, one obtains insteadTherefore, a measurement of Higgs pair production in the final state at the HL-LHC might be even above the threshold for discovery, provided the effects due to mis-identification of light and charm jets as b-jets can be reduced.

Conclusions and outlook

In this work we have presented a feasibility study for the measurement of Higgs pair production in the final state at the LHC. Our strategy is based on the combination of traditional cut-based analysis with state-of-the-art multivariate techniques. We take into account all relevant backgrounds, in particular the irreducible 4b and the reducible 2b2j and 4j QCD multijets. We have illustrated how the 2b2j component leads to a contribution comparable to that of QCD 4b production, due to a combination of parton shower effects, b-quark pair radiation, and selection requirements. We have also demonstrated the robustness of our analysis strategy under the addition of significant PU. In particular, we have explored two scenarios, and , and we found a comparable overall signal significance in the two cases.

Combining the contributions from the resolved, intermediate and boosted categories, we find that, for ab, the signal significance for the production of Higgs pairs turns out to be . This indicates that, already from the final state alone, it should be possible to claim observation of Higgs pair production at the HL-LHC. Our study also suggests possible avenues that the LHC experiments could explore to further improve this signal significance. One handle would be to reduce the contribution from light and charm jet mis-identification, ensuring that the irreducible 4b background dominates over the 2b2j component. This would allow one to enhance almost to the discovery level; see Table 9. It would also be advantageous to improve the b-tagging efficiency, allowing to achieve higher signal yields. Another possibility would be to improve the mass resolution of the Higgs reconstruction in high-PU environments, and, more generally, to optimise the PU subtraction strategy in order to reduce the impact of PU in the modeling of kinematic variables and the associated degradation in the MVA discrimination.

Another challenging aspect of the measurement of Higgs pairs in the final state is achieving an efficient triggering strategy. In order to reduce the rate from background QCD processes sufficiently, while being able to access the relevant regimes, (multi-)jet triggers using b-quark tagging information online for one or more jets are likely to be necessary. The additional rejection provided by these triggers could enable events to be selected efficiently, with four jets down to  GeV in the resolved category, and boosted Higgs decays in large-R jets down to jet transverse momenta of  GeV. In addition, good control of the multijet backgrounds and the experimental systematics of the MVA inputs will be important to achieve these sensitivities.

Our strategy relies on the modeling of the kinematic distributions of signal and background events, since these provide the inputs to the MVA discriminant. In this respect, it would be important, having established the key relevance of the channel for the study of Higgs pair production, to revisit and improve the theoretical modeling of our signal and background simulation, in particular using NLO calculations matched to parton showers both for signal [17, 35] and for backgrounds [63, 76].

One important implication of this work is that it should be possible to significantly improve the accuracy on the extraction of the Higgs trilinear coupling from a measurement of the cross section, as compared to existing estimates. A determination of in our approach is however rather non-trivial, involving not only generating signal samples for a wide range of values of , but also repeating the analysis optimisation, including the MVA training, for each of these values. This study is left to a future publication, where we will also compare the precision from the final state with the corresponding precision that has been reported from other final states such as and . It will also be interesting to perform this exercise for a 100 TeV hadron collider [11-14]. While at 100 TeV the signal yields would be increased, also the (gluon-driven) QCD multijet background would grow strongly. Revisiting the present analysis, including the MVA optimisation, at 100 TeV would also allow us to assess the accuracy of an extraction of the trilinear coupling from the final state at 100 TeV.

In this work we have considered only the SM production mechanism, but many BSM scenarios predict deviations in Higgs pair production, both at the level of total rates and of differential distributions. In the absence of new explicit degrees of freedom, deviations from the SM can be parametrised in the EFT framework using higher-order operators [14, 48]. Therefore, we plan to study the constraints on the coefficients of these effective operators that can be obtained from measurements of various kinematic distributions in the process. Note that the higher rates of the final state as compared to other final states, such as , allow for better constraints upon operators that modify the high-energy behaviour of the theory, for instance, it would become possible to access the tail of the distribution.

As in the case of the extraction of the Higgs trilinear coupling , such a study would be a computationally intensive task, since BSM dynamics will modify the shapes of the kinematic distributions and thus in principle each point in the EFT parameter space would require a re-optimisation with a newly trained MVA. In order to explore efficiently the BSM parameters without having to repeat the full analysis for each point, modern statistical techniques such as the Cluster Analysis method proposed in Ref. [46] might be helpful.

Table 10

LO and NLO cross sections at the generation level for the single-Higgs background processes listed above, computed using MadGraph5_aMC@NLO. The subsequent decays and the corresponding branching fractions are not included in these generation-level cross sections

Sample	LO	NLO	K-factor
Zh (13 TeV)	6.5×10^-1 pb	7.7×10^-1 pb	1.19
tt¯h (13 TeV)	3.8×10^-1 pb	4.6×10^-1 pb	1.29
bb¯h (13 TeV)	4.9×10^-1 pb	6.1×10^-1 pb	1.22

Table 11

The values of the branching fractions applied to the single-Higgs background processes from Table 10, corresponding to the most updated PDG values

Sample	Decay	Branching fraction
Zh	(Z→bb¯)(h→bb¯)	0.086
tt¯h	(W→qq¯)^2(h→bb¯)	0.26
bh¯h	h→bb¯	0.57

Table 12

Signal and background cross sections at the end of the cut-based analysis (before the MVA is applied), in the case without PU. We separate the results into the three exclusive categories used in our analysis

	Sample	Pre-MVA cross section (fb)
		Boosted	Intermediate	Resolved
Signal	hh→bb¯bb¯	3.5×10^-1	2.2×10^-1	1.2×10^0
Backgrounds	QCD multijet	2.5×10^+2	1.8×10^+2	4.9×10^+3
	Z(→bb¯)h(→bb¯)	2.0×10^-2	1.2×10^-1	7.5×10^-1
	tt¯h(→bb¯)	5.1×10^-2	6.3×10^-3	4.0×10^-1
	bb¯h(→bb¯)	2.3×10^-3	5.5×10^-3	2.6×10^-1