The analytic leading color contribution to the Higgs-gluon form factor in QCD at NNLO

We present an analytic three-loop result for the leading color contribution to the Higgs-gluon form factor in QCD. The leading color contribution is given at next-to-next-to-leading order by the $N_c^2$-term in QCD with $N_c$ colors. The main focus of this article lies on the evaluation of the relevant Feynman integrals with a special emphasis on the elliptic sector.


Introduction
A precise knowledge of the gluon-gluon-Higgs (ggH) form factor with arbitrary quark masses and virtualities is essential for many Higgs physics applications. For example, it is a main ingredient in the total cross sections for single and double Higgs production at the Large Hadron Collider (LHC). Therefore, the evaluation of the ggH form factor in QCD has a long history: The initial groundwork was done by completing the exact two-loop result for arbitrary quark masses in Refs. [1][2][3][4]. First computations at three loop level go back to large-mass expansions 1 in the top quark mass [5,6]. More recent results exploit the partial knowledge of the threshold behavior of the form factors [7] and apply the techniques of Padé approximations [8], which in principle is sufficient for phenomenology. Subsequently, the applicability of the Padé approximations for top quark induced loops was proven shortly after in Ref. [9]. Furthermore, using the latter result the b-quark mass effects could be included exactly for the first time.
In contrast to the purely numerical result of Ref. [9], we choose a different approach in this article. We want to continue our work started in Ref. [10] towards a full analytic result for the ggH form factor in QCD with a single massive quark by completing the calculation for the leading color contribution.
Nowadays, the integration-by-parts reductions for three loop form factors are no longer a bottleneck and are manageable on a laptop within a few weeks. The following evaluation of master integrals can be grouped into two categories. The master integrals that can be expressed in terms of multiple polylogarithms are by now very well understood and their evaluation can be automated to a high degree.
The master integrals of the leading color contribution to the ggH form factor contain only a single elliptic sector. In order to give analytic results for this sector a new class of iterated integrals has to be introduced with integration kernels depending on complete elliptic integrals of the first kind. The results in this article can be evaluated numerically via an included implementation of these iterated integrals.
In section 2, the computational setup and toolchain of the calculation is described. Section 3 is devoted to the evaluation of the various types of master integals. In section 4, we present the final result of this paper. Finally, our conclusion and an outlook to possible future work is given in section 5.

Computational Setup
The computational setup of the leading color contribution to the ggH form factor consists of multiple steps. At first, the complete set of Feynman diagrams in QCD with two incoming gluons and one outgoing Higgs boson is generated with the computer program qgraf [60]. Sample diagrams are shown in fig. 1. The output files of qgraf are further processed by q2e [61,62] in order to insert Feynman rules in Feynman gauge. In a next step, the Feynman diagrams are mapped to a predefined set of topologies via the tool exp [61,62], which generates output files for the computer algebra system FORM [63].
The tensor structure of the Feynman amplitude M ab;µν is restricted by Bose symmetry and Ward identities to where q µ 1 and q ν 2 are the components of the momenta of the external gluons. Only the form factor C is relevant for physical quantities. It can be projected out via C = δ ab [(q 1 · q 2 )g µν − q 2µ q 1ν ] N A (q 1 · q 2 ) 2 (2 − d) M ab;µν , (2.2) where N A is the number of gauge generators and d = 4 − 2 the number of space-time dimensions.
A custom FORM-code implements eq. (2.2) and expresses C in terms of scalar Feynman integrals. The individual group theory factors of the Feynman diagrams are evaluated with the FORM package color [64], which expresses the color factors in terms of the Casimir invariants C F , C A and d abc F d abc F . In the QCD gauge group SU (N c ) these Casimir invariants have the values with N A = N 2 c − 1 and the normalization T F = 1/2. The scalar Feynman integrals are then reduced to 409 master integrals via integration-byparts relations [65,66] using Laporta's algorithm [67] implemented in the C++ program Kira [68,69].
In the following, we restrict ourselves to the leading color contribution that is given by the N 2 c -term of the form factor C after inserting eq. (2.3). Only 196 master integrals contribute to this term.

Master integrals
In order to study the kinematic dependency of master integrals the method of differential equations [70][71][72][73] turned out to be an indispensable tool. The coupled system of differential equations can be simplified significantly if it is expressed in terms of an appropriate kinematic variable. For our purposes, we choose in accordance with Ref. [10] the variable x related to the variable z = M 2 H /(4m 2 q ) + i0 by, with the mass M H of the Higgs boson and the mass m q of the massive quark 1 . Arranging all master integrals in a vector F (x, ), where the integrals are ordered by an increasing number of propagators in such a way that integrals of the same sector stay together, leads to a system of differential equations whereM (x, ) is a block-triangular matrix. The differential equations of a certain sector can therefore be written as where f (x, ) contains only the master integrals of the sector under consideration and g(x, ) the master integrals with a smaller number of propagators. Eq. (3.3) is the starting point for all the calculation methods described in the following subsections. Solving the master integrals in F (x, ) sector by sector always allows us to consider g(x, ) in eq. (3.3) as already known.
The system of differential equations in eq. (3.3) fixes the master integrals up to some boundary conditions. A suitable boundary condition for all master integrals under consideration is an expansion around x ≈ 1 corresponding to a large quark mass m q compared 1 In Ref. [10] the variable z was named τ .
The analytic results for the master integrals derived in the context of this article have all been checked against numerical results obtained by sector decomposition [78][79][80] at an euclidean point x = 0.01 via the tool FIESTA [81].
The method treated in section 3.1 is applicable to the vast majority of master integrals under consideration, i.e. the master integrals expressible in terms of generalized polylogarithms (GPLs) with argument x. The two sectors described in section 3.2 can also be expressed in terms of GPLs but require an argument of √ x. Section 3.3 covers the only elliptic sector present in the leading color contribution. Section 3.4 deals with two further sectors, where the homogenous parts of eq. (3.3) are of the type described in Section 3.1 but the inhomogeneities g(x, ) contain integrals of the two other types.

Canonical sectors in x
178 of the 196 master integrals can be expressed completely in terms of GPLs of argument x. Samples of master integrals falling in this category are shown in fig. 2. The differential equations, eq. (3.3), of these so-called canonical master integrals can be cast into an -form [82] d via a rational transformation, where the matrix M (x) does not depend on anymore. A basis transformation that transforms a system of differential equations of the form eq. (3.3) into an -form of eq. (3.4) can be constructed with the tool epsilon [83], which implements an algorithm presented in Ref. [84].
The matrix M (x) of eq. (3.4) is fuchsian, i.e. it has only simple poles in x and thus can be written as The master integrals under consideration have singularities at with a sixth root of unity r = e iπ/3 = (1 + i √ 3)/2 and the golden ratio ϕ = (1 + √ 5)/2.
In order to simplify the calculation even further a choice of basis has been made in which the matrix B(x, ) in eq. (3.4) has only simple poles in x as well.
The evaluation of the master integrals in this basis can be done in analogy to the master integrals of Ref. [10]. We introduce an evolution operator U (x, x 0 ; ) which fulfills Making an ansatz leads to a recursive solution where we used eq. (3.5). Comparing this to the definition of GPLs [85], reveals that the integrals in eq. (3.9b) evaluate to GPLs of arguments x and x 0 and indices in S.  Figure 3: Canonical sectors in √ x. The notation is the same as for fig. 2.
The full solution of the integrals in f (x, ) are now given by where f (x 0 , ) are the boundary conditions around x 0 = 1, i.e. the asymptotic heavyquark expansions of the master integrals.
Since the integrals in g(x , ) are expressed in terms of GPLs of argument x (without rational function prefactors) and B(x , ) has only simple poles in x , the integrals evaluate to linear combinations of GPLs of argument x as well. Hence, the solutions f (x, ) are suitable to serve as input into the inhomogeneous parts of higher canonical sectors.

Canonical sectors in √ x
The differential equations in x of the two sectors shown in fig. 3 cannot be transformed into -form. The sector in fig. 3(a) consists of three master integrals and the sector in fig. 3(b) of only a single master integral. Transforming the homogeneous parts of the differential equations of those two sectors to fuchsian form using epsilon reveals halfinteger eigenvalues of the residue matrices at 0 and ∞ indicating that a change of variable is required 2 .
The change of variable x = y 2 transforms eq.
This reparameterization effectively multiplies the residues at 0 and ∞ by a factor of two, eliminating the half-integer eigenvalues, while all other singular points of M (x, ) are split into two distinct points due to the partial fraction decomposition 1 (3.13) The differential equations in eq. 3.12 are now in a form applicable to the tool epsilon in order to construct a basis transformation into -form.
At this point, the master integrals can basically already be solved in terms of GPLs of argument y = √ x in an analogous manner to section 3.1 but there is a technical problem. The master integrals of lower sectors entering eq. 3.12 through g(y, ) are expressed in terms of GPLs of argument x and have to be reparameterized by the variable y first. It is straightforward to rewrite GPLs of argument x = y 2 into GPLs of argument y by applying recusively, where we used the substitution t = s 2 and a partial fraction decomposition. Unfortunately, this transformation leads to a huge increase of the expression size. Since every GPL generates 2 m terms, where m is the number of a i = 0, available computing resources quickly become insufficient.
A solution to this problem is to introduce a new type of multiple polylogarithms, defined by 3 P 1,··· ,1 0,··· ,0 n times (y) = 2 n ln n (y) n! , P n 1 ,··· ,n N a 1 ,··· ,a N (y) = y 0 dt f n 1 a 1 (t)P n 2 ,··· ,n N a 2 ,··· ,a N (t) , On the one hand, this definition allows rewriting every GPL in x as a single P -function in y via G a 1 ,··· ,a N (x) = P 1,··· ,1 a 1 ,··· ,a N (y) . On the other hand, GPL-like integration kernels appearing in eq. (3.12) can always be 3 Note that the kernel f 0 0 (t) = 2/t 2 has a second order pole and thus should be excluded. Figure 4: Elliptic sector. The notation is the same as for fig. 2. expressed in terms of the kernels defined in eq. (3.15b) using the equation Eq. (3.9) and eq. (3.11) can now be adapted to the integration kernels of eq. (3.15b) and allow solving the master integrals in terms of the P -functions defined in eq. (3.15).
In order to use existing implementations for numerical evaluation (e.g. Refs. [87][88][89]), one can easily rewrite the P -functions in terms of GPLs via recursive application of (3.18)

Elliptic sector
The most complicated sector of the calculation is shown in fig. 4. The master integrals of this sector cannot be expressed solely in terms of GPLs but involve elliptic integrals.
As a starting point for this elliptic sector, we define the basis integrals and a normalization factor Since the internal mass has been set to unity in eq. (3.20), the scalar products of the external momenta are given by with z defined above eq. (3.1). Even though the integrals I 1 (x, ), · · · , I 6 (x, ) are already finite for → 0, it is more convenient to proceed in a different basis where the differential equations are fuchsian and the eigenvalues of the residues are free from resonances. Such a basis is given by  . . . where It turns out that after expanding the functions a j (x, ) in , only the seven coefficients a (0) 6 (x) and a 3 (x) contribute to the leading color form factor. Since the homogeneous part of the differential equation for a 3 (x, ) vanishes for → 0, the coefficient a where the matrix M (x) and the inhomogeneity b(x) can be found by expanding the differential equations for the a j (x, ) in .
For later convenience, we introduce yet another variable t by 4 and define once again a new basis by After transforming the differential equation to this basis, it is straightforward to read off the third-order differential equation fulfilled by c 1 (t),    3 (t), vanishes completely which makes it almost trivial to solve for these coefficients after c 1 (t) has been found.
In Section 3.3.1 we describe how to solve the homogeneous part of eq. (3.29), while in Section 3.3.2 we include the inhomogeneity k(t) into the calculation which requires the introduction of a new class of iterated integrals.

Homogeneous solution
As a first step in order to evaluate the master integrals of the elliptic sector, the homogenous differential equation, eq. (3.29) with k(t) = 0, has to be solved. This third-order differential equation is a symmetric square of a second-order differential equation. Thus, a complete set of three independent homogeneous solutions of eq. (3.29) is given by 1,h (t) = ψ 1 (t)ψ 2 (t) , c where ψ 1,2 (t) solve which can be easily checked with the MAPLE 5 package ReduceOrder [90].
The differential equation (3.31) has solutions in terms of Gauss' hypergeometric functions, found with another MAPLE package hypergeometricsols [91]. The solutions eq. (3.32) are related to the complete elliptic integral of the first kind K(z) via Unfortunately, the function ψ 1 (t) in eq. (3.32a) turns out to be discontinuous at arguments that map onto the branch cut of the 2 F 1 -function. The position of this discontinuity (expressed through the variable x related to t by eq. (3.27)) is depicted in fig. 5, where we restrict the discussion to the physically relevant upper complex half-plane. A continuous solution can be defined by ; inner region , ψ 1 (t) + 2ψ 2 (t) ; outer region . (3.34) An easy to implement criterion for the inner region is given by  with α = t(4 + t) 5 (t 2 + 8t + 20)(t 2 + 6t + 4) 2 .

Inhomogeneity
In this section, we discuss the solution to eq. (3.29) including the inhomogeneity k(t).
Since at the boundary x 0 = 1 (or equivalently t 0 = −4 − 2i), the function c 1 (t) has to vanish, the solution for c 1 (t) becomes which can be obtained by a variation of constants. The function w(t) is the Wronskian of the second-order differential equation eq. (3.31) and is given by .
When solving the integrals in eq. (3.37), one encounters a new class of iterated integrals, defined by The integration kernels f j (x) that appear in this article are 6 .
Note that the integration in eq. (3.39) is defined over the variable x. The variable t ≡ t(x) in eq. (3.40) is given by eq. (3.27b).
The notation [x 0 ] at the lower limit of eq. (3.39) indicates that a regularization procedure might be required. In general the integral has an expansion around a point x 1 in powers of √ x − x 1 and ln(x − x 1 ). Hence, the integrals are regularized in such a way that the constant term of the expansion around x 0 has to vanish. In order to not further complicate the calculation, we set x 0 = 1 for the remainder of the paper and neglect it in the notation, but we note that for x 0 = 0, this regularization method is compatible with the standard procedure of GPLs with trailing zeros. Thus defined, the functions also obey the usual shuffle relation of iterated integrals. Details to a numerical evaluation of this class of iterated integrals are deferred to appendix A.
(a) (b) Figure 6: Mixed sectors. The notation is the same as for fig. 2.
After changing the variable in eq. (3.37) from t to x, one can insert eq. (3.41) and rewrite the HPLs in k(t) in terms of the new iterated integrals (i.e. with a lower limit x 0 = 1). The shuffle relation of iterated integrals leads then to integrals of the form of eq. (3.39) with integration kernels from eq. (3.40).
The coefficients a 3 in eq. (3.28a) require an additional integration, but in doing so no further kernels have to be introduced.

Mixed sectors
The differential equations of the sectors depicted in fig. 6 have homogeneous parts which can be transformed into -form as described in section 3.1. Unfortunately, the inhomogeneities of both sectors depend on master integrals of the canonical sectors in √ x and the sector in fig. 6(b) also depends on the elliptic master integrals.
The differential equations can therefore be written as where g x (x, ) only contains the canonical master integrals in x (section 3.1), g y (x, ) the canonical master integrals in y = √ x (section 3.2), and g E (x, ) the elliptic master integrals (section 3.3).
We can now split the master integrals into three parts, so that each part obeys a differential equation where the inhomogeneity only depends on one type of master integrals, i.e.
The boundary conditions could be assigned completely to f x (x, ), This setup allows to solve eq. (3.44a) as described in section 3.1, eq. (3.44b) as in section 3.2 and eq. (3.44c) as in section 3.3.
We want to point out, that solving for f E (x, ) in terms of the iterated integrals, eq. (3.39), does not lead to additional integration kernels besides the ones defined in eq. (3.40).

Result
In this section, we fix the notation of the final result and present it in numerical form.
For the sake of legibility, the lengthy analytic result has been deferred to appendix C. The calculation has been done for a SU (N c ) gauge group; QCD is reached for N c = 3.
The form factor C, defined in eq. (2.2), has a power series expansion in the MS renormalized strong coupling constant α s , with the vacuum expectation value v. The quark mass has been renormalized in the on-shell scheme.
Starting at two-loop, the contributions are still infrared divergent. The universal structure of these divergencies can be subtracted via [92] C 0 = C 0 , (4.2a) where the factors I (1) g and I (2) g are given by [92,93] and We expand the finite contributionsC i in powers of N c as In order to fix the notation, we provide the one-loop result, The leading color contributionC (2) 2 is shown in fig. 7. The analytic result will be presented, along with the two-loop results, in appendix C and in electronic form in an ancillary file, see appendix B.1.
We checked analytically that our result agrees with the heavy-quark expansion of Refs. [6,8,94] . The non-analytic terms of the threshold expansion, provided in Ref. [7], could also be reproduced.

Conclusions and outlook
In this paper, we presented for the first time an analytic result for the leading color contribution to the three-loop ggH form factor in QCD with a finite mass of the mediating quark. Previously at three-loop, analytic expressions were only available for contributions involving one light-fermion loop [10]. A full result that includes also sub-leading color terms is still only known in numerical form [9].
We showed that our result requires the introduction of a new class of iterated integrals, eqs. (3.39, 3.40), with integration kernels involving elliptic integrals. These iterated integrals can be evaluated numerically to high precision, see Appendix A.
For the sub-leading color contributions, we expect more elliptic sectors to appear. If these sectors have homogeneous solutions that can be solved as easily as the sector discussed in this paper, it should be possible to provide analytic results for these contributions as well.
This would most certainly require the introduction of a larger set of integration kernels for the iterated integrals.

Acknowledgments
We

A Numerical evaluation of iterated integrals
In this appendix, we discuss the numerical evaluation of the iterated integrals defined in eqs. (3.39, 3.40).
It is possible to derive a recursive algorithm for a series expansions for an iterated integral around an arbitrary point x 1 in powers of √ x − x 1 and ln(x − x 1 ), where almost all coefficients are expressed through coefficients of an iterated integral with one integration less. Only the constant term in this series expansion cannot be determined that way.
The strategy for a numerical evaluation is to start at x 0 , where the constant terms vanish by definition, and fit the constant term of an expansion around a point x 1 , not too far away, by evaluating both series at a point in between that lies within the convergence radii. By repeating this procedure for points along a path, every point of interest can be reached.
Section A.1 explains the basics for a series expansion of iterated integrals. Section A.2 describes the expansion of the integration kernels in eq. (3.40).
A C++ implementation is included in the supplementary material and will be discussed in Section B.2.

A.1 Series expansion of iterated integrals
As an ansatz, we assume that eq. (3.39) has a series expansion around x = x 1 of the form while the integration kernels expand to Shifting the integration variable in eq. (3.39) yields Inserting eqs. (A.1, A.2) leads to The integral is solved by The constant in eq. (A.5) depends on the lower limit and contributes only to the coefficient c 1) that cannot be expressed through lower weight coefficients, anyway. Hence, eq. (A.4) becomes This equation allows us to read off the coefficients of eq. (A.1), The initial conditions for the recursion are given by the expansion coefficients of the weight zero function I(x 0 , x) = 1,

A.2 Series expansion of integration kernels
Since the expansion coefficients of the integration kernels, eq. (3.40), enter eq. (A.7), it is disirable to find efficient approaches to determine these quantities. The GPL-like kernels, eq. (3.40a), are trivial to expand. Unfortunately, this is not true for all the other integration kernels as they involve the squared homogeneous solution ψ 2 (t) given in eq. (3.34).
Instead of ψ 2 (t), we consider the function which fulfills a homogeneous third-order differential equation in the variable x, This differential equation can be solved with a power series ansatz, that is sufficient for all points x 1 .
Furthermore, the integration kernels involve rational functions in t, which have to be expanded in x. We perform these expansions in three seperate steps. The expansion of such a function in t as well as the expansion of t in x can be performed in a straightforward manner. In a third step, these expansions have to be composed using Faà di Bruno's formula [102].
Via the Cauchy product, the expansion of φ(x) can be combined with the x-expansion of a rational function in t in order to form the integration kernels in eq. (3.40d) and eq. (3.40c). The kernel τ (x), eq. (3.40b), requires in addition the construction of a reciprocal power series.

B Supplementary material
As supplementary material to this article, we provide the main result in electronic form, see section B.1. In addition, we present in section B.2 an implementation of the algorithm discussed in appendix A for a numerical evaluation of the iterated integrals defined in eqs. (3.39, 3.40).

B.1 ggh-nc2.m
The ancillary file ggh-nc2.m contains the results of this article in a Mathematica 7 readable format using the following notation:

B.2.1 Installation
The libiint package provides a CMake 9 build system to compile and install the library and the Mathematica interface.
The following dependencies have to be properly installed before the compilation can be started: • FLINT: Fast Library for Number Theory 10 , • Arb -a C library for arbitrary-precision ball arithmetic 11 , • yaml-cpp: A YAML parser and emitter in C++ 12 .
Obviously, Wolfram Mathematica is required in addition to build the Mathematica interface. where the second line may require root privileges depending on the installation paths. You might also have to adjust Mathematica's $Path variable so that it includes the installation directory of the Mathematica interface.

B.2.2 Mathematica interface
The library libiint obtains its usefulness only through the Mathematica interface, which is why we refrain from discussing direct access to the library in this article.
After a successfull installation, you should be able to load the package via The package offers a few commands to the user: • IIntInit[prec] Initialize libiint and set the working precision to prec bits. Verbose output can be enabled with the option Verbose->True.

• IIntCreate[expression]
Register all IInts in expression and replace them by IIntObj objects.
• IIntMatchEuclidean[q,a] Match the constant terms of all registered IInts along a path in the euclidean region from 1 to 0. For a point x 1 where the constant term is already known, the next point along the path x 2 is chosen so that the matching point x fulfills |x−x 1 | ≤ min(qr 1 , a) and |x − x 2 | ≤ min(qr 2 , a), where r 1 (r 2 ) is the convergence radius for an expansion around x 1 (x 2 ).
• IIntMatchPhysical[q,a] Same as IIntMatchEuclidean but for a path in the physical region.

• IIntSave[filename]
Save all constant terms to filename in the YAML file format.

• IIntLoad[filename]
Load all constant terms from filename previously saved by IIntSave.

• IIntEvaluate[expression]
Evaluate all IIntObjs in expression with numeric arguments x. x must lie within the convergence radius of the closest point where a constant term is available. of eq. (C.5) at various kinematic points. The first two scripts are used to generate the YAML files containing the constant terms. If one wants to evaluate the form factor only at the predefined points, this step that usually takes a few hours can be skipped since the package already includes the files euclidean.yaml and physical.yaml with all points irrelevant for the evaluation removed in order to save disk space.
The GPLs in eq. (C.5) are evaluated via PolyLogTools [103] and GiNaC [104], which have to be installed as well. In order to run the scripts in the eval/ directory, the file ggh-nc2.m (see section B.1) has to be copied (or symlinked) to the eval/ directory. The eval.m script is run with 13 math -script eval.m and should produce an output of the form: The evaluation point x is given at the far left, while the result evaluated at x is given after the arrow.

C Analytic result
In this appendix, the analytic result for the coefficientC (2) 2 is presented along with the two-loop contribution.
In order to express the result in a more compact form, we introduce a few short-hand notations. For the three types of iterated integrals discussed in this article, we neglect the arguments, i.e.
Furthermore, we add another integration kernel, to the set of standard GPL integration kernels, where r denotes a sixth-root of unity r = e iπ/3 ; i.e. GPLs with some indices beingr can be interpreted as the sum of all possible GPLs, wherer = r orr = r * , e.g.