On the impact of non-factorisable corrections in VBF single and double Higgs production

We study the non-factorisable QCD corrections, computed in the eikonal approximation, to Vector-Boson Fusion single and double Higgs production and show the combined factorisable and non-factorisable corrections for both processes at Oαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({\alpha}_s^2\right) $$\end{document}. We investigate the validity of the eikonal approximation with and without selection cuts, and carry out an in-depth study of the relative size of the non-factorisable next-to-next-to-leading order corrections compared to the factorisable ones. In the case of single Higgs production, after selection cuts are applied, the non-factorisable corrections are found to be mostly contained within the factorisable scale uncertainty bands. When no cuts are applied, instead, the non-factorisable corrections are slightly outside the scale uncertainty band. Interestingly, for double Higgs production, we find that both before and after applying cuts, non-factorisable corrections are enhanced compared to the single Higgs case. We trace this enhancement to the existence of delicate cancellations between the various leading-order Feynman diagrams, which are partly spoiled by radiative corrections. All contributions studied here have been implemented in proVBFH v1.2.0 and proVBFHH v1.1.0.


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 and in particular its couplings to itself and to the other Standard Model particles [3].One of the key channels for studying the Higgs boson is the Vector-Boson Fusion (VBF) production mode, where the Higgs boson is produced together with two (typically) hard and forward jets.This process has been the focus of several recent fixed order theoretical calculations [4][5][6][7][8][9][10].
A common point between all these calculations is that they are performed in the factorised approximation, which corresponds to the limit where partons from the two colliding protons are treated as coming from two identical copies of QCD that interact exclusively through the electroweak sector.When all emissions are integrated over, this approximation is referred to as the structure function approach [11].Due to colour conservation, this approach is exact up to NLO, but starts to be violated from NNLO onwards, where colour-singlet two-gluon exchanges between the incoming partons are neglected.Since these non-factorisable contributions are colour suppressed compared to their factorisable counterparts, it has generally been assumed that they can also safely be neglected [5].
Recently it has been shown that the impact of the nonfactorisable corrections at NNLO can be estimated in the so-called eikonal approximation [12].Although this calculation confirms that their impact is moderate, it was found that these contributions also receive a π 2 -enhancement due to their connection to the Glauber scattering phase, which partially overcomes the effects of colour suppression.
Given these findings, the purpose of this paper is twofold.Firstly, we investigate the validity of the approximation employed in Ref. [12] for single Higgs production outside of tight VBF cuts, in order to estimate the leading non-factorisable corrections on the inclusive VBF cross section.We then conduct an in-depth phenomenological study of the factorisable and non-factorisable corrections, and establish the relative impact of the latter for a range of selection cuts and observables.Secondly, we extend the calculation of Ref. [12] to study the impact of non-factorisable corrections to the production of a pair of Higgs bosons in VBF.In this case, contrary to single Higgs, it is well known that the rather small LO cross-section is the result of delicate cancellations of more than one order of magnitude between the different Feynman diagrams that contribute to the process, shown in figure 5.While QCD radiative corrections in the factorisable approximation affect equally all Born diagrams and are not expected to spoil this cancellation, the same cannot be expected a priori for the non-factorisable ones.As we will demonstrate in this paper, unitarity ensures that non-factorisable radiative corrections preserve the delicate pattern of cancellations observed at leading order, leading to overall suppressed contributions, both at the inclusive and at the differential level.
In figure 1, we provide a summary of the impact of O(α 2 s ) corrections to single Higgs VBF production as a function of the selection cuts on the rapidity separation ∆y jj and the invariant mass m jj .Figure 1a shows the ratio of the factorisable corrections to the LO cross section.In general they are suppressed by an order of magnitude compared to the factorisable corrections.In figure 2 we show the same comparison but for di-Higgs production.As can be seen in figure 2a the factorisable corrections have a more complicated dependence on the cut values compared to single Higgs VBF production, first decreasing with an increase in both cuts and then finally increasing in size as both cuts become large.However, the non-factorisable corrections shown in figure 2b have a similar kinematical dependence to the single-Higgs case, although they are somewhat suppressed as is the case for their factorisable analogue.
We note that in addition to the non-factorisable corrections studied in this paper, a number of known perturbative corrections to VBF Higgs production are usually neglected.These include t/u-channel interference and s-channel contributions [13], single-quark line contributions [14], and loop induced interferences between VBF and gluon-fusion Higgs production [15].These corrections are small within typical VBF cuts and we do not consider them here.The NLO corrections in the electroweak coupling have also been studied in Ref. [13].
The rest of the paper is structured as follows: in section 2 we provide a review of the known QCD corrections to VBF single Higgs production, and describe how to perform a similar estimate of the non-factorisable corrections to di-Higgs production in the eikonal approximation.In section 3 we compare factorisable and non-factorisable corrections for VBF single Higgs production in a realistic setup.In section 4 we discuss the impact of the non-factorisable corrections to VBF di-Higgs production.In section 5 we give our conclusions.

QCD corrections in VBF Higgs production 2.1 Factorisable corrections
Both in single and double Higgs production via VBF, the Higgs bosons are emitted by the electroweak vector bosons exchanged between the two scattering partons.Schematically, the Born process for the emission of an arbitrary number of Higgs bosons can be depicted as in figure 3a.
In the factorised approximation, the VBF cross section is then expressed as a double deep inelastic scattering (DIS) process, see Fig. 3b, for which the cross section is given by [11] Here V = W ± , Z corresponds to the mediating boson with mass m V and squared propagator is the hadronic tensor and dΩ VBF is the VBF phase space.
The matrix element of the vector-boson fusion sub-process is denoted as M V,µν .The hadronic tensor can be expressed as Pi,µ Pi,ν where we have defined Pi,µ = P i,µ − Pi•qi q 2 i q i,µ and F V i (x, Q 2 ) are the standard DIS structure functions with i = 1, 2, 3.
For single Higgs production, given by the diagram T in Fig. 4, M V,µν can be written as By using the known DIS coefficient functions up to order α 3 s [16][17][18][19][20], this can be used to evaluate the inclusive VBF cross section to single Higgs production up to N 3 LO in the factorised approximation.By combining an inclusive NNLO calculation with the corresponding fully differential NLO prediction for electroweak Higgs production in association with three jets [6], one can obtain fully differential results at NNLO through the projection-to-Born method [7] or the antenna subtraction method [8].
The factorisable QCD corrections to the di-Higgs process can be calculated in the same way as for its single Higgs analog, expressing the cross section in the form of eq. ( 1), but with M now referring to the di-Higgs matrix element.The Higgs pair production process differs from the single Higgs case only at the interaction between the vector and Higgs bosons, where an additional Higgs can arise from an intermediate vector or Higgs boson, or from the hhV V quartic coupling.The V V → hh sub-process at LO can be expressed as [21] where we have defined the propagators and k 1 , k 2 are the momenta of the final state Higgs bosons and λ and ν are the trilinear Higgs self-coupling and the vacuum expectation value of the Higgs field respectively.The matrix element arises from the four Feynman diagrams shown in Figure 5, which we label T 1 , T 2 , B 1 and B 2 .We stress here that, while we are including the bosons' widths for completeness, they play no role for the estimation of QCD corrections to Higgs production in VBF.

Non-factorisable corrections
The factorisable approach described above, which includes diagrams such as the one represented in figure 3b, is exact up to NLO due to colour conservation.At NNLO this is no longer true, as in particular two gluons in a colour singlet state can be emitted between the two quark lines, as shown in figure 6.As the gluons have to be in a colour singlet state, these diagrams will be colour suppressed compared to their factorisable counterparts.For this reason it has long been argued that they can be neglected when considering NNLO corrections to VBF [5].
Due to the complexity involved in computing the two-loop non-factorisable corrections, very little has been known about them beyond the fact that they are colour suppressed.However, very recently [12] significant progress was made, when it was shown that the corrections can be estimated within the eikonal approximation [22][23][24][25].This calculation exploits the fact that when typical VBF cuts are applied, the VBF cross section can be expanded in the ratio of the leading jet transverse momentum over the total partonic centre-of-mass energy In this kinematical configuration, the authors of Ref. [12] conclude that the non-factorisable corrections receive a π 2enhancement connected to the presence of a Glauber phase, which can partially compensate their colour suppression.Indeed, it turns out that for VBF single Higgs production, the non-factorisable corrections can contribute up to 1% in certain regions of phase space, making them larger than the factorisable N 3 LO corrections.In what follows we will use the same approximation to estimate the impact of non-factorisable corrections for the case of double Higgs production as well.
In order to see how the NNLO non-factorisable corrections can be estimated in the eikonal approximation both for single and double Higgs production, let us consider a generic VBF Born diagram, which we will call D, for the production of an in principle arbitrary number of Higgs bosons, see Fig. 3a.In what follows this diagram will represent either the Born diagram for VBF single Higgs production T of Fig. 4, or any of the Born diagrams for double Higgs production T 1 , T 2 , B 1 or B 2 in Fig. 5.
It is important to stress here that, somewhat counterintuitively, we will be considering QCD corrections on each single diagram separately, and not on the full Born matrix element.Since we are interested in computing the NNLO QCD corrections to this class of processes, we imagine dressing the diagram D with 1-loop or 2-loop QCD corrections, as depicted in Fig. 6, where we provide two representative diagrams for illustration only.It turns out that, at least up to two loops in QCD, we can limit ourselves to diagrams where the gluons are in a colour-singlet configuration, i.e. exchanged between the two quark lines.All other configurations do not contribute to the cross-section due to colour conservation.Therefore, the calculation of the one-and two-loop QCD corrections in the eikonal approximation reduces effectively to the corresponding calculation in QED, with the colour-averaged effective coupling Following Ref. [12], let us consider the process q(p 1 ) + q(p 2 ) → q(p 3 ) + q(p 4 ) + X(P ) where X(P ) can represent one or multiple Higgs bosons produced in vector-boson fusion.At leading order, we call the momenta flowing in the two vector bosons respectively The leading term in the eikonal approximation can then easily be obtained by employing light-cone coordinates, which make transparent the separation between the dynamics in the plane spanned by the momenta of the incoming quarks and the plane transverse to them [22][23][24][25].For a momentum k µ we indicate by k ± the light-cone coordinates and by k those in the transverse plane, i.e. we write and we choose a reference frame such that the incoming quark momenta have each one light-cone component different from zero It turns out that both at one and two loops, at leading order in the eikonal approximation, the quark propagators coupled to the soft gluons simplify and, after summing over all permutations of the gluons and the vector bosons, the quark propagators recombine in terms of delta functions of the light-cone components of the loop momenta.This allows one to effectively decouple the light-cone dynamics from the one in the two-dimensional plane transverse to the momenta of the incoming quarks and one is left with the calculation of the effective two-dimensional loop diagrams shown schematically in Fig. 7. Extra care has to be taken when considering diagrams of type (c) and (d) in Fig. 5, where even in the eikonal approximation the lightcone components cannot be neglected in the propagator of the central vector boson.Their effect can nevertheless be effectively included by modifying the mass of the central vector boson to an effective mass M V → µ V , whose value depends on the light-cone components of the scattering quarks and of the Higgs boson. 1 We provide its explicit value later on, see Eq. ( 16).With this, one can easily write the one-and two-loop QCD corrections in the eikonal approximation in a rather compact form.
By calling q 1 and q 2 the transverse components of the momenta q 1 and q 2 in (9), and indicating schematically with {q, µ V } the set of transverse momenta of the Higgs bosons produced and the effective mass of the intermediate vector boson, as described above, we can write for the generic Born diagram D M (1)

M
(2) where M (n) D are the corrections to the Born diagram D coming from the exchange of n gluons, χ are functions which depend on the (transverse) kinematics of the corresponding Born diagram and on the effective mass µ V , and the effective coupling α s was defined in eq. ( 7).Finally, the factor 1/2! comes from the symmetrisation of the two identical gluons [12].We stress once more that, if we are interested in double Higgs production, this happens separately for each of the Born diagrams in Fig. 5.We also remind the reader that this is true as long as we limit ourselves to colour-singlet gluon exchange.Given the considerations above, it is easy to see that QCD corrections to the Born diagram of single Higgs production T , or to T 1,2 for double Higgs, reduce to the computation of a two-dimensional one-or two-loop trianglelike integral, while the corrections to B 1,2 , involve the computation of more complicated box-like loop integrals.Moreover, it should also be clear that for T , T 1 and T 2 , the QCD corrections only depend on the momenta q 1 , q 2 and are therefore equal in all three cases.Putting everything together, we find similarly to Ref. [12] χ (1) where we defined k 12 = k 1 + k 2 and have introduced a fictitious gluon mass λ to regulate the residual IR divergences.Also, we have removed the dependence on the momenta {q, µ V } since for these diagrams q 1 + q 2 = q and there is no vector boson propagator that depends on µ V .
Let us consider now the two box-like topologies, which have a non-trivial dependence on the momenta of the two Higgs bosons.Calling their momenta q 3 and q 4 and using q 4 = −q 1 − q 2 − q 3 (all momenta are incoming), we find where we put q ij = q i + q j and defined in addition the "transverse-plane" Mandelstam variables s = (q 1 + q 1 ) 2 , t = (q 1 + q 3 ) 2 , u = (q 1 + q 4 ) 2 , as well as Similarly to the previous case, we regulated the residual IR divergences with a gluon mass λ.Notice that in defining χ B2 we need to swap the entire four-momenta q µ 3 ↔ q µ 4 , which also affects the value of µ V as defined in Eq. ( 16).
The integrals above can be computed in many different ways, most notably making use of the reduction of all 1-loop and certain 2-loop n-point functions, with n ≥ 3, in d = 2 space-time dimension, to lower-point topologies.Also a direct computation of the integrals using their Feynman parameter representation can be attempted, which turns out to be particularly simple for the triangle integrals χ (1) T (q 1 , q 2 ) and χ (2) T (q 1 , q 2 ), see Ref. [12].While the analytic computation is conceptually straightforward, the result, in particular for what concerns the box-type integrals, can become very cumbersome due to their dependence on a large number of scales and are not particularly illuminating.
Nevertheless, since we are dealing with two-dimensional euclidean integrals, it turns out to be entirely straightfor-ward to produce very compact one-fold integral representations for them by extracting the logarithmic divergences as λ → 0 and integrating directly on the 2-dimensional loop momenta in polar coordinates.This remains true at two loops, where one can first integrate out the gluonic one-loop sub-bubble, and then proceed in the very same way as for the one-loop integrals.This allows us to get all results as one-fold integrals over simple algebraic functions and at most powers of logarithms.
We write down the results for the one-and two-loop triangles as and similarly for the boxes where the function f T and f B depend on the corresponding transverse momenta and χ (j) B2 can be obtained from χ (j) B1 by swapping q µ 3 ↔ q µ 4 as in Eq. ( 15) and Eq. ( 16).In order to write their analytic expression, we start off by parametrising the kinematics in the two-dimensional transverse plane as q 1 = (q 1x , 0) , q 2 = (q 2x , q 2y ) , q 3 = (q 3x , q 3y ) (19) with q 4 = −q 1 − q 2 − q 3 , and we introduce the shorthand notation The functions can then be written as follows f (1) with r ij = r i − r j , rj = −r j /M 2 V and the six roots read where r * j indicates complex conjugation and The complex arguments and the three square-roots above might seem somewhat unappealing, in particular because lengthy but fully analytic representation can be obtained for all these functions in terms of polylogarithms.Nevertheless, our results involve only integrals of logarithms and exhibit a very high degree of symmetry, both moving from one to two loops and going from 3-to 4point functions.Moreover it is straightforward to rewrite the integrals to make them explicitly real, at the price of introducing inverse trigonometric functions.Finally, as a curiosity, it turns out that performing the calculation in this way the results can be effortlessly generalised to higher-point integrals, i.e. for an arbitrary number of Higgs bosons in the final state.
With the definitions above, the non-factorisable QCD corrections to the total amplitude for single and double Higgs production can be written, respectively, as where for single Higgs we have simply while for double Higgs we find which of course implies a much richer interference pattern.
More explicitly, we find for the cross-section for single Higgs production where dσ LO is the leading-order cross section given in (1), α s is the effective coupling in eq. ( 7), and the NNLO nonfactorisable contributions only depend on the functions f T .
As an illustration, and in order to compare this case to di-Higgs production, it is useful to compute the corrections in the limit where all transverse scales become small compared to the vector-boson mass, i.e. q 2 1,2 M 2 V .In that limit, all integrals become trivial and we find [12] In the case of double Higgs production, the form of the corrections is rather cumbersome but still entirely straightforward and we prefer to avoid writing down the formulas explicitly.On the other hand, if we consider the same limit as above, ie q 2 1,2 ∼ q 2 3,4 M 2 V , with the additional assumption that the Higgs bosons are produced with zero rapidity, the formulae simplify considerably.In order to present the result, we divide the LO cross-section in three contributions as where dσ LO T T is the contributions stemming solely from diagrams T 1 and T 2 , σ LO BB from B 1 and B 2 and σ LO T B from the interference of the two classes of diagrams, see Fig. 5.With this, we find that the nonfactorisable corrections at NNLO take the suggestive form dσ NNLO  HH,nf ∼ α2 where in this approximation H , and we have defined y = 1+z 2 z 2 and z = M V M H .For the numerical evaluation we have used M V = M W = 80.398 GeV and M H = 125 GeV.Eq. (33) shows that the three contributions to the Born cross-section for di-Higgs production can receive radiative corrections which are different at the 10% level.The cross-section for HH production at LO is the result of delicate cancellations of more than one order of magnitude between the three different contributions in eq. ( 32), as can be seen in table 1.These cancellations are a well known manifestation of the role that the Higgs boson has in restoring unitarity in the Standard Model.Since we are working in the eikonal approximation, one could wonder whether this approximation could spoil these cancellations and induce in this way artificially large NNLO QCD corrections on the di-Higgs cross-section.As a matter of fact, as we will see in detail below, the corrections conspire so to preserve the unitarity induced cancellations among the various components and produce rather small effects both at the inclusive and at the differential level.
In figure 8, we show the α2 s coefficient of each of the T T , T B and BB contributions as a function of a cut on the maximum transverse momentum, max(p t,j1 , p t,H1 ) < p t,max , with the only requirement that |y H,1 , y H2 | < 1.For technical reasons this plot only includes the charged-current sub-process.At small values of p t,max we can observe a convergence of these coefficients to the analytic expression given in eq.(33).Above p t,max ∼ M W there is a transition to different numerical values of the coefficients, leading to a potential spoliation of the cancellation present in the LO cross section.
Another interesting limit to study is the case where the Higgs bosons' transverse momenta are small, i.e. q 2 3 , q 2 4 M 2 V which implies that q 2 1 ∼ q 2 2 .In this limit, assuming again that the Higgs bosons are produced with zero rapidity, it is easy to show that the cross section becomes where the coefficients C T T , C T B and C BB now depend on the variable x = q 2 1 /M 2 V = q 2 2 /M 2 V , in addition to y defined above.As their analytical expression is rather lengthy, we prefer to provide them in appendix A.Here we notice that  Fig. 8: α2 s coefficient of the T T , T B and BB chargedcurrent contributions to the NNLO non-factorisable corrections, normalised to the corresponding LO term, as a function of a cut on the maximum transverse momentum, max(p t,j1 , p t,H1 ) < p t,max .this result reproduces eq. ( 33) when x → 0, as expected.On the other hand, as x grows the radiative corrections to the three contributions differ widely and in particular we see that as x → ∞ the cross section takes the form (35) In the limit where q 2 1 ∼ q 2 2 ∼ M 2 V , ie when x = 1 we obtain the numerical expression (36) This leads us to conclude that the non-factorisable corrections can grow very rapidly with x, ie with the transverse momentum of the jets.We note that the eikonal approximation is strictly speaking only valid when x ∼ 0 and that the rapid growth could be seen (at least in part) as a consequence of the breakdown of the approximation.
Before concluding this section, it is worth noting that at O(α 2 s ) there are both loop-induced and real emission di-agrams that contribute to the non-factorisable corrections discussed above.Nevertheless, it is well known that real emission diagrams do not contribute to leading order in the eikonal approximation, and the whole cross-section in this limit stems from the virtual contributions only.This is also demonstrated by the fact that IR divergences cancel between the two-loop and the one-loop squared amplitudes.We stress here that the real emission diagrams have been computed for the single Higgs case in [26], and could be used to compute non-factorisable corrections beyond the leading eikonal approximation, once the full two-loop amplitudes become available.
3 Results for single Higgs VBF production

Setup
In order to investigate the size of the various QCD corrections, we study 13 TeV proton-proton collisions, in a setup identical to Ref. [7].We use a diagonal CKM matrix, full Breit-Wigners for the W , Z and the narrowwidth approximation for the Higgs boson.We take the NNPDF 3.0 parton distribution functions at NNLO with α s (M Z ) = 0.118 (NNPDF30 nnlo as 0118) [27], as implemented in LHAPDF-6.1.6[28].We consider five light flavours and ignore contributions with top quarks in the final state or internal lines.We set the Higgs mass to M H = 125 GeV, compatible with the experimentally measured value [29].Electroweak parameters are set according to known experimental values and tree-level electroweak relations.As inputs we use M W = 80.398 GeV, M Z = 91.1876GeV and G F = 1.16637 × 10 −5 GeV −2 .For the widths of the vector bosons we use Γ W = 2.141 GeV and Γ Z = 2.4952 GeV.The central factorisation, µ F , and renormalisation, µ R , scales are set to when computing factorisable corrections.We compute the residual scale uncertainties by varying this scale up and down by a factor 2 keeping µ R = µ F , which was shown in Ref. [7] to encompass almost the same scale uncertainty bands as a full 7-point scale variation (i.e.where µ R and µ F are varied independently by a factor 2 with 1 2 ≤ µR µF ≤ 2).For the purpose of comparing these effects, we compute the non-factorisable corrections using the same central scale, which differs from the renormalisation scale choice µ R = √ p t,j1 p t,j2 in Ref. [12].The residual scale uncertainties for these last predictions have been obtained using the full 7-point scale variation.
In the following we will discuss results both fully inclusively in the VBF jets, and under a set of representative VBF selection cuts.To pass our VBF selection cuts, events should have at least two jets with transverse momentum p t > 25 GeV; the two hardest (i.e.highest p t ) jets should have absolute rapidity |y| < 4.5, be separated by a rapidity ∆y j1,j2 > 4.5, have a dijet invariant mass m j1,j2 > 600 GeV and be in opposite hemispheres (y j1 y j2 < 0).We define jets using the anti-k t algorithm [30], as implemented in FastJet v3.1.2[31], with radius parameter R = 0.4.

Validity of the eikonal approximation
The calculation of the non-factorisable NNLO corrections in VBF given in Ref. [12] is carried out as an expansion in ξ truncated to lowest order in ξ, see eq. ( 6).The authors argue that this ratio is typically of the order 1  6 , based on experimental measurements of the p t,j1 and m jj spectra [35,36], and hence that the relative error associated with truncating the power expansion at the leading order is roughly 1  36 .This analysis is performed under the VBF cuts given in sec.3.1 which guarantee large √ s because of the requirement on the invariant mass of the di-jet system.
In this section we investigate in some detail how robust this approximation remains when no cuts are applied to the jets.Although such an inclusive setup is not of much phenomenological interest, it is of theoretical interest, given that not only the factorisable NNLO corrections are known fully inclusively, but also the N 3 LO ones [9].
In the left panel of figure 9 we show the normalised VBF cross section integrated in ξ, defined as We show the cross section fully inclusively and under VBF cuts.Under VBF cuts the cross section clearly lives below ξ ∼ 0.2, whereas the fully inclusive cross section receives contributions all the way up to ξ ∼ 0.4.However almost 85% of the events have ξ < 0.2 implying that the approximation used in Ref. [12] is valid in a large region of the inclusive VBF phase space.In fact, the average value of ξ is below 0.12 for all values of the rapidity of the Higgs Boson and for moderate transverse momenta, p t,H , as can be seen in the two right panels of figure 9.At large p t,H , the average value of ξ increases almost linearly with p t,H and hence the eikonal approximation starts to break down.This is not unexpected as p t,H is balanced by the jet transverse momenta, which by definition have to take moderate values in order to keep ξ small.In the region of phase space where the eikonal approximation breaks down, i.e. when ξ becomes large, ( 29) is no longer valid.However, in this region the non-factorisable corrections are not expected to receive a Glauber phase enhancement partially mitigating the colour suppression, as this enhancement arises only in the eikonal limit.

Fiducial results
In this section, we provide results for both factorisable and non-factorisable corrections on differential and fiducial cross sections.Although they have been presented separately in Refs.[7,12], we show here for the first time the combined factorisable and non-factorisable NNLO prediction to VBF Higgs production.

With VBF cuts
In figures 10 and 11 we compare the size of the factorisable and non-factorisable corrections to VBF Higgs production under the selection cuts of section 3.1.In the upper panels we show the NLO prediction.The lower panels show various predictions normalised to the NLO prediction.In blue we show the factorisable NNLO prediction with its associated scale uncertainty band.The red curve shows the combined NNLO factorisable and non-factorisable prediction.In the bulk of the phase space, the non-factorisable corrections are small and within the scale uncertainty bands.However, it is interesting to observe that for large p t,j2 and p t,H the corrections can in certain regions become larger than the factorisable scale uncertainty.This makes the nonfactorisable corrections of potential relevance in boosted Higgs boson searches.On the other hand it is clear from figure 9 that the eikonal approximation is not reliable at very high values of p t,H , and the corrections should therefore be applied with care.A summary of the impact of O(α 2 s ) corrections on the fiducial cross section is shown in figure 1a as a function of the ∆y jj and m jj selection cuts, requiring also two R = 0.4 anti-k t jets with p t > 25 GeV and |y j | < 4.5.The corrections have only a mild dependence on the cuts, decreasing from roughly −4% at low cuts to around −3% at larger cut values.The non-factorisable corrections shown in figure 1b on the other hand show a stronger dependence on the cuts.In general they are suppressed by an order of magnitude compared to the factorisable corrections.They increase in size with an increase in the ∆y jj cut, and decrease as the m jj cut increases.The first effect is related to the Glauber enhancement which grows with the separation of the jets.The decrease of the non-factorisable corrections as the m jj cut increases is consistent with figure 11, where we observe that these corrections change sign around 2.5 TeV.It is important to keep in mind that these results will strongly depend on the choice of jet radius.Beyond LO the VBF cross section is well-known to be affected by real radiation escaping the tagging jets [37], while the NNLO non-factorisable corrections are independent of the jet radius.Therefore, one should compare the non-factorisable contribution not only to their factorisable counter-part, but also to the size of the scale uncertainty bands, particularly for large R values when the NNLO factorisable corrections become numerically small but their scale uncertainty remains large.

Without selection cuts
As discussed in section 3.2 the eikonal approximation is only formally valid when considering fiducial VBF production.It can however still provide a useful estimate of the size of the non-factorisable corrections in inclusive production.This is of particular interest, as the inclusive factorisable N 3 LO corrections are available for comparison in this regime.
In figure 12 we show in the lower panels for p t,H and |y H | in blue the factorisable NNLO prediction, and in green the factorisable N 3 LO prediction, both normalised to the NLO curve shown in the upper panel.We also show the combined factorisable and non-factorisable NNLO prediction in two setups: the non-factorisable corrections computed according to eq. ( 29) everywhere in phase space  (red) and the non-factorisable corrections computed according to eq. ( 29) when ξ < 0.2 and set to 0 otherwise (dashed-orange).This last procedure is used to verify that differential observables do not receive significant contributions from the large ξ region.We observe that in the bulk of the phase space the numerical difference between the red and the dashed-orange curves is very small -of the order of a few percent.It is therefore clear that the bulk of the non-factorisable corrections in the eikonal approximation come from the ξ < 0.2 region, even without VBF cuts.This is consistent with figure 9 which shows that the mean value of ξ is typically below 0.2.We observe that the non-factorisable NNLO corrections are typically larger than the factorisable N 3 LO ones, and that they are not covered by the NNLO scale variation uncertainties.
In fact, the non-factorisable NNLO corrections are almost O(40%) of the factorisable ones at this order.However we stress again that the non-factorisable corrections computed in the eikonal approximation do not necessarily provide reliable predictions in the full VBF phase space, as subleading terms can become relevant.It should also be noted that this large effect stems not from an enhancement of the non-factorisable effects, but rather from an order of magnitude decrease in the factorisable corrections when no cuts are applied.

Results for di-Higgs VBF production
We will now investigate the impact of non-factorisable contributions to the VBF Higgs pair production process.The electroweak parameters are set identically to the previous single Higgs study detailed in 3.1, with a width Γ H = 4.030 MeV for the internal Higgs propagator, while the final state Higgs bosons remain in the narrow-width approximation.For the factorisable corrections the central renormalisation and factorisation scales are now set to and uncertainties from missing higher orders are again estimated by varying the scales symmetrically up and down by a factor two, as was discussed in section 3.1.For the non-factorisable corrections we pick the same central scales, but when showing the residual scale uncertainty envelope, we perform the full 7-point scale variation, i.e. varying independently µ R and µ F by a factor 2 but keeping the ratio to the interval 1 2 ≤ µR µF ≤ 2. The NNLO corrections are calculated with proVBFHH v1.2.0 [38,39].

Validity of the eikonal approximation
Similarly to what we did for single Higgs production, we start by examining the validity of the eikonal approxima- tion.We expect the eikonal approximation to be valid when all transverse scales are small compared to the total centreof-mass energy.To test this statement quantitatively, we define where t and u are defined as below eq.( 15).In figure 13, we show in the left panel the normalised di-Higgs VBF cross section integrated in ξ HH , both fully inclusively and under VBF cuts.Here we see that compared to the single Higgs process, the ξ HH distribution with no cuts is contained to lower values below ξ HH ∼ 0.25.With VBF cuts the two distributions are very similar, and we therefore expect the eikonal approximation to be valid also in the di-Higgs process.In particular, from the right panel of figure 13, it is clear that the approximation only starts to break for very large transverse momentum values of the Higgs pair.

Fiducial results
In this section we discuss the impact of the non-factorisable NNLO corrections to di-Higgs VBF production computed in section 2.2.As was discussed there, the non-factorisable corrections are characterised by an interesting interference  Table 2: Fiducial cross sections for single and double Higgs VBF production at NLO, along with the corresponding NNLO corrections.The quoted uncertainties correspond to scale dependence, while statistical errors at NNLO are about 0.5 .For details on the scale variation see sec 3.1.
their interference, σ T B , c.f. figure 5.As one can see, the di-Higgs cross section at LO is the result of cancellations spanning several orders of magnitude.For the individual sets of diagrams, the non-factorisable corrections are below 1% and one can show that the combined 1-and 2-loop contribution in each case is always negative.It is interesting to note that the relative correction to σ T T of −0.327% is very close to the correction found in the single Higgs process of −0.32% under identical cuts (c.f.figure 1b).
To put the size of the NNLO non-factorisable corrections into context, we compare them to the factorisable NNLO corrections.In figure 14 we show the corrections to the two hardest jets, normalised to the NLO cross section.For low to moderate jet transverse momenta the non-factorisable corrections are at the few permille level and typically smaller in size than their factorisable counterparts.For p t,j2 they can grow to 1−2% when the transverse momentum becomes large, see figure 14.We also note here that a similar growth of the NNLO non-factorisable corrections in the jet p t distributions can be observed also in single Higgs production, see figure 10.We should also stress here that, as the jet transverse momenta grow, the eikonal approximation becomes less reliable and hence one should use the results in this region with caution.
In figure 15 we show the transverse momentum distribution of the Higgs with larger transverse momentum and of the di-Higgs system.The corrections remain moderate, and tend to be substantially smaller than the factorisable corrections.
Finally, in figure 16 we show the dijet rapidity separation and invariant mass.The non-factorisable corrections tend again at the few permille level in both observables, and smaller than the factorisable corrections.

Conclusions
In this paper we have studied the relative sizes of factorisable and non-factorisable QCD corrections to both single and double VBF Higgs production.A summary of the results is given in table 2, which shows the NLO fiducial cross section of single and di-Higgs production and the corresponding NNLO corrections.This study was made possible by recent advances in estimating the non-factorisable terms contributing to the NNLO cross section [12], which we extended to the di-Higgs process.We have presented the combined factorisable and non-factorisable NNLO corrections, as implemented in the public code proVBFH v1.2.0 for single Higgs VBF production.We find that for typical selection cuts the non-factorisable NNLO corrections are small and mostly contained within the factorisable scale uncertainty bands.For large jet and Higgs transverse momenta, the non-factorisable corrections can become comparable to the factorisable ones.In this region, it however not clear that the eikonal approximation used to estimate the corrections remains valid.
We also showed that the corrections computed in Ref. [12] can be used to provide an estimate for the nonfactorisable corrections for the fully inclusive VBF phase space.In this case we find that the non-factorisable NNLO corrections are of the same order as the NNLO factorisable corrections, and moderately larger than the factorisable N 3 LO corrections.This is in contrast with the usual statement that non-factorisable corrections can be neglected at this order [5] for inclusive quantities.We stress that this estimate comes from an extrapolation of the eikonal approximation into a regime beyond where it is expected to remain valid, and should therefore only be taken as an estimation of the true size of non-factorisable NNLO corrections to fully inclusive VBF.
Finally, we have implemented the non-factorisable correction to the Higgs pair production process in VBF, which is available in proVBFHH v1.2.0.We find that the nonfactorisable corrections to double-Higgs production follow similar patterns as were found in the single-Higgs case.There is however a suppression of the non-factorisable corrections which comes from a delicate cancellation of the various Born diagrams.
Public versions of the codes used in this article are available online [40].These results pave the way for precision measurements of the Higgs sector at the LHC and HL-LHC, as well as for further studies of non-factorisable effects and their interplay with the choice of jet radius.

- 3 .Fig. 1 :Fig. 2 :
Fig. 1: Single Higgs VBF production: Ratio of the factorisable (a) and non-factorisable (b) NNLO corrections relative to LO for fiducial cross sections with two R = 0.4 anti-k t jets satisfying p t > 25 GeV and |y j | < 4.5, as a function of the m jj and ∆y jj selection cut.

Fig. 3 :
Fig. 3: Born diagram for the production of n Higgs bosons in VBF (a) and representative 2-loop factorisable corrections (b).

Fig. 7 :
Fig. 7: Non-factorisable 1-loop (a) and 2-loop (b) corrections in the eikonal approximation.Notice that these are two-dimensional euclidean diagrams in the plane transverse to the incoming quark momenta.

Fig. 9 :
Fig. 9: (left): The normalised integrated cross section as a function of ξ fully inclusively (purple) and under the VBF cuts of sec.3.1 (green) for single Higgs production through VBF.(right): The average of ξ as a function of y H and p t,H .

Fig. 10 :
Fig. 10: Upper panel: NLO prediction for VBF production with cuts for the transverse momentum of the two leading jets and the Higgs boson.Lower panel: Ratio of the factorisable NNLO prediction to NLO (blue) and of the full NNLO prediction to NLO (red).The blue bands represent the scale uncertainty of the NNLO factorisable prediction.

Fig. 11 :
Fig. 11: Upper panel: NLO prediction for VBF production with cuts for the invariant mass and rapidity separation of the tagging jets, as well as for the rapidity of the Higgs boson.Lower panel: Ratio of the factorisable NNLO prediction to NLO (blue) and of the full NNLO prediction to NLO (red).The blue bands represent the scale uncertainty of the NNLO factorisable prediction.

Fig. 12 :Fig. 13 :
Fig. 12: panel: NLO prediction for inclusive VBF production for the transverse momentum and rapidity of the Higgs boson.Lower panel: Ratio of the factorisable NNLO prediction to NLO (blue) and of the full NNLO prediction to NLO (red).The blue bands represent the scale uncertainty of the NNLO factorisable prediction.The orange dashed curve shows the full NNLOprediction after applying a cut ξ < 0.2 on the non-factorisable component.The factorisable N 3 LO prediction is shown in green.

Fig. 15 :
Fig. 15: Same as figure 14 for (a) transverse momentum of the hardest Higgs boson (b) transverse momentum of the di-Higgs system.

Table 1 :
Fiducial cross section for vector-boson fusion di-Higgs production at 13 TeV under the cuts of sec.3.1.The first row indicates the Born contribution in femtobarn of the triangle diagrams, box diagrams and their interference.The second row shows the 1-loop squared non-factorisable (1-loop NF) correction in percent of the Born results for the same three contributions.The third row shows the same but for the 2-loop times tree-level non-factorisable (2-loop NF) contribution.The last row shows the same breakdown but for the sum of both contributions.The last column shows the sum of the contributions across each row.pattern which is not present in the single Higgs process.In table 1 we exemplify this by showing the LO fiducial cross section under the cuts of section.3.1 and their NNLO non-factorisable corrections.We split the cross section into the contribution coming from only the T 1 and T 2 topologies, σ T T , only the B 1 and B 2 topologies, σ BB and