Two-loop planar master integrals for Higgs$\to 3$ partons with full heavy-quark mass dependence

We present the analytic computation of all the planar master integrals which contribute to the two-loop scattering amplitudes for Higgs$\to 3$ partons, with full heavy-quark mass dependence. These are relevant for the NNLO corrections to fully inclusive Higgs production and to the NLO corrections to Higgs production in association with a jet, in the full theory. The computation is performed using the differential equations method. Whenever possible, a basis of master integrals that are pure functions of uniform weight is used. The result is expressed in terms of one-fold integrals of polylogarithms and elementary functions up to transcendental weight four. Two integral sectors are expressed in terms of elliptic functions. We show that by introducing a one-dimensional parametrization of the integrals the relevant second order differential equation can be readily solved, and the solution can be expressed to all orders of the dimensional regularization parameter in terms of iterated integrals over elliptic kernels. We express the result for the elliptic sectors in terms of two and three-fold iterated integrals, which we find suitable for numerical evaluations. This is the first time that four-point multiscale Feynman integrals have been computed in a fully analytic way in terms of elliptic functions.


Introduction
At the Large Hadron Collider (LHC), the main production mode of the Standard Model (SM) Higgs boson is via gluon-gluon fusion. The Higgs boson does not couple directly to the gluons, the interaction being mediated by a heavy-quark loop. That makes the evaluation of the radiative corrections to Higgs boson production via gluon-gluon fusion challenging, since the Born process is computed through one-loop diagrams, the next-to-leading order (NLO) QCD corrections involve the computation of two-loop diagrams, the next-to-nextto-leading order (NNLO) corrections the computation of three-loop diagrams, and so on. In fact, fully inclusive Higgs production is known up to NLO [1,2], while Higgs production in association with one jet [3] and the Higgs p T distribution [4] are known only at leading order.
The evaluation of the radiative corrections simplifies considerably in the Higgs effective field theory (HEFT), where the heavy quark is integrated out and the Higgs boson couples directly to the gluons, effectively reducing the computation by one loop. For fully inclusive Higgs production, the HEFT is valid when the Higgs mass is smaller than the heavy-quark mass, m H m Q . Thus it is expected to be a good approximation to the full theory (FT), which gets corrections from the top-mass contribution and from the top-bottom interference. In fact, using the FT NLO computation as a benchmark, one can see that the HEFT NLO computation approximates very well the FT NLO computation, since the top-bottom interference and the top-mass corrections are about the same size although with opposite sign [5]. At NNLO, the FT mass corrections are expected to be in the percent range, which is though competitive with the precision of the HEFT computation at next-to-next-to-next-to-leading order (N 3 LO) [6,7].
For Higgs production in association with one jet or for the Higgs p T distribution, using the leading-order results [3,4] as a benchmark one can show that the HEFT is valid when m H m Q and the jet or Higgs transverse momenta are smaller than the heavy-quark mass, p T m Q [8,9]. In the HEFT, Higgs production in association with one jet [10,11] and the Higgs p T distribution [12] are known at NNLO. No complete FT results are known beyond the leading order. Approximate NLO top-mass effects have been computed, and shown to be small and to agree well with the HEFT for p T m top [13][14][15] and up to p T ∼ 300 GeV [16]. However, they are expected to be non-negligible in the high p T tail. Finally, it is worth noting that in many New Physics (NP) models, the high p T tail of the Higgs p T distribution is sensitive to modifications of the Higgs-top coupling [17][18][19].
In this paper, we report on the analytic computation of all the planar master integrals which are needed to compute the two-loop scattering amplitudes for Higgs→ 3 partons, with full heavy-quark mass dependence. These are relevant to compute the FT NNLO corrections to fully inclusive Higgs production and the FT NLO corrections to Higgs production in association with one jet or to the Higgs p T distribution.
The differential equations method [20][21][22][23][24] has proven to be one of the most powerful tools to compute (dimensionally regularized) loop Feynman integrals. In particular, the reduction of the Feynman integrals to a set of linearly independent integrals, dubbed master integrals [25][26][27][28], through integration-by-parts identities, the exploration of new classes of special functions such as multiple polylogarithms [29,30], and a better understanding of their functional properties [31][32][33], have made the technique increasingly efficient. However, until recently the method was mostly applied in relatively simple kinematic situations, with the Feynman integrals depending on few scales, while complicated integrals needed a caseby-case analysis.
A major breakthrough was made in [34], where a canonical form of the differential equations for Feynman integrals was proposed. A key idea is that the canonical basis can be found by inspecting the singularity structure of the loop integrand. More precisely, one computes the leading singularities, i.e. maximal multidimensional residues of the loop integrand [35,36]. The fact that this can be done before the differential equations are set up renders this technique extremely efficient 1 . When considering differential equations for a set of integrals defined to be pure functions of uniform weight, all relevant information about the analytic properties of the result is manifest at the level of the equations. Moreover, it is possible to find an analytic expression for the master integrals in terms of iterated integrals over algebraic kernels in a fully algorithmic way, up to any order of the dimensional regularization parameter (see [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] for many applications of these ideas). It is important to note that these ideas also streamline calculations whose output cannot be immediately written in terms of multiple polylogarithms, but where Chen iterated integrals [56] are the appropriate special functions, see e.g. [51,57]. This class of functions will also be important in this paper.
Beyond Chen iterated integrals, there are cases where elliptic functions appear. This is typically related to several equations being coupled in four dimensions, see e.g. [57,58]. The appearance of elliptic functions can be also anticipated by inspecting the maximal cuts of the corresponding loop integrands [59]. In this case the precise form of the canonical basis is not yet known, and presumably finding it will involve a generalization of the concept of leading singularities.
Over the last two decades a lot of effort has been made to understand the analytic properties of Feynman integrals which go beyond the multiple polylogarithms case, mostly related to the so-called sunrise diagram [38,[60][61][62][63][64][65][66][67][68][69]. However, to the best of our knowledge, such a generalized class of Feynman integrals has not been used so far in a fully analytic computation of a four-point multiscale scattering amplitude. In this paper, we compute in the Euclidean region all the planar master integrals relevant for Higgs→ 3 partons, retaining the full heavy-quark mass dependence, which include two elliptic integral sectors.
We write down the differential equations following the approach of [34]. We find that most integrals can be expressed in terms of Chen iterated integrals [56]. The corresponding function alphabet depends on three dimensionless variables and contains 49 letters, underlining the complexity of the problem. Having a fast and reliable numerical evaluation in mind, we derive a representation of all functions up to weight two in terms of logarithms and dilogarithms. Following [57], this allows us to write the weight three and four functions in terms of one-fold integral representations. We find the latter suitable for numerical evaluation. We show that the two remaining integral sectors involve elliptic functions. We analyze the corresponding system of coupled equations, and solve them in a suitable variable. An important tool is to reduce the problem to a one-variable problem (a similar strategy has been used in [70] to effectively rationalize the alphabet of multiscale processes). The solution at any order in ǫ can be expressed in terms of iterated integrals involving elliptic kernels. We then show that using auxiliary bases and basis shifts, the result for the elliptic sectors can be expressed in terms of two and three-fold iterated integrals, which we find suitable for numerical evaluation.
The outline of the paper is as follows. In section 2 we briefly discuss the reduction to the master integrals and the kinematics of the processes under consideration. In section 3 we review the differential equations method in the context of pure functions of uniform weight, i.e. the canonical basis approach. In section 3.2 we show that when a canonical basis exists the solution can be expressed to all orders of the dimensional regularization parameter in terms of multiple polylogarithms, also when a rational parametrization of the alphabet is not possible. We derive a one-fold integral representation of the result up to weight four which is suitable for fast and reliable numerical evaluation. In section 4 we discuss in detail how to analytically solve the elliptic sectors in terms of iterated integrals over elliptic kernels. In section 5 we discuss the class of functions used to represent the elliptic sectors. In section 6 we conclude and discuss future directions. We also provide six appendices in which we collect more details about the calculation. In appendix A we write the explicit expressions for the canonical form of the master integrals, or conversely for the basis choice in the elliptic case. In appendix B we show the 125 master integrals in the pre-canonical form. In appendix C, we give the alphabet for the master integrals. In appendix D we list the dilogarithms we used to express the master integrals at weight two.
In appendix E we give more details about the one-fold integral representation in terms of which we express the master integrals not depending on elliptic functions. Finally in appendix F we show that the maximal cut of the six-denominator elliptic sector provides useful information about the class of functions which characterise the sector.

Notations and conventions
The leading order QCD contribution to Higgs decay to three partons, or alternatively to Higgs production in hadronic collisions, is a process mediated by a loop of heavy quarks. This is due to the fact that the SM Higgs boson does not couple directly to massless particles. The decay channels are H → ggg and H → gqq; the production channels are gg → gH, gq → qH and qq → gH. The one-loop Feynman diagrams for all these processes can be described using the four-denominator topology 2 (and subtopologies) depicted in fig. 1. Figure 2: Planar seven-denominator topologies for the NLO contribution to the cross section of Higgs boson production in association with a jet in proton collisions, with full heavy-quark mass dependence.
At NLO in α S , Feynman diagrams with up to seven propagators contribute to the processes above. They can all be described using the eight different planar seven-propagator topologies (and their subtopologies) depicted in fig. 2. We parametrized all eight topologies into nine-propagator integral families and we reduced the corresponding dimensionally regularized integrals to a minimal set of independent integrals, dubbed master integrals, using the computer program FIRE [71][72][73] combined with LiteRed [74]. The list of denominators defining the integral families and additional details about this part of the calculation are provided in appendix A.
The most general integral is defined in D = 4 − 2ǫ space-time dimensions as, I i a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 ,a 8 ,a 9 = where i is the family index, and a i are integers. The reduction process leads to a set of 125 master integrals, shown in fig. 3, that may be of relevance to more than one physical process. We shall focus here on a Higgs boson decaying to three partons and on Higgs+jet production. These processes differ by the physical phase-space region. Defining, where p 2 1 = p 2 2 = p 2 3 = 0, the relevant physical regions are both with the internal quark mass m 2 > 0. The integrals are functions of three dimensionless invariants, Figure 3: Master integrals in pre-canonical form. Internal plain thin lines represent massless propagators, while thick lines represent the top propagator. External plain thin lines represent massless particles on their mass-shell. External dashed thin lines represent the dependence on s, t, or m 2 H . The external dashed thick line represents the Higgs on its mass-shell. The squared momentum p 2 can assume the values p 2 = s, t, m 2 H . The squared momentum q 2 can assume the values q 2 = s, m 2 H . The squared momentum r 2 can assume the values r 2 = s, t.
In this paper we evaluate the integrals in the Euclidean region where no branch cuts are present, or rather in the subset there-of which has, It is then possible to analytically continue the result to the physical region using the Feynman prescription, by assigning a positive infinitesimal imaginary part to the external invariants and a negative infinitesimal imaginary part to the internal masses. The analytic continuation of the master integrals will be provided elsewhere.
The full basis of master integrals we evaluated in this paper is listed in appendix A. The explicit results for the master integrals require about 200 MB to be stored in electronic form, and can be obtained upon request to the authors.

Differential equations
In order to analytically compute the master integrals we rely on the differential equations method [20][21][22][23][24]. All the integrals discussed in this paper can be expressed in terms of multiple polylogarithms except eight of them, which involve elliptic functions. In the polylogarithmic case we find a modified basis of integrals that are pure functions of uniform weight [34]. In this basis the differential equations take a canonical form and can be readily solved. This basis is found by choosing integrals with constant leading singularities. In the case of elliptic functions the appropriate generalization of the notion of leading singularity has not yet been worked out. It is nevertheless possible to choose a basis where the elliptic nature of the integrals is manifest and the problem can be reduced to the solution of second order differential equations, as we discuss in section 4.

General features of differential equations for Feynman integrals
Denoting a set of N basis integrals by f , the set of kinematical variables by x, and working in D = 4 − 2ǫ dimensions, it is possible to define a system of first order linear differential equations for the integrals, that can be written in total generality as, where we used the shorthand ∂ m = ∂/∂x m , and A m (x, ǫ) is an N × N matrix with rational entries of its variables. The matrix A m (x, ǫ) satisfies the integrability condition, The choice of the basis is not unique. Performing a basis change f → Bf the system of differential equations transforms according to In [34] it was conjectured that performing a basis change with algebraic coefficients, for integral sectors expressible in terms of multiple polylogarithms, it is possible to factorize out the ǫ dependence of the differential equations, Such a system of differential equations is said to be in canonical form. In order to discuss the properties of the solution it is convenient to write the differential equations in differential form, whereÃ is a matrix such that, The matrix elements ofÃ(x) are Q-linear combinations of logarithms. The arguments of the logarithms are known as letters, while the set of linearly independent letters is known as alphabet. The main virtue of the canonical system of differential equations is that its solution is elementary, and it can be written for general ǫ in terms of a path-ordered exponential, where P is the path ordering operator along the integration path C, connecting the boundary point to x, while f (0, ǫ) are boundary conditions for f (x, ǫ). In practice it is convenient to express the solution as a power series around ǫ = 0. Denoting with f (i) (x) the coefficient of ǫ i , we have, 8) and the different orders of the solution are related by the following recursive relation, The previous relation shows that the solution is expressed to all orders of ǫ in terms of Chen iterated integrals [56]. The solution is a pure function of uniform weight corresponding to the order of the ǫ expansion. The specific choice of the integral basis leading to the canonical form was achieved using the ideas outlined in [34]. In particular, it is expected that integrals with constant leading singularities [36] satisfy canonical differential equations. Using generalized cuts we look for combinations of integrals with simple leading singularities, that can be normalized to unity rescaling the candidate integrals. This typically leads to a form close to the canonical form. The remaining unwanted terms can be then algorithmically removed from the differential equations shifting the integral basis [42,44,57].

Polylogarithmic representation for algebraic alphabets
The alphabet (see appendix C for the explicit alphabet of the integral families) of the canonical integrals discussed in this paper contains 8 independent square roots that cannot be simultaneously rationalized via a variable change. This means that it is not possible to integrate (3.9) directly in terms of multiple polylogarithms [29].
However we can find an expression in terms of these functions by making a suitable ansatz in terms of polylogarithms of a given weight. The main task is to find suitable function arguments, as we discuss presently. This strategy is streamlined using the concept of symbol [29,31,75] of an iterated integral. The symbol corresponds to the integration kernels defining the iterated integrals. Since the integral basis is chosen to be of uniform weight, the symbol of the solution is completely manifest in our differential equations approach. Denoting by f The corresponding polylogarithmic functions can be found proceeding in the following algorithmic steps (see also [31,32]). First, one generates a list of function arguments as monomials in the letters appearing in the alphabet. For the classical polylogarithms Li n (x), one requires that 1 − x factorizes over the letters appearing in the alphabet 3 (a caveat is that in principle spurious letters might be needed [32]). For Li 2,2 (x, y), the condition is that 1 − x, 1 − y, 1 − xy factorize over the alphabet. Similar factorization properties are required for higher weight functions. Second, for each weight i, one chooses a maximal set of linearly independent functions for the alphabet. The linear independence can be verified using the symbol. By construction, we can solve the differential equations at every order in terms of this set of functions. Third, we determine the terms in the kernel of the symbol at weight i by writing the most general ansatz in terms of the lower weight functions, and solving the differential equations at O(ǫ i ) for the free coefficients of the ansatz. Finally, we recover transcendental additive constants imposing boundary conditions. Note that no assumptions were made on the rationality of the alphabet letters, so that the above steps generalize the algorithm of [32] to algebraic cases. Note also that, as opposed to a purely symbol-based approach, using the knowledge of the differential equations and of the boundary conditions, the solution is fully determined. In practice the alphabet under consideration is quite large, and a reasonably fast computer implementation of the algorithm above up to weight four is challenging. We can nevertheless use the algorithm to reconstruct polylogarithmic functions up to weight two, 3 When square roots are present it might be difficult to directly check factorization over the alphabet. In practice we can proceed as follows. We consider the logarithm of the function whose factorization we want to check, and we equate it to a generic linear combination of the logarithms of the alphabet (ansatz). Since additive constants are irrelevant at the symbol level, we derive the identity with respect to each variable. We then specialize the resulting system of equations for the (rational) free coefficients of the ansatz to many numeric values of the variables. If a solution exists the argument factorizes as desired over the alphabet and the solution defines the factorized form.
for which the alphabet letters contributing to the result are a relatively small subset of the full alphabet. The full set of linearly independent dilogarithms for the four families is listed in appendix D.
Having a representation of the weight two functions in terms of classical polylogarithms at hand is in fact very useful. As was shown in ref. [57], this can be used to write down useful one-dimensional integral representations for the remaining weight-three and weightfour functions.
Following [57], we use the Chen integral representation of the solution to write down a one-fold integral representation at weight three and four. Parametrizing the integration path C with α ∈ [0, 1], (3.9) translates to an iterated integral, In this language when the weight-two functions are known analytically, the weight-three functions are one-fold integrals. Initially, the weight-four functions are two-fold iterated integrals of differentials of logarithms, and they can be converted to one-fold integrals integrating by parts (see appendix E for a detailed discussion). The boundary conditions required to fix the solution are determined using the regularity of the pre-canonical integrals and the behavior of the algebraic factors defining the canonical basis in the boundary point. We find it convenient to use the boundary point The values of our integrals at this point correspond to the large heavy-quark limit so that one can apply the corresponding well-known graph theoretical prescriptions [76][77][78]. In the limit all the canonical integrals vanish except those that factor into products of one-loop integrals of which one is massless and thus diverges in the limit. These are however known analytically to all orders [24,79]. With this choice of the boundary point we can parametrize the integration path as, We have validated the analytic expressions performing numerical checks against the computer program FIESTA [80][81][82] for randomly selected points in the Euclidean region (2.6).

Elliptic integral sectors
The last two integral sectors of Family A (see appendix A), integrals f A 66 − f A 73 , turn out to be expressed in terms of elliptic functions. Using the language of the differential equations, the homogeneous part for sector I A 1,1,0,1,1,1,1,0,0 is not cast in canonical form, as the solution is expressed in terms of complete elliptic integrals. In appendix F we show that these properties can be verified a priori analyzing the maximal cut of the integrals. In section 4.1 we show that we can reduce the problem to the solution of a second order differential equation. In section 4.2 we show that using a proper unidimensional parametrization of the integrals the relevant second order differential equation can be solved with elementary (k 2 +p 1 ) 2 Figure 4: The four master integrals of the elliptic sector I A 1,1,0,1,1,1,1,0,0 . techniques. In section 4.3 we show that employing two auxiliary bases we obtain a two-fold iterated integral representation of the integral sector. The highest sector of Family A is I A 1,1,1,1,1,1,1,0,0 . In this case the homogeneous part of the differential equations can be cast in canonical form, however they depend via inhomogeneous terms on the lower elliptic sector. In section 4.4 we write the result as a three-fold integral. We found these integral representations to be suitable for precise and reliable numerical evaluations. When implemented in Mathematica the evaluation of both the elliptic sectors in one Euclidean point takes about 10 minutes using one processor, with about eight-digit accuracy.

Sector I
The integral sector I A 1,1,0,1,1,1,1,0,0 has four master integrals, shown in fig. 4, which are expressed in terms of elliptic functions, although its subtopologies do not involve them. We start by considering the following basis of finite integrals, We parametrize the integrals through the linear parametrization (3.12), and we define the differential equations with respect to the new parameter using the chain rule, where h is a vector, whose components are given in eq. (4.1). The differential equations have the following form, where g(α, ǫ) is the vector of the subtopologies, C (0) (α) and C (1) (α) are 4 × 4 matrices and D (1) (α) is a 4 × 65 matrix. In particular, the matrix C (0) (α) has the form, The last two integrals are decoupled from each other, but this is not required for the applicability of the method described here. It is manifest that the equations for the first two integrals are coupled.
We look for a solution in power series around ǫ = 0, The coefficients of the power series satisfy the following first order differential equations, where h (i) (α) is the unknown and the other terms define the inhomogeneous part. A two-bytwo system of first order differential equations for the first two components of h(α) defines a second order differential equation for the first component, where p 1 (α) and q 1 (α) depend on the matrix elements of C (0) (α), and are the same for every i, while r (i) 1 (α) is a function of the inhomogeneous part of (4.6). Once two homogeneous solutions of (4.7), y 1 (α) and y 2 (α), have been found, a particular solution can be determined using the method of the variation of constants. In general we get, where the arbitrary constants c i are fixed by the boundary conditions, and where w(α) is the Wronskian of the homogeneous solutions, Once h     1 (α). In order to solve the last two integrals we solve the respective first order differential equations, which depend on h (i) 1 (α) and h (i) 2 (α) via the inhomogeneous terms. This shows that, when computed in this way, h In order to optimize the numerical evaluation it is important to get rid of the extra integration. Furthermore, since at O(ǫ 4 ) these integrals would be expressed in terms of five iterated integrations, one integration must be spurious. In the non-elliptic case one is able to remove extra integrations using integration by parts. However in the elliptic case in order to perform an integration by parts one needs to integrate over elliptic functions, which is in general not possible analytically. We show how this is done in section 4.3.

Solution of the second order differential equation
The possibility of solving algorithmically a second order differential equation is related to the number of its singular points, including the point at infinity. If there are up to three singular points the equation can be cast in the form of the hypergeometric equation and two linearly independent solutions can be expressed in terms of hypergeometric functions [83]. Similar algorithms exist when four singular points are present. On the other hand if more than four singular points are present the solution requires a case by case analysis.
After differentiating with respect to the Mandelstam variables, the second order differential equation for I A 1,1,0,1,1,1,1,0,0 has six singular points. We show that using the parametrization (3.12) the solution can be reduced to the three singular point case.
Once h 1 (x(α), ǫ) is made explicit as in (4.1), the coefficients of the second order differential equation (4.7) are, and, where, We see that after using parametrization (3.12) we are left with three singular points, which are the two roots of d 1 (α) = 0 and the point at infinity. The homogeneous solutions of (4.7) can be then readily found 4 to be where the function k(z) is, and K(z) is the complete elliptic integral of the first kind 5 , The complete elliptic integral of the second kind is defined as, We have the following relations for the derivatives of the complete elliptic integrals, and, 1 (α) and its first derivative, it is expressed in terms of complete elliptic integrals of the first and second kind, of the same arguments as in (4.13). The Wronskian of the two homogeneous solutions is defined in terms of the derivatives above. Its expression is a rational function of the integration variable α, and in our case it reads, This property can be proven by using the Legendre identity, (4.20) Thanks to the overall normalization factor we chose for h 1 (x, ǫ), it is elementary to determine boundary conditions and use them to fix the free constants of the general solution (4.8). Integral I A 1,1,0,1,1,1,1,0,0 is regular for α = 0, so that h 1 (0, ǫ) = ∂ α h 1 (0, ǫ) = 0 and c 1 = c 2 = 0.

Auxiliary bases and solution in terms of two-fold iterated integrals
Since we need to evaluate the components of h (4.1) through O(ǫ 4 ), all the I integrals of eq. (4.1) need to be computed through O(ǫ 0 ), except I A 1,1,0,1,1,1,2,0,0 which must be evaluated through O(ǫ). Higher orders are irrelevant for two-loop processes. In general, the result for a master integral at O(ǫ i ) is obtained integrating over subtopologies through O(ǫ i−1 ) and, if coupled to them, over integrals of the same topology at O(ǫ i ). In section 3.2 we saw that weight-two functions can be expressed in terms of logarithms and dilogarithms, and weight-three functions can be reduced to one-fold integrals. This implies that, because of the general form of (4.3), h 3 (α) is expressed in terms of one-fold integrals, while h (4) and h (4) 2 (α) are expressed in terms of up to two-fold integrals. On the other hand h (4) 3 (α) and h (4) 4 (α) involve three-fold iterated integrals. In order to avoid considering more than two iterated integrations we introduce two auxiliary bases. We look for bases with differential equations of the form of (4.3), where the first integral is h 1 (α, ǫ) and the second integral is linearly independent of h 1 (α, ǫ) and h 2 (α, ǫ), the auxiliary integrals being independent of each other. In this way we compute the two remaining integrals as linear combinations of h 1 (α, ǫ) and its first derivative, generating at most two-fold iterated integrals. In practice we found two auxiliary bases equal to basis (4.1) modulo replacing, in turn, h 2 (α, ǫ) with ǫ 4 I A 1,2,0,1,1,1,1,0,0 and ǫ 4 I A 1,1,0,1,1,1,2,−1,0 . The full (finite) basis for the integral sector is then chosen to be, Interestingly, if we consider the differential equations for f A 66 −f A 69 , they are fully coupled and cannot be solved directly. We could nevertheless solve them with the help of auxiliary bases. A   1,1,1,1,1,1,1,0,0 The highest elliptic sector is I A 1,1,1,1,1,1,1,0,0 . It has four master integrals, shown in fig. 5, and it depends on the elliptic subsector I A 1,1,0,1,1,1,1,0,0 via inhomogeneous terms in the differential equations. Using the criteria outlined in [34] we can find a basis satisfying,

Sector I
(4.22) v(α, ǫ) is a four-dimensional basis vector for the highest elliptic sector, g(α, ǫ) is the vector of the subtopologies, F (1) (α) is a 4 × 4 matrix, G (0) (α) and G (1) (α) are 4 × 69 matrices. The homogeneous part is in canonical form, while this is not the case for the subtopologies. When solving the above equation for a given power of ǫ, we have to integrate over subsectors of the same order due to the G (0) (α) matrix. For numerical optimization it is convenient to get rid of such integrals. Matrix elements of G (0) (α) corresponding to non-elliptic subsectors are removed with a basis shift, as described in [44,57]. In order to remove G (0) (α) entries corresponding to elliptic subsectors we proceed as follows. Let us consider the i th component of v, which fulfills the equation, where k ij (α), with i = 1, . . . , 4 and j = 1, 2, are known algebraic functions and e 1 , e 2 are two coupled integrals of an elliptic subsector, satisfying, where b ij (α) are functions to be determined. After the basis shift the equation for v i reads, In order to remove terms proportional to e 1 (α, ǫ) and e 2 (α, ǫ), their coefficients must vanish, i.e. b ij (α) must fulfill the equations, with j = 1, 2. For fixed i, the above equation is a two-by-two system of first order differential equations. The matrix defining the system is the transpose of the matrix defining (4.24). This implies that if y 1 (α) and y 2 (α) are the homogeneous solutions of (4.24) and w(α) is their Wronskian, the solutions of (4.27) are, where c is an overall constant. Their Wronskian is c 2 /w(α). Therefore with the method of the variation of constants the full expression for b i1 (α) reads, where L i (α) are functions of k i1 (α) and k i2 (α), and where two arbitrary integration constants have been set to zero. In addition, we set the lower integration bound to 1 but we have the freedom to choose a different value. Usually this is dictated by the properties of the integrand, that might have non-integrable singularities for specific integration bounds. Once b i1 (α) is known it is elementary to obtain b i2 (α) using the same differential equations. For sector I A 1,1,1,1,1,1,1,0,0 the integrals that need to be shifted are f A 71 and f A 73 , as e 1 and e 2 defined via (4.25) are equal to f A 66 and f A 67 respectively. y 1 (α) and y 2 (α) are the same as those of (4.13) and, while L 1 and L 3 vanish.
In general the integrals of (4.29) are not known analytically in closed form. Since after the basis shift they will contribute to the matrix elements of the differential equations, one might wonder if such a basis change is convenient in practice, as our main goal was to get rid of one integration. In practice, because of the simple form of (4.30), its numerical evaluation takes O(10 −3 ) sec. In this form the result for the elliptic sector at O(ǫ 4 ) is in terms of three-fold integrals, while their numerical performance is comparable to the one of two-fold integrals. Alternatively, it is possible to series expand the complete elliptic integrals of eq. (4.29) and then perform the integrations analytically 6 . In this way the result for the integral sector can be expressed in terms of two-fold integrals 7 .

The class of functions
In order to discuss the general structure of the solution of sector I A 1,1,0,1,1,1,1,0,0 let us introduce the following shorthands for the complete elliptic integrals defined in section 4.2, are expressed as linear combinations of the class of functions, where E (σ) can be one of the following complete elliptic integrals, where σ ∈ {−1, 1}. F(t) denotes a linear combination of pure weight-two and weight-three functions, belonging to the subtopologies, multiplied by either derivatives of logarithms or derivatives of algebraic functions, with respect to α 8 . Interestingly, weight-three functions are never multiplied by derivatives of logarithms, but only by the following simple inverse square roots (modulo functions depending only on rescaled Mandelstam invariants), The same class of functions has been found in [84] for the massive crossed triangle. See [62,63,67,85] for results in terms of elliptic polylogarithms [86], and [64-66, 69] for a related class of functions.
In order to decouple integral sector I A 1,1,1,1,1,1,1,0,0 from sector I A 1,1,0,1,1,1,1,0,0 , in section 4.4 we performed a non-algebraic basis shift of f A 71 and f A 73 , involving integrals of complete elliptic integrals, that we denote here with the following shorthands, G(t) has the same properties as F(t) described above, but the prefactors of pure weight-three functions are any of the algebraic functions,

Conclusion and perspectives
In this paper we presented the analytic computation of all the planar master integrals which are necessary to evaluate the two-loop amplitudes for Higgs → 3 partons, with the full heavy-quark mass dependence. They occur in the NNLO corrections to fully inclusive Higgs production and in the NLO corrections to Higgs plus one jet production in hadron collisions. The result is expressed in terms of iterated integrals over both algebraic and elliptic kernels. This is the first time that Feynman integrals for four-point multiscale amplitudes involving elliptic functions are computed in a fully analytic way. While it was generally believed that the analytic computation of multiscale loop integrals with many internal massive lines was out of reach with present analytic tools, this work shows that new ideas involving the proper parametrization of the integrals, an optimal basis choice, and the subsequent solution with the differential equations method in terms of elliptic iterated integrals, are effective to treat such problems. The computation of the non-elliptic integral sectors has been performed with the differential equations method applied to a set of basis integrals defined to be pure functions of uniform weight. The presence of many square roots that cannot be simultaneously rationalized makes the direct solution of these equations in terms of multiple polylogarithms not possible. We have shown that the Chen iterated integral representation plus the knowledge of the boundary conditions provide the information needed to integrate the system in terms of a minimal polylogarithmic basis, circumventing in this way the necessity to rationalize the square roots of the alphabet. To do so we used an algorithm for the integration of symbols with general algebraic alphabets, generalizing well established algorithms for the rational case.
We have seen that the crucial point for the computation of the elliptic sectors is the solution of the associated homogeneous second order differential equation. We noticed that a very simple univariate reparametrization of the integrals makes the equation elementary and standard tools are sufficient to solve it. The central point is that the fewer singular points are present in higher-order differential equations, the simpler is their solution. It will be important to further investigate and develop the idea of what is the proper parametrization of the integrals yielding the simplest singular structure of the equations. The univariate parametrization has also the benefit that only one set of differential equations has to be solved, while in the traditional approach one has to iteratively solve multiple sets of equations, one for each variable, which might be highly non-trivial when elliptic functions are involved.
In contrast to the non-elliptic sectors, we did not use the notion of canonical basis for the elliptic sectors. Instead, we showed that the problem can be completely solved in total generality, once the relevant higher order homogeneous equations have been solved. However it will be important to extend the notion of canonical basis to elliptic cases. First, this will clarify the class of functions needed to represent the answer -in our case we used a rather general class that might still contain spurious information. Second, it is natural to expect that the explicit results for canonical integrals will be relatively compact. In order to define a canonical basis in the elliptic case, the notion of leading singularity has to be generalized, which is beyond the scope of the present paper (see appendix F for a discussion about the maximal cut of those integrals, which would be the starting point for defining a generalization of leading singularity in the elliptic case). In particular, we know [37][38][39] that it is possible to obtain a form of the differential equations with only Fuchsian singularities and linear in ǫ. This is valid for any Feynman integral and it is another natural starting point for finding a canonical basis.
We showed that for the sake of stable and precise numerical evaluations we can express elliptic iterated integrals through O(ǫ 4 ) in terms of one and two-fold iterated integrals for the non-elliptic and elliptic sectors, respectively. We found these representations suitable for numerical evaluation. In principle, as the integrands are known functions, it should be possible to achieve a series representation of the solution, though we did not attempt it as the integral representation already showed satisfying performance. It will be important to develop general purpose numerical routines for elliptic iterated integrals, so that one can take advantage of such analytic expressions also when higher loop orders are considered, i.e. when more iterated integrals are needed.
PITN-GA-2012-316704. JMH is supported in part by a GFK fellowship and by the PRISMA cluster of excellence at Mainz university.

A Integral basis
In this appendix, we provide the explicit form of the integral families we used to parametrize the integrals defined in eq. (2.1). We call them: family A, B, C, and D.
For each family, we perform an independent reduction to the master integrals. Then we perform a change of basis that maps the master integrals into the canonical form. We For each family of integrals we define the corresponding system of differential equations, that we then solve as discussed in sections 3 and 4.
Note that, in general, there is an overlap among the master integrals of the different families. Making the appropriate correspondences, we can reduce the process to the computation of 125 master integrals. In the next appendix, we draw these 125 (pre-canonical) master integrals and we link them to the corresponding canonical form.
We label with p 1 , p 2 , and p 3 the momenta of the massless partons, and with p 4 = p 1 + p 2 + p 3 the momentum of the Higgs. The loop momenta are labeled with k 1 and k 2 . Finally, we use the shorthand p ij = p i + p j .
Family A is defined by the nine propagators, with the extra restriction that a 8 and a 9 are non-positive. The family contains 73 master integrals. Below, we give the basis transformation between pre-canonical and canonical forms.

C Alphabet
In this appendix we list the alphabet for the four integral families defined in section 2. We introduce the following shorthands for the set of 13 square roots, They appear in the alphabet in the following 8 linearly independent combinations, Referring to the matrixÃ defined in (3.4) the alphabets of the four families can be written in terms of the following linearly independent 49 letters,

D Weight-two functions
In section 3.2 we described how to express the non-elliptic master integrals in terms of a minimal set of logarithms and dilogarithms up to weight two, while the weight three and four components are expressed as one-fold integrals over linear combinations of weight-one and weight-two functions with algebraic coefficients. In this appendix we list the basis choice we made for the set of linearly independent dilogarithms required to express the master integrals of each family at weight two. They are chosen to be single-valued in the Euclidean region x 3 < x 2 < x 1 < 0.

E One-fold integral representations
We consider a system of differential equations for a set of integrals f (x, ǫ) in canonical form [34] defined by a matrixÃ(x), . If the boundary value f (0, ǫ) is singular, the integral above cannot be computed numerically since the integrand has non-integrable singularities in α = 0. However in our case all the divergent integrals are factorisable into products of one-loop integrals, which are known analytically [24,79]. In such cases we need to define the integrals that, viaÃ(x), depend on the singular ones.
Assume that integral f k (x, ǫ) has a singular boundary condition f k (0, ǫ), and that it is known analytically to all orders of ǫ. Consider an integral f n (x, ǫ), with n = k, with a regular boundary condition f n (0, ǫ). Using eq. (E.2) we can write it as, where we made explicit the dependence on the singular integral f (i) k (x). Since by assumption f (i) k (x) is known analytically, we can directly compute the second integral on the right hand side. Note that the fact that f n (0, ǫ) is regular ensures that the second integral is convergent even if f (i) k (0) is singular. Finally we can perform an integration by parts and reduce the other integrals to the form of (E.3).

F Maximal cut of the elliptic sectors
We show that the maximal cut [35,89] of I A 1,1,0,1,1,1,1,0,0 provides useful information about the class of functions needed to represent the result. We cut the six visible propagators. We parametrize the two loop momenta using the spinor-helicity formalism [90] (see [91,92] for a different formalism), (F.1) We get the following two-fold integral result for the maximal cut, The integrand is the square root of a quartic polynomial in z 8 , with four different roots. This means that the integrand has two genuine branch cuts that cannot be removed by any change of variables, yielding an elliptic function upon integration [59,93].
(F. 5) We note that in the limit m 2 → 0 the Jacobian reduces to, J(k 2 )| m 2 →0 = s 3 12 (k 2 + p 1 ) 2 , (F. 6) reproducing the well known result for the cut of the massless case. Localizing the contour onto the two genuine propagators of (F.4), will yield an expression similar to (F.2) -an inverse square root of a quartic polynomial with no repeated roots.