Higgs pair production in vector-boson fusion at the LHC and beyond.

The production of pairs of Higgs bosons at hadron colliders provides unique information on the Higgs sector and on the mechanism underlying electroweak symmetry breaking (EWSB). Most studies have concentrated on the gluon-fusion production mode which has the largest cross section. However, despite its small production rate, the vector-boson fusion channel can also be relevant since even small modifications of the Higgs couplings to vector bosons induce a striking increase of the cross section as a function of the invariant mass of the Higgs boson pair. In this work we exploit this unique signature to propose a strategy to extract the hhVV quartic coupling and provide model-independent constraints on theories where EWSB is driven by new strong interactions. We take advantage of the higher signal yield of the [Formula: see text] final state and make extensive use of jet-substructure techniques to reconstruct signal events with a boosted topology, characteristic of large partonic energies, where each Higgs boson decays to a single collimated jet. Our results demonstrate that the hhVV coupling can be measured with 45% (20%) precision at the LHC for [Formula: see text] (3000) fb[Formula: see text], while a 1% precision can be achieved at a 100 TeV collider.

Introduction

Following the discovery of the Higgs boson in 2012 [1, 2], the measurement of its couplings to the other standard model (SM) particles has become one of the main goals of the LHC programme. In this respect double Higgs production provides a unique handle, in particular since it allows the extraction of the trilinear Higgs self-coupling . In addition to constraining , in the vector-boson fusion (VBF) channel double Higgs production also probes the strength of the Higgs non-linear interactions with vector bosons at high energies. This process can thus help establish the nature of the Higgs boson, whether it is a composite or elementary state or whether or not it emerges as a Nambu–Goldstone boson (NGB) of some new dynamics at the TeV scale [3-5].

Many scenarios of new physics beyond the SM (BSM) generically predict enhanced cross sections for Higgs pair production with or without the resonant production of new intermediate states; see for example Refs. [6-20]. For this reason, searches for Higgs pair production at the LHC by ATLAS and CMS have already started at 8 TeV [21-25] and TeV [26-29], and will continue during Runs II and III, as well as at the high-luminosity LHC (HL-LHC) upgrade with of integrated luminosity. On the other hand, in the SM the small production rates make a measurement of Higgs pair production extremely challenging even at the HL-LHC, and the ultimate accuracy could only be achieved at a future TeV hadron collider [17, 30–34].

Similarly to single Higgs production, the dominant mechanism for Higgs pair production is the gluon-fusion mode [35]. This channel has been extensively studied in the literature and several final states have been considered, including , , and (for a list of feasibility studies, see for example Refs. [17, 30, 36–43]). Working in the infinite top mass approximation, the gluon-fusion di-Higgs production cross section was calculated at NLO in [44] and NNLO in [45]. The resummation of soft-gluon emissions was performed at NNLL in [46, 47]. Beyond the limit, the impact of top-quark mass effects on NLO QCD corrections was first determined in [48] through a reweighting technique based on an approximate two-loop matrix element and by [49, 50] in a expansion. Recently, the full NLO calculation was performed by [51]. Matching the fixed order computations to a parton shower was done at LO in [52] and at NLO in [53].

Recent studies indicate that Higgs pair production in gluon fusion at the HL-LHC will allow the extraction of the Higgs self-coupling with accuracy, with details varying with the analysis and the specific final state; see Refs. [54-58] for the latest ATLAS and CMS estimates, as well as [17, 30, 38, 59, 60].

Higgs pairs can also be produced in the VBF channel [3, 4, 61–64] where a soft emission of two vector bosons from the incoming protons is followed by the hard scattering, with . In the SM, the VBF inclusive cross section at 14 TeV is around 2 fb – more than one order of magnitude smaller than in gluon fusion. QCD corrections give a 10% increase and have been computed at NLO in Refs. [53, 65] and at NNLO in Ref. [63]. Production in association with W or Z bosons, known as the Higgsstrahlung process [65-67], or with top quark pairs [68], exhibit even smaller cross sections.

Despite its small rate, Higgs pair production via VBF is quite interesting since even small modifications of the SM couplings can induce a striking increase of the cross section as a function of the di-Higgs mass. Specific models leading to this behaviour are, for instance, those where the Higgs is a composite pseudo-NGB (pNGB) of new strong dynamics at the TeV scale [69]. In these theories, the Higgs anomalous couplings imply a growth of the cross section with the partonic centre-of-mass energy, , where f is the pNGB decay constant [3]. This enhanced sensitivity to the underlying strength of the Higgs interactions makes double Higgs production via VBF a key process to test the nature of the electroweak symmetry breaking dynamics and to constrain the hhVV quartic coupling. A first study of double Higgs production via VBF at the LHC was performed in Ref. [4], for a mass GeV, by focusing on the 4W final state. Following the discovery of the Higgs boson, more studies of the hhjj process at the LHC were presented in Refs. [61, 62, 64, 70].

In this work, we revisit the feasibility of VBF Higgs pair production at the LHC and focus on the final state. While this final state benefits from increased signal yields due to the large branching fraction of Higgs bosons to bottom quarks, in the SM [35], it also suffers from overwhelming large QCD multijet backgrounds. In this respect, the remarkable VBF topology, characterized by two forward jets well separated in rapidity and with a large invariant mass, together with a reduced hadronic activity in the central region, provides an essential handle to disentangle signal events from the QCD background. Additionally, the di-Higgs system will acquire a substantial boost in the presence of BSM dynamics. It is thus advantageous to resort to jet-substructure techniques [71] in order to fully exploit the high-energy limit and optimize the signal significance.

We will thus focus on the kinematic region where the invariant mass of the Higgs pair, , is large because modifications of the couplings between the Higgs and vector bosons cause the tail of this distribution to become harder in the signal whereas the background is not modified. Therefore, this region exhibits the highest sensitivity to the modified Higgs couplings and in particular to the deviations in the hhVV quartic coupling . Given that for large the Higgs bosons can be produced boosted, improved discrimination can be achieved using jet substructure, and to this end we use scale-invariant tagging [10, 43] to smoothly combine the resolved, intermediate and boosted topologies.

Our analysis takes into account all the main reducible and irreducible backgrounds: QCD multijet production, Higgs production via gluon fusion (where additional radiation can mimic the VBF topology), and top-quark pair production. We pay special attention to the role of light and charm jets being misidentified as b-jets which can contribute sizeably to the total background yield. For instance, in the channel, the 2b2j background is comparable to the 4b component [43].

We quantify the constraints on the Higgs quartic coupling that can be obtained from VBF di-Higgs production at the LHC 14 TeV with and as well as at a future circular collider (FCC) with a centre-of-mass energy of 100 TeV and a total luminosity of . We find that, despite the smallness of the production cross sections, the LHC with can already constrain the hhVV coupling with an accuracy of around its SM value at the 1- level, which is further reduced to at the HL-LHC and down to the 1% level at the FCC. Our results strongly motivate that searches for VBF Higgs pair production at the LHC should already start during Run II.

The structure of this paper is as follows. In Sect. 2 we present the general parametrization of the Higgs couplings which we adopt and review its impact on VBF Higgs pair production. Then in Sect. 3, we discuss the analysis strategy used to disentangle the signal from the background events in the final state. Our main results are presented in Sect. 4, where we quantify the potential of the VBF di-Higgs process to measure the hhVV coupling at various colliders and discuss the validity of the effective field theory expansion. Finally in Sect. 5, we conclude and discuss how our analysis strategy could be applied to related processes. Technical details are collected in three appendices which describe the Monte Carlo event generation of signal and background events (Appendix A), the fits to the tail of the distribution for backgrounds (Appendix B), and the validation studies of the QCD multijet event generation (Appendix C).

Higgs pair production via vector-boson fusion at hadron colliders

We begin by reviewing the theoretical framework for Higgs pair production via vector-boson fusion in hadronic collisions. First, we introduce a general parametrization of the Higgs couplings in the effective field theory (EFT) framework. Then we consider the values that these couplings take in specific models. Finally, we briefly discuss the validity of the EFT approximation and the possible contribution of heavy resonances to this process.

General parametrization of Higgs couplings

Following Ref. [4], we introduce a general parametrization of the couplings of a light Higgs-like scalar h to the SM vector bosons and fermions. At energies much lower than the mass scale of any new resonance, the theory is described by an effective Lagrangian obtained by making a derivative expansion. Under the request of custodial symmetry the three NGBs associated with electroweak symmetry breaking parametrize the coset SO(4) / SO(3). They can be fitted into a matrixwith GeV the Higgs vacuum expectation value. Assuming that the couplings of the Higgs boson to SM fermions scale with their masses and do not violate flavour, the resulting effective Lagrangian in [4] can be parametrized aswhere V(h) denotes the Higgs potential,The parameters , , , , and are in general arbitrary coefficients, normalized so that they equal 1 in the SM. The Higgs mass is fixed to be GeV [72].

As the notation in Eq. (2) indicates, the coefficients , , and control the strength of the hVV, hhVV and hhh couplings, respectively. The coefficients and instead modify the Higgs coupling to fermions and quartic self interaction. Thus, they do not affect the double-Higgs production cross section in the VBF channel. In Fig. 1, we show the tree-level Feynman diagrams, in the unitary gauge, that contribute to Higgs pair production in the vector-boson fusion channel at hadron colliders. In terms of the general parametrization of Eq. (2), the left, middle, and right diagrams scale with , , and , respectively.

Fig. 1

Tree-level Feynman diagrams contributing to Higgs pair production via VBF. In terms of Eq. (2), the left, middle, and right diagrams scale with , , and , respectively

Table 2

Acceptance and VBF selection cuts applied to signal and background events after jet clustering and b-tagging. The central jet veto is applied on jets with pseudo-rapidity in the interval , where and are the pseudo-rapidities of the VBF tagging jets

	14 TeV	100 TeV
Acceptance cuts
p_T_j(GeV)≥	25	40
p_T_b(GeV)≥	25	35
|η_j|≤	4.5	6.5
|η_b|≤	2.5	3.0
VBF cuts
|Δy_jj|≥	5.0	5.0
m_jj(GeV)≥	700	1000
Central jet veto: p_T_j_3(GeV)≤	45	65

In the SM, a cancellation dictated by perturbative unitarity occurs between the first and second diagrams. This is best understood by describing the process as a slow emission of the vector bosons by the protons followed by their hard scattering into a pair of Higgs bosons [73]. For generic values of and , the amplitude of the partonic scattering grows with the energy until the contribution from the new states at the cutoff scale unitarizes it. The leading contribution in the energy range comes from the scattering of longitudinal vector bosons and is given byup to and corrections. In scenarios with , the growth of the partonic cross section with thus provides a smoking-gun signature for the presence of BSM dynamics [3].

In the parametrization of Eq. (2), the amplitude for the process can be decomposed as follows:where , , and are numerical coefficients. In the present work, we will focus on the quartic coupling and set and to their SM values. This is justified for since the ATLAS and CMS measurements of Higgs production cross sections, when analyzed in the context of a global fit of Higgs properties [74-76] typically set bounds on at the level of 10–20%, depending on the specific assumptions made – see for example [77-79] and the references therein. Tighter limits on can be derived from electroweak precision tests in the absence of additional BSM contributions [80].

On the other hand, the trilinear Higgs coupling  (where ) only has loose experimental constraints so far. As an illustration, a recent ATLAS search for non-resonant Higgs pair production at 13 TeV in the final state [27] translates into the bound at the 95% confidence level. Achieving precision in the measurement of will thus most likely require the full HL-LHC statistics. Focusing on VBF production, as anticipated and further discussed in the following, gaining sensitivity to is achieved by reconstructing events with large values of . In this kinematic region, it turns out that the sensitivity to is reduced, indicating that our analysis is not optimal to probe the Higgs trilinear coupling. For these reasons, setting is a good approximation in the context of the present analysis. We can then defineand this way the total cross section will be parametrized asHowever, while setting is a very good approximation, fixing is not as equally well justified. In particular, it would be more prudent to treat as a Gaussian distributed nuisance parameter centered around its SM value with a width corresponding to the current experimental precision. To do this, a similar expression to Eq. (7) above can be derived by neglecting the subleading effects involving . In this case, Eq. (7) is replaced byWe will use this expression to evaluate the impact of on the derived bounds on at the end of Sect. 4. The values of the SM cross section and of the parameters A, B are reported in Table 1 for and 100 TeV, both after acceptance cuts and after applying all the analysis cuts as discussed in Sect. 3 – see Appendix D for the values of the parameters in bins of . We will make extensive use of this parametrization in Sect. 4 where we present our results in terms of the sensitivity on . Note that the value of A and B increase after imposing all cuts precisely because they have been optimized to enhance the sensitivity on .

Table 1

Coefficients of Eqs. (7) and (9) as obtained through a fit of Montecarlo points. The cuts are listed in Table 2 and Eqs. (15)–(17)

s	Cuts	σ_sm (fb)	A	B	C	D
14TeV	Acceptance	0.010	-5.19	29.5	-0.939	0.854
	All	0.0018	-8.18	67.5	-0.699	0.325
100TeV	Acceptance	0.20	-9.18	306	-0.699	0.584
	All	0.030	-20.7	1080	-0.516	0.251

Although we do not attempt to extract with our analysis, it is still interesting to discuss the dependence of the total cross section on this parameter. By fixing and defining , the cross section can now be parametrized asThe coefficients C and D are also reported in Table 1. As opposed to the previous case, now their values decrease after applying the full set of cuts, reflecting that the sensitivity on is suppressed by our analysis which aimed at measuring . Extracting would require retaining the events close to the hh threshold but this kinematic region is totally dominated by the background and turns out to be of little use. A measurement of the Higgs trilinear coupling in the VBF channel using the final state thus does not seem feasible even at the FCC. However, other final states might exhibit better prospects at 100 TeV.

Models

The Lagrangian of Eq. (2), with arbitrary values of the coefficients , , , and , describes a generic light scalar, singlet of the custodial symmetry, independently of its role in the electroweak symmetry breaking. In specific UV models, however, the coefficients are generally related to each other and their values are subject to constraints which depend on whether the Higgs-like boson h is part of an doublet. For example, in the SM, all the parameters in Eq. (2) are equal to 1 and terms denoted by the ellipses vanish. In this case the scalar h and the three NGBs combine to form a doublet of which is realized linearly at high energies.

Composite Higgs theories are another example where the electroweak symmetry is realized linearly in the UV, though in this case non-linearities in the Higgs interactions can be large and are controlled by the ratio , where f is the pNGB decay constant. For instance, minimal SO(5) / SO(4) models [81, 82] predictOn the other hand, the value of the Higgs trilinear coupling is not determined by the coset structure alone, and depends on how the Higgs potential is specifically generated. For instance, in the MCHM5 model with fermions transforming as vector representations of SO(5) [82], the Higgs potential is entirely generated by loops of SM fields and the Higgs trilinear coupling is predicted to beA precision model-independent determination of would thus provide stringent constraints on a number of BSM scenarios. To begin with, if the Higgs boson belongs to an electroweak doublet, as suggested by the LHC data, and the modifications to its couplings are small, then the values of and are in general predicted to be correlated [3, 5]:where . This follows because there is a single dimension-6 effective operator ( in the basis of Ref. [3]) which controls the shift in both couplings. Therefore, a high-precision measurement of can test whether the Higgs boson belongs to a doublet in case a deviation is observed in  [5].

Another interesting case is the scenario where the Higgs-like boson is not part of a doublet, and in fact does not play any role in the electroweak symmetry breaking mechanism, known as the light dilaton scenario [83-88]. In this model, invariance under dilatations implies , a condition that can be tested if the two couplings in Eq. (12) can be measured with comparable precision. Moreover, comparing and can also provide information on the coset structure in the case of a composite NGB Higgs [5].

A large variety of BSM scenarios also exists where the Higgs trilinear coupling receives large modifications while the value of the other couplings are close to the SM prediction. Higgs portal models fall in this class; see for example the discussion in [17, 36] and also [89, 90]. However, since our analysis is not sensitive to , we will not consider these scenarios any further and will always assume .

In the following, we will take as a representative benchmark scenario a model with corresponding to , with the other couplings set to their SM values, namely .

Validity range of the effective theory

As discussed above and shown by Eq. (4), the amplitude of the partonic scattering grows with the energy in the EFT described by Eq. (2). However, this behaviour holds only below the typical mass scale of new BSM states, i.e. below the cutoff scale of the effective theory. When the invariant mass of the di-Higgs system becomes large enough, the EFT approximation breaks down and it becomes necessary to take into full account the contribution from the exchange of new states, such as vector and scalar resonances. These resonances eventually tame the growth of the scattering amplitude at large energies to be consistent with perturbative unitarity bounds. In the context of composite Higgs scenarios the impact of resonances on VV scattering has been explored for instance in [91-96].

While we do not include the effect of such resonances in this work, we will report the sensitivity on as a function of the maximum value of the invariant mass of the di-Higgs system (see Fig. 13 in Sect. 4), as suggested by Ref. [97]. This comparison allows one to assess the validity of the EFT description once an estimate of is provided in terms of the masses and couplings of the UV dynamics. The result of this analysis – which is discussed in detail at the end of Sect. 4 – confirms that the EFT is valid over the full range of that is used to derive limits on . The explicit inclusion of scalar and vector resonances and their phenomenological implications for Higgs pair production via VBF is left for future work.

Fig. 13

Left the expected precision for a measurement of at the 68% CL as a function of for the LHC and the FCC. The grey area indicates the region where , obtained by setting . Right the 95% CL exclusion limits in the plane, assuming Eq. (20), where again the grey area corresponds to the non-perturbative regime, defined by . The transition between solid and dashed curves occurs at the last bin with at least one event

Analysis strategy

In this section we present our analysis of double Higgs production via VBF. First, we discuss how signal and background events are reconstructed and classified. This includes a description of the jet reconstruction techniques adopted, the b-tagging strategy, and the event categorization in terms of jet substructure. Then we illustrate the various selection cuts imposed to maximize the signal significance, in particular the VBF cuts, as well as the method used to identify the Higgs boson candidates. Finally, we present the signal and background event rates for the various steps of the analysis, and discuss how the signal cross sections are modified when is varied as compared to its SM value.

Event reconstruction and classification

Signal and background events are simulated at leading order (LO) by means of matrix-element generators and then processed through a parton shower (PS). The detailed description of the event generation of signal and backgrounds can be found in Appendix A. The dominant background is given by QCD multijet production, while other backgrounds, such as top-quark pair production and Higgs pair production via gluon fusion, are much smaller. After the parton shower, events are clustered with FastJet v3.0.1 [98] using the anti- algorithm [99] with a jet radius .

The resulting jets are processed through a b-tagging algorithm, where a jet is tagged as b-jet with probability if it contains a b-quark with GeV. In order to account for b-jet misidentification (fakes), jets which do not meet this requirement are also tagged as b-jets with probability or depending on whether they contain a c-quark or not. Only events with four or more jets, of which at least two must be b-tagged, are retained at this stage.

Fully exploiting the final state requires efficient b-tagging capabilities in both the resolved and boosted regimes as well as a good rejection of fakes. Both ATLAS and CMS have presented recent studies of their capabilities in terms of b-tagging and light-jet fake rejection for both topologies; see Refs. [100-104] and the references therein. In the present study, we have considered two representative b-tagging working points:The first point is consistent with the current performance of ATLAS and CMS, and is the one adopted as baseline in this paper. The second working point is more optimistic and is intended to assess how much one could gain with a more efficient b-tagger. As our results will show (see Fig. 14), using WP2 leads to a marginal improvement in our analysis. For simplicity, we applied the efficiencies in Eq. (13) to a jet based on its constituents. This is sufficient for the purpose of the current analysis, which is namely to demonstrate the sensitivity of double Higgs production via VBF to . Accordingly, we leave a detailed study of b-tagging including hadronization effects and dependence to future studies.

Fig. 14

Signal significance in the scenario, as a function of the di-Higgs invariant mass at the HL-LHC and the FCC, for the two b-tagging working points of Eq. (13)

Subsequently to b-tagging, events are classified through a scale-invariant tagging procedure [10, 43]. This step is crucial to efficiently reconstruct the Higgs boson candidates and suppress the otherwise overwhelming QCD backgrounds while at the same time taking into account all the relevant final-state topologies. The basic idea of this method is to robustly merge three event topologies – boosted, intermediate and resolved – into a common analysis. This is particularly relevant for our study, since, as discussed in Sect. 2, the degree of boost of the di-Higgs system strongly depends on the deviations of from its SM value.

This scale-invariant tagging strategy is schematically represented in Fig. 2. First of all, the b-tagged jets are ordered in and the constituents of the hardest two jets are then re-clustered using the Cambridge/Aachen (C/A) algorithm [105] with . Each C/A jet is processed with the BDRS mass-drop (MD) tagger [106]. This jet-substructure tagger has two parameters: and . And, in this work, we set and as in the original BDRS study. To determine if a given jet arises from the decay of a massive object, the last step of the clustering for jet j is undone, giving two subjets and which are ordered such that . Then, if the two subjets satisfy the conditionswhere is the angular separation between the two subjets, j is tagged as a jet with a mass drop. Else, the procedure is applied recursively to until a mass drop is found or the C/A jet is fully unclustered.

Fig. 2

Schematic representation of the analysis strategy adopted in this work

Jets are mass-drop tagged only if they satisfy the following additional requirement: at least two b-quarks must be contained within the jet, each of which with GeV, and with a minimal angular separation . The request of a second b-quark completes our b-tagging algorithm in the case of boosted jets. Other more sophisticated approaches to b-tagging could have been considered – e.g., using ghost-association between large-R MD-tagged jets and small-R b-tagged jets [23, 43], or accounting for an efficiency which depends on the jet  [107]. The approach followed here is at the same time simple yet realistic enough for a first feasibility study with the caveat that a more complete analysis should treat b-tagging more in line with the actual performance of the ATLAS and CMS detectors (and in particular it should include a full detector simulation).

The use of the BDRS mass-drop tagger allows us to classify a given signal or background event under one of the three categories: boosted, if two mass-drop tags are present; intermediate, for an event with a single mass-drop tag; and resolved, if the event has no mass-drop tags. In the resolved category, events are only retained if they contain at least four b-tagged jets, while at least two b-tagged jets, in addition to the MD-tagged jet, are required in the intermediate one. By construction, this classification is exclusive, i.e., each event is unambiguously assigned to one of the three categories. This exclusivity allows the consistent combination of the signal significance from the three separate categories.

Following the event categorization, acceptance cuts to match detector coverage are applied to signal and background events. These cuts are listed in the upper part of Table 2, and have been separately optimized for the LHC TeV and the FCC 100 TeV. We require the of the light (b-tagged) jets to be larger than 25 GeV (25 GeV) at 14 TeV and than 40 GeV (35 GeV) at 100 TeV, respectively. Concerning the pseudo-rapidities of light and b-tagged jets, and , at the LHC the former is limited by the coverage of the forward calorimeters, while the latter is constrained by the tracking region where b-tagging can be applied. At 100 TeV, we assume a detector with extended coverage of the forward region up to of 6.5 [108].

VBF selection cuts

Subsequently to the acceptance cuts, we impose a set of selection cuts tailored to the VBF topology which is characterized by two forward and very energetic jets with little hadronic activity between them. In particular, we cut on the rapidity separation and the invariant mass of the two VBF tagging jets, and impose a central jet veto (CJV) on the hardest non-VBF light jet in the central region. The VBF tagging jets are defined as the pair of light jets satisfying the acceptance cuts of Table 2 with the largest invariant mass . This definition is robust with respect to soft contamination from the underlying event (UE) and pile-up (PU) and to the contribution of b-jets mistagged as light jets.

Figure 3 shows the distribution of the rapidity separation and invariant mass of the VBF tagging jets at 14 and 100 TeV after the acceptance cuts. In each case, we show the results for the signal (SM and benchmark) and for the total background. The signal distributions exhibit the distinctive VBF topology, with two VBF tagging jets widely separated in rapidity and with a large invariant mass. This is in contrast with the backgrounds where both the and distributions peak at zero. In Fig. 3, as well as in the subsequent figures, kinematic distributions have been area-normalized and then rescaled by a common factor such that the largest bin in the plot is of unit height.

Fig. 3

Distribution of the rapidity separation (upper) and the invariant mass (bottom panels) of the VBF tagging jets at 14 TeV (left) and 100 TeV (right panels), for signal (SM and ) and background events after the acceptance cuts The vertical line indicates the value of the corresponding cut from Table 2. The distributions have been area-normalized and rescaled by a common factor

Based on the distributions of Fig. 3 we identified appropriate values of the VBF cuts, listed in Table 2 and represented in each panel by a vertical dash-dotted line. It is important to tailor these cuts to the specific centre-of-mass energy, 14 and 100 TeV, to avoid losing a substantial fraction of the signal events. One should also take into account that the large rapidity separation between the VBF tagging jets in signal events results from jets pairs with a large invariant mass, given that these two variables are strongly correlated [4]. This large separation in rapidity is especially useful in the final state to trigger on signal events, providing a significant improvement compared to the same final state produced in gluon fusion where triggering issues are more severe [42, 43].

Figure 3 clearly highlights that in order to maximize the acceptance of events with VBF topology – the detectors must have a good coverage of the forward region. This issue is particularly relevant at 100 TeV as illustrated in Fig. 4 which shows the pseudo-rapidity distribution of the most forward light jet. At 100 TeV, this peaks at around 5, so a detector instrumented only up to would lose more than 50% of signal events. The discontinuity in the FCC case delineates the edge of the b-jet acceptance region, , above which no b-tagging is attempted and b-jets contribute to the light-jet yield.

Fig. 4

Distribution of the pseudo-rapidity of the most forward light jet at 14 and 100 TeV. Both curves have a discontinuity at the edge of the corresponding b-tagging region which is delineated by the solid (dashed) grey vertical lines in the case of 14 (100) TeV. This discontinuity is more clearly visible in the 100 TeV curve and is purely due to combinatorics

Turning to the transverse momentum of the light jets, Fig. 5 shows the distributions of the three hardest jets at 14 and 100 TeV for SM signal events. One can see that while the leading jet is typically quite hard, the subleading ones are rather soft. It is thus important to avoid imposing a too stringent cut in , in order not to suppress the signal. Fortunately, in contrast from the gluon-fusion process, adopting a soft cut is not a problem since triggering can be performed based on the VBF topology. Comparing the distributions at 14 and 100 TeV, their shapes turn out to be rather similar, shifted towards larger values at TeV. This justifies the harder cut in this case (see Table 2), also required to reduce the contamination from UE and PU.

Fig. 5

Distributions of the of the leading, subleading, and third light jet at 14 TeV (left panel) and 100 TeV (right panel) for the SM signal. Distributions are area-normalized as in Fig. 3

As mentioned above, another characteristic feature of VBF production is a reduced hadronic activity in the central region between the two VBF tagging jets. This follows because the latter are not colour-connected since the production of the central system only involves electroweak bosons. For this reason, a CJV cut is commonly imposed in VBF analyzes. This cut vetoes light jets, with pseudo-rapidity , lying between those of the VBF tagging jets, , above a given threshold.

The effect of the CJV is illustrated in Fig. 6 where we show the distribution of the of the third light jet, , for the SM signal and the total background. Although the latter has a harder spectrum than the signal, imposing too stringent a veto is not advantageous. This is because the final state leads to a non-negligible amount of hadronic activity in the central region, for instance due to gluon radiation from the b quarks and to b-jet misidentification. Based on these results, in our analysis, we impose a CJV with the threshold value reported in Table 2 and shown in the plots by the dot-dashed line.

Fig. 6

Distribution of the of the third light jet at 14 TeV (left panel) and 100 TeV (right panel) for the SM signal and the total background, including only events where this jet lies within the pseudo-rapidity region between the VBF jets. The vertical line indicates the CJV cut

Higgs reconstruction

The next step in our analysis is the reconstruction of the Higgs boson candidates. This is done separately for each of the three event categories. In the resolved category, starting with the six hardest b-jets in the event,1 we reconstruct the first Higgs boson candidate by identifying it with the pair of b-jets whose invariant mass is closest to the Higgs mass, GeV. Out of the remaining b-jet pairs, the one with an invariant mass closest to is then assigned to be the second Higgs boson candidate, . In the case of the intermediate and boosted categories, each of the mass-drop tagged jets is identified with a Higgs candidate. The second Higgs candidate in the intermediate category is then formed by considering the five hardest b-jets in the event and selecting the pair whose mass is closest to .

The invariant mass distributions of the Higgs candidates for the signal (SM and ) and the total background are shown in Fig. 7. The peak around GeV is clearly visible for signal events, especially in the case of . The smearing of the signal distribution of the second Higgs candidate arises from out-of-cone radiation effects which reduce the reconstructed mass. It is largest in the SM, while it is reduced in the scenario and in particular at TeV, due to the larger boost of the Higgs bosons. The small peak in the background distributions for is artificially sculpted by the analysis selection cuts. The fact that the efficiency for the reconstruction of the Higgs bosons is similar in the SM and for the benchmark is another validation of the scale-invariant tagging, since while in the SM most events lead to resolved topologies, the scenario is dominated by the boosted category (see Fig. 8 below).

Fig. 7

Invariant mass distribution of the leading () and subleading () Higgs candidates for signal (SM and ) and background events at 14 TeV (left) and 100 TeV (right)

Fig. 8

Invariant mass distribution of the di-Higgs system at 14 TeV (left) and 100 TeV (right) after all analysis cuts, for the signal (SM and ) and the total background. We show the contribution from resolved and boosted events as well as the sum of the three categorieslabelfig

After reconstructing the Higgs candidates, we require that their invariant masses, and , are reasonably close to the nominal mass. In the resolved category, these conditions are The mass window in Eq. (16) is centered around the mass of the first candidate rather than the nominal Higgs mass in order to make the cut robust against UE and PU effects. For the intermediate and boosted categories, the invariant mass of all Higgs candidates is required to satisfy Eq. (15). The cuts of Eqs. (15) and (16) are especially effective in suppressing the QCD backgrounds which have almost featureless distributions in the Higgs mass regions. Finally, we impose an additional cut on the invariant mass of the di-Higgs system:This condition greatly reduces the background rates while leaving the interesting kinematic region at large – where deviations from the SM signal mostly appear – unaffected.

Figure 8 shows the distribution, after all the cuts listed in Table 2 and Eqs. (15)–(17), for signal (SM and ) and the total background at 14 and 100 TeV. In each case, we show both the sum of the three event categories and the individual contributions from resolved and boosted events. The intermediate category, which contributes very little in all cases, is not shown. This comparison helps illustrate the relative weight of the boosted and resolved categories. For signal events in the SM, the vast majority are classified in the resolved category as expected since in this case the boost of the di-Higgs system is small except at 100 TeV and for large values. On the other hand, in the case of , the energy growth of the partonic cross section induces a much harder spectrum. This implies that, already at 14 TeV, a substantial fraction of events falls in the boosted category which becomes the dominant one at 100 TeV. For , the crossover between the resolved and boosted categories takes place at TeV for both colliders, although this specific value depends on the choice of the jet radius R [10]. Unsurprisingly, background events are always dominated by the resolved topology.

Signal and background event rates

Now that we have presented our analysis strategy we can discuss the actual impact on the cross sections and event rates of the various steps of the cut flow. In Table 3 we report the cross sections at 14 and 100 TeV after acceptance, VBF, Higgs reconstruction, and cuts of Table 2 and Eqs. (15)–(17), respectively, for both the signal (SM and ) and for the total background.

Table 3

Cross sections, in fb, at TeV (upper table) and TeV (lower table) after the successive application of the acceptance and VBF cuts (Table 2) and of the Higgs reconstruction cuts (Eqs. (15)–(17)), for signal events (SM and ) and for the total background

	Cross sections (fb)
	Acceptance	VBF	Higgs reco.	m_hh cut
14 TeV
Signal SM	0.011	0.0061	0.0039	0.0020
Signal c_2V=0.8	0.035	0.020	0.017	0.011
Bkgd (total)	1.3×10^5	4.9×10^3	569	47
100 TeV
Signal SM	0.22	0.15	0.11	0.033
Signal c_2V=0.8	3.4	2.7	1.9	1.6
Bkgd (total)	1.9×10^6	1.9×10^5	9.5×10^3	212

At TeV, we find that the VBF di-Higgs signal in the SM is rather small already after the basic acceptance cuts. On the other hand, the signal event yield is substantially increased for as illustrated by the benchmark value of leading to more than a factor 3 (5) enhancement compared to the SM after the acceptance (all analysis) cuts. The fact that this cross-section enhancement for the scenario is more marked at the end of the analysis is not a coincidence: our selection cuts have been designed so as to improve the sensitivity to by increasing the signal significance in the large- region.

From Table 3 we also find that a similar qualitative picture holds at TeV with the important difference that, in this case, the event rate is already substantial in the SM which yields events after the acceptance cuts with ab. The cross section enhancement at TeV as compared to TeV is driven by the larger centre-of-mass energy and leads to a signal rate greater by a factor 20 (17) after the acceptance (all analysis) cuts in the SM, and by a factor in the scenario. At TeV, the ratio of signal cross sections in the and SM scenarios is after acceptance (all analysis) cuts. Note however, that at both 14 and 100 TeV, even after all analysis cuts the background is still much larger than the signal (either SM or ) at the level of inclusive rates. It is only by exploiting the large- region that the former can be made small enough to achieve high signal significances.

Table 3 is also useful to assess the relative impact on the signal and the total background of each of the cuts imposed. In the case of the VBF cuts, we find that the background is drastically reduced, by more than one order of magnitude, at the cost of a moderate decrease of the signal cross sections. The Higgs mass window requirement is also instrumental to further suppress the backgrounds, especially the QCD multijets which are featureless in , while leaving the signal mostly unaffected. A final reduction of the background, by around another order of magnitude, is achieved through the cut. The relative impact of each cut on signal events is similar in the SM and for .

Figure 9 graphically illustrates the dependence of the di-Higgs production cross section on the couplings and . The left panel shows the cross section in SM units as a function of and after applying the acceptance cuts of Table 2 (dashed curves) and after all the analysis cuts (solid curves). The sensitivity on is particularly striking, for example the cross section for is enhanced by a factor compared to its SM value after all cuts. This sensitivity is the key ingredient for measuring with good precision, even though the SM cross section itself cannot be extracted with comparable accuracy. In particular, as we will show in Sect. 4, the sensitivity to derives mainly from the tail of the distribution. This observation elucidates the enhancement (suppression) in the sensitivity to in Fig. 9 after the application of all the cuts which, for instance, remove the threshold region up to GeV at the LHC(FCC).

Fig. 9

Left panel VBF di-Higgs cross section, in units of the SM value, as a function of (thick blue) and (thin red), after acceptance cuts (solid) and all analysis cuts (dashed). Right panel ratio of VBF di-Higgs cross section between 100 and 14 TeV as a function of

The right panel of Fig. 9 shows, instead, the ratio between the VBF di-Higgs cross sections at and 14 TeV. Given the larger centre-of-mass energy of the FCC, it is expected that this ratio grows rapidly for and, indeed, it can reach values as high as 300 for . As will be demonstrated in Sect. 4, this effect allows for much more precise measurements of at the FCC, with uncertainties reduced by a factor 20 as compared to the HL-LHC. The results of Fig. 9 are of course consistent with the findings of Table 3.

From Fig. 9, we also observe that the sensitivity of the signal on the Higgs trilinear coupling is relatively weak even for large variations and it is reduced by our analysis strategy. As already mentioned, the last feature is expected because the sensitivity to  comes from events near the di-Higgs threshold, , which are removed by our cuts due to the overwhelming backgrounds in that region. This weak dependence of the VBF di-Higgs cross section on , together with the large event rates for background after all the analysis cuts (see Table 3), suggest that the VBF process is not suitable to extract the Higgs self-coupling.

Let us now discuss the background processes. As mentioned above and discussed in Appendix A, there are two types of processes that contribute to the final-state signature under consideration. The first type are QCD processes and in particular multijet and top-quark pair production in association with additional hard radiation. The second is Higgs pair production in the gluon–gluon fusion channel in association with additional jets, where the latter can mimic the VBF topology, as in single-Higgs production.

In the case of QCD multijet processes, it is important to account for the effects of both the 4b and the 2b2j backgrounds (where we label each process by its matrix-element level content; as explained in Appendix A, additional jets are generated by the parton shower). The latter process can lead to events being classified as signal when light jets are misidentified as b-jets or when a gluon splits into a pair during the parton shower. Even with a small light-jet mistag rate of , it can have a contribution to the total background comparable to or bigger than the 4b process. Details of the generation of the QCD backgrounds and on the associated validation tests are presented in Appendix A and Appendix C.

Concerning gluon-fusion Higgs pair production in association with additional hard jets, similarly to single-Higgs VBF production there will be certain configurations that mimic the VBF topology as emphasized, for example, in Ref. [64]. In contrast to the VBF channel, however, the Higgs pair production in gluon fusion does not exhibit any enhancement in the tail of the distribution. This substantially reduces its contamination to the region with the highest sensitivity to in our analysis. Note also that a harder distribution could be generated by higher-order EFT operators, for instance those leading to a contact interaction of the form gghh, as in Ref. [17]. The investigation of this scenario is, however, left for future work.

In Table 4, following the structure of Table 3, we give the cross sections at 14 and 100 TeV for the individual background processes (and their sum) after the acceptance, VBF, Higgs reconstruction, and cuts. We find that in all steps in the cut-flow the dominant background component is QCD multijet production, both at 14 and 100 TeV. After all analysis cuts, the 2b2j component is a factor 10 larger than 4b at 14 TeV while they are of similar size at 100 TeV.

Table 4

Same as Table 3, now listing separately each background process

	Acceptance	VBF	Higgs reco.	m_hh cut
LHC 14 TeV
4b	1.18×10^4	613	54	4.45
2b2j	1.14×10^5	4.31×10^3	514	42.6
tt¯jj	150	4.75	0.732	0.0706
gg→hh	0.98	0.0388	0.0223	0.00857
Total	1.3×10^5	4.9×10^3	569	47
FCC 100 TeV
4b	3.93×10^5	4.59×10^4	2.61×10^3	106
2b2j	1.52×10^6	1.46×10^5	6.88×10^3	104
tt¯jj	9.76×10^3	832	55	1.47
gg→hh	24.8	2.48	1.31	0.0892
Total	1.9×10^6	1.9×10^5	9.5×10^3	212

Other backgrounds, including gluon-fusion di-Higgs production, are much smaller than QCD multijets. Note however, that the former is actually larger than the VBF signal for SM couplings with a cross section at 14 TeV of fb after acceptance (all) cuts, compared to fb for the VBF case. On the other hand, this fact does not affect the measurement , since, as we show next, the sensitivity comes from the large tail where the gluon-fusion component is heavily suppressed.

The decomposition of the total background in terms of individual processes as a function of is shown in Fig. 10 where the components are stacked on top of each other so that the content of each bin matches the total background cross section. In the TeV case, the 4b background dominates for large while the 2b2j one is instead the most important for small . The TeV case is similar with one exception, namely that the gluon-fusion di-Higgs background becomes the dominant one for very high invariant masses, TeV. Such an extreme region, however, is phenomenologically irrelevant due to the very small rates of both signal and background even at a 100 TeV collider.

Fig. 10

Decomposition of the total background into individual processes as a function the di-Higgs invariant mass after all analysis cuts have been imposed, except for the cut

Results

The last column of Table 3 indicates the cross sections for the signal and total background after imposing all analysis cuts. We observe that the background, dominated by QCD multijets, still has a much larger cross section than the signal, both in the SM and in the benchmark scenario. As anticipated, the additional handle which we can now exploit to increase the signal significance is the different behaviour of the distribution for the signal and the background, in particular when departs from its SM value. The latter has a sharp fall-off at large values, while, instead, the signal exhibits a much harder spectrum for . This cross-section growth implies that, for sufficiently large, there will be a crossover value of where the signal overcomes the background.

This behaviour is illustrated in Fig. 11 where we show the invariant mass distribution of the Higgs pairs after all analysis cuts, at 14 and 100 TeV, for the signal (SM and ) and the total background. In the case of the benchmark scenario with , the crossover between signal and background is located at TeV (TeV) at TeV (TeV). We also observe that, for invariant masses above this crossover, the ratio between the signal and the backgrounds keeps increasing steeply.

Fig. 11

The di-Higgs distribution at 14 TeV (left) and 100 TeV (right) after all analysis cuts showing the results for the signal (SM and ) and for the total background

With the final results of our analysis in hand, we can now estimate the expected sensitivity to deviations in the hhVV coupling, parametrized as , by exploiting the information contained in the full differential distribution (as opposed to using only the total number of events satisfying all cuts from Table 3). To achieve this, we first bin our results in and then follow a Bayesian approach [109] to construct a posterior probability density function. We include two nuisance parameters, and , to account for the uncertainty associated with the background and signal event rate, respectively. The parameter encodes the theoretical uncertainties on the di-Higgs cross section and the branching fraction . We conservatively assume a uncertainty uncorrelated in each bin.

Concerning , we expect that an actual experimental analysis of di-Higgs production via VBF would estimate the overall normalization of the different background components by means of data-driven techniques. We assume a 15% uncertainty arising from the measurement and subsequent extrapolation of the dominant QCD multijet background; see for example a recent ATLAS measurement of dijet cross sections [110]. The background nuisance parameter, , is conservatively also assumed to be uncorrelated among bins. In addition, while we already rescale the background cross sections to match existing NLO and NNLO results (see Appendix A), there still remains a sizeable uncertainty in their overall normalization from missing higher orders, in particular for the QCD multijet components. For this reason, below, we explore the robustness of our results upon an overall rescaling of all the background cross sections by a fixed factor.

The posterior probability function constructed in this way readswith and denoting, respectively, the number of predicted (for a generic value of ) and observed (assuming SM couplings) events for a given integrated luminosity in the ith bin of the di-Higgs invariant mass distribution , given by2:In Eq. (19), and indicate the signal (for a given value of ) and total background cross sections, respectively, for the i-th bin of the distribution. The functional form of is given by Eq. (7) and the value of the coefficients in bin i are given in Appendix D. We denote by the prior probability distribution of the coupling.

As justified above, in the evaluation of Eq. (18) we set and assume that the two nuisance parameters are normally distributed. We have verified that assuming instead a log normal distribution leads to similar results. In addition, we take a Poissonian likelihood in each bin and assume the prior probability to be uniform. The resulting posterior probabilities are shown in Fig. 12 for the LHC with fb (LHC) and ab (HL-LHC), and for the FCC with ab. To produce this figure, as well as to determine the values reported in Tables 5 and 6, we included all bins with at least one event.

Fig. 12

Posterior probability densities for at the LHC for fb (LHC) and ab (HL-LHC) and for the FCC with ab. Note the different scales of the axes in the two panels

Table 5

Expected precision (at 68% probability level) for the measurement of at the LHC and the FCC, assuming SM values of the Higgs couplings. We show results both for the nominal background cross section , and for the case in which this value is rescaled by a factor 3

	68% probability interval on δ_c_2V
	1×σ_bkg	3×σ_bkg
LHC_14	[-0.37, 0.45]	[-0.43, 0.48]
HL-LHC	[-0.15, 0.19]	[-0.18, 0.20]
FCC_100	[0, 0.01]	[-0.01, 0.01]

Table 6

95% probability upper limits on the fiducial signal strength,

	95% probability upper limit on μ
	1×σ_bkg	3×σ_bkg
LHC_14	109	210
HL-LHC	49	108
FCC_100	12	23

From Fig. 12, we can determine the expected precision for a measurement of at the LHC and the FCC in the case of SM values of the Higgs couplings. The 68% probability intervals for the determination of at the LHC and the FCC are listed in Table 5. This is the central result of this work. To assess its robustness with respect to our estimate of the background cross sections, we also provide the same intervals in the case of an overall rescaling of the total background by a factor 3. Furthermore, we can also assess the effect of varying on the bound on by treating as a nuisance parameter and marginalizing over it. The leading effect of varying comes from the term at the amplitude level – see Eq. (4) – and can be included using the parametrization of Eq. (8). The neglected dependence is subleading and arises from the interference of diagrams proportional to and . We take to be Gaussian distributed with a mean equal to 1 (i.e., its SM value) and a width equal to 4.3, 3.3, and 2% at the LHC Run II, HL-LHC, and FCC respectively. In case of the LHC (both Run II and HL), the width of the Gaussian corresponds to the projected sensitivity from the two parameter fit by ATLAS [111]. The effect of marginalizing over is subleading in both LHC scenarios and weakens the bound on . We find that the results of Table 5 change by 2% for LHC and 7% for HL-LHC. The effect at the FCC is much larger causing the bound on to be rather than 0.01. This is not surprising and indicates that a joint likelihood would be required at the FCC.

From Table 5 we find that the coupling, for which there are currently no direct experimental constraints, can already be measured at the LHC with with a reasonably good accuracy: with 68% probability. This accuracy is only marginally degraded if the background is increased by a factor 3. A better precision, of the order of , is expected at the HL-LHC with . This estimate is robust against an overall rescaling of the background cross section. Finally, we find a very significant improvement at the FCC with , where a measurement at the 1% level could be achieved, providing an unprecedented test for our understanding of the Higgs sector.

It is interesting to compare these results with the experimental precision expected on the fiducial VBF di-Higgs production cross section after all analysis cuts, expressed in terms of the signal strength parameter normalized to the SM result, . Table 6 shows the 95% upper limits on for the nominal background cross section and after rescaling the latter by a factor 3. The comparison with Table 5 clearly shows that the high precision expected on can obtained despite the rather loose constraints that can be obtained on the VBF di-Higgs cross section even at TeV. As already discussed, this behaviour follows from the strong dependence of the signal cross section on ; see Fig. 9.

The results of Tables 5 and 6 have been obtained by making full use of the information contained on the di-Higgs invariant mass distribution . However, the EFT expansion might break down at large enough values of , corresponding to large partonic centre-of-mass energies, and some assessment on the validity of our procedure is thus required. In particular, results can be consistently derived within the EFT framework only if the new physics scale is smaller than the largest value of included in the analysis.

As stressed in Ref. [97], constraining requires making assumptions on the structure of the UV dynamics extending the SM. For example, for the case where the new physics is characterized by a single coupling strength  and mass scale  [3], one naively expectsTherefore, for maximally strongly coupled UV completions (with ) it is possible to derive the following upper limit,which makes explicit the connection between the value of and the new physics scale . The validity of the EFT can thus be monitored by introducing a restriction on the bins used in the construction of the posterior probability Eq. (18), so that , and then determining how the sensitivity on varies as a function of  [97].

The precision on , defined though the symmetrized 68% probability interval, is shown in Fig. 13 as a function of for the LHC and the FCC. As expected, increasing , i.e. making the cut less stringent, leads to stronger constraints. Eventually, flattens out when is large enough that all the bins which contain at least one event are included in the posterior probability of Eq. (18). Bins which contain at least one event are depicted in Fig. 13 with a solid curve, while those containing less than one event are depicted by a dashed curve. The grey area in Fig. 13 corresponds to the non-perturbative region where , obtained by setting in Eq. (21), the most optimistic assumption compatible with the validity of the EFT expansion.

As an additional way to quantify the validity of the EFT approach in our analysis, we derive the region of exclusion in the plane  [97], corresponding to the limits on derived as a function of . This is shown in the left panel of Fig. 13, making use of Eq. (20) and then setting . The grey area in the upper part of the plot indicates the non-perturbative region, defined by . We find that the dominant constraints on arise from a region in the parameter space where the EFT expansion is valid, both at the LHC and at the FCC. The results from the two panels of Fig. 13 indicate that our EFT analysis is robust and that it can be used to derive stringent bounds on in the absence of new explicit degrees of freedom.

Finally, in Fig. 14 we show the signal significance, , in the = 0.8 scenario as a function of the di-Higgs invariant mass at the HL-LHC and the FCC. The results are presented for the two b-tagging working points defined in Eq. (13). As already discussed, these have been chosen so that the first (WP1) is consistent with the current ATLAS and CMS performances, while the second (WP2) assumes an improved detector performance. One can observe that the signal significance of each individual bin is at most at the HL-LHC (though the precise numbers depend on the specific choice of binning), while at the FCC one finds much higher signal significances, with already for TeV and then increasing very rapidly for higher values of . Figure 14 clearly shows that the signal significance depends very mildly on the specific details of the b-tagging performance and that operating at WP2 instead of WP1 implies only a minor improvement.

To summarize, we have demonstrated how Higgs boson pair production via VBF can be used to provide the first direct constraints on the coupling already at the LHC with fb (Table 5), which at a TeV collider would become a high-precision measurement with potentially sub-percent accuracy. We have also assessed (Fig. 13) the robustness of our strategy and the validity of the underlying EFT expansion. Our analysis clearly highlights the unique physics potential of extending current di-Higgs searches at the LHC to the vector-boson fusion channel.

Conclusions and outlook

The measurement and study of Higgs pair production is one of the cornerstones of the LHC program as well as of any future hadron collider. It provides unique information on the Higgs sector and on the mechanism underlying electroweak symmetry breaking, and allows a direct test of the strength of the Higgs boson self-interactions. On the other hand, it is a challenging measurement and the low production rates require large integrated luminosities to achieve reasonable signal significances. While most studies of Higgs pair production so far have concentrated on the gluon-fusion channel which has the largest cross section, we have shown in this work how the vector-boson fusion channel can impose stringent constraints on Higgs couplings that are not directly accessible by other means, in particular on the hhVV quartic coupling .

Exploiting the high signal yield of the final state, we have presented a detailed feasibility study of the measurement of Higgs boson pairs in the vector-boson fusion channel at the LHC and at a future 100 TeV hadron collider. A key ingredient of our strategy is provided by the fact that deviations of the Higgs couplings to vector bosons from the parabola significantly harden the distribution resulting in a large fraction of events with a boosted Higgs pair. The subsequent decays into pairs can then be reconstructed by means of jet substructure techniques. While QCD backgrounds are very large, we have shown how the combination of selection cuts exploiting the VBF topology and the growth of the distribution when the Higgs couplings depart from their SM values leads to a remarkable model-independent sensitivity to the coupling.

Our results demonstrate that at the LHC with an integrated luminosity of (3000) fb the hhVV coupling can be measured with () precision at the 68% probability level, reaching 1% accuracy at a TeV collider. Therefore, stringent constraints on this so far unknown coupling can be obtained already before the start of the HL-LHC data taking. Our analysis provides strong motivation for the ATLAS and CMS collaborations to extend their searches for Higgs pair production to the VBF channel already during Runs II and III. On the other hand, we also find that the VBF channel is clearly inferior to the gluon-fusion channel for a measurement of the Higgs self-coupling , at least for the final state studied here, since the sensitivity to arises from the threshold region where QCD backgrounds overwhelm the signal even for sizeable modifications of the Higgs couplings with respect to their SM values.

There are several possible avenues for future work. On one hand, it might be interesting to study the possibility to enhance the sensitivity to by means of a multivariate analysis (MVA), such as those used in [43], in order to dynamically determine the optimal set of selection cuts and optimize the discrimination between signal and background events. Further, it should be possible to quantify the constraints on additional EFT operators that can contribute to the di-Higgs VBF signal yield and that have not been considered in this work. Finally, a complete analysis should include a full detector simulation, especially for the reconstruction of the forward jets and of the Higgs boson candidates, and a b-tagging strategy able to reproduce more closely the one adopted by the LHC experiments.

Table 7

Parton-level generation cuts for the VBF di-Higgs signal and the QCD background samples at and 100 TeV. In this table, j refers to light quarks and gluons, b to bottom quarks, and is the invariant mass of the two light partons in the event with largest

	14 TeV		100 TeV
	Signal	QCD bkg.	Signal	QCD bkg.
p_T_j(GeV)≥	25	20	25	20
p_T_b(GeV)≥	25	20	25	20
|η_j|≤	4.5	3	10	3
|η_b|≤	2.5	3	3	3
ΔR_bb≥	0.1	0.1	0.1	0.1
ΔR_jb≥	0.4	0.1	0.4	0.1
ΔR_jj≥	4	0.1	4	0.1
m_jj(GeV)≥	600	–	600	–

Table 8

Parton-level cross sections for VBF Higgs pair production in the final state at 14 and 100 TeV (including branching ratios) after the acceptance cuts of Table 7.We show the cross sections for different values of the couplings normalized to the SM values. We also provide the number of events at the HL-LHC for and at the FCC for

		Signal: VBF hh→bb¯bb¯
		LHC 14 TeV		FCC 100 TeV
{c_V,c_2V,c_3}		σ (fb)	N_ev(L=3ab^-1)	σ (fb)	N_ev(L=10ab^-1)
{1,1,1}	SM	0.26	780	15	1.5×10^5
{1,0,1}		4.4	1.3×10^4	593	5.9×10^6
{1,2,1}		2.5	7.5×10^3	471	4.7×10^6
{1,0,0}		5.8	1.7×10^4	656	6.6×10^6
{1,0,-1}		7.5	2.3×10^4	731	7.3×10^6
{1,1,0}		0.64	1.9×10^3	30	3.0×10^5
{1,0.8,1}	Benchmark	0.58	1740	48	4.8×10^5

Table 9

Parton-level cross sections for the background processes considered in this work after the acceptance cuts of Table 7. We indicate the programs used, the numbers of events generated and whether they are weighted [W] or unweighted [UW], the LO cross sections and the corresponding K-factors and the relevant reference. In the case of the 4b2j and 2b4j backgrounds, the same K-factors as those for the 4b and 2b2j, respectively, were applied. For the 4b2j sample, the 6M(2M) events correspond to those generated at the LHC(FCC) centre-of-mass energies, respectively. The generation cuts used for the ALPGEN  samples are different; see text

Background processes
			σ_LO (fb)		K-factor		Refs.
Process	Program	N_ev	LHC14	FCC100	LHC14	FCC100
4b	Sherpa2.2	50M [W]	1.1×10^6	1.6×10^7	1.7	1.7	[112]
2b2j	Sherpa2.2	50M [W]	2.6×10^8	3.8×10^9	1.3	1.3	[112]
tt¯jj	Sherpa2.2	10M [W]	1.9×10^4	1.6×10^6	1.6	1.6	[122]
4b2j	ALPGEN	6M(2M) [UW]	5.4×10^4	2.4×10^6	1.7	1.7	–
2b4j	ALPGEN	260k [UW]	10^7	5.2×10^8	1.3	1.3	–
gg→hh→bb¯bb¯	aMC@NLO	1M [UW]	6.2	272	2.4	2.2	[47]

Table 10

The fit coefficients in Eq. (B.1) for each background process at 14 and 100 TeV

s (TeV)	Background	a	b	c
14	4b	0.492	14.2	0.512
	2b2j	-2.65	24.9	0.198
	ttjj	-18.3	48	-3.63
	ggf	-8.27	30.1	-1.82
100	4b	-4.82	32.3	-0.539
	2b2j	-2.97	36.8	-0.118
	ttjj	-3.31	38.4	-0.265
	gg→hh	3.61	15.9	0.695

Table 11

The bin by bin fit coefficients obtained by fitting MC events to Eq. (7) for the LHC with TeV. The first column labelled ‘Bin’ gives the bin number in question. The bin definition is given in the text; see also footnote 2. The last three columns labelled , , and give the coefficients of the correlation matrix among the three fit parameters

Bin	σ_sm (fb)	A	B	ρ_σA	ρ_σB	ρ_AB
1	(3.08±0.05)×10^-4	-3.61±0.0957	6.89±0.217	0.3	-0.62	-0.66
2	(5.84±0.0691)×10^-4	-3.96±0.0706	7.27±0.161	0.32	-0.62	-0.68
3	(7.17±0.0654)×10^-4	-4.3±0.0573	9.71±0.148	0.32	-0.68	-0.63
4	(7.31±0.049)×10^-4	-4.96±0.0462	14.9±0.146	0.35	-0.77	-0.59
5	(5.98±0.0507)×10^-4	-6.39±0.0677	26.7±0.287	0.42	-0.86	-0.59
6	(4.19±0.0677)×10^-4	-8.28±0.157	50.1±0.93	0.5	-0.92	-0.6
7	(2.38±0.0571)×10^-4	-11.9±0.302	103±2.65	0.61	-0.96	-0.67
8	(1.15±0.0447)×10^-4	-15.1±0.597	174±7.05	0.67	-0.98	-0.71
9	(4.93±0.176)×10^-5	-19.4±0.712	301±11	0.72	-0.99	-0.75
10	(1.68±0.148)×10^-5	-29.1±2.65	758±67.8	0.8	-1	-0.81
11	(5.41±1.02)×10^-6	-42.3±8.58	1.69×10^3±320	0.83	-1	-0.83
12	(1.28±0.506)×10^-6	-47.4±22.6	(3.88±1.53)×10^3	0.76	-1	-0.77
13	(8.55±13.4)×10^-7	-16.5±29.7	(2.27±3.57)×10^3	0.76	-1	-0.77
14	(3.5±3.57)×10^-7	-0.00901±26.4	(1.28±1.31)×10^3	-0.064	-1	0.055
15	(6.24±3.38)×10^-7	-3.59±5.3	105±60.7	0.19	-0.98	-0.24
16	(6.33±1.97)×10^-7	-2.88±1.57	22.1±8.59	0.15	-0.9	-0.36

Table 12

Same as Table 11 but for the FCC with TeV

Bin	σ_sm (fb)	A	B	ρ_σA	ρ_σB	ρ_AB
1	(5.98±0.152)×10^-3	-2.87±0.191	7.8±0.513	0.14	-0.57	-0.35
2	(1.43±0.0185)×10^-2	-4.29±0.106	11.8±0.317	0.23	-0.62	-0.46
3	(2.25±0.0246)×10^-2	-5.63±0.0999	21.9±0.395	0.29	-0.74	-0.46
4	(2.36±0.045)×10^-2	-7.95±0.212	50.8±1.27	0.37	-0.86	-0.47
5	(1.88±0.0277)×10^-2	-11.7±0.203	109±1.86	0.5	-0.93	-0.56
6	(1.14±0.0226)×10^-2	-14.6±0.339	210±4.5	0.56	-0.97	-0.6
7	(5.68±0.196)×10^-3	-24.9±0.985	729±25.8	0.69	-0.99	-0.7
8	(2.73±0.224)×10^-3	-43.2±4.07	2.05×10^3±170	0.77	-1	-0.78
9	(1.1±0.104)×10^-3	-61.3±7.08	5.54×10^3±523	0.77	-1	-0.77
10	(2.97±0.806)×10^-4	-135±43.2	(1.96±0.532)×10^4	0.83	-1	-0.83
11	(4.83±2.99)×10^-5	-66.2±124	(9.7±6.01)×10^4	0.32	-1	-0.32
12	(1.71±1.57)×10^-5	-114±276	(1.88±1.73)×10^5	0.37	-1	-0.37
13	(1.39±31.2)×10^-5	-163±3.64×10^3	(1.21±27.3)×10^5	1	-1	-1
14	(4.32±219)×10^-6	-1.59±689	(1.35±68.1)×10^5	-0.012	-1	0.012
15	(2.66±6.97)×10^-5	-0.0026±32.2	(3.53±9.29)×10^3	-0.091	-1	0.09
16	(2.6±1.8)×10^-5	-0.0018±7.38	144±109	-0.38	-0.98	0.36