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.

JHEP05(2020)149 In the present publication, we present an exact result for the form factor in QCD with a single massive quark. In particular, we compute the diagrams of figure 1 with both quark loops with the same flavour, as well as the complete set of diagrams with only one massivequark loop. A result in QCD with several massive quarks would still require a calculation of the diagrams figure 1 with massive quarks of different flavour.
Our results are certainly necessary to answer the question whether Padé approximants are indeed sufficient phenomenologically as claimed in ref. [9]. Independently, the knowledge of exact quark mass dependence of the form factor opens the possibility of including b-quark mass effects exactly.
The paper is organised as follows. In the next section, we introduce our conventions and define finite remainders of the form factor after infrared renormalisation. We use this opportunity to provide explicit formulae for the scale dependence of the form factor as well. We subsequently describe the methodology that has allowed us to obtain not only a -2 -JHEP05(2020)149 high precision numerical result but also high-order expansions around the three physical singularities: infinite quark mass (large-mass expansion), intermediate-quark production threshold (threshold expansion) and vanishing quark mass (high-energy expansion). Finally, we present our results and compare them to previous work, in particular, the Padé approximants of ref. [9]. This main text is closed with conclusions and outlook. The three expansions are reproduced in separate appendices. The last appendix presents the contents of the supplementary material that contains our results in electronic form.

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: 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:

JHEP05(2020)149
In the limit M → ∞: 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 Higgsgluon 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:

JHEP05(2020)149
The dependence of the finite remainder on the scale logarithm, L µ , is thus given by the β-function only: 1 C (1) (2.15) 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: 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µ.

JHEP05(2020)149
Explicitly: For instance, this result allows to obtain eqs. (2.13) and the scale dependence of C Z after using eqs. (2.15).
Finally, let us note that our results can be used to obtain the three-loop form factor before factorisation of the infrared divergences with the help of the two-loop result provided to O 2 in ref. [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 [16] with the help of a C++ implementation [17] of the Laporta algorithm [18]. The same reduction has also been exploited to construct a system of first-order homogeneous linear differential equations [19,20]: 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):  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.
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. [21] and first applied to a physical problem in ref. [22]. To this end, we have used the In order to obtain these, we have used a high-order large-mass expansion, see e.g. [24], 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 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 figure 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 JHEP05(2020)149 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. (3.3) unique. Since eqs. (3.3) are linear, there is: 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. [25]. 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.

Results
Since the scale logarithms of the three-loop coefficient of the finite remainder are entirely determined from the analytically known lower order results, see eqs. (2.15), we only present 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 figure 4. 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 figure 5. An improvement of the Padé approximants has recently appeared in the proceedings [26]. The respective plots are also shown in figures 4 and 5. Clearly, the agreement with the exact result is worse for n l = 5 and better for n l = 0.
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 figure 3. Assuming an off-shell Higgs-boson with a partonic center-ofmass 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-quark loops. On the other hand, in the case of the production of an on-shell Higgs boson through a b-quark loop, ∆ (2,0) ≈ 10%. Furthermore, the difference grows rapidly with the Higgs-boson offshellness, √ s. Hence, the approximation is rather poor for b-quark loops. In the same figure, we also show ∆ (2,0) using the improved Padé approximant of ref. [26]. We note that   Figure 4. 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] (left panel) and improved to [7,1] in ref. [26] (right panel), 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é approximations as estimated in refs. [9] and [26]. 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 . the approximation is now better for b-quarks. Nevertheless, ∆ (2,0) > 10% for an off-shell Higgs boson of 400 GeV.
Our exact result is a sample of C The exact result for C 3/4 ≤ ρ < 1 -high-energy expansion, appendix C and figure 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 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.     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. achieved thanks to the knowledge of one-and two-loop results in analytic form. This translation is independent of infrared renormalisation. 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.  The exact expansion coefficients are provided in the supplementary material. We agree with refs. [6,7] up to O z 4 and with ref. [9] up to O z 6 .
B Threshold expansion