Precise predictions for double-Higgs production via vector-boson fusion

Theoretical predictions with next-to-next-to-leading order (NNLO) QCD accuracy combined with the next-to-leading order (NLO) electroweak (EW) corrections are presented for differential observables of the double-Higgs production process via vector-boson fusion. While the QCD corrections were previously known, the EW ones are computed here for the first time. The numerical results are obtained for a realistic experimental set-up at the LHC and are presented in the form of fiducial cross sections and differential distributions. Within this setup we find that the VBF approximation employed in the NNLO QCD correction is accurate at the sub-percent level. We find that the NLO EW corrections within the fiducial volume are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\,6.1\%$$\end{document}-6.1%, making them of almost the same order as the NLO QCD corrections. In some kinematic regions they can grow as large as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\,30\%$$\end{document}-30% making them the dominant radiative corrections. When the EW corrections are combined with the NNLO QCD corrections we find a total correction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\,14.8\%$$\end{document}-14.8%. The results presented here thus comprise the state-of-the-art theoretical predicition for the double-Higgs production via vector-boson fusion, which will be of value to the high-luminosity programme at the LHC.


Introduction
The discovery of the Higgs boson at the Large Hadron Collider (LHC) [1,2] opened the door to a new era in high-energy physics, dominated by the quest for precision, and which will culminate in the next decade with the ambitious highluminosity programme at the LHC [3]. The determination of a e-mail: frederic.dreyer@physics.ox.ac.uk (corresponding author) b e-mail: alexander.karlberg@physics.ox.ac.uk c e-mail: jlang@physik.uzh.ch d e-mail: mpellen@hep.phy.cam.ac.uk all the fundamental parameters of the Higgs sector will be at the heart of this programme. Beyond the mass and the width of the Higgs boson, its interactions with the other particles of the Standard Model (SM) will be precisely scrutinised. One of the most far-reaching tasks will be the determination of the Higgs self-couplings. These are fundamental parameters of the SM Lagrangian that determine the shape of the Higgs potential. The investigation of these properties will rely notably on the investigation of processes involving a pair of Higgs bosons in the final state [4][5][6][7][8][9][10][11][12][13][14][15][16].
The production of two Higgs bosons via vector-boson fusion (VBF) i.e. pp → HHjj is a particularly important process for the determination of the triple-Higgs coupling [17]. While gluon fusion is the dominating double-Higgs production mode [18][19][20][21][22], the VBF channel offers unique opportunities for measurements of Higgs pair production. In VBF, the Higgs bosons are produced at leading order from heavy gauge bosons that are themselves radiated off two quark lines. These two quarks offer a useful handle as tagging jets for its experimental measurement, providing a promising channel for studies of the trilinear and quartic Higgs couplings at the LHC [17,23].
While single-Higgs production via VBF has been measured [24], its double-Higgs counterpart is not yet accessible with the current experimental data set. The reason is that the cross section at 14 TeV is of the order of 1 fb as it is a purely electroweak (EW) process and requires very exclusive event selections in order to single it out from its background. Nonetheless, this process is already used to look for new physics [25]. The high-luminosity programme of the LHC is aiming at collecting about 3000 fb −1 in the next two decades. This should allow the observation of the double-Higgs production via VBF. To that end, precise and reliable theoretical predictions are required.
In this article we provide for the first time next-to-leading order (NLO) EW corrections and combine them with the existing next-to-next-to-leading order (NNLO) QCD corrections. Both types of corrections are very different in magnitudes and shapes as they account for different physical effects. The QCD ones are usually larger and their missing higher-order terms can be well estimated by a variation of the renormalisation and factorisation scales. In Refs. [26,27], it has been shown that only NNLO QCD corrections provide reliable predictions at the differential level. Recently, the N 3 LO corrections to the inclusive cross section have also been computed, and were shown to be at the few permille level [28]. The EW corrections on the other hand are typically rather suppressed and appear in the high-energy tail of differential distributions where the effect of Sudakov logarithms become large.
The present article is organised as follow: in the first part, the details of the computation are provided. In particular the numerical inputs and the event selection are explained. In the second part, the results are given in the form of cross sections and differential distributions. Finally, the conclusion contains a short summary of the article as well as concluding remarks.

Description of the process
The process of interest is the double-Higgs production via VBF, which can be expressed as pp → HHjj. (1) It is a purely EW process which is defined at tree level at the order O α 4 . At this order, it contains the VBF topology which is defined as two quark lines that radiate heavy gauge bosons that fuse to give rise to two Higgs bosons, as shown in Fig. 1a. Alternatively, the two Higgs bosons can be radiated from a heavy gauge boson decaying subsequently into two quarks. The latter topology is usually referred to as Higgs Strahlung. While these two topologies are intimately related, they entail rather different physical effects. It has lead, on the experimental side, to consider the two processes separately by imposing rather stringent cuts on the invariant mass and the rapidity difference of the two tagging jets. This has motivated the use of the VBF approximation for theoretical predictions, which neglects s-channel contributions as well as t/u interferences in a gauge-invariant way. Additionally the two quark lines are not allowed to exchange virtual or real gluons. One advantage of this approach is that it greatly simplifies the computations. In particular, it has allowed the computation of QCD corrections up to N 3 LO of single and double Higgs production [26,[28][29][30][31][32][33][34]. The differential cross section of the di-Higgs process was also com-puted at NLO with matching to parton showers in this approximation [35].
In the present article, we have used the full computations at LO, NLO QCD, and NLO EW while we have utilised the VBF approximation at NNLO. In order to justify their combination, we have used a rather exclusive experimental set-up where the VBF approximation holds at the per-mille level. In addition, we have added a correction factor to the NNLO corrections to account for the differences between the full and the VBF computation. More precisely we have checked that the difference between the full computation and the VBF-approximate one is identical at LO and NLO QCD in all differential distributions. Recently, extensive studies investigating the quality of such approximation at NLO in QCD for vector-boson scattering (VBS) [36] and the EW production of a Higgs boson in association with three jets [37] have been performed. We have refrained from presenting similar results for the process at hand as it is very close to processes mentioned above. Instead we have displayed the correction factor used in all differential distributions which indicates the difference between the full and the approximate computation.

• QCD corrections
In the full computation at NLO, the real corrections consist in all the contributions of the type pp → HHjjj at order O α s α 4 . The virtual corrections entail all the possible gluon insertions on a given quark line, with gluon exchanges between the two quark lines vanishing for colour reasons. At NNLO, for which an example diagram is shown in Fig. 1b, heavy-quark loop-induced diagrams are neglected, as well as t-/u-channel interferences, single-quark line contributions and loop induced interferences between VBF and gluon-fusion production. These contributions have been shown to be suppressed to less than a percent in the single-Higgs case [38][39][40], but are not known for Higgs-pair production. Although they could be sizeable there is no apriori reason to expect them to be enhanced in di-Higgs production. A dedicated future study is however necessary to confirm whether or not these effects can be neglected. The non-factorisable diagrams involving the exchange of two gluons between the quark lines, have recently been estimated in Ref. [27]. Unlike their counterparts in single Higgs VBF production [41], they have been shown to be of the same order as the corrections arising from gluon exchanges limited to one quark line, the so-called factorisable corrections shown in Fig. 1b. In this work we therefore also provide an estimate of the non-factorisable corrections, but will show them separately from the factorisable corrections. Unless explicitly specified, when referring to NNLO QCD corrections, we will always mean the factorisable ones. We compute the factorisable NNLO QCD corrections using the projection-to-Born method as detailed in Ref. [26].

• EW corrections
For the EW corrections the real radiations are made of the pp → jjHHγ channels at order O α 5 . At the same order, the virtual corrections are obtained by inserting EW particles anywhere possible in the tree-level topologies, an example of which is shown in Fig. 1c. Note that at the order O α 5 , photon-induced contributions also arise. These have been neglected in the present work as these have been shown to be rather small for similar processes [38,42,43]. 1 Note that EW corrections to single-Higgs production have been computed for the first time in Refs. [38,46] and are available in HAWK [44]. Later they have also been obtained in VBFNLO [47,48].
As mentioned previously, all LO and NLO predictions are based on the full computation, i.e. without employing the VBF approximation. These have been obtained from the Monte Carlo MoCaNLO, which has already been used for a variety of processes and in particular VBS ones [42,43,49] at NLO EW and NLO QCD. The matrix elements are provided by Recola [50][51][52] which internally uses the Collier library [53,54] to evaluate tensor integrals.
On the other hand, the NNLO QCD corrections have been obtained from proVBFHH v1.1.0 [26,28] which uses the projection-to-Born method [32] to compute the fully differential NNLO corrections in the VBF approximation. In order to correct for the mismatch between this computation and the full computation used for the LO and NLO computations, we compute a differential correction factor and obtain the NNLO cross section provided below in the following way where the full quantities refer to the computations with no approximations and the VBF one to the relative NNLO corrections in the VBF approximation. At the differential level, the NNLO predictions are obtained in the same way. We have checked numerically that if one were to obtain K full/VBF using instead the NLO cross section, the results do not change within statistical uncertainties. This is only true under the very stringent cuts that we will introduce below. The non-factorisable NNLO QCD corrections, δ NF NNLO QCD , have also been obtained using proVBFHH v1.1.0 as computed in Ref. [27]. Since they are very different in nature from the factorisable corrections, we do not correct them by K full/VBF , and show them separately from all other results.
Finally, to combine the NNLO QCD prediction and the NLO EW ones, we have followed the Higgs cross section working group recommendation [55] for VBF which reads with σ full NLO EW = σ full LO +δ full NLO EW . The same formula has been applied differentially. With such a prescription we provide thus best predictions at fixed order in the SM for double-Higgs production in VBF. We stress that we do not include δ NF NNLO QCD in this quantity, but rather quote it separately.

Numerical inputs
Our predictions are obtained for proton-proton collisions at the LHC running with a centre-of-mass energy of √ s = 14 TeV. The 5-flavour scheme is used throughout the computation i.e. the bottom quarks are considered massless. Bottom quarks are also included in the jet definition. For the parton distribution functions (PDF) and α s , the NNPDF31_nnlo_as_0118_luxqed set [45] has been used for all computations i.e. at LO, NLO, and NNLO as implemented in LHAPDF6 [56]. The choice of the central renormalisation and factorisation scale is the same as the one used in Ref. [26] and reads with M H the mass of the Higgs boson and p T,HH the transverse momentum of the di-Higgs system. In order to estimate the impact of missing higher-order QCD corrections, we perform a 3-point scale variation defined by μ ren = μ fac = {1/2, 1, 2} × μ. Note that this method is applied for the LO predictions and the factorisable QCD corrections only. For the EW corrections, such a method does not provide a good estimate of missing higher-order and thus NLO EW corrections are simply applied as an overall factor. The masses and widths used for the numerical simulations read The mass of the bottom quark is taken to be zero in accordance with the 5-flavour scheme. The width of the top quark is also taken to be zero as no top quarks are produced resonantly. The value taken for the Higgs-boson mass is taken from the report of the Higgs cross section working group [55]. Note that the width has been taken to zero as the Higgs bosons are on-shell external states. They can nonetheless appear as internal propagator in the splitting H * → HH.
As the invariant mass of the off-shell Higgs boson tends to be far from the on-shell mass and since the value of the width, H = 4.07 × 10 −3 GeV, is small, setting the width to zero has no numerical impact. 2 The pole masses and widths used for the simulations are obtained from the measured on-shell 2 For the NNLO corrections the width has been used for the internal Higgs propagators.
(OS) values [57] for the W and Z bosons according to with V = W, Z.
The G μ scheme [58] is used for all computations and is translated to α via In all the LO and NLO computations, the intermediate W/Z-boson resonances are treated in the complex-mass scheme [59][60][61] to ensure gauge independence of all amplitudes.

Event selection
The experimental cuts are rather generic and typically intend to single out VBF EW contributions from their QCD background. To that end, high invariant-mass and large rapidity difference between the two tagging jets are required. The tagging jets are defined by requiring that each jet fulfils the following condition: p T,j > 25 GeV and |y j | < 4.5. (9) The jets are clustered using the anti-k t algorithm [62] with R = 0.4, using FastJet v3.3.0 [63] in the case of proVBFHH. The two hardest jets in the transverse momentum fulfilling these requirements are required to obey the VBF-selection cuts which read m j 1 j 2 > 600 GeV and |y j 1 − y j 2 | > 4.5. (10) Note that we do not apply a cut ensuring that both tagging jets are in opposite hemispheres as in Ref. [26]. Also, the event selection is completely inclusive in the final state Higgs bosons and no cuts of any kind are applied to them.

Results
In this section we show numerical results for cross sections and differential distributions. Given that the QCD corrections have already been presented in Ref. [26,27], the discussion focuses more on their combination with the EW ones which are presented here for the first time. Nonetheless, some distributions were not shown in Ref. [26] and are therefore discussed here in more details. In the following, the NNLO QCD × NLO EW predictions are sometimes referred to as state-of-the-art predictions. Finally, the differences between the full and the VBF computation are highlighted.
In Table 1, fiducial cross sections and higher-order corrections are displayed for the event selection presented in Sect. 2. They are expressed both in femto barn and in per cent. The numbers in per cent are with respect to the LO cross section. The numbers in parenthesis indicate the statistical error while the additional information on σ full LO and σ NNLO QCD×NLO EW gives the scale variation estimate. Note that the total statistical uncertainty is not obtained by adding the individual statistical uncertainties in quadrature, as these are all correlated.
One of the main messages of Table 1 is that the QCD corrections are negative as for similar signatures such as single Higgs-production via VBF or VBS at the LHC. In addition, the higher-order QCD corrections dramatically reduce the uncertainty associated with missing QCD higher orders. In particular, it goes from [+ 10.5%, − 8.8%] at LO to [+ 0.3%, − 0.06%] at NNLO in QCD. We note that the nonfactorisable NNLO QCD corrections are the only positive corrections, and that their contribution almost exactly cancels the factorisable NNLO QCD corrections. This is a coincidence of the particular cuts used here.
The second important point is the size of the EW corrections. It has recently been found (and further confirmed in Refs. [43,64]) that large EW corrections are an intrinsic feature of VBS at the LHC [49]. It originates from the quantum numbers of the particles involved in the process as well as the large scale induced by the massive t-channel exchange [65]. For such processes, the corrections reach about − 15 to − 20% of the LO prediction. On the other hand, for single-Higgs production via VBF, EW corrections have been found to be around − 5% [38,46]. It is thus interesting to observe that, despite having a higher typical scale, the magnitude of the EW corrections for double-Higgs production via VBF is very close to the single-Higgs one. In particular, in VBS the typical scale (the invariant mass of the four leptons) is m 4 ∼ 390 GeV while the VBF case it is even larger with m HH ∼ 610 GeV. In the same way as in Ref. [49], one can derive a leading-logarithmic approximation using Ref. [66]. Because the quantum numbers of the Higgs boson, such as the effective EW Casimir operator (see Eq. (B.10) in Ref. [66]), are significantly smaller than the ones of the Z or W gauge bosons, the logarithm coefficients are reduced with respect to the VBS case. For example, the coefficient of the double logarithms, which is directly proportional to the effective EW Casimir operator, is smaller by about a factor two. This implies, that VBF does not feature intrinsic large EW corrections as VBS.
The QCD corrections on the other hand tend to be somewhat smaller for double-Higgs production compared to single Higgs. This is due to the larger energy transfer in the t-channel which leads to harder jets and a higher dijet invariant mass. This in turn means that fewer events are lost due to QCD radiation. Overall, the state-of-the-art prediction displays a correction of about − 15% with respect to the LO prediction. Finally, the numerical value of the correction factor is K full/VBF = 0.99220(11) at the level of the fiducial cross section. It means that for the (rather exclusive) fiducial volume chosen here, the VBF approximation is reliable below the per-cent level. As shown later, this correction factor is not constant over the kinematic range and thus motivates its incorporation in our final predictions.
In Fig. 2, several transverse-momentum distributions are shown. The first two are the transverse momentum of the hardest jet and second hardest jet. Comparing the EW corrections to the QCD ones, we observe that their shapes are rather similar. These corrections are driven by Sudakov logarithms that grow negatively large towards the high-energy region. In this case, high energy refers to the high transverse momenta. For the hardest jet, it goes to about − 25% at 400 GeV, while for the second hardest jet, the value − 25% is reached already at 200 GeV. Concerning the transverse momenta of the Higgs bosons, ordered by their transverse momentum, the overall behaviour is similar and the corrections grow smoothly towards higher momenta. It is interesting to notice that in this case the difference between the full and the VBF computation is of 10-20% at 400 GeV. This is the most striking effect of the VBF approximation that we have observed in the present study. But it happens in a rather suppressed part of the phase space. For other distributions, the results are very much in line with Ref. [55] where the difference between the full and the VBF computation have been found to be small.
The rapidity distributions of the hardest jet and Higgs boson are shown in Fig. 3. For the rapidity distribution of the hardest jet, the NLO QCD corrections go from − 20% at 0 rapidity to about + 4% at 4.5 rapidity. The NNLO QCD corrections are smaller and vary within 4% across the phase space. The NLO EW corrections, on the other hand, hardly change over the displayed rapidity range and essentially inherit the overall normalisation. The rapidity distribution of the hardest Higgs boson displays even more stable corrections. The NLO QCD corrections fluctuate by less than 5% over the whole spectrum and the NNLO QCD corrections are also very stable (few per cent variation). The NLO EW corrections show a small shape distortion at the per-cent level. In both cases, the difference between the full and the VBF computation is minimal and do not exceed few per cent. Note also that for these distributions, the NNLO non-factorisable Table 1 The fiducial cross section for the process pp → HHjj, expressed in fb and in per cent, computed according to Eq. (4) at 14 TeV and under the selection cuts given in Sect. 2. The numbers in per cent are with respect to the LO cross section. The errors given in parenthesis are purely statistical whereas the additional uncertainties quoted for σ full LO and σ NNLO QCD×NLO EW are the QCD scale variations. We also show δ NF NNLO QCD separately. The value of the correction factor to go from the VBF approximation to the full computation is K full/VBF = 0.99220(11) Finally, in Fig. 4 the invariant-mass and transversemomentum distributions of the di-Higgs system are shown. For the invariant mass, the EW corrections are positive (above 5%) at 200 GeV. In the second and third bin, the maximum of the distribution is reached. The corrections become then negative about − 10% at 2 TeV. The QCD corrections on the other hand are rather flat and largely inherit the overall normalisation. In the tail of the distribution at about 2 TeV, the VBF computation start to diverge from the full one at a level of 5%. The transverse momentum of the two Higgs bosons displays much larger shape distortions. The NLO QCD corrections are maximal at low and high transverse momentum and are minimal around 200 GeV. This originates from kinematic constraints at low transverse momentum i.e. in the first three bins (the cut on the jet transverse momenta being 25 GeV). With higher-order QCD radiation, this constraint is relaxed allowing a relative increase in the cross section. At high transverse momentum, the running of the strong coupling is then taking place. Such an effect has already been observed for single Higgs production [32] and is also visible (even if less pronounced) in the transverse-momentum distributions of the Higgs boson in Fig. 2c and d. The NNLO QCD corrections slowly increase towards high transverse momentum to + 5% at 600 GeV. The EW corrections on the other hand displays a typical Sudakov behaviour. The corrections become rather large ( −25%) at 600 GeV. But such an energy corresponds to a rather extreme part of the phase space which is suppressed by more than three orders of magnitude with respect to the maximum of the distribution. Finally, the invariant mass and the rapidity difference of the two tagging jets are also displayed. These are typical observables used by experimental collaborations to enhance the VBF signal. Going to high invariant mass, the EW corrections slightly increase to reach about − 10%. This value is also obtained for larger rapidity difference of the two tagging jets. Such a kinematic (large invariant mass and large rapidity difference) is the typical VBF kinematic meaning that making the phase space cuts even more exclusive would increase the EW corrections but not dramatically. Overall, both at the level of the cross section and for differential distribution, the findings regarding the EW corrections are very similar to the ones for single Higgs production [38].
Finally, a few remarks regarding the non-factorisable NNLO QCD corrections. For most observable shown here the non-factorisable NNLO QCD corrections are of roughly the same size as the factorisable NNLO QCD corrections, although shapes differ for most observables. This warrants their inclusion when high-precision predictions are necessary. However, one has to be careful as the non-factorisable corrections are here computed in the eikonal approximation, which is only valid whenever all transverse momentum scales are much smaller than the partonic centre-of-mass energy. In particular this means that the approximation breaks down whenever the transverse momentum of a jet becomes large. As can be seen in Fig. 2a eikonal approximation is expected to be valid as discussed in Ref. [27].

Conclusion
The high-luminosity phase of the LHC will allow to probe rare processes, giving insight into fundamental interactions at higher energies. One of these rare processes is the production of two Higgs bosons through VBF. To that end, SM predictions are critical for the measurements of the process and associated searches for new physics phenomena.
In the present article we have provided state-of-the-art predictions within the SM in a realistic experimental setup at the LHC at 14 TeV. They feature the full LO predictions wit NLO QCD and EW corrections without relying on any approximation. These have then been combined with the existing NNLO QCD corrections to obtain the first predictions at NNLO QCD × NLO EW. The NLO EW corrections are presented for the first time here as well as the correspond-ing state-of-the-art predictions. Also, some of the differential distributions are shown at NNLO QCD here for the first time.
We have found that the EW corrections display the typical Sudakov behaviour in the high-energy limits. The corrections are of − 6% for the fiducial cross section while they can typically grow up to − 10 to − 20% in differential distributions. In general, the corrections are rather similar to the ones for single-Higgs productions and do not display very large EW corrections. The situation is thus rather close to the one of QCD corrections that largely share similarities between single-and double-Higgs production.
For the event selection chosen here, we predict a fiducial cross section of 0.67 fb. This amounts to a correction factor of − 14.8% with respect to the LO predictions. At the level of the differential distributions, the corrections can be larger and give rise to shape distortions of the order of 40%/50%. We note that for the High Luminosity LHC with an expected integrated luminosity of 3000 fb −1 this amounts to more than 2000 events in the rather strict fiducial volume used here.
Finally, we have also analysed the differences between the full computation including s-channel contributions and the VBF approximated one. This is also the first time that such results are shown in that respect. They are in rather good agreement with previous findings for single-Higgs production. In line with previous studies for similar processes, we have found that the VBF approximation becomes unreliable in rather inclusive set-ups or in extreme regions of the phase space.
The present results provide detailed information regarding higher-order corrections in double-Higgs production via VBF. It should thus be used as a guideline by the experimental collaborations in their quests for the measurement of this process during the high-luminosity phase of the LHC. It can also serve as a reference in the corresponding newphysics searches for future collider experiments that will help us unravel the mysteries of fundamental physics. edges support from the Swiss National Science Foundation (SNF) under contract BSCGI0-157722. The research of MP has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (Grant agreement no. 683211).

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Predictions derived from the calculation presented in this article can nonetheless be obtained from the authors upon request].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .