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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b\bar{b} b\bar{b}$$\end{document}bb¯bb¯ 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathscr {L}=300$$\end{document}L=300 (3000) fb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}-1, 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 a e-mail: fady.bishara@physics.ox.ac.uk b e-mail: roberto.contino@sns.it c e-mail: j.rojo@vu.nl 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][4][5].
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 bbγ γ , bbτ + τ − , bbW + W − and bbbb (for a list of feasibility studies, see for example Refs. [17,30,[36][37][38][39][40][41][42][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 m t → ∞ 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 1/m t 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 O(1) accuracy, with details varying with the analysis and the specific final state; see Refs. [54][55][56][57][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][62][63][64] where a soft emission of two vector bosons from the incoming protons is followed by the hard V V → hh scattering, with V = W, Z . 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][66][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 V V → hh cross section with the partonic centre-of-mass energy,σ ∝ŝ/ f 4 , 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 hhV V quartic coupling. A first study of double Higgs production via VBF at the LHC was performed in Ref. [4], for a mass m h = 180 GeV, by focusing on the 4W final state. Following the discovery of the Higgs boson, more studies of the hh j j 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 hh → bbbb final state. While this final state benefits from increased signal yields due to the large branching fraction of Higgs bosons to bottom quarks, BR(H → bb) = 0.582 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 jetsubstructure techniques [71] in order to fully exploit the highenergy limit and optimize the signal significance.
We will thus focus on the kinematic region where the invariant mass of the Higgs pair, m hh , 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 hhV V quartic coupling c 2V . Given that for large m hh 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 gg → hh → bbbb channel, the 2b2 j background is comparable to the 4b component [43].
We quantify the constraints on the Higgs quartic coupling c 2V that can be obtained from VBF di-Higgs production at the LHC 14 TeV with L = 300 and 3000 fb −1 as well as at a future circular collider (FCC) with a centre-of-mass energy of 100 TeV and a total luminosity of 10 ab −1 . We find that, despite the smallness of the production cross sections, the LHC with 300 fb −1 can already constrain the hhV V coupling with an accuracy of +45% −37% around its SM value at the 1-σ level, which is further reduced to +19% −15% 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 bbbb final state. Our main results are presented in Sect. 4, where we quantify the potential of the VBF di-Higgs process to measure the hhV V 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 m hh 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 with v = 246 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 as where V (h) denotes the Higgs potential, The parameters c V , c 2V , c ψ , c 3 , and c 4 are in general arbitrary coefficients, normalized so that they equal 1 in the SM. The Higgs mass is fixed to be m h = 125 GeV [72]. As the notation in Eq. (2) indicates, the coefficients c V , c 2V , and c 3 control the strength of the hV V , hhV V and hhh couplings, respectively. The coefficients c ψ and c 4 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 treelevel 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 c 2V , c 2 V , and c V c 3 , respectively.
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 c V and c 2V , the amplitude of the partonic scattering V V → hh grows with the energy √ŝ until the contribution from the new states at the cutoff scale unitarizes it. The leading contribution in the energy range m W √ŝ ≡ m hh comes from the scattering of longitudinal vector bosons and is given by  thus provides a smoking-gun signature for the presence of BSM dynamics [3].
In the parametrization of Eq. (2), the amplitude for the process pp → hh j j can be decomposed as follows: where A, B, and C are numerical coefficients. In the present work, we will focus on the quartic coupling c 2V and set c V and c 3 to their SM values. This is justified for c V since the ATLAS and CMS measurements of Higgs production cross sections, when analyzed in the context of a global fit of Higgs properties [74][75][76] typically set bounds on c V − 1 at the level of 10-20%, depending on the specific assumptions madesee for example [77][78][79] and the references therein. Tighter limits on c V can be derived from electroweak precision tests in the absence of additional BSM contributions [80].
On the other hand, the trilinear Higgs coupling c 3 (where c 3 = λ/λ sm ) 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 bbbb final state [27] translates into the bound σ (hh)/σ sm (hh) 27 at the 95% confidence level. Achieving O(1) precision in the measurement of c 3 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 c 2V is achieved by reconstructing events with large values of m hh . In this kinematic region, it turns out that the sensitivity to c 3 is reduced, indicating that our analysis is not optimal to probe the Higgs trilinear coupling. For these reasons, setting c V = c 3 = 1 is a good approximation in the context of the present analysis. We can then define and this way the total cross section will be parametrized as However, while setting c 3 = 1 is a very good approximation, fixing c V = 1 is not as equally well justified. In particular, it would be more prudent to treat c V 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 c 3 . In this case, Eq. (7) is replaced by We will use this expression to evaluate the impact of c V on the derived bounds on δ c 2V at the end of Sect. 4. The values of the SM cross section σ sm and of the parameters A, B are reported in Table 1 for √ s = 14 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 m hh . We will make extensive use of this parametrization in Sect. 4 where we present our results in terms of the sensitivity on δ c 2V . Note that the value of A and B increase after imposing all cuts precisely because they have been optimized to enhance the sensitivity on c 2V .
Although we do not attempt to extract c 3 with our analysis, it is still interesting to discuss the dependence of the total cross section on this parameter. By fixing c V = c 2V = 1 and defining δ c 3 ≡ c 3 − 1, the cross section can now be parametrized as The 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 c 3 is suppressed by our analysis which aimed at measuring c 2V . Extracting c 3 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 bbbb 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 c V , c 2V , c ψ , c 3 and c 4 , 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 c i 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 SU (2) L 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 SU (2) L 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 ξ ≡ v 2 / f 2 , where f is the pNGB decay constant. For instance, minimal SO(5)/SO(4) models [81,82] predict On 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 be A precision model-independent determination of c 2V 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 c 2V and c 2V are in general predicted to be correlated [3,5]: where δ c 2 V ≡ c 2 V − 1. This follows because there is a single dimension-6 effective operator (O H in the basis of Ref. [3]) which controls the shift in both couplings. Therefore, a highprecision measurement of c 2V can test whether the Higgs boson belongs to a doublet in case a deviation is observed in c V [5].
Another interesting case is the scenario where the Higgslike 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][84][85][86][87][88]. In this model, invariance under dilatations implies δ c 2V = δ c 2 V , a condition that can be tested if the two couplings in Eq. (12) can be measured with comparable precision. Moreover, comparing c 2V and c V 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 c 3 , we will not consider these scenarios any further and will always assume δ c 3 = 0.
In the following, we will take as a representative benchmark scenario a model with c 2V = 0.8 corresponding to δ c 2V = −0.2, with the other couplings set to their SM values, namely δ c V = δ c 3 = 0.

Validity range of the effective theory
As discussed above and shown by Eq. (4), the amplitude of the partonic scattering V V → hh 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 V V scattering has been explored for instance in [91][92][93][94][95][96].
While we do not include the effect of such resonances in this work, we will report the sensitivity on δ c 2V 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 c 2V 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 m hh that is used to derive limits on δ c 2V . The explicit inclusion of scalar and vector resonances and their phenomenological implications for Higgs pair production via VBF is left for future work.

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 c 2V 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 antik t algorithm [99] with a jet radius R = 0.4.
The resulting jets are processed through a b-tagging algorithm, where a jet is tagged as b-jet with probability ε(b-tag) if it contains a b-quark with p b T > 15 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 ε(c-mistag) or ε(q, g-mistag) 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 bbbb final state requires efficient btagging 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][101][102][103][104] and the references therein. In the present study, we have considered two representative btagging 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 δ c 2V . Accordingly, we leave a detailed study of b-tagging including hadronization effects and p T dependence to future studies. 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 topologiesboosted, 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 c 2V 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 p T and the constituents of the hardest two jets are then re-clustered using the Cambridge/Aachen (C/A) algorithm [105] with R C/A = 1.2. Each C/A jet is processed with the BDRS mass-drop (MD) tagger [106]. This jetsubstructure tagger has two parameters: μ and y cut . And, in this work, we set μ = 0.67 and y cut = 0.09 as in the original BDRS study. To determine if a given jet arises from Boosted cat.

MD tags and at least 4 b-jets
Reject event Intermediate cat.
Resolved cat.
yes yes yes no no no Fig. 2 Schematic representation of the analysis strategy adopted in this work the decay of a massive object, the last step of the clustering for jet j is undone, giving two subjets j 1 and j 2 which are ordered such that m j 1 > m j 2 . Then, if the two subjets satisfy the conditions where R j 1 , j 2 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 j 1 until a mass drop is found or the C/A jet is fully unclustered.
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 p T b ≥ 15 GeV, and with a minimal angular separation R bb ≥ 0.1. 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 ghostassociation between large-R MD-tagged jets and small-R b-tagged jets [23,43], or accounting for an efficiency which depends on the jet p T [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 14 TeV and the FCC 100 TeV. We require the p T 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 pseudorapidities of light and b-tagged jets, η j and η b , 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]. 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 η j3 in the interval η min and η min j are the pseudo-rapidities of the VBF tagging jets 14 TeV 100 TeV Acceptance cuts Central jet veto: p T j 3 (GeV) ≤ 45 65

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 y j j ≡ |y lead j − y sublead j | and the invariant mass m j j 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 m j j . This definition is robust with respect to soft contamination from the underlying event (UE) and pileup (PU) and to the contribution of b-jets mistagged as light jets. Figure 3 shows the distribution of the rapidity separation | y j j | and invariant mass m j j 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 c 2V = 0.8 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 y j j and m j j distributions peak at zero. In Fig. 3, as well as in the subsequent figures, kinematic distributions have been areanormalized and then rescaled by a common factor such that the largest bin in the plot is of unit height.
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 LHC 14 TeV  Table 2. The distributions have been area-normalized and rescaled by a common factor these two variables are strongly correlated [4]. This large separation in rapidity is especially useful in the bbbb 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 |η| = 4.5 would lose more than 50% of signal events. The discontinuity in the FCC case delineates the edge of the b-jet acceptance region, |η| ≤ 3, above which no b-tagging is attempted and b-jets contribute to the light-jet yield.
Turning to the transverse momentum of the light jets, Fig. 5 shows the p T 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 p T , in order not to suppress the signal. Fortunately, in contrast from the gluon-fusion process, adopting a soft p T cut is not a problem since triggering can be performed based on the VBF topology. Comparing the p T distributions at 14 and 100 TeV, their shapes turn out to be rather similar, shifted towards larger values at 100 TeV.  (100) TeV. This discontinuity is more clearly visible in the 100 TeV curve and is purely due to combinatorics This justifies the harder p T cut in this case (see Table 2), also required to reduce the contamination from UE and PU.
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 η j 3 , lying between those of the VBF tagging jets, η max j > η j 3 > η min j , above a given p T threshold. The effect of the CJV is illustrated in Fig. 6 where we show the distribution of the p T of the third light jet, p T j 3 , 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 bbbb 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.

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 h 1 by identifying it with the pair of b-jets whose invariant mass is closest to the Higgs mass, m h = 125 GeV. Out of the remaining b-jet pairs, the one with an invariant mass closest to m h 1 is then assigned to be the second Higgs boson candidate, h 2 . 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 m h .
The invariant mass distributions of the Higgs candidates for the signal (SM and c 2V = 0.8) and the total background are shown in Fig. 7. The peak around m h = 125 GeV is clearly visible for signal events, especially in the case of h 1 . The smearing of the signal distribution of the second Higgs candidate h 2 arises from out-of-cone radiation effects which reduce the reconstructed mass. It is largest in the SM, while it is reduced in the c 2V = 0.8 scenario and in particular at 100 TeV, due to the larger boost of the Higgs bosons. The small peak in the background distributions for h 1 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 c 2V = 0.8 benchmark is another validation of the scale-invariant tagging, since while in the SM most events lead to resolved topologies, the c 2V = 0.8  scenario is dominated by the boosted category (see Fig. 8 below). After reconstructing the Higgs candidates, we require that their invariant masses, m h 1 and m h 2 , are reasonably close to the nominal mass. In the resolved category, these conditions are The mass window in Eq. (16) This condition greatly reduces the background rates while leaving the interesting kinematic region at large m hh -where deviations from the SM signal mostly appear -unaffected. Figure 8 shows the m hh distribution, after all the cuts listed in Table 2

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 m hh cuts of Table 2 and Eqs. (15)- (17), respectively, for both the signal (SM and c 2V = 0.8) and for the total background. At 14 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 c 2V = 1 as illustrated by the benchmark value of c 2V = 0.8 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 c 2V = 0.8 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 c 2V by increasing the signal significance in the large-m hh region.
From Table 3 we also find that a similar qualitative picture holds at 100 TeV with the important difference that, in this case, the event rate is already substantial in the SM which yields 2000 events after the acceptance cuts with L = 10 ab −1 . The cross section enhancement at 100 TeV as compared to 14 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 100 (150) in the c 2V = 0.8 scenario. At 100 TeV, the ratio of signal cross sections in the c 2V = 0.8 and SM scenarios is ∼15 (50) 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 c 2V = 0.8) at the level of inclusive rates. It is only by exploiting the large-m hh 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 m h , while leaving the signal mostly unaffected. A final reduction of the background, by around another order of magnitude, is achieved through the   Figure 9 graphically illustrates the dependence of the di-Higgs production cross section on the couplings c 2V and c 3 . The left panel shows the cross section in SM units as a function of δ c 2V = c 2V − 1 and δ c 3 = c 3 − 1 after applying the acceptance cuts of Table 2 (dashed curves) and after all the analysis cuts (solid curves). The sensitivity on c 2V is particularly striking, for example the cross section for |δ c 2V | 1 is enhanced by a factor ∼50 compared to its SM value after all cuts. This sensitivity is the key ingredient for measuring c 2V 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 δ c 2V derives mainly from the tail of the m hh distribution. This observation elucidates the enhancement (suppression) in the sensitivity to δ c 2V (δ c 3 ) in Fig. 9 after the application of all the cuts which, for instance, remove the threshold region up to m hh = 500 (1000) GeV at the LHC(FCC).
The right panel of Fig. 9 shows, instead, the ratio between the VBF di-Higgs cross sections at √ s = 100 and 14 TeV. Given the larger centre-of-mass energy of the FCC, it is expected that this ratio grows rapidly for δ c 2V = 0 and, indeed, it can reach values as high as 300 for δ c 2V 1. As will be demonstrated in Sect. 4, this effect allows for much more precise measurements of c 2V 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 c 3 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 c 3 comes from events near the di-Higgs threshold, m hh 2m h , 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 c 3 , 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 2b2 j 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 bb pair during the parton shower. Even with a small light-jet mistag rate of O(1%), 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 m hh distribution. This substantially reduces its contamination to the region with the highest sensitivity to c 2V in our analysis. Note also that a harder m hh 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 m hh 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 2b2 j component is a factor 10 larger than 4b at 14 TeV while they are of similar size at 100 TeV.
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 0.98 (0.009) fb after acceptance (all) cuts, compared to 0.11 (0.002) fb for the VBF case. On the other hand, this fact does not affect the measurement c 2V , since, as we show next, the sensitivity comes from the large m hh tail where the gluon-fusion component is heavily suppressed.
The decomposition of the total background in terms of individual processes as a function of m hh 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 14 TeV case, the 4b background dominates for large m hh while the 2b2 j one is instead the most important for small m hh . The 100 TeV case is similar with one exception, namely that the gluon-fusion di-Higgs background becomes the dominant one for very high invariant masses, m hh 10 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.

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 c 2V = 0.8 benchmark scenario. As anticipated, the additional handle which we can now exploit to increase the signal significance is the different behaviour of the m hh distribution for the signal and the background, in particular when c 2V departs from its SM value. The latter has a sharp fall-off at large m hh values, while, instead, the signal exhibits a much harder spectrum for c 2V = 1. This crosssection growth implies that, for |δ c 2V | sufficiently large, there will be a crossover value of m hh 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 c 2V = 0.8) and the total background. In the case of the benchmark scenario with c 2V = 0.8, the crossover between signal and background is located at m hh 2 TeV (4 TeV) at 14 TeV (100 TeV). We also observe that, for invariant masses m hh above this crossover, the ratio between the signal and the backgrounds keeps increasing steeply.
With the final results of our analysis in hand, we can now estimate the expected sensitivity to deviations in the hhV V  coupling, parametrized as δ c 2V = c 2V − 1, by exploiting the information contained in the full m hh 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 m hh and then follow a Bayesian approach [109] to construct a posterior probability density function. We include two nuisance parameters, θ B and θ S , to account for the uncertainty associated with the background and signal event rate, respectively. The parameter θ S encodes the theoretical uncertainties on the di-Higgs cross section and the branching fraction BR(h → bb). We conservatively assume a 10% uncertainty uncorrelated in each m hh bin.
Concerning θ B , 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 bb cross sections [110]. The background nuisance parameter, θ B , is conservatively also assumed to be uncorrelated among m hh 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 reads with N i (θ i B , θ i S ) and N i obs denoting, respectively, the number of predicted (for a generic value of c 2V ) and observed (assuming SM couplings) events for a given integrated luminosity L in the ith bin of the di-Higgs invariant mass distribution m hh , given by 2 : In Eq. (19), σ i sig (c 2V ) and σ i bkg indicate the signal (for a given value of c 2V ) and total background cross sections, respectively, for the i-th bin of the m hh distribution. The functional   12 23 form of σ i sig (c 2V ) is given by Eq. (7) and the value of the coefficients in bin i are given in Appendix D. We denote by π(c 2V ) the prior probability distribution of the c 2V coupling.
As justified above, in the evaluation of Eq. (18) we set δ B(S) = 0.15 (0.1) 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 L(N i |N i obs ) in each bin and assume the prior probability π(c 2V ) to be uniform. The resulting posterior probabilities are shown in Fig. 12 for the LHC with L = 300 fb −1 (LHC 14 ) and L = 3 ab −1 (HL-LHC), and for the FCC with L = 10 ab −1 . 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.
From Fig. 12, we can determine the expected precision for a measurement of δ c 2V at the LHC and the FCC in the case of SM values of the Higgs couplings. The 68% probability intervals for the determination of c 2V 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 c V on the bound on δ c 2V by treating c V as a nuisance parameter and marginalizing over it. The leading effect of varying c V comes from the (c 2V − c 2 V ) 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 c 2 V and c V c 3 . We take c V 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 c V is subleading in both LHC scenarios and weakens the bound on δ c 2V . We find that the results of Table 5 change by 2% for LHC 14 and 7% for HL-LHC. The effect at the FCC is much larger causing the bound on δ c 2V to be O(0.04) 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 c 2V coupling, for which there are currently no direct experimental constraints, can already be measured at the LHC with 300 fb −1 with a reasonably good accuracy: +45% −37% 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 +19% −15% , is expected at the HL-LHC with 3 ab −1 . This estimate is robust against an overall rescaling of the background cross section. Finally, we find a very significant improvement at the FCC with 10 ab −1 , 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  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 c 2V can obtained despite the rather loose constraints that can be obtained on the VBF di-Higgs cross section even at 100 TeV. As already discussed, this behaviour follows from the strong dependence of the signal cross section on c 2V ; 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 m hh . However, the EFT expansion might break down at large enough values of m hh , 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 m hh 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 g * and mass scale [3], one naively expects Therefore, for maximally strongly coupled UV completions (with g * 4π ) it is possible to derive the following upper limit, which makes explicit the connection between the value of δ c 2V and the new physics scale . The validity of the EFT can thus be monitored by introducing a restriction on the m hh bins used in the construction of the posterior probability Eq. (18), so that m hh ≤ m max hh , and then determining how the sensitivity on δ c 2V varies as a function of m max hh [97]. The precision on δ c 2V , defined though the symmetrized 68% probability interval, is shown in Fig. 13 as a function of m max hh for the LHC and the FCC. As expected, increasing m max hh , i.e. making the cut less stringent, leads to stronger constraints. Eventually, δ c 2V flattens out when m max hh is large enough that all the m hh 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 δ c 2V > δ max c 2V , obtained by setting = m max hh 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 ( , g * ) [97], corresponding to the limits on δ c 2V derived as a function of m max hh . This is shown in the left panel of Fig. 13, making use of Eq. (20) and then setting = m max hh . The grey area in the upper part of the plot indicates the non-perturbative region, defined by g * ≥ 4π . We find that the dominant constraints on δ c 2V 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 δ c 2V in the absence of new explicit degrees of freedom.
Finally, in Fig. 14 we show the signal significance, S/ √ B, in the c 2V = 0.8 scenario as a function of the di-Higgs invariant mass m hh 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 S/ √ B 2 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 S/ √ B 5 already for m hh 1.5 TeV and then increasing very rapidly for higher values of m hh . 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 c 2V coupling already at the LHC with L = 300 fb −1 (Table 5), which at a 100 TeV collider would become a high-precision measurement with potentially subpercent 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 hhV V quartic coupling c 2V .
Exploiting the high signal yield of the bbbb 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 c 2V = c 2 V significantly harden the m hh distribution resulting in a large fraction of events with a boosted Higgs pair. The subsequent decays into bb 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 m hh distribution when the Higgs couplings depart from their SM values leads to a remarkable model-independent sensitivity to the c 2V coupling.
Our results demonstrate that at the LHC with an integrated luminosity of L = 300 (3000) fb −1 the hhV V coupling can be measured with +45% −37% ( +19% −15% ) precision at the 68% probability level, reaching 1% accuracy at a 100 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 bbbb final state studied here, since the sensitivity to λ arises from the threshold region m hh 2m h 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 c 2V 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.
Matter: One Solution for Three Mysteries (DaMeSyFla). J. R. is supported by the ERC Starting Grant "PDF4BSM".
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: Monte Carlo event generation
In this appendix we discuss the event generation of the signal and background process used in this analysis. In each case, we discuss the programs used, the input parameters and theoretical settings, and the cross-checks that have been performed.
Signal events for Higgs pair production via VBF have been generated at leading order (LO) with MadGraph5 aMC@NLO [112] using a tailored UFO [113] model that allows for generic values of the c V , c 2V , and c 3 couplings in the Lagrangian Eq. (2) (see [114,115] for other UFO models of Higgs EFTs). We generated 1M unweighted events for each value of c 2V and for centre-of-mass energies of 14 and 100 TeV. The size of the signal sample is dictated by the condition of achieving an adequate coverage of the large m hh region. Even so, for small deviations of c 2V with respect to its SM value, this region is very difficult to populate. We overcome this limitation by performing a fit to the cross section in each m hh bin as a function of c 2V using the general parametrization of Eq. (7).
The MadGraph5_aMC@NLO generation of signal events uses the NNPDF2.3LO [116] set with α S (m Z ) = 0.119 interfaced via LHAPDF6 [117]. The calculation is performed in the n f = 4 scheme and the factorization and renormalization scales are taken to be the W boson mass, μ F = μ R = M W . To account for higher-order effects, we apply an NNLO/LO K -factor 1.1 on the inclusive cross section [63] both at 14 and 100 TeV. The decays of the Higgs bosons into bb pairs are performed within MadGraph5_aMC@NLO with adjusted parameters to ensure that the branching ratio corresponds to Higgs Cross-Section Working Group value [35] of BR(h → bb) = 0.582. Since in this work we only modify the c 2V coupling, assuming the SM value of the total Higgs width is a very good approximation.
With this setup, we have generated signal events in the SM and for a wide range of values of δ c 2V . At the matrix-element level, we applied the generation cuts listed in Table 7. The corresponding cross sections are summarized in Table 8, where we provide the SM results as well as the predictions for various BSM scenarios defined by the couplings {c V , c 2V , c 3 }, including the benchmark scenario c 2V = 0.8. In addition, the number of expected events at the HL-LHC, assuming an integrated luminosity of L = 3 ab −1 , and at the FCC 100 TeV for L = 10 ab −1 , are also listed. Following the parton-level event generation, the resulting Les Houches event files for signal events are showered using the Pythia8 Monte Carlo generator [118] v8.212 using the Monash 2013 Tune [119]. No hadronization or underlying event (UE) effects are included for simplicity. Although we have not attempted to simulate the effects of pile-up (PU) in our analysis, recent studies [43] indicate that modern PU mitigation techniques [120] should be efficient enough to minimize the PU contamination even for the bbbb final state.
Concerning the generation of the QCD background processes, their parton-level cross sections are summarized in Table 9. In each case we indicate the programs used for the event generation, the number of MC events N ev that have been generated, whether the events are weighted or unweighted, and the LO cross sections along with the corresponding higher-order K -factors at 14 and 100 TeV. As discussed in Sect. 4, our results include a 15% systematic uncertainty on the total background normalization (assuming a data-driven determination in an experimental analysis), and we also assessed the robustness of our estimate for the measurement of δ c 2V in case of an overall upwards shift of the backgrounds by a factor 3.
For QCD multijet production we have explored two complementary approaches for event generation. First of all, we used ALPGEN [121] to generate a large sample of unweighted events at the matrix-element level and showered them with Pythia8. This approach, however, presents a difficulty, because the differential cross section with respect to the di-Higgs invariant mass m hh falls very rapidly for increasing m hh . We thus find that an unrealistically large sample of unweighted events would be needed in order to adequately populate the tail of the m hh distribution.  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 N ev and whether they are weighted [W] [47] To bypass this limitation, we generated 50M weighted 4b and 2b2 j events with Sherpa v2.2 [123] and then processed them with the built-in shower. The main advantage of this approach is that events with small weights are not discarded by the unweighting and, therefore, enough events with large m hh (and thus small weight) are still kept. One possible drawback is that, for the same number of events, the unweighted sample provides a better estimate of the total cross section than the weighted one. In our case this is not an issue, since with our weighted sample we achieve statistical uncertainties of order 2% (to be compared with 0.01% for a same-size unweighted sample), which is more than sufficient for our purposes. On the other hand, to generate a large enough number of events in a reasonable amount of CPU time we are forced to use a lower multiplicity final state, and thus we rely on the parton shower to generate the additional light partons (see Table 9). To validate this procedure, we explicitly verified that, in the region where the two approaches lead to sufficient statistics, the Sherpa calculation based on 4b (2b2 j) matrix elements and the ALPGEN one, based on 4b2 j (2b4 j) matrix elements, lead to comparable results (see Appendix C). This agreement indicates that the parton shower does a reasonable job in modelling additional hard radiation in the phase-space region of interest.
The Sherpa samples in Table 9 have been generated using the NNPDF3.0 NNLO set [124] with strong coupling α S (m 2 Z ) = 0.118 and with n f = 4 active quark flavours, and use as factorization and renormalization scales μ F = μ R = H T /2, with where n is the final-state matrix-element multiplicity and m t,i and p t T,i are the transverse mass and transverse momentum of the i-th final-state parton. The generation-level acceptance cuts applied to these samples are listed in Table 7. As for the other background samples, LO cross sections have been rescaled by the best available higher-order K -factors.
Higgs pair production via gluon fusion is simulated at LO using MadGraph5_aMC@NLO for loop-induced processes [125]. The cross section is rescaled to match the inclusive NNLO+NNLL calculation [47] yielding a K -factor, σ NNLO+NNLL /σ LO = 2.4 (2.2) at 14 TeV (100 TeV). Partonlevel events are then showered with Pythia8 using the same settings as for the signal VBF di-Higgs samples. No generation cuts are applied, and the resulting cross sections (including the branching fraction into bbbb) are listed in Table 9. While the hard-scattering process that is generated, gg → hh, does not include additional jets, initial state radiation from the gluon legs can give a VBF-like topology with two forward jets characterized by a large invariant mass, and thus contribute to the total hh j j yield. We have verified that our calculation provides a reasonable description of the relevant kinematical distributions, such as m hh , by comparing with the gg → hh j j process computed in the EFT approximation. As shown in Sect. 4, while the gluon-fusion contamination to the VBF signal is substantial close to threshold, it is marginal in the large m hh region where the sensitivity to c 2V is the highest.

Appendix B: Fitting the tail of the m hh distribution for background processes
The background processes considered in this work (see Appendix A) exhibit a steep fall-off for large values of m hh , the invariant mass distribution of the reconstructed Higgs pair. Even with the use of weighted events, it is difficult to adequately populate this region. To obtain a reliable estimate of the cross section there, it is thus necessary to introduce a fitting procedure. In this appendix we discuss how the fits to the m hh distributions for the background processes, and the associated validation tests, have been performed.
We found that, far from the di-Higgs production threshold, the following functional form provides a reasonable description of the m hh distribution: where √ s is the collider centre-of-mass energy. This specific functional form is a modified version of one of those suggested in Ref. [126] for the fit of the di-photon invariant mass distribution at the LHC. We explored other choices, finding a comparable fit quality.
In Eq. (B.1), the fit parameters {a, b, c} are determined by minimizing the χ 2 , where n is the number of bins in the m hh distribution used as input to the fit, σ (th) i is the theoretical prediction for the cross section in the bin m (i) hh , and is the statistical uncertainty associated with the N i weighted Monte Carlo events that populate bin i with weights {w (i) k }. Note that when the MC events are unweighted, δ i ∝ √ N i , as expected. The resulting fit coefficients are shown in Table 10 for both centre-of-mass energies and for the four background processes considered. We have also verified that the results for the fit parameters are stable with respect to variations of the binning in m hh .
To validate the fitting procedure, we show in Fig. 15 the m hh distribution of the reconstructed di-Higgs system for the total background, where each individual process has been fitted separately and then added up to construct the solid histograms. The fit results are compared to the cross sections from the Monte Carlo generation (indicated by filled circles) with their corresponding statistical uncertainty, Eq. (B.3). In the lower panels of Fig. 15 we show the fit residuals, defined as the difference between MC and fit divided by the statistical uncertainty in each bin. The fact that the parametrized fit agrees with the MC calculations at the 2-σ level or better (in units of the uncertainty of the latter) for a wide range of m hh demonstrate the goodness of these fits.
Note that, in the case of the QCD multijet backgrounds, we exclude the first bin to ensure that the fit is not affected by artificial features in the m hh distribution induced by the analysis selection cuts. Furthermore, in the case of the gluon-fusion di-Higgs background, it is necessary to exclude the first few bins of the m hh distribution from the fit. The reason is that in this case there is a production threshold at 2m h , and, as a result, the distribution does not decrease monotonically with m hh unless one is far enough from threshold. By excluding these bins, we avoid biasing the resulting fit in the tail of the m hh distribution, the region where a functional form such as Eq. (B.1) does provide an equally satisfactory description as for the rest of background processes.

Appendix C: Validation of the QCD multijet generation
As discussed in Appendix A, the modelling of the tail of the m hh distribution for background processes is particularly challenging. One reason is because all the backgrounds considered are characterized by a steep power-like fall-off with m hh above the di-Higgs production threshold, and therefore it becomes necessary to introduce a cross-section parametrization (see Appendix B) to be able to cover this region. In addition, in the case of the QCD multijet backgrounds, there are different options available for the modelling of the m hh distribution, in particular the multiplicity of the matrixelement calculation prior to the parton shower. Ideally, one should generate all relevant final-state partonic multiplicities and merge them to avoid double counting, either at LO [127][128][129] or at NLO [130,131]. For the purposes of this feasibility study, however, we found it sufficient to generate LO samples using 4b and 2b2 j matrix elements with Sherpa, with additional hard radiation provided by the parton shower. Note that our approach could introduce some double counting from gluon splittings into bb, which if anything would increase the background cross section and make our estimates of δ c 2V more conservative. Moreover, let us recall that, as discussed in Sect. 4, an actual experimental analysis would estimate the overall background normalization by means of data-driven techniques.
To validate the robustness of our simulation of the m hh distributions for the QCD multijet background, we have compared the Sherpa calculation, based on 4b and 2b2 j matrix elements and weighted events, with the ALPGEN simulation, based on 4b2 j and 2b4 j matrix elements and unweighted events showered with Pythia8. In principle, ALPGEN should provide a better description of the kinematics of the VBF jets, since it is based on higher-multiplicity matrix elements but it has the drawback that populating the tail of the m hh with unweighted events is very CPUtime intensive. On the other hand, it turns out that the two approaches give comparable results for the m hh distribution in the region where both approaches lead to sufficient MC statistics, validating the use of the Sherpa for the calculation of background cross sections.
In Fig. 16 we show the m hh distribution for the 4b background at 14 and 100 TeV, comparing the ALPGEN and SHERPA calculations. We find good agreement for the entire m hh range, indicating that the SHERPA event generation does a reasonable job in modelling the VBF tagging jets, and demonstrating that it can be reliably used to populate the large m hh region for the multijet backgrounds with high efficiency. The differences between the two calculations are at most a factor two, typically less, well within the typical size of the theoretical uncertainties for LO multijet calculations. Figure 17 shows a similar comparison for the invariant mass of the two VBF tagging jets, m j j , and for the transverse momentum of the hardest light jet in the event, p T j 1 . For the m j j distribution, the agreement between the ALPGEN and SHERPA calculations is also good for the entire kinematical range. On the other hand, for the transverse momentum of the leading jet the ALPGEN calculation leads to harder (softer) spectra than the SHERPA one at high (low) values of p T j 1 . This can be understood from the fact that j 1 will be typically generated by the parton shower in the latter case, and by the matrix element in the former. These differences are however inconsequential for our analysis, since the cuts imposed on  the p T of the VBF tagging jets are relatively mild, see Table 2, and thus the event selection will not be affected. The validation studies discussed in this appendix demonstrate that, for the purposes of the present analysis, our approach to event generation based on SHERPA for the simulation of the QCD multijet backgrounds is justified. On the other hand, they also highlight that future studies aiming to enhance the separation between signal and background events from shape comparisons of kinematical distributions, such as multivariate analysis [43], would require an improved modelling of the QCD multijet backgrounds. This could be achieved by using merging techniques to combine QCD jet samples of different multiplicities.

Appendix D: Coefficients of the δ c 2V fit
The dependence of the signal cross section on δ c 2V , as parametrized in Eq. (7), is required in order to construct the likelihood function. In this appendix, we list the coefficients of Eq. (7) in each m hh bin. The coefficients are extracted by fitting MC events after all cuts have been applied as discussed in Sect. 3 with the exception of the m hh cut. We use 15 equally spaced bins on a log scale starting from 250 GeV up to 6(30) TeV for √ s = 14(100) TeV. In addition, we define an overflow bin up to the centre-of-mass energy. The results are shown in Tables 11 and 12. The fit error on each coefficient is also provided along with the off-diagonal entries of the correlation matrix ρ i j where i, j ∈ {σ, A, B}. Table 11 The bin by bin fit coefficients obtained by fitting MC events to Eq. (7) for the LHC with √ s = 14 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 ρ 0 A , ρ 0B , and ρ AB give the coefficients of the correlation matrix among the three fit parameters