Exact quark-mass dependence of the Higgs-gluon form factor at three loops in QCD

We determine the three-loop form factor parameterising the amplitude for the production of an off-shell Higgs boson in gluon fusion in QCD with a single massive quark. The result is obtained via a numerical solution of a system of differential equation for the occurring master integrals. The solution is also used to determine the high-energy and threshold expansions of the form factor. Our findings may be used for the evaluation of virtual corrections generated by top-quark and b-quark loops in Higgs boson hadroproduction cross sections at next-to-next-to-leading order.

: Complete set of Feynman diagrams with two fermion loops contributing to the Higgs-gluon form factor at three-loop order. The fermion loop connected to the Higgsboson line corresponds to a massive quark. The quark of the second fermion loop may be either massive or massless.

Finite remainders
Consider the amplitude for the fusion of two gluons of momenta p 1,2 , helicities λ 1,2 and adjoint-representation colors a 1,2 , followed by the production of one, possibly off-shell, Higgs boson: − iM g(p 1 , λ 1 , a 1 ) + g(p 2 , λ 2 , a 2 ) → H ≡ iδ a 1 a 2 ( 1 · p 2 ) ( 2 · p 1 ) − ( 1 · 2 ) (p 2 · p 1 ) 1 v α s π C . (2.1) Here, v is the Higgs-doublet Vacuum Expectation Value. The coupling of a single quark field, Q, of mass M = 0 to the Higgs-boson field, H, is given by the tree-level Lagrangian term −MQQH/v. Finally, the gluon polarisation vectors are normalised as follows: The Form Factor C is expanded in the strong coupling constant, α s , and the number of massless quark flavors, n l : The strong coupling is defined in the MS scheme with massive-quark decoupling. Its dependence on the renormalisation scale µ is given by the β-function for n l massless quarks, α s ≡ α (n l ) s (µ). Contributions C (n,n) = 0, n > 0 are only due to coupling constant renormalisation. The massive-quark mass, M , is defined in the on-shell scheme implying the same for the Yukawa coupling. The dimensionless form-factor expansion coefficients depend on two variables only: The leading contribution is: (2. 6) In the limit M → ∞: C (0) z = 0 = 1 3 . (2.7) Hence, the amplitude Eq. (2.1) may be obtained at M → ∞ from the Higgs-Effective-Theory tree-level Lagrangian: where G a µν is the standard QCD field-strength tensor, L QCD = −1/4 G a µν G a µν + L matter . Beyond leading order, the form factor is infrared divergent after renormalisation. The results presented in this publication correspond to Conventional Dimensional Regularisation with space-time dimension d = 4 − 2 . The infrared divergences may be factorised yielding the Finite Remainder, C I , of the form factor: where the two-loop I-operator of Catani [12] (see Ref. [13] for the specific case of the Higgs-gluon form factor) is given by: with the first two coefficients of the QCD β-function: and: (2.12) In general, C (n,n) I = 0, n > 0. However: Just as the form factor itself, the I-operator, Eq. (2.10), is independent of the scale µ (up to two-loop order of course). In consequence: 14) The dependence of the finite remainder on the scale logarithm, L µ , is thus given by the β-function only 1 : (2.15) 1 Notice that the I-operator of Ref. [11] (see Eq. (3.7b) of that publication) is missing a scale-dependent factor in the Hg-term (compare to Eq. (4.38) of Ref. [13]). With this difference, the I-operator of Ref. [11] is not scale invariant and C (2) I contains an additional contribution to the single scale-logarithm term, Hg/4 C (0) Lµ.
A different finite remainder, C Z , is obtained if the factorisation of infrared divergences is performed in the MS scheme [14]. Define: The solution at two-loops is: with the anomalous dimensions: (2.20) Since the dependence on the highest-power of n l in Eq. (2.3) is only due to the pure poles in the minimal ultraviolet renormalisation constant Z αs , it must be cancelled by the, equally minimal, constant Z. Thus: The scale dependence of C Z , on the other hand, is non-trivial: The conversion between the two infrared schemes is achieved with the help of: Explicitly: (2.24) -6 -For instance, this result allows to obtain Eqs. (2.13) and the scale dependence of C Z after using Eqs. (2.15).

Technicalities
The three-loop diagrams corresponding to the amplitude Eq. (2.1) have been reduced to a set of (master) integrals, M i (z, ), via Integration-By-Parts identities [15] with the help of a C++ implementation [16] of the Laporta algorithm [17]. The same reduction has also been exploited to construct a system of first-order homogeneous linear differential equations [18,19]: where the coefficients A ij (z, ) are rational functions in z and . Truncated -expansions have been subsequently substituted to represent the master integrals. A large-mass expansion (see below) of each M i has been used to determine the lowest power of , n i , with non-vanishing coefficient, while the amplitude and the differential equations have been used to determine the highest power of , n i , necessary to obtain the amplitude at O 0 . Let the coefficients of the truncated -expansions be denoted with I k (z): where k i have been chosen to avoid overlap of the k-indices of the expansion coefficients I k of different master integrals. The coefficients I k satisfy a system of first-order homogeneous linear differential equations derived from Eqs. (3.1): where the coefficients B kl (z) are rational functions in z. Instead of seeking an analytic solution of Eqs. (3.3), we have solved the system numerically as proposed originally in Ref. [20] and first applied to a physical problem in Ref. [21]. To this end, we have used the Boost [22] library odeint. In particular, we have chosen the Bulirsch-Stoer algorithm, bulirsch_stoer_dense_out. In order to keep the numerical precision of the results under control, we have used the Boost library multiprecision with a gmp/mpc backend. The floating point containers were requested to represent 100 decimal digits. A local error of 10 −40 has been requested from the differential equation solution.
The numerical solution of Eqs. (3.3) requires a boundary value for each I k . In order to obtain these, we have used a high-order large-mass expansion, see e.g. [23], around z = 0. The expansion must have unit radius of convergence 2 in z, since the nearest singularity of the master integrals is at z = 1. The expansion has been obtained using diagrammatic methods for the first few coefficients. It has been subsequently extended with the help of  Figure 2: Contours for the numerical solution of the differential equations for the master integrals. The points on the abscissa correspond to singularities of the differential equations. Every time a contour reaches the real axis, the interval between singularities is explored in both directions.
the differential equations. As boundary point, we have chosen z = 1/4(1 + i), well within the radius of convergence. Because of the presence of singularities in the coefficients B kl , we have used evolution contours shown in Fig. 2. An additional solution has also been obtained starting from z = 1/4(0.7 + 0.7i) in order to control the error of the final result.
Having high-precision values of the master integrals allows to obtain expansions around arbitrary points, even around singularities. In the course of the present work, we have obtained threshold and high-energy expansions. They are necessary to evaluate the threeloop coefficient of the form factor in the vicinity of z = 1 and 1/z = 0 respectively. In general, expansions of I k are of power-log type, since an expansion in of the master integrals has already been performed: where l k , m k , m k ∈ Z, and y = √ 1 − z for the threshold expansion, while y = 1/z for the high-energy expansion. In practice, the expansions are truncated at an affordable order considering the available computing ressources. For each I k , only one c k ≡ c klm for some l and m, is necessary to make the solution of Eqs.
In order to obtain c klm and thus also F kl (y), we have used an efficient C++ software that was originally developed for Ref. [24]. Upon choosing a suitable y point where the threshold or the high-energy expansion has excellent convergence, we were able to obtain c k with high precision.  Figure 3: Comparison of the three-loop coefficient of the finite remainder, Eq. (2.9), at n l = 5, L µ = 0 (five massless quarks, renormalisation scale µ 2 = −s), with the default Padé approximation, [6,1], constructed in Ref. [9], as function of z = s/4M 2 with √ s the center-of-mass energy of the Higgs boson and M the mass of the single massive quark. The bands correspond to the uncertainty of the Padé approximation as estimated in Ref. [9]. The lower plot shows the absolute difference between the approximation and the exact result. Also shown is the large-mass expansion (LME) of the three-loop coefficient of the finite remainder truncated at O z 2 , O z 4 and O z 100 . our findings at L µ = 0.
We first note that our result for C (2,1) I agrees perfectly with Ref. [11]. Remains to compare with the Padé approximants of Ref. [9] for C (2) . A comparison for the case of five massless quarks is presented in Fig. 3. We observe that the uncertainty estimates of the approximants are reliable over most of the range of z. Slightly larger deviations are observed for the n l = 0 case as demonstrated in Fig. 4. In order to understand the phenomenological relevance of the difference between the exact result and its Padé approximation for n l = 0, we consider the quantity: as a proxy for the error induced on the partonic cross section. We acknowledge the limitations of ∆ (2,0) in this respect due to the size of the real-radiation corrections to the cross section at higher orders. We expect that the actual effect is about 1/2 of ∆ (2,0) , at least for a top-quark loop. For simplicity, we fix the value of the strong coupling at α s = 0.1. ∆ (2,0) is plotted in Fig. 5. Assuming an off-shell Higgs-boson with a partonic center-of-mass energy, √ s, of up to 1 TeV produced through a top-quark loop, there is ∆ (2,0) < 1%. Hence, the Padé approximant provides an excellent approximation for top-  The domain of physical z values may be compactified with the following mapping: The exact result for C 3/4 ≤ ρ < 1 -high-energy expansion, Appendix C and Fig. 8.

Conclusions and outlook
With the results presented in this work, the Higgs-gluon form factor is known exactly at three loops in QCD with a single massive quark. This is sufficient for applications to  Table 1: Numerical values of the three-loop coefficient of the finite remainder C (2) I at n l = 0, L µ = 0, for 1/4 ≤ ρ ≡ z/(4 + z) < 1/2.
Higgs-boson hadroproduction in the five-flavour scheme, where the massive quark is the top. In this case, we have confirmed that an approach based on Padé approximants [9] is sufficient to obtain sub-percent precision for physical observables. On the other hand, our result removes any uncertainties on the value of the form factor present in Ref. [9]. Once b-quark loops are considered at non-vanishing b-quark mass, our result becomes indispensable, since Padé approximants potentially induce errors on physical predictions in the ten-precent range.
For the presentation of our results, we have used two different infrared-renormalisation schemes. On the other hand, we have chosen to renormalise the Yukawa coupling in the on-shell scheme. Fortunately, a translation to any other scheme, e.g. MS, can be easily achieved thanks to the knowledge of one-and two-loop results in analytic form. This translation is independent of infrared renormalisation.   Table 2: Numerical values of the three-loop coefficient of the finite remainder C (2) I at n l = 0, L µ = 0, for 1/2 ≤ ρ ≡ z/(4 + z) ≤ 3/4.
In principle, our calculation can also be used to obtain the form factor for the process H → γγ, as well as processes involving pseudo-scalars instead of a scalar. We intend to provide these results in forthcoming publications.
Finally, we stress that a complete knowledge of the form factor at three loops in the most general case requires the evaluation of diagrams with two different massive quarks. This can be achieved with numerical methods presented here, for example by fixing the ratio of the b-and top-quark masses. We leave this problem to future work.
Our results are available in computer readable form, see Appendix D.
A Large-mass expansion (a n,0 + a n,1 L s ) z n , The exact expansion coefficients are provided in Ref.
-18 - B Threshold expansion We agree with Ref. [10] for the coefficients of the first three non-analytic terms: