Non-factorizable virtual corrections to Higgs boson production in weak boson fusion beyond the eikonal approximation

Non-factorizable virtual corrections to Higgs boson production in weak boson fusion at next-to-next-to-leading order in QCD were estimated in the eikonal approximation [1]. This approximation corresponds to the expansion of relevant amplitudes around the forward limit. In this paper we compute the leading power correction to the eikonal limit and show that it is proportional to first power of the Higgs boson transverse momentum or the Higgs boson mass over partonic center-of-mass energy. Moreover, this correction can be significantly enhanced by the rapidity of the Higgs boson. For realistic weak boson fusion cuts, the next-to-eikonal correction reduces the estimate of non-factorizable contributions to fiducial cross section by O(30) percent.


Introduction
At the Large Hadron Collider (LHC), Higgs bosons are frequently produced in weak boson fusion (WBF).This process has a recognizable signature, characterized by two energetic low-p ⊥ jets in the opposite hemispheres.Higgs boson production in WBF has been measured by CMS [2,3] and ATLAS [4,5] collaborations.The measured WBF cross section agrees with the Standard Model prediction to within 20%.Further improvements in the experimental exploration of Higgs boson production in weak boson fusion are expected during the Run III and the high-luminosity phase of the LHC.
Theoretical understanding of Higgs boson production in WBF is very advanced.It is based on the knowledge of next-to-leading order (NLO) [6,7] and next-to-next-to-leading order (NNLO) [8][9][10][11] QCD corections, as well as mixed QCD-EW [12] corrections to this process.N 3 LO QCD corrections have also been computed [13]; they change the leadingorder cross section by just about one permille.
It is to be noted, however, that all these studies were performed in the so-called factorization approximation where contributions due to gluon exchanges between two incoming fermion lines are neglected.These effects, that we will refer to as non-factorizable corrections, are color-suppressed and, for this reason, are expected to be smaller than the factorizable ones [8,9].However, virtual non-factorizable corrections, which start contributing to the WBF cross section at NNLO QCD, exhibit a peculiar enhancement by two powers of π.This enhancement was first observed when the two-loop non-factorizable amplitude was computed in the leading eikonal approximation [1].To better understand these two-loop effects and to establish the validity of the eikonal approximation for phenomenological analyses of Higgs boson production in weak boson fusion, it is essential to go beyond the leading term in the eikonal expansion.
Since the calculation of exact non-factorizable contributions, which requires the twoloop five-point amplitude with five independent kinematic variables and two masses, is currently not possible, it is reasonable to explore the possibility to extend the eikonal expansion beyond the forward limit.In this paper we make the first step in that direction and compute the leading power correction to the eikonal limit of non-factorizable five-point WBF amplitude.
The remainder of this paper is organized as follows.In the next section, we describe kinematics of weak boson fusion and explain how we use it to set up an expansion around the eikonal limit.In Sections 3 and 4, we derive integral representations for one-and two-loop amplitudes which contribute to non-factorizable corrections to WBF; these representations retain the next-to-eikonal accuracy.In Section 5, we explain how the infra-red finite, twoloop non-factorizable correction can be derived from these integral representations.In Section 6 we analyze the numerical impact of the computed next-to-eikonal corrections and show that they change the current estimate of the non-factorizable contribution to the WBF cross section by about O(30) percent.We conclude in Section 7. Discussion of the analytic computation of one-and two-loop non-factorizable amplitudes is relegated to appendix.The analytic results for the amplitudes can be found in an ancillary file provided with this submission.

Kinematics of Higgs production in weak boson fusion
We begin with the discussion of the kinematics of Higgs production in the WBF process q(p 1 ) + q(p 2 ) → q(p 3 ) + q(p 4 ) + H(p H ) . (2.1) We perform the Sudakov decomposition of the four-momenta of the outgoing quarks and write Employing the on-shell conditions p 2 3 = 0, p 2 4 = 0, we find1 where s = 2p 1 • p 2 is the partonic center-of-mass energy squared.The WBF events are selected by requiring that two tagging jets with a relatively small transverse momentum are present in opposite hemispheres; this ensures that α 3 ∼ β 4 ∼ 1 and that p 2 3,⊥ ∼ p 2 4,⊥ s.
We define two auxiliary vectors q 1 and q 2 which describe momentum transfers from the quark lines to the Higgs boson.They read where δ 3 = 1 − α 3 and δ 4 = 1 − β 4 .It follows from the momentum conservation condition that p H = q 1 + q 2 . (2.5) Upon squaring the two sides of this equation and some rearrangements, we find We can use Eq.(2.6) to fully specify the relevant aspects of WBF kinematics around the forward limit.Indeed, given the proximity of the Higgs boson mass and electroweak boson masses, and the fact that the important contribution to WBF cross section comes from kinematical configurations where the transverse momenta of tagging jets are comparable to m H and m W,Z , the above equation implies Note that we introduced a parameter λ to indicate the smallness of various ratios in the above equation.We consider central production of Higgs bosons so that neither forward nor backward direction is preferred.Then δ 3 ∼ δ 4 and We note that, with the required accuracy, the two parameters δ 3,4 can be written as follows where p H,⊥ is the transverse momentum and y H is the rapidity of the Higgs boson in the partonic center-of-mass frame.We will use the above relations between kinematic parameters to construct the expansion of one-and two-loop non-factorizable WBF amplitudes in the following sections.

One-loop non-factorizable contributions to WBF
We consider the one-loop non-factorizable QCD corrections to Higgs boson production in WBF.To avoid confusion, we note that they do not contribute to the WBF cross section at NLO since their interference with the leading order amplitude vanishes because of color conservation.Nevertheless, since the one-loop amplitude is needed for the construction of the NNLO QCD corrections, we need to discuss it.
To write the non-factorizable amplitude in a convenient way, we assume that the coupling of the vector boson V to the Higgs boson is given by ig V V H g µν and that the coupling of the massive vector boson to quarks is vector-like, −ig W γ µ .Since we work with massless quarks, their helicities are conserved and we can reconstruct non-factorizable contributions for V = Z and V = W from the results that are reported below.
We write the one-loop non-factorizable amplitude as follows where T a ij denote the generators of the SU (3) color group and A 1 stands for the colorstripped one-loop amplitude2 In Eq. (3.2), we used the notation to define propagators of virtual bosons.In addition, following the conventions in Fig. 1, we introduced two quark currents where we assumed that the incoming fermions are left-handed.In writing Eq. (3.4) we employed the quantities ρ i (k), i = 1, 2, 3, 4 to describe quark propagators; they read We would like to construct an expansion of the amplitude in Eq. (3.2) in powers of λ.To understand how to do tatt, we introduce the Sudakov parametrization of the loop momentum k 1 and write The integration measure in Eq. (3.2) becomes The various propagators in Eq. (3.2) are linear polynomials in α 1 and β 1 .Hence, integration over either one of these two variables can be easily performed using the residue Region α 1 Table 1: Kinematic regions relevant for one-loop non-factorizable contributions.Symmetric regions are not shown.
The one-loop amplitude, shown on the left, can be constructed by contracting the currents for the upper and lower fermion lines.The current for the upper fermion line ) is shown on the right.
theorem.The resulting integrand is a product of (at most) quadratic polynomials in the other variable so that the structure of singularities can be easily analyzed.Performing this analysis and assuming that the transverse loop momentum can either be of the same order as the transverse momenta of the outgoing jets or of the same order as the center-of-mass energy, we come to the conclusion that the following loop-momenta regions, 3 shown in Table 1, need to be considered.The first region is the so-called Glauber region; the second one is "Glauber-soft", the third one is soft, the fourth is collinear and the last one is hard.Using the scaling of the loop-momentum components as indicated in Table 1, we estimate the contributions of the various regions to the one-loop amplitude.We find We note that the leading order WBF amplitude scales as λ −2 and that, as follows from Eq. (3.8), the expansion of the one-loop amplitude proceeds in powers of √ λ.To compute O( √ λ) correction to the virtual amplitude, we need to account for the contributions of regions a), b) and c) to first subleading power and the contribution of region d) to leading power in the expansion in λ.
We begin with the discussion of region a).Using momentum scaling in Table 1, we simplify the various propagators that appear in the integrand in Eq. (3.4).To present the result in a compact way, we introduce the following quantities (3.9) In region a), all inverse propagators scale as O(λ).To compute the first subleading correction we need to keep all terms that scale as λ 3/2 and neglect all terms that scale as λ 2 .We find If we use the simplified propagators shown in Eq. (3.10) to compute the amplitude A 1 , we observe that integrations over α 1 and β 1 factorize.We then write where In Eq. (3.11) σ is a cut-off parameter that forces β 1 and α 1 to stay in the region since this choice will allow us to use the same cut-off σ to study the Glauber-soft region.
We note that we replaced k1 with k1,⊥ in the currents when writing Eq. (3.11); this is justified since α 1 and β 1 terms in the Sudakov expansion of k provide O(λ) and not O( √ λ) corrections in region a).Hence, if we aim at computing the non-factorizable amplitude with O( √ λ) relative accuracy, we can discard them.In fact, to compute the amplitude with O( √ λ) relative accuracy, terms with k1,⊥ in Eq. (3.11) can be dropped altogether.Indeed, since k 1,⊥ ∼ √ λ, if we retain it in one of the terms that appear either in Φ µν or in Φµν , the other current should be computed at leading λ-power.However, in this case and terms with k1,⊥ lead to the vanishing contributions pi k1,⊥ pi = 0, i = 1, 2 , (3.16) since p 2 1,2 = 0 and p 1,2 • k 1,⊥ = 0. Furthermore, in region a) we can expand the remnants of weak boson propagators that appear in Eq. (3.11).Keeping terms that provide O( √ λ) corrections, we find (3.17) Focusing on Φ µν , we simplify the expression for the current, use Eq.(3.17) and obtain4 where To compute Φ, we use valid for z a ∈ [−σ, σ].Furthermore, we need Neglecting the σ-dependent terms that will cancel with the contribution from the Glaubersoft region, we obtain A similar computation for Φµν gives where Combining these results for Φ and Φ and neglecting all terms beyond desired O( √ λ) corrections, we obtain the following contribution to the one-loop amplitude from the Glauber region We then proceed with the discussion of the contribution of region b) with the mixed scaling α 1 ∼ λ and β 1 ∼ √ λ.According to Eq. (3.8), we require the contribution of this region through first subleading terms.However, it is easy to see that, in actuality, the contribution of region b) starts at O(λ −3/2 ) and, therefore, should be computed at leading power only.
To understand why this is the case, we first discuss the currents J µν and Jµν and, in particular, the numerators of the contributing terms.Since we work with O( √ λ) accuracy, in region b) we should replace k 1 with k 1 → β 1 p 2 + k 1,⊥ in both currents.Suppose we do this replacement in J µν .Since these terms already provide an O( √ λ) correction, the current Jµν should be taken at leading power.Since at leading power Jµν ∼ p µ 2 p ν 2 , it is easy to see that all contributions of vector k 1 drop from the current J µν once the Lorentz indices are contracted.
However, if we account for k 1 in the current Jµν , the situation is different.In this case, since i) k 1 is independent of α 1 , ii) it appears with different signs in the two terms in Jµν , and iii) Jµν is contracted with J µν computed at leading power, the corresponding contribution vanishes after integration over α 1 .
Having concluded that, similar to the Glauber region, we can drop k 1 from the fermion currents, we note that the current J µν (k 1 , −q 1 − k 1 ) in region b) can be further simplified.Indeed, using the fact that β 1 ∆ 1 /s, Θ 3,1 /s, we expand the current and obtain This equation implies that in region b) the current scales as O(1) and not as O(λ −1/2 ) as a naive estimate suggests.This suppression occurs because of the cancellation between two terms in brackets in Eq. (3.26).This means that the contribution of the region b) starts at λ −3/2 , so that all ingredients needed to compute the amplitude in region b), except the current J µν (k 1 , −q 1 − k 1 ), are to be taken at leading power in λ.Hence, we find where Φ is still given by Eq. (3.24) and Calculation of this integral is straightforward.We obtain5 Performing a similar computation for a symmetric region β ∼ λ, α ∼ √ λ, we obtain Combining the contributions of regions a) and b), we find We turn our attention to region c) which corresponds to the soft scaling According to Eq. (3.8) we require the contribution of this region through first subleading power.However, a more careful analysis shows that the contribution of this region is suppressed stronger than originally expected.To see this we note that in the soft region, to leading power, the currents vanish.For example, the expression for and we have set it to zero because poles of the fermion propagators have already been accounted for when the Glauber region was analyzed.Hence, to obtain a non-vanishing contribution from the soft region, subleading terms in both currents J µν and Jµν are needed.
The subleading contributions to the currents scale as O(1) and not as 1/ √ λ as a naive estimate for the currents' scaling would suggest.This implies that at variance with the original estimate M (c) ∼ λ −2 in Eq. (3.8), the contribution of the soft region is suppressed by an additional power of λ.For this reason, the soft region is not needed for computing the two-loop non-factorizable amplitude with O( √ λ) accuracy.The contribution of the collinear region can be analyzed in the same way.Since, in this case, the amplitude scales as M d ∼ λ −3/2 , both currents need to be taken at leading power.We find It is clear that the contraction of the two currents in Eq. (3.33) vanishes.Hence, we conclude that collinear regions do not provide the O( √ λ) corrections to the leading term in the eikonal expansion.Since, obviously, the hard region is not relevant as well, we conclude that, with O( √ λ) accuracy, the one-loop non-factorizable contribution is given by the sum of the Glauber and Glauber-soft contributions in Eq. (3.31).
Having performed this analysis, we note that the final result for the two regions a) and b) can be obtained by simply computing the functions Φ and Φ from the following unexpanded expressions It is straightforward to integrate over β 1 and α 1 in Eq. (3.34).Indeed, focusing on the function Φ, we note that, if we close the integration contour in the upper half plane, only the residue at β 1 = Θ 3,1 /(sα 3 ) contributes.We then find

.35)
Expanding this result in δ 3 , performing a similar computation for Φ, and keeping only the relevant terms in the product of Φ and Φ, we obtain Eq. (3.31).
Finally, it is convenient to write the one-loop non-factorizable amplitude by extracting exact (i.e.not expanded in powers of λ) Born amplitude.The latter reads Using it, we write The function C 1 reads (3.38) We note that the above expression includes both the leading and the first subleading terms in the expansion of the one-loop amplitude in powers of √ λ.The function C 1 can be computed analytically and expressed through logarithmic and dilogarithmic functions; the corresponding discussion can be found in appendix.

Two-loop non-factorizable contributions to WBF
We continue with the computation of two-loop non-factorizable QCD corrections to Higgs boson production in weak boson fusion.The two-loop non-factorizable amplitude is written as where 2) The overall factor 1/2! comes from the symmetrization of two identical gluons and are bosonic propagators. 6Similarly to the one-loop case, in Eq. (4.2), we defined two quark currents; the conventions are explained in Fig. 2. The currents read and (4.5) The two-loop amplitude, shown on the left, is thought in terms of the currents that make it up.On the right, we define of the generalized upper current J µνα (k 1 , k 2 , −k 12 − q 1 ) used in the calculation of the two-loop amplitude.
To integrate over the loop momenta k 1,2 , for each of them (and also for their linear combinations) we need to consider regions shown in Table 1.We write The leading contribution comes from the Glauber region where Similar to the one-loop case, the leading correction arises from the mixed region where some of the α-or β-components scale as √ λ.Both in the Glauber region and in the mixed region, the loop momenta in the numerators of both currents J µνα and Jµνα can be discarded.The reason for this is the same as in the one-loop case and we do not repeat this analysis here.
Building on the experience with the one-loop calculation reported in the previous section, we can make the following observation.To obtain O( √ λ) correction, we only need to consider the cases where i) one or both β 1,2 components of the loop momenta scale as √ λ and both α 1,2 scale as λ, or ii) the other way around.If one of the two α's and one of the two β's scale as √ λ, then α 12 and β 12 also scale as √ λ.As the result, both currents in Eq. (4.4) and Eq.(4.5) are suppressed by O( √ λ).Thus, the contribution of this region is suppressed by O(λ) and can be discarded.We conclude that, if we want to construct an integrand which is valid both in the Glauber region and in the mixed region, we need to write an expression that incorporates √ λ corrections to one of the currents and that the other current should be taken at leading order.
These considerations also guide the expansion of the propagators in powers of λ to make them valid in both the Glauber region and in the mixed region.To write the approximate expressions, we define for i ∈ {1, 2, 12}, where α 12 = α 1 + α 2 , β 12 = β 1 + β 2 etc. and obtain We emphasize that the above expressions for propagators are valid both in the Glauber region and in the mixed region.Because of that, we can use them to compute the two-loop non-factorizable amplitude with O( √ λ) accuracy in the same way as Eq.(3.34) was used to do that in the one-loop case.
Using the expanded propagators, we simplify the currents in Eq. (4.5) and write the amplitude as where and To integrate over β 1,2 and α 1,2 it is useful to rearrange terms in the curly brackets in Eqs.(4.10, 4.11).Focusing on the integrand in Eq. (4.10), we rewrite it as follows We use the above representation to compute the function Φ in Eq. (4.10).We note that the first term in Eq. (4.12) can be discarded, because of the location of its poles.Indeed, To compute the contribution of the second term in Eq. ( 4.12), we close the integration contours in the lower half-planes for both integration variables.We obtain To compute the contribution of the third term in Eq. ( 4.12), we close the integration contours for both β 1 and β 2 in the upper half-planes.The result reads (4.15) Adding up Φ 1,2,3 , we find the following expression for the function Φ which provides the combined contribution of both the Glauber region and the mixed region The calculation for Φ proceeds in an identical way.We obtain (4.17) Finally, putting everything together and retaining terms that provide O( √ λ) corrections, we find the following result for the two-loop non-factorizable amplitude It remains to analyze the contributions of the other regions to the two-loop nonfactorizable amplitude.This analysis proceeds along the lines of the discussion of the one-loop case.It relies on the fact that for soft and collinear gluons, fermion currents simplify dramatically.Consider, for example, the case where k 1 is Glauber and k 2 is soft.Naively, this region would contribute at O(λ −2 ) so that we need to account for subleading contributions from this region.In practice, the contribution is O(λ) suppressed compared to a naive estimate.
Indeed, if k 2 is soft and k 1 is Glauber, then k 12 is also soft.To understand how the currents simplify in this case, consider Eq. (4.12).Since β 12 ∼ β 1 ∼ √ λ λ, the leading contribution in the last line of Eq. (4.12) vanishes; we then find that the current in Eq. (4.12) scales as λ −1 , at variance with the naive scaling λ −3/2 .We note that we ignore the pole at β 1,2 = 0 for the same reason as in the one-loop case, see Eq. (3.32).Since both currents exhibit this behavior, we conclude that the contribution of this region to the amplitude scales as O(λ −1 ) and not as O(λ −2 ) as naively expected.For this reason, it is not relevant for the calculation of the two-loop amplitude with the O( √ λ) accuracy.Similar to the one-loop case, we write the two-loop amplitude as where M 0 is defined in Eq. (3.36) and the function C 2 reads This function looks analogous to the one-loop function C 1 , c.f. Eq. (3.38).It is relatively straightforward to compute C 2 analytically; the corresponding discussion can be found in appendix.

Infrared pole cancellation and the finite remainder function
To compute the double-virtual non-factorizable contribution to the differential WBF cross section, we square the one-loop amplitude in Eq. (3.37) and calculate the interference of the two-loop amplitude in Eq. (4.19) with the Born amplitude.Summing over spins and colours, we find where α s = g 2 s /4π is the strong coupling constant, 7 dσ LO is the exact Born differential cross section for Higgs boson production in WBF and C nf characterizes the non-factorizable corrections.The function C nf reads and all terms that are suppressed stronger than O( √ λ) are supposed to be discarded when computing it.
We note that functions C 1 and C 2 are infra-red divergent; these divergences arise when the loop momenta k i,⊥ , i = 1, 2, vanish.Computing these functions and expanding in , we find Using these results in Eq. ( 5.2), we obtain which is infra-red finite and can be computed for ε = 0.The fact that the double-virtual contribution to non-factorizable corrections in WBF is finite through O( √ λ) is in accord with Catani's formula for infra-red divergences of generic two-loop amplitudes applied to the WBF process [17].Analytic results for the function C nf can be found in the ancillary file provided with this submission.

Numerical results and phenomenology
It is instructive to study the results of the calculation in several ways.First, we compare the analytic results for the function C nf at leading order in the λ-expansion against numerical results8 reported in Ref. [1] and find good agreement.Second, to explore the accuracy of our result in a realistic setting, we compare the one-loop amplitude including leading and first sub-leading terms in the λ-expansion, with the exact one-loop non-factorizable amplitude A 1 .To this end, we generate events that pass the WBF cuts [19], use them to evaluate both amplitudes, and compute the following quantity In Eq. (6.1),A 1 is the exact amplitude, A 1 is the leading eikonal amplitude and A a&b 1 is given in Eq. (3.31).We expect that in WBF kinematics X δ ∼ O( √ λ) and we would like to check if this is indeed the case.
WBF events are required to contain at least two jets with transverse momenta p ⊥,j > 25 GeV and rapidities |y j | < 4.5.The two jets must have well-separated rapidities, |y j 1 − Figure 4: Eikonal and next-to-eikonal contributions to the transverse momentum and rapidity distributions of the leading jet.In the upper pane, leading eikonal contribution is plotted with a red, dashed line and the next-to-eikonal one with a green, solid line.In the lower pane, we show the ratio of next-to-eikonal to eikonal contributions.We note that in the upper left pane, absolute values are shown.See text for further details.
use dynamical renormalization and factorization scales9 We set the mass of the W boson to m W = 80.398 GeV, the mass of the Z boson to m Z = 91.1876GeV, and the mass of the Higgs boson to m H = 125 GeV.The Fermi constant is taken to be G F = 1.16637 × 10 −5 GeV −2 .For 13 TeV proton-proton collisions, we find that the non-factorizable, double-virtual contribution to Higgs boson production in WBF evaluates to σ V V = (−3.1 + 0.53) fb , ( where we display contributions of leading and next-to-leading terms in the λ-expansion. We emphasise that the next-to-eikonal correction is calculated by excluding kinematic configurations where |y H | > 1 in the partonic center-of-mass frame, in addition to conventional WBF cuts that we listed earlier.It follows from Eq. (6.5) that the correction to the leading eikonal approximation amounts to O(17%).
We now turn to the discussion of kinematic distributions.In Fig. 4, we display nonfactorizable corrections to transverse momentum and rapidity distributions of the leading jet.The comparison of leading and next-to-leading eikonal contributions in lower panes shows that next-to-leading eikonal corrections range from ten to fifty percent.They appear to modify the leading order eikonal contribution by O(50%) for higher values of p ⊥,j 1 .This enhancement is partially related to the fact that the leading eikonal contribution changes sign at around p ⊥,j 1 ∼ 2m W , which is the reason for rapidly changing ratio of eikonal factors shown in the lower pane.
Figure 5: Eikonal and next-to-eikonal contributions to the transverse momentum and rapidity distributions of the Higgs boson.In the upper pane, leading eikonal contribution is plotted with a red, dashed line and the next-to-eikonal one with a green, solid line.In the lower pane, we show the ratio of next-to-eikonal to eikonal contributions.
The non-factorizable contributions to Higgs boson transverse momentum and rapidity distributions are shown in Fig. 5.The relation between eikonal and next-to-eikonal contributions are similar to what was observed for the fiducial cross section as well as p ⊥ and rapidity distributions of the leading jet.

Conclusion
We computed the two-loop virtual non-factorizable QCD corrections to Higgs boson production in weak boson fusion through next-to-leading order in the eikonal expansion.We found that such an expansion proceeds in powers of p ⊥,H / √ s ∼ m H / √ s and explained how to simplify the integrand of the two-loop amplitude to calculate both the leading and the next-to-leading terms in such an expansion.We observed that combining individual diagrams before integrating over loop momenta leads to significant simplifications in the calculation.This happens because contributions of some of the virtual-momenta regions, that are relevant for computing next-to-eikonal corrections in individual Feynman diagrams, receive additional suppression in the full amplitude and start contributing only at next-to-next-to-leading power.
We have derived compact integral representations for the double-virtual non-factorizable amplitude at both leading and next-to-leading power in the eikonal expansion.We have also explained how to compute the two-loop amplitude analytically and provided the analytic results in the ancillary file.
The numerical impact of next-to-eikonal corrections is significant although, given the overal smallness of non-factorizable contributions, they do not change the original conclusions of Refs.[1,21].Nevertheless, we find that, typically, the next-to-eikonal corrections change the estimate of the non-factorizable contributions based on the leading term in the eikonal expansion by O(20) percent.
As a final comment, we note that other sources of non-factorizable contributions to WBF cross sections, including double-real emission and the real-virtual corrections, were recently studied in Ref. [22].It was found that, thanks to the WBF cuts, all the contributions beyond the double-virtual ones are tiny and cannot impact the phenomenological studies of Higgs production in WBF in any way.The results reported in this reference allow us to estimate the contribution of the non-factorizable double-virtual corrections to the WBF cross section with a precision that is likely better than O(10) percent.Since the non-factorizable contribution itself is just O(1) percent of the total WBF cross section, the remaining uncertainties stemming from the imprecise knowledge of the two-loop virtual amplitude are irrelevant.We conclude that the current understanding of non-factorizable effects is sufficient for phenomenological studies of Higgs production in weak boson fusion envisaged for the Run III and the high-luminosity phase of the LHC.The master integrals are displayed in Fig. 6.Although we need these integrals at d = 2, we find it more convenient to study them first in four dimensions.In particular, at d = 4, we easily obtain the canonical basis [28] using the Magnus series expansion method [29].We then transform the integrals to d = 2 using the dimensional recurrence relations [30].In four dimensions, the canonical basis reads Note that all the g's are normalized to be dimensionless and can be regarded as functions of x, y and z only.The canonical basis vector g = (g 1 , g 2 , g 3 , g 4 , g 5 , g 6 ) T satisfies a differential equation in the dlog form, Eq. (A.6) can be recursively solved order-by-order in and the solutions are expressed in terms of Chen's iterated integrals [31] with some boundary constants that cannot be determined from the differential equations alone.For the integrals g 1,..,6 these constants can be computed with a relative ease since all canonical integrals, except g 2 and g 4 , vanish when x = y = z = 0.The non-vanishing integrals g 2 and g 4 at this kinematic point evaluate to As the result, the solutions of the system Eq.(A.6) can be expressed in terms of multiple polylogarithms.In fact, since we need g only through O( 2 ), relevant expressions for integrals involve logarithms and dilogarithms of u, v, w.To express them in terms of x, y, z, we use the following formulas u, v = x − y + z ± r 2 2z , w = 2 + z − r 1 2 .(A.12)

6 Figure 6 :
Figure6: Two-dimensional two-loop master integrals.The thick and thin internal lines represent massive and massless propagators, respectively.Red lines have mass m V .The black thick lines correspond to external "massive" legs.A dot on the internal line means raising the power of corresponding propagator by one.