The complete set of two-loop master integrals for Higgs + jet production in QCD

In this paper we complete the computation of the two-loop non-planar master integrals relevant for Higgs plus one jet production initiated in arXiv:1907.13156. We compute the integrals by defining differential equations along contours in the kinematic space, and by solving them in terms of one-dimensional generalized power series. This method allows for the efficient evaluation of the integrals in all kinematic regions, with high numerical precision. We show the generality of our approach by considering both the top- and the bottom-quark contributions. This work along with arXiv:1907.13156, arXiv:1609.06685, arXiv:1907.13234 provides the full set of master integrals relevant for the NLO corrections to Higgs plus one jet production, and for the real-virtual contributions to the NNLO corrections to inclusive Higgs production in QCD in the full theory.


Introduction
The main production mode of the Higgs boson at the Large Hadron Collider (LHC) is via gluon fusion. In perturbative Quantum Chromo Dynamics (QCD) the production is mediated by a quark loop that couples to the final-state Higgs. The quark-Higgs coupling is proportional to the quark mass, hence the largest contribution is given by corrections involving a top-quark. Being mediated by a quark loop, the leading-order (LO) corrections require the computation of one-loop amplitudes while the next-to-leading-order (NLO) corrections require the computation of two-loop amplitudes and so on. The inclusive LO corrections to Higgs production have been computed in the full theory at LO in [4] and at NLO in [5,6]. On the other hand, the computation of the higher order corrections is much more challenging, and results in the full theory are not available yet. The computation can be considerably simplified in the limit where the top quark is assumed to be infinitely heavy, while the other quarks are assumed to be massless. This limit is known as the Higgs Effective Field Theory (HEFT). The next-to-next-to-leading-order (NNLO) QCD corrections have been computed in the HEFT in [7][8][9] while, more recently, the corresponding next-to-next-to-next-to-leading-order (N 3 LO) corrections have been computed in [10,11].
In addition to inclusive cross-sections, differential cross sections play an important role in the study of the properties of the Higgs boson. In particular, the Higgs may couple to particles not predicted by the Standard Model, and many of these effects are best studied by observing the high transverse momentum (p T ) distribution of the Higgs [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. In the full theory the Higgs plus jet production cross section and the p T distribution are only known at LO. At NLO, the top-quark contributions have been computed in [27], while the top-bottom interference was computed in [28] by combining the HEFT with an asymptotic expansion around small bottom-mass. At higher perturbative order no result is available in the full theory, and only partial results are known in the HEFT. More specifically, the NNLO corrections to the Higgs plus one jet production and the Higgs p T distribution are known in the HEFT. However, while the HEFT approximation works well for inclusive observables, it diverges very rapidly for high-energy differential observables, such as the high p T distribution of the Higgs (see e.g. [29] and references therein).
To this date no complete result for the Higgs plus jet amplitudes at NLO is available in the full theory. The first step in this direction has been taken in [2] and, more recently, in [1,3], where the planar master integrals and one of the two non-planar families of master integrals at two loops have been computed in terms of one-dimensional generalized power series. This technique is not constrained in any way to a particular kinematic region or a specific configuration of the relevant masses, and allows to efficiently compute the master integrals by assuming the full dependence on all the mass scales. In this paper, we apply this technique to compute the remaining family of non-planar master integrals. Besides the NLO QCD corrections to Higgs plus jet production, these master integrals are a fundamental building block for the NNLO inclusive corrections to Higgs production in the full theory, where the Higgs plus jet amplitudes appear as the single real radiation contribution.
The paper is organised as follows. In section 2 we define the non-planar integral family computed in this paper. In section 3 we review the differential equations method for dimensionally regulated Feynman integrals, and we discuss the structure of the differential equations of our integral family. In section 4 we describe our solution strategy, i.e. we solve differential equations along contours in the space of kinematic invariants in terms of generalized power series. In section 5 we show how the expansion strategy is used to evaluate the master integrals in a very large set sample of points of the physical region, for both the top-and bottom-quark contributions. In section 6 we draw our conclusions and we discuss directions for future work.

The integral family
As discussed in ref. [1], six seven-propagator integral families contribute to the two-loop QCD contribution to H + jet production. Of these families, four are planar and have been computed in ref. [2], and of the remaining, non-planar, one was computed in ref. [1] and the remaining one, denoted family G, will be the topic of the present paper.
That integral family is defined by I a 1 a 2 a 3 a 4 a 5 a 6 a 7 ;a 8 a 9 = 9 P a 1 1 P a 2 2 P a 3 3 P a 4 4 P a 5 5 P a 6 6 P a 7

7
(2.1) with Only P 1 -P 7 can appear as genuine propagators, so we have a 8 and a 9 restricted to the non-positive integers. The kinematics is p 2 1 = p 2 2 = p 2 3 = 0 and additionally where m 2 denotes the squared mass of the quark that couples to the Higgs, and p 2 4 the squared mass of the Higgs. By using integration-by-parts (IBP) [30][31][32][33] reduction methods [34,35], we identify a set of 84 master integrals for this family, whose diagrams are shown in Fig. 2. With those master integrals we defined a basis of Feynman integrals which is presented in Appendix A.

Differential equations for the integral family
Given a basis of N master integrals I( , s), where is the dimensional regulator defined by D = 4 − 2 and s = {s 1 , . . . , s n } is a set of n Lorentz invariants, it is possible to define a closed system of linear, first order differential equations [36][37][38][39][40] for I( , s i ) that in full generality reads, where ∂ s i ≡ ∂ ∂s i and M s i is a set of N × N matrices. The choice of the basis integrals is not unique, and by performing a basis change B = T I the differential equations transform according to,  Figure 2: The 84 master integrals. Shown on the figure is the sector, i.e. the set of propagators, to which the master integrals belong. Higher powers of propagators, numerators, or prefactors are not shown. External momenta are labelled using p ij = p i +p j and p 4 = p 1 +p 2 +p 3 . Masses (internal as well as external) are indicated with a thicker line.
In Ref. [41] it was conjectured that with a proper basis choice it is possible to cast the differential equations for Feynman integrals in the following simplified form, where the dependence on is factorised and the matrices A s i ( s) depend only on the invariants s. Such a system of differential equations is said to be in canonical form, and the basis B is referred to as the canonical basis. A canonical system of differential equations is equivalent to the following equation in differential form, d B( , s) = dÃ( s) B( , s) , (3.4) where, by construction, the matrixÃ satisfies The differential equation (3.4) can be formally solved in terms of a path-ordered exponential B( , s) = P exp γ dÃ( s) B( , s 0 ) , (3.6) where γ is a contour connecting a boundary point s 0 to s, and P is the path-ordering operator. In dimensional regularization we are generally interested in a solution around = 0. By performing the expansion for small , the path-ordered exponential translates to iterated integrals over the entries ofÃ( s). Specifically, the solution to all orders of is So far we made no assumptions on the class of functions arising from the iterated integrals of Eq. (3.7). A large class of master integrals admits a canonical basis such that the matrixÃ( s) is a Q-linear combination of logarithms of rational or algebraic arguments. This form also implies that the transformation T( , s) to the canonical basis is rational or algebraic. In the rational case, the solution can be directly expressed in terms of multiple polylogarithms [42], defined recursively as, G(a 2 , . . . , a n , t), (3.8) with G(, x) ≡ 1 and G( 0 n , x) ≡ log(x) n n! . On the other hand, the general algebraic case is much less understood. In some cases it is possible to rationalize the algebraic functions by a suitable reparametrization of the invariants, reducing in this way the problem to a rational one. If a rational parametrization is not available, it is possible, in some case, to define an ansatz for the solution in terms of polylogarithms of suitably chosen (algebraic) arguments. The unknown parameters of the ansatz are then fixed by solving the differential equations and by imposing boundary conditions (see e.g. [1,2,[43][44][45][46][47]). Nonetheless, in the general algebraic case, it is not known whether the differential equations always admit a solution in terms of multiple polylogarithms.
More recently, a lot of progress has been made in the study of Feynman integrals which evaluate to elliptic generalizations of multiple polylogarithms (eMPLs) [2,11,. However, while in some cases it is possible to define a basis that casts the differential equations in canonical form (see e.g. [75,76,80,81]), little is known on their general analytic properties and how to systematically solve them in terms of elliptic multiple polylogarithms.
As we will see in the next sections, the differential equations for the integral family considered in this paper depend on complicated algebraic functions. Moreover some equations are coupled, and their solution involves functions of elliptic type. In this case finding a closed form solution for the integrals seems to be out of reach with current technology. Nonetheless, having phenomenological applications in mind, we follow a different approach, based on the series solution of the differential equations along contours in the space of the kinematic variables [3].

Canonical integrals
We denote the set of Lorentz invariants as s = {s 1 , s 2 , s 3 , s 4 } = {s, t, p 2 4 , m 2 }. The first 71 master integrals of the basis chosen for family G are such that the system of differential equations are in canonical form. Namely, they satisfy Eq. (3.4). The matrixÃ can be constructed by using the following iterative definitions, and takingÃ For this integral family, the matrixÃ( s) is a Q-linear combination of 76 logarithms (the so-called letters) depending on 11 different square roots. The full set of letters and square roots is presented in Appendix B. It will be interesting to study the problem of finding a polylogarithmic solution to these equations by using the methods recently put forward in [47]. The first two of these evaluate to elliptic integrals of the first kind, while the latter evaluates to a combination of logarithms. This corresponds to the two elliptic curves,

Elliptic integrals
being present in the results. The integrals B 72 −B 75 are planar, and indeed that sector is equivalent to the sector of the integral A 66 discussed in refs [2,3]. Likewise the integrals B 76 −B 79 are merely a crossing thereof with p 1 ↔ p 2 , corresponding to t ↔ u. The fact that eq. (3.13) does not evaluate to elliptic integrals does not mean that such structures are absent in the un-cut integrals, as elliptic curves would appear at the sub-maximal cuts corresponding to the sectors B 72 −B 75 and B 76 −B 79 . The appearance of functions of elliptic type is also manifest by analyzing the relevant system of differential equations. More specifically, integrals B 72 −B 84 satisfy, where the vector G 72−84 ( s, ) depends on the canonical integrals B 1−71 ( s, ), and the homogeneous matrix has the schematic form, where the lines separate the 3 elliptic sectors. We see that integrals 72,75 and 76,79 are coupled. The homogeneous solution of integrals 72 and 76 are linear combinations of the elliptic integrals (3.11) and (3.12) respectively [60], evaluated along all possible pairs of branching points of the respective elliptic curve [63]. By using the variation of the constants method, it is possible to solve integrals 72 and 76 in terms of iterated integrals over (products of) the elliptic integrals (3.11) and (3.12). The remaining integrals 73-75, 77-84 are coupled to 72 and 76 via inhomogeneous terms of the differential equations, and they can be solved in terms of iterated integrals over (products of) the same elliptic integrals (3.11) and (3.12). We remark that the presence of multiple elliptic curves renders the functional form of the solution an open problem (but see e.g. [75] for progress in this direction).

Series expansion along contours
We consider the series expansion strategy outlined in [1,3] (see also [11,[83][84][85][86][87][88] for the application of expansion methods to single scale processes, and [89][90][91][92][93][94] for expansion methods applied to multiscale problems in particular kinematic limits). The strategy relies on parametrizing the integrals along straight line segments, for which we solve the corresponding differential equations in terms of one-dimensional generalized series. We briefly review the strategy here, and highlight aspects that are specific to the integral family under consideration. We start from the system of differential equations of the basis defined in Appendix A, which has the form, where M = s i M s i ( , s)ds i , and where we otherwise suppress variable dependence in the notation. For convenience, we put m 2 = 1 without loss of generality. We consider a generic line parametrized as, The differential equations along this line take the form, ∂λ . Next, we expand the differential equations in to obtain a system for each order in . In particular, we let, Note that both expansions start at finite order for our choice of basis. The system of differential equations now takes the following form, order-by-order in , where we separated out the homogeneous part M Let us consider first the polylogarithmic sectors and review the series solution strategy for those. The system of differential equations becomes simply, where M (1) λ = ∂Ã/∂λ , andÃ was defined in Eq. (3.4). Hence, the work for the polylogarithmic sectors amounts to solving a sequence of first order differential equations without homogeneous part. The general solution is easily found from a single integration, Importantly, we solve the integration by performing series expansions in λ. The only algebraic terms in our basis and in the matrices are square roots, so that the expansions of the matrix elements are in terms of integer and half-integer powers of λ. After integrating, the series expansions will also contain logarithmic terms. Therefore, in general each integration in Eq. (4.7) is of the form, where q is an integer or half-integer, and p is a non-negative integer. It may easily be verified that such integrals evaluate to sums of terms of the same type, by using integration-by-parts identities to reduce the power of the logarithm inside the integral. As shown in Eq. (3.16), we simplified our basis in such a way that in each elliptic sector at most 2 integrals are coupled together, in particular the pairs of integrals 72, 75 and 76, 79. These integrals can be solved by combining their first order differential equations into second order differential equations. Integrals 74 and 78 can be solved from their first order differential equation, but, in contrast to the polylogarithmic sectors, their differential equations have a homogeneous component. Lastly, integrals 73, 77, 80, 81, and 82 satisfy a first order differential equation without a homogeneous component, and can be solved in the same manner as the polylogarithmic integrals. For completeness, we discuss solving these two cases next.
Consider a first order differential equation with homogeneous component of the form, The solution to the homogeneous part is easily found to be, up to an arbitrary multiplicative constant. The full solution to Eq. (4.9) is then given in terms of µ(λ) by, Now, consider a second order differential equation of the form, Given two solutions µ 1 (λ) and µ 2 (λ) to the homogeneous part of the differential equation, the general solution can be written using the method of variation of parameters as where d 1 and d 2 are complex constants to be fixed from boundary conditions. The remaining challenge is to find two distinct homogeneous solutions µ 1 (λ) and µ 2 (λ) that are not related by a rescaling. From the well-known Frobenius method (see e.g. [95] for an extensive review of the method), it follows that we may always find one series solution of the form µ 1 (λ) = λ r + λ r ∞ k=1 µ 1,k λ k . The values for r and µ 1,k may be found up to the desired order in λ by plugging µ 1 (λ) into the homogeneous differential equation as an ansatz, and solving order-by-order in λ for the unknowns. The lowest order in λ gives a quadratic equation in r called the indicial equation. By picking r to be the largest root of the indicial equation, it is guaranteed that we may solve for the remaining unknowns µ 1,k with k ≥ 1.
It remains to find a second homogeneous solution. This may be done in the following way. First we write the second homogeneous solution µ 2 (λ) as µ 2 (λ) = µ 1 (λ)h(λ). Plugging this expression into the homogeneous part of Eq. (4.12), we find a new equation, which we recognize as a first order homogeneous differential equation for h (λ), which we know how to solve. This way, we obtain the second homogeneous solution µ 2 (λ). Thus we may now use Eq. (4.13) to compute the full solution to Eq. (4.12).

Boundary conditions
To fix our system of differential equations we need a suitable boundary point. Similar to [1], we work in the heavy mass limit parametrized by, where λ is a line parameter that goes to zero. Using the method of asymptotic expansions in the parametric representation [96][97][98][99][100][101], we may obtain values of our basis integrals in the heavy mass limit. The final result turns out to be very simple, We note that the homogeneous solution of the differential equation for B 78 along γ hm (λ) is proportional to λ, and hence we are not able to determine the boundary constant for B 78 directly from Eq. (4.16). It may be verified that B 78 is also zero at order λ 1 in the heavy mass limit, and hence the constant multiplying the homogeneous solution may be put to zero for this integral.

Convergence of the series
A trait of the expansion strategy is that each expansion at a given point along a line has a limited range of convergence. Namely, each expansion at a given point is valid up to the distance of the point to the nearest singularity. Thus, to obtain results along a given line, numerous expansions along segments of the line have to be patched together in order to reach a given point in phase space. In particular, to cross a singularity we may perform an expansion at the singularity, and fix its boundary conditions from an expansion at a neighbouring point along the line. We employ the following strategy for deciding along which line segments to expand: • First we create a list A of all singularities of the matrix elements of M λ on the line γ(λ) along which we seek to integrate. By singular point we mean any non-analytic point of the differential equations. In our case, these are the zeros of the denominators of the matrix elements, and the zeros of the square roots.
• Some of the singularities may be complex. We replace each complex singularity λ sing = λ sing re + iλ sing im in the list A by three real points: λ sing re − λ sing im , λ sing re and λ sing re + λ sing im .
• Next, we consider a Möbius transformation λ = g(λ ) for each triplet (a, b, c) of neighbouring points in A, such that g −1 ({a, b, c}) = {−1, 0, 1}. Note that a series (in λ ) centered at λ = 0 will have a radius of convergence greater than or equal to 1.
• To obtain results along γ(λ) from λ 0 to λ 1 , we have to match expansions along neighbouring line segments, which are expressed in terms of Möbius transformed line parameters, say λ and λ . We may find a matching point between two neighbouring expansions by solving λ = −λ , assuming λ corresponds to the line segments lying on the right.
• In general, one may find that this condition picks λ and λ to be very close to 1 and -1, respectively, where both series may be very slowly converging. This can be solved by adding additional expansion points along the line segments. In particular, we may consider new expansion points between -1 and 1, such that upon matching neighbouring expansions, neither gets evaluated further than a certain fraction of the distance to the nearest singularity. We will refer to the inverse of this fraction by the parameter k. Hence, with k = 2, the expansion points are chosen such that no series is evaluated beyond half its radius of convergence. The situation is illustrated in Figure 3.
We note that in general we may encounter both spurious, physical, and non-physical singularities. The spurious singularities are singularities that only appear in the elements of M, but which are not singularities of the basis integrals themselves. The physical singularities are threshold singularities, in our case s = 4m 2 and p 2 4 = 4m 2 . For those, it is important to cross the singularity according to Feynman prescription, which tells us to interpret s and p 2 4 as having an infinitesimally small positive imaginary part. Furthermore, we should make sure to assign the same imaginary part to the square roots in our basis that are associated with physical singularities. Specifically, some of our basis integrals have the prefactors 4m 2 − p 2 4 and √ 4m 2 − s, which are analytically continued as 4m 2 − p 2 4 − iδ and √ 4m 2 − s − iδ for an infinitesimally small δ > 0. Lastly, there are also non-physical singularities, which can arise from rational prefactors in the basis, or from square roots in the basis that do not correspond to physical singularities. Since these singularities are introduced by the basis choice, we are free to assign every non-physical root in the basis the standard branch, i.e. we consider the argument to carry the imaginary part +iδ.
To improve the convergence of our series solutions, we compute their diagonal Padé approximants and evaluate those instead at each (matching) point. Since we are dealing with generalized series that may in general include powers of logarithms, we collect first on powers of logarithms and compute the Padé approximant for each series that multiplies a given power.

Results for top and bottom quarks
In this section we present explicit results that were obtained using the expansion method described in the previous section. Specifically, we used our method to compute the integrals in 10000 points covering the physical region, for both the top-and bottom-quark corrections, and we present plots thereof. The relevant physical region for Higgs plus jet production is given by, We may map the physical region to the unit square by using the parametrization, Since we chose to work with m 2 = 1, the value for p 2 4 is given by m 2 H /m 2 q where m H denotes the mass of the Higgs particle, and m q denotes the mass of the internal quark. For the top quark we approximate the ratio by p 2 4 = 13/25, while for the bottom quark we consider the ratio p 2 4 = 323761/361. For the case of the top quark, the particle production threshold s = 4m 2 corresponds to z = 13/100. For the sake of the presentation of the plots, we use a Möbius transformation to map z = 13/100 to 1/2, while keeping z = 0 and z = 1 fixed. Thus, we consider the following parametrizations of the physical regions of the top and bottom quark contributions, To produce plots in these regions we seek to compute n 2 evenly spaced points on the unit square for all basis integrals, and in particular we let n = 100, so that we obtain 10000 points in total. We explain next how we obtained results in these points. For convenience we use the notation a → b to denote a line, we denote coordinates in the physical regions by pairs (l, z), and we denote the heavy mass limit by 0. The following discussion applies to both the top and bottom region, given their respective set of (l, z)-coordinates.
We computed the boundary data for the horizontal lines of the top and bottom physical regions on a laptop using a single core. We then computed all horizontal lines on a cluster with 48 cores. The run for the horizontal lines of the top quark and the run for the horizontal lines of the bottom quark, both took a few hours to complete on the cluster.

Cross-checks of the expansions
We have performed several checks of our results. The first class of cross-checks was performed by evaluating multiple points by reaching them with different contours, finding full agreement. These cross-checks are summarised in Table 1.
For the top-quark integrals we compared our results against FIESTA [102] for multiple points of the physical region finding full agreement. For the physical region of the bottom quark we checked most of the integrals against FIESTA and SecDec [103] finding full agreement. However, for some of the integrals, these programs encounter numerical instabilities. We then perform different checks. Firstly we cross checked our results against FIESTA in the point (s = 53, t = −11, p 2 4 = 23, m 2 = 1) finding full agreement. This provides a direct check of the analytic continuation past the thresholds s = 4m 2 and p 2 4 = 4m 2 . In addition, we have performed a numerical cross-check against a private Line(s).
Evaluated at #Segments (k = 2) Max relative error 0 →  Table 1: Internal cross-checks of our results. The cross-check lines are compared to the data that we generated for the plots, which was computed in the manner illustrated in Fig. 4. We compare the results by computing the relative error B code [104,105], which exploits the loop tree duality [106][107][108][109] for the numerical evaluation of Feynman integrals. In particular, we compared integrals B 72

Conclusion
In this paper we computed a family of two-loop non-planar master integrals relevant for the QCD corrections to Higgs plus one jet production in the full theory. Our result, together with [1][2][3], provide the full set of master integrals required for the computation of the NLO corrections to Higgs plus one jet production, and the NLO corrections to the p T distribution of the Higgs. Moreover, our results provide the full set of master integrals relevant for the single-real radiation contributions to the NNLO corrections to inclusive Higgs production.
The computation was performed by using the differential equations method. More specifically, we defined an integral basis such that most of the integrals satisfy differential equations in canonical form. Three integral sectors are coupled, and their solution involve functions of elliptic type. Having phenomenological applications in mind, we solved the differential equations along contours in the space of kinematic invariants, in terms of onedimensional generalized power series. More specifically, given a boundary point where the value of the integrals is known, we defined the differential equations along a contour connecting a boundary point to a new point of the kinematic regions. In this way the problem was effectively reduced to a single-scale one, and finding the series solution was algorithmic. We showed that this method is efficient, and can be repeated in order to compute the integrals in any point of the kinematic regions. The analytic continuation of the series solution across the physical thresholds is straightforward, as it requires only the analytic continuation of logarithms and square roots.
In order to show the generality of our approach, we computed the master integrals for both the top-and bottom-quark mass. Moreover, we explicitly obtained results for a large set of points covering the relevant physical regions. The typical evaluation time is of the order of 1 second per integral, with a relative accuracy of order 10 −24 , on a single CPU core. If needed, the numerical precision can be made arbitrarily high by increasing the truncation order of the power series. These features render our methods well suited for Monte Carlo phase-space integrations.
We remark that the applicability of our methods does not rely on the number of physical scales, specific kinematic configurations, or a particular form of the differential equations. For this reasons, we believe that our approach will be relevant for the computation of several processes of phenomenological interest.

A Canonical basis and basis for elliptic sectors
In this section we provide the set of 84 basis integrals used in this paper, written in terms of the set of master integrals depicted in Figure 2 and defined as in Eq. (2.1).