Vector-Boson Fusion Higgs Pair Production at N$^3$LO

We calculate the next-to-next-to-next-to-leading order (N$^3$LO) QCD corrections to vector-boson fusion (VBF) Higgs pair production in the limit in which there is no partonic exchange between the two protons. We show that the inclusive cross section receives negligible corrections at this order, while the scale variation uncertainties are reduced by a factor four. We present differential distributions for the transverse momentum and rapidity of the final state Higgs bosons, and show that there is almost no kinematic dependence to the third order corrections. Finally we study the impact of deviations from the Standard Model in the trilinear Higgs coupling, and show that the structure of the higher order corrections does not depend on the self-coupling. These results are implemented in the latest release of the proVBFH-incl program.


I. INTRODUCTION
Following the discovery of the Higgs boson in 2012 [1,2], it has become a primary focus of the experimental program of the Large Hadron Collider (LHC) to measure its properties. In particular, the measurement of the selfcoupling of the Higgs boson will be crucial both to further our understanding of the electroweak symmetry breaking mechanism, and to constrain possible new physics beyond the Standard Model (SM).
The simplest process with sensitivity to the trilinear Higgs coupling at hadron colliders is the Higgs pair production process, which has already been the focus of significant experimental studies [3][4][5][6][7][8][9][10][11][12][13][14][15]. Due to the low cross sections, production rates at the LHC are very small. For this reason, processes with two final state Higgs bosons are posed to play a key role at the high energy LHC (HE-LHC) and a future 100 TeV circular collider (FCC) in probing the Higgs sector. It is therefore important to have precise theoretical predictions for the dominant channels.
As for single-Higgs production, the leading contribution at the (HE-)LHC comes from gluon-gluon fusion [16]. This has been calculated up to next-to-next-to-leading order (NNLO) [17,18] matched to threshold resummation at next-to-next-leading logarithmic (NNLL) accuracy [19], and including finite top mass effects [20].
In this article, we focus on the vector-boson fusion (VBF) Higgs pair production channel, shown at leading order in figure 1. While it is only the second largest channel after gluon-gluon fusion, the VBF production mode is one of particular interest for several reasons: the presence of two tagging jets allows for a significant reduction of the large backgrounds through an appropriate choice of cuts; it is particularly sensitive to deviations from the SM in the trilinear Higgs coupling [21]; and is also a promising channel for measurements of the hhV V quartic coupling at the LHC [22].
Because of the important role that double Higgs pro-duction via VBF will play at the LHC and beyond, substantial efforts have been made to calculate its cross section to high accuracy. The differential cross section has been calculated up to next-to-leading order (NLO) [21,23] with matching to parton shower [24], and up to next-to-next-to-leading order when integrating out all hadronic final states [25]. We present here the calculation of di-Higgs production up to next-to-next-to-next-to-leading order (N 3 LO), which is also the first calculation at this order of a 2 → 4 process. Together with a companion paper presenting the fully differential NNLO calculation [26], this brings the VBF double Higgs channel to the same theoretical accuracy as single-Higgs VBF production [27][28][29][30]. These results are obtained using the structure function approach [31], which is the limit in which there is no partonic exchange between the two protons, and in which all radiation has been integrated over. Since the single-gluon exchange is zero for color reasons, this approximation is exact at NLO, while it has been shown to be accurate to more than 1% at NNLO for the single-Higgs process [32][33][34]. Because the presence of an additional Higgs boson does not impact the color flow between the hadrons, this limit is expected to be just as valid for Higgs pair production.
This paper is structured in the following way: In section II we present the details of our calculation, in section III we present results for the inclusive cross section, while differential distributions are given in section IV. We give our conclusions in section V.

II. HIGGS PAIR PRODUCTION IN VBF
We start by setting up the formalism needed to calculate the inclusive cross section up to third order in the expansion in the strong coupling constant, which is analogous to the single-Higgs one.
lated as a double deep inelastic scattering (DIS) process, and can be written as [31] Here G F is Fermi's constant, m V and ∆ 2 V are the mass and squared propagators of the mediating W or Z bosons, and √ s is the collider center-of-mass energy. We defined Q 2 i = −q 2 i and x i = Q 2 i /(2P i · q i ) as the usual DIS variables, where q i is the four-momentum of the vector boson V i and P i that of the initial proton. Finally W V µν is the hadronic tensor and dΩ VBF is the four particle VBF phase space. The matrix element of the V V → hh sub-process is expressed as [35] where k 1 , k 2 are the final state Higgs momenta, which satisfy k 1 + k 2 = q 1 + q 2 , λ is the trilinear Higgs selfcoupling and ν is the vacuum expectation value of the Higgs field. (1) is given by where the F V i (x, Q 2 ) functions are the standard DIS structure functions with i = 1, 2, 3, which can be expressed as a convolution of the parton distribution functions (PDF) with the short distance coefficient functions To evaluate equation (4), it is useful to define the singlet and non-singlet distributions q S , q NS,i , as well as the nonsinglet valence distribution q V NS and the asymmetry δq ± We can then decompose the quark coefficient functions into non-singlet and pure-singlet parts, and define the valence coefficient function The neutral current structure functions can now be ex- pressed as . The vector and axial-vector coupling constants v Z i and a Z i are given by For the charged current case, the structure functions can be written as where we have again , and the couplings are simply a W j = v W j = 1 √ 2 . We can calculate corrections up to N 3 LO by making use of the known three-loop coefficient functions [36][37][38][39][40], whose parameterized expressions have been implemented in HOPPET v1.2.0-struct-func-devel [41]. Examples of three-loop diagrams included in this calculation are shown in figure 2. PDF4LHC14_nnlo_mc To calculate the variation of the cross section with different choices of factorization and renormalization scales, we compute the scale dependence to third order in the coefficient functions as well as in the PDFs.
We start by evaluating the running coupling for α s where we defined β 0 = (33 − 2n f )/12π and β 1 = (153 − 19n f )/24π 2 . The coefficient functions can then easily be expressed as an expansion in α s (µ R ). To evaluate the dependence of the PDFs on the factorization scale µ F , we can integrate the DGLAP equation, using Expressing the PDF in terms of an expansion in α s (µ R ) evaluated at µ F , it then straightforward to evaluate equation (4) for any choice of the renormalization and factorization scales up to N 3 LO.
To estimate the theoretical uncertainty due to missing higher order corrections, we calculate the envelope of seven different scale choices, taking where we keep 1 2 ≤ µ R /µ F ≤ 2 and µ 0 is the central scale choice. We set the central renormalization and factorization scales to the vector boson virtuality of the corresponding sector, Q 1 or Q 2 . PDF4LHC15_nnlo_mc  For the numerical integration, we use the phase space parameterization of VBFNLO [42]. Unless otherwise specified the center-of-mass energy is set to the expected energy of the HE-LHC, which is 27 TeV. For all simulations, we use the PDF4LHC15 nnlo mc set [43] with a four-loop evolution of the strong coupling, starting from an initial condition α s (M Z ) = 0.118. We set the mass of the Higgs boson to m H = 125 GeV. The electroweak parameters are set to the PDG values [44], with m W = 80.379 GeV, m Z = 91.1876 GeV and G F = 1.16637 · 10 −5 GeV −2 . The narrow-width approximation is used for the final state Higgs bosons, while Breit-Wigner distributions are used for internal bosons, taking Γ W = 2.085 GeV, Γ Z = 2.4952 GeV, and Γ H = 4.030 · 10 −3 GeV.

III. TOTAL CROSS SECTION
We start by providing results for the inclusive cross section.
In figure 3, we show the dependence of the total cross section on the renormalization and factorization scales for each order in QCD. One can clearly observe the conver-σ (14 TeV) [fb] σ (27 TeV)   We show the dependence of the total cross section as a function of center-of-mass energy in figure 4. Here we see that at even very high energies, the third order corrections are fully contained within the NNLO scale variation bands, with an almost constant K-factor. One should note that this is somewhat dependent on the choice of central scale, and less dynamical scales such as an m h or p t,HH based prescription will lead to third order corrections that can deviate from the NNLO uncertainty bands in certain kinematic regions or at sufficiently high energies.
We detail the precise value of the cross section and its scale variation uncertainties in table I. Values are given for three reference center-of-mass energies: the 14 TeV LHC, the 27 TeV HE-LHC and the 100 TeV FCC. For each of these energies, we provide inclusive cross sections at each order in perturbative QCD, along with the corresponding scale variation envelope. We can observe that while the corrections are at the level of a few permille only, the scale uncertainty bands are reduced by more than a factor four when going from NNLO to N 3 LO.
A comment is due on the impact of contributions beyond those included in the DIS limit. There are a number of corrections to the Born diagrams shown in figure 1 beyond those due to the radiation of additional partons. These should be included where possible for precise phenomenological predictions.
In particular, the s-channel production mode, while suppressed to a few permille after VBF cuts, contributes to about 16% to the total cross section for 27 TeV collisions, and can therefore not be neglected. It can be calculated to NLO using the MadGraph5 aMC@NLO framework [45] and can be straightforwardly included.
Furthermore, NLO electroweak corrections are currently unknown and expected to be sizeable. They can be estimated from dominant light quark induced channels using Recola(Collier)+MoCaNLO [46][47][48][49] for the di-Higgs and single-Higgs VBF process, comparing the lat- ter to HAWK [50]. For VBF Higgs pair production the EW corrections to the inclusive cross section lie between −5% and −7%. Compared to the single-Higgs VBF correction of roughly −5% (using the same set-up, i.e. excluding photonic and b-quark channels), the double Higgs VBF process thus does not seem to receive large VBS-like corrections [51,52]. One can therefore expect the full electroweak corrections to be at least at the same level as the NNLO QCD corrections, and significantly larger than the N 3 LO corrections.
There are also a number of α 2 s and α 3 s contributions that are neglected in the structure function approximation, notably: the double and triple gluon exchange between the two quark lines; heavy-quark loop induced production; t-/u-channel interferences; single-quark line contributions; loop induced interferences between VBF and gluon-fusion Higgs production. These have been shown to contribute at the few permille level in single Higgs VBF production [32][33][34]53], and we therefore expect that they can be neglected.
The impact on the cross section of PDF and α s uncertainties can be evaluated using the PDF4LHC15 nnlo mc pdfas set, and is of about 2.1%. Finally, there is also a theoretical PDF uncertainty, due to missing higher orders in the determination of the PDFs. In this paper we use an NNLO pdf set to evaluate an N 3 LO cross section, since N 3 LO sets are currently unavailable. The uncertainty due to these missing higher order terms come from two sources. They are dominated by missing third order corrections to the coefficient functions relating physical observables to PDFs, and can be estimated to about 8 using the method presented in [29]. The second source of corrections is due to unknown four-loop splitting functions [54] appearing in the DGLAP evolution, which have been estimated to be negligible [29].

IV. DIFFERENTIAL DISTRIBUTIONS
The calculation described in section II is inclusive over final state QCD radiation. One can thus not obtain differential predictions with respect to the jet kinematics without using the projection-to-Born method [28] and combining it with a higher multiplicity NNLO prediction. However, we have full access to the kinematics of the Higgs bosons, and it is therefore straightforward to compute differential observables with respect to the their momenta. Let us now focus on several differential distributions of particular interest. We will again consider here a 27 TeV proton-proton collider except where otherwise specified.
In figure 5, we show the transverse momentum p t,HH , rapidity y HH and invariant mass m HH distributions of the Higgs pair for each order in QCD. The latter is of particular interest, as the Higgs pair invariant mass can  be notably sensitive to deviations due to physics beyond the SM [22,55,56]. We see that once we get to the third order, there is almost no kinematic dependence to the K-factor, except at very high rapidities, where the N 3 LO corrections can bring changes to the central value of about one percent. The N 3 LO scale variation bands are always fully contained within the NNLO scale uncertainties, but are about four times thinner.
We order the Higgs bosons according to their transverse momentum. Figure 6 provides the transverse momentum distribution of both the harder (p t,H1 ) and the softer (p t,H2 ) Higgs. The third order corrections to these observables are negligible, with again a large reduction in scale uncertainties. The corrections to the rapidity distributions of the two Higgs bosons are shown in figure 7. We note here that the α 3 s contribution has almost no kinematic dependence up to rapidities of |y H | ∼ 4, with the scale variation bands being again fully contained by the theoretical uncertainties of the previous order.
Finally, let us study the impact of the trilinear Higgs self coupling, λ, by varying the corresponding factor in equation (2). Constraining the trilinear coupling is of particular interest, since many scenarios of new physics beyond the SM predict significant deviations of this value. Examples of such models are SO(5)/SO(4) minimal composite Higgs [57,58] and dilaton models [59].
To study the impact of deviations of this type, we define λ = κλ SM , where λ SM = m 2 h /2v, and consider a range of values for κ.
The total cross section up to N 3 LO as a function of κ is given in figure 8, both for a 27 TeV HE-LHC and for a 100 TeV FCC. One can observe that while the inclusive cross section changes by several orders of magnitude, there is as expected almost no dependence of the higher order corrections on κ beyond leading order.
In figure 9, we show the kinematics of the Higgs pair for several values of κ. We see that even very small deviations in the trilinear Higgs coupling have a substantial impact on the cross sections, both in the normalization and shape of the distributions. In particular, the rapidity and invariant mass of the Higgs pair are particularly sensitive to changes in λ.
The lower panels in figure 9 show the ratio to the central value obtained with κ = 1, with the leading order predictions shown as dashed lines. One can see that the changes to the cross sections from variations in κ can be substantial. However there is essentially no change in the N 3 LO/LO K-factor.

V. CONCLUSIONS
In this article, we have completed the first N 3 LO calculation of a 2 → 4 process, namely the production of a Higgs boson pair through VBF. This calculation was made possible by the factorizable nature of the higher order QCD corrections. Together with the fully differential NNLO calculation presented in a companion paper [26], this brings the di-Higgs production channel to the same theoretical accuracy as has been achieved for the single-Higgs process, opening up the prospect of precision studies of the Higgs sector through Higgs pair production at the HE-LHC and at future hadron colliders.
We have presented differential distributions of the Higgs pair transverse momentum, rapidity and invariant mass for the 27 TeV HE-LHC. The corrections are at the few permille level, however the calculation of the third order leads to a substantial reduction in scale uncertainties. The convergence of the perturbative series is very stable at this order, with almost no kinematic dependence to the N 3 LO corrections, except at very high rapidities. The N 3 LO scale variation bands are always fully contained within the second order scale uncertainties, but are over four times thinner.
Finally we studied the impact of deviations from the SM in the trilinear Higgs self-coupling on the N 3 LO distributions. Small deviations of this constant can substantially change both the total cross section and the shape of the distributions. However the structure of the higher order QCD corrections is unaffected by variations in the coupling, with the N 3 LO/LO K-factor staying constant over a range broad range of λ values. The results presented here have been implemented in the version 2.0.0 of proVBFH-incl [60] which provides predictions for both single and double Higgs inclusive cross sections up to N 3 LO in QCD.
This article provides also the first element for a fully differential N 3 LO calculation of VBF Higgs pair production. This could be achieved by combining the present inclusive calculation with a differential NNLO computation of the electroweak production of two Higgs bosons in association with three jets.