Scale dependence of non-factorizable virtual corrections to Higgs boson production in weak boson fusion

The renormalization-scale dependence of the non-factorizable virtual corrections to Higgs boson production in weak boson fusion at next-to-next-to-leading order in perturbative QCD is unusually strong, due to the peculiar nature of these corrections. To address this problem, we compute the three-loop non-factorizable contribution to this process which accounts for the running of the strong coupling constant, and show that it stabilizes the theoretical prediction.


Introduction
The non-factorizable corrections to Higgs boson production in weak boson fusion (WBF) are color-suppressed and, for this reason, are expected to be smaller than the factorizable ones. 1 However, virtual non-factorizable corrections, which start contributing to the WBF cross section at next-to-next-to-leading order (NNLO) in perturbative QCD, exhibit a peculiar enhancement by two powers of π.This enhancement was first observed when the two-loop amplitude was computed in the leading eikonal approximation in Ref. [11].Recently, the calculation reported in Ref. [11] was extended in two important ways.First, in Ref. [12] the calculation of real-virtual and double-real non-factorizable corrections to Higgs boson production in weak boson fusion was performed.It was shown that these contributions are negligible in comparison with O(π2 )-enhanced eikonal virtual corrections.Second, in Ref. [13] it was found that the leading power correction to the eikonal contribution reduces the prediction based on the eikonal approximation by about 20 percent.
An important observation that emerged from the studies described above is that the dependence of these predictions on the renormalization scale is significant and can easily reach 20 − 30 percent. 2 This feature is the direct consequence of the fact that the nonfactorizable corrections appear for the very first time at NNLO in the perturbative QCD expansion.For this reason the mechanism responsible for compensating the dependences of physical quantities on the renormalization scale is not yet present in the results of Refs.[11][12][13], in spite of the fact that they are part of the NNLO QCD corrections to Higgs boson production in WBF.
To address this problem, one should calculate the N3LO QCD non-factorizable corrections.This calculation is currently unfeasible.A much simpler possibility is to compute the corrections associated with the running of the QCD coupling constant, starting from the fermion bubble insertions into the gluon propagators, and then re-writing such contributions through the full QCD β-function in the spirit of the BLM approach [14].Since, as follows from Ref. [11], the leading contribution to the non-factorizable corrections to WBF is related to the potential scattering, it appears quite natural to include an improved description of the interaction potential between two quarks in QCD as the first step towards a more reliable theoretical prediction.This is what we do in this paper.
The rest of the paper is organized as follows.In Section 2 we explain how to accommodate the running coupling constant into the computation of non-factorizable corrections and derive a convenient one-dimensional representation for the three-loop O(β 0 α3 s ) corrections.In Section 3 we discuss the analytic computation of these contributions.In Section 4 we show that the computed corrections strongly reduce the renormalization scale dependence of the theoretical predictions for the non-factorizable contributions to Higgs boson production in WBF, making them more precise and reliable.We conclude in Section 5. Analytic formulas for O(β 0 α 3 s ) corrections to Higgs boson production in WBF can be found in an ancillary file provided with this submission.

Fermion-bubble corrections
According to the analyses of Refs.[11,13] the leading O(α 2 s ) QCD non-factorizable contributions to Higgs boson production in weak boson fusion can be written in the following way The function C nf , which describes the effect of the non-factorizable corrections, reads where . (2.3) In the above equation, we have used as well as Furthermore, m V is the mass of the electroweak gauge boson 3 and we employ boldface fonts to denote two-dimensional Euclidean vectors.For example, p 3,4 are the transverse momenta of the tagging jets in the WBF process pp → 2j + H.
An important feature of non-factorizable corrections is that the function C nf is infrared finite although the functions C 1,2 are infrared divergent.To see this, it is simplest to use the integral representations for C 1,2 and write C nf as . (2.6) Using explicit expressions for the ∆-functions, it is easy to show that where the rank-2 tensor T ij is non-singular in the limit of vanishing k 1,2 .This fact ensures that C nf is not infrared divergent and, as a consequence, in Eq. (2.6) we can replace the space-time dimensionality d with 4, so that d − 2 → 2.
It is now straightforward to introduce the running coupling constant into Eq.(2.6) since all we need to do is to modify the gluon propagators ∆ 1 and ∆ 2 in Eq. (2.6).We find4 where ∆i = ∆ i 1 + µ 2 e 5/3 , (2.9) is the number of massless quark flavors and µ is the renormalization scale of the coupling α s in Eq. (2.1).We note that by using Eq.(2.9) in Eq. (2.8) we describe the all-order impact of the running coupling constant on C nf but, since we are interested in computing O(β 0 α 3 s ) corrections to Higgs boson production in WBF, we only use this formula for bookkeeping purposes.As we explain below, eventually, we expand Eq. (2.8) in powers of α s to the required order.
Although one can work with C nf in Eq. (2.8), it is more convenient to compute C 1 and C 2 separately.However, C 1,2 are infrared divergent so that extracting them from Eq. (2.8) and calculating them separately requires an infrared regulator.Since this regulator should work efficiently for integrands with 1/k 2 and 1/k 2 ln k 2 , it seems that the best option is an analytic regulator where each gluon propagator is raised to the power 1 + ν.Hence, we define and we will use these two auxiliary functions to compute C nf through O(β 0 α s ).We begin with the calculation of C 1 (ν).Introducing Feynman parameters and integrating over the transverse momentum, we obtain the following representation for this function . (2.11) To write Eq. (2.11) we introduced the following functions and the following dimensionless quantities that we will use throughout the calculation.
With this result at hand, it is straightforward to compute C 2 (ν 1 , ν 2 ).Indeed, in this case we can first integrate over one of the two loop momenta, keeping their sum fixed.We obtain where ν 12 = ν 1 + ν 2 .Using this expression in Eq. (2.10), we find that the remaining integration over k 12 is identical to the one-loop case provided that we replace ν with ν 12 .Hence, we find It is now easy to see that to compute C nf through O(β 0 α s ) we need to take ν 1 = ν 2 = ν, expand the quantity through first order in ν and identify ν with α s β 0 /(2π).Writing the expansion of C 1 (ν) in powers of ν as follows5 we obtain where We can easily derive convenient one-dimensional integral representations for the coefficients 1 , i = 0, 1, 2, by expanding Eq. (2.11) in powers of ν [15,16].We find (2.20) We will discuss below how to use these representations to complete the analytic computation of C (0),( 1),( 2) but we emphasize that the integral representations in Eq. (2.20) are quite convenient for numerical calculations.
Before we discuss the analytic computation in full generality, we study C nf in some kinematic limits where compact formulas can be derived.The simplest case to consider is when all transverse momenta are small |p 3 | ∼ |p 4 | ≪ m V .Then, r 1 → 0 and r 12 → 1.As the result, the expression for C 1 (ν) simplifies and we obtain lim We then find Following Ref. [14], we define the appropriate renormalization scale µ * as the scale for which corrections proportional to β 0 vanish.From the above equation it follows that at small transverse momenta, the appropriate renormalization scale for the non-factorizable corrections in WBF is (2.24) As the result, the dependence on t in the integrand in Eq. (2.11) disappears and we obtain It is then straightforward to compute C nf also in this case by expanding the hypergeometric function in powers of ν.We do not provide the result of such a calculation here and do not discuss the corresponding "optimal" scale choice because C nf in this case is not positive definite, which leads to pathological results for µ * .However, we that we can investigate the leading term in the expansion of C nf in the limit |p 3 | ≫ m V .In this case, the leading asymptotic behavior is described by the following formula It follows that the appropriate renormalization scale in this case is In practice, the transverse momenta of the outgoing jets relevant for Higgs production in weak boson fusion are below 200 GeV.For such transverse momenta, the formula for the renormalization scale in Eq. (2.27) leads to µ * ≤ m V ∼ m H . Therefore, since for the smallest transverse momenta µ * ∼ 0.4m V and for the highest relevant transverse momentum µ * ∼ m H , it appears that the traditional choice of the scale variation interval used in fixed-order computations of Higgs production in WBF, roughly covers the two extreme choices of the renormalization scale discussed above.In Section 4 we will illustrate that this is indeed the case by an explicit numerical computation for the fiducial cross section and many kinematic distributions.

Analytic computation of the function C 1
In this section, we discuss the analytic computation of the function C 1 (ν).In principle, such a computation requires a straightforward application of partial fractioning and integration by parts.However, one should perform these operations in the right order to minimize the size of intermediate expressions to the extent possible.Furthermore, it is important to use symmetries to write the result as compactly as possible.To this end, we note that the expression for C 1 (ν) in Eq. (2.11) is invariant under the interchange of x and y, and we will use this symmetry when writing the result of the integration.We find it convenient to integrate by parts in Eq. (2.11) first.We obtain where It is straightforward to compute I 1 .We find where the symmetry under the interchange of x and y is made manifest.Expanding the above formula in powers of ν, we obtain where x + 1 x − y y + (y − 1) ln(y + 1) + 2 ln 2 (y + 1) + Li 2 (−y) + π 2 6 + (x ↔ y), I x + 1 x − y − y + (1 − y) ln(y + 1) + (1 − y) ln 2 (y + 1) + ln(y) ln 2 (y + 1) To compute I 2 , we expand the integral in Eq. (3.2) in powers of ν and obtain the following expression with We note that two quadratic polynomials r 1,2 appear in the denominators in the above integrands.To continue with the t-integration, it is important to factorize them.We find where u, v are given by the following expression In the physical region (z − x − y) 2 − 4xy < 0 so that u and v are complex conjugates of each other.It is easy to see that under the transformation For r 2 , we find where It follows that the two roots t ± are outside of the integration region t ∈ [0, 1].Repeatedly applying partial fractioning and integration by parts, we recognize that all integrals that appear in Eq. (3.2) are of the following form where t a,..,d are elements of the following set Such integrals can be easily written as linear combinations of polylogarithmic functions through weight 3.
We emphasize one more time that, although the integration over t is straightforward, writing the result in an optimal way requires some effort.To this end, we will use the variables x, y, z where possible and only employ u and v where necessary.We will also make the x ↔ y symmetry manifest.Finally, it turns out to be convenient to introduce a new variable w defined as and to use it when writing the result.We find which is explicitly symmetric under x ↔ y.For C 1 , we find where + (x + 1)(y + 1)(w(−x + y + z + 2) − 2) (w + 1)(x − y)((x + 1) − w(y + 1)) ln(w) ln(x + 1) It is easy to check that A 1 is real and that Hence, the final result for C The function C (2) 1 admits a representation similar to C (2) The functions A 2,1 , A 2,2 , B 2,1 , B 2,2 are, unfortunately, very lengthy and we do not present them here.However, they can be found in the ancillary file provided with this submission.

Numerical results
We now discuss the numerical results of the calculation.We consider proton collisions at the LHC running at an energy of 13.6 TeV and adopt parameters and selection criteria from Ref. [12].The Higgs boson is considered stable with a mass of m H = 125 GeV.The vector boson masses are set to m W = 80.398 GeV and m Z = 91.1876GeV with widths of Γ W = 2.105 GeV and Γ Z = 2.4952 GeV respectively.The Fermi constant value of G F = 1.16639 × 10 −5 GeV −2 is used to derive the weak couplings and the CKM matrix is set equal to the identity matrix.
We use NNPDF31-nnlo-as-118 parton distribution functions [17] and α s (m Z ) = 0.118 for all calculations reported below.The evolution of both parton distribution functions and the strong coupling constant is obtained directly from LHAPDF [18].We fix the factorization scale µ F = m H throughout this calculation.
We require that a WBF event contains at least two jets with transverse momenta higher than 25 GeV.The tagging jets are required to have rapidities between −4.5 < y < 4.5 and should be separated by a large rapidity interval |y j 1 − y j 2 | > 4.5.Their invariant mass should be larger than 600 GeV.The two jets should be in opposite hemispheres in the laboratory frame; this is enforced by requiring the product of their rapidities to be negative, y j 1 y j 2 < 0. In principle, jet identification requires a particular jet algorithm but this is not relevant with only two quarks in the final state.For the results reported below, we employ the two-loop amplitude in the leading eikonal approximation as summarized in Eq. (2.8), and we do not include the next-to-eikonal power correction computed in Ref. [13].
The main conclusion of our analysis is that the O(β 0 α 3 s ) corrections significantly reduce the scale dependence of the non-factorizable contributions.To illustrate this point, we first show results for fiducial non-factorizable contributions to the WBF cross section at leading O(α 2 s ) order, and then compare them with the results at next-to-leading order where we only include the O(β 0 α 3 s ) correction.We find where we have used the renormalization scale µ = m H to obtain the central values, and µ = m H /2 and µ = 2m H to obtain values described by superscripts and subscripts in Eq. (4.1), respectively.It follows from Eq. (4.1) that the scale variation is reduced very significantly once O(β 0 α 3 s ) contributions are included.The same statement applies to kinematic distributions.We illustrate this in Figures 1  and 2 where examples of transverse momenta, rapidity and two-jet invariant mass distributions are shown.In all plots, the upper pane displays the leading order (tree-level) distribution.In the lower pane ratios of non-factorizable corrections to leading order distributions are shown.Bands correspond to the range of theoretical predictions obtained with the renormalization scales from the interval m H /2 ≤ µ ≤ 2m H at leading O(α 2 s ) and nextto-leading O(β 0 α 3 s ) orders, respectively.We observe that accounting for NLO O(β 0 α 3 s ) corrections stabilizes theoretical predictions by strongly reducing their dependence on the renormalization scale.We note that the O(β 0 α 3 s ) results are sometimes outside the O(α 2 s ) scale-variation bands; this mostly happens at large(r) values of the transverse momenta and invariant masses.

Conclusions
In this paper we computed O(β 0 α 3 s ) corrections to the non-factorizable contribution to Higgs boson production in weak boson fusion.These corrections reflect the impact of the running of the QCD coupling constant on the non-factorizable contribution to the WBF cross section and are responsible for stabilizing the dependence of the theoretical prediction on the renormalization scale.Indeed, we find that after including these O(β 0 α 3 s ) corrections, the dependence of the cross section on the renormalization scale reduces from about O(20) percent to below O(5) percent.Similar reductions of the scale dependence are observed in theoretical predictions for major kinematic distributions including transverse momenta and rapidity distributions of the tagging jets and the Higgs boson.
We provided a simple one-dimensional integral representation of the O(β 0 α 3 s ) nonfactorizable corrections as well as the analytic formulas for these corrections.Although the analytic results are complex, they can easily be implemented into partonic Monte Carlo and used to obtain phenomenological predictions.In fact, we have used the onedimensional integral representation of these corrections to cross check the results of the analytic computation.Under realistic running conditions, analytic formulas provide a significant speed-up whereas a one-dimensional integral representation is a slow but robust way to compute cross sections and distributions.

3 sFigure 1 :
Figure 1: Transverse momenta distributions of hardest and next-to-hardest jets in Higgs boson production in weak boson fusion.The upper panes show the LO (tree-level) distributions and the lower panes show the ratio of non-factorizable contributions to LO for corrections of O(α 2 s ) (blue) and O(β 0 α 3 s ) (red).The factorization scale, µ F , is kept fixed and only the renormalization scale, µ, is varied.See text for further details.

3 sFigure 2 :
Figure 2: Non-factorizable corrections to Higgs boson rapidity distribution (left panes) and the invariant mass distribution of the two final state jets (right panes).The upper panes show the LO (tree-level) distributions and the lower panes show the ratio of nonfactorizable contributions to LO for corrections of O(α 2 s ) (blue) and O(β 0 α 3 s ) (red).The factorization scale, µ F , is kept fixed and only the renormalization scale, µ, is varied.See text for further details.