Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops

We obtain generalised power series expansions for a family of planar two-loop master integrals relevant for the QCD corrections to Higgs + jet production, with physical heavy-quark mass dependence. This is achieved by defining differential equations along contours connecting two fixed points, and by solving them in terms of one-dimensional generalised power series. The procedure is efficient and can be repeated in order to reach any point of the kinematic regions. The analytic continuation of the series is straightforward and we present new results below and above the physical thresholds. The method we use allows to compute the integrals in all kinematic regions with high precision. Performing a series expansion on a typical contour above the physical threshold takes on average $\mathcal{O}(1 \text{ second})$ per integral with worst relative error of $\mathcal{O}(10^{-32})$, on a single CPU core. After the series is found the numerical evaluation of the integrals in any point of the contour is virtually instant. Our approach is general and can be applied to Feynman integrals provided that a set of differential equations is available.


Introduction
The computation of Feynman integrals is a central ingredient for the prediction of collider experiments. In the past decades, we have seen an enormous progress in our capabilities to efficiently compute Feynman integrals in closed analytic form or in a purely numerical way. From the analytic side, several techniques are available. Some of the most effective techniques are the differential equations method [1][2][3][4][5] and the direct integration methods [6,7]. In dimensional regularisation, one is able to reduce a given (generally large) set of scalar Feynman integrals to a minimal set of linearly independent integrals, called master integrals (MIs), by using integration-by-parts identities (IBP) [8][9][10][11]. Once a basis is identified, it is possible to define a closed system of first order linear differential equations, that can be solved in terms of iterated integrals [12]. Our understanding of differential equations for Feynman integrals has been further refined by the identification of canonical bases of integrals [13], which make the solution of the equations in terms of iterated integrals completely algorithmic. In some cases Feynman integrals can be computed in terms of special functions known as multiple polylogarithms (MPLs) [14] or, more recently, in terms of their elliptic generalisation, elliptic multiple polylogarithms (eMPLs) [15][16][17][18] (for analytic results involving functions of elliptic type see e.g. [17,). Even though having a representation in terms of known functions is important from the conceptual and practical side, in recent years this approach has become challenging. For state-of-the-art computations one usually encounters multi-loop integrals depending on several mass scales. In this case the differential equations exhibit complicated analytic structures, and their solution in terms of known special functions is not well understood yet. Similar challenges are encountered when solving Feynman integrals by direct integration, e.g., in Feynman parameter space. From the purely numerical side, several methods are available to compute Feynman integrals by using Monte Carlo integration techniques [50,51]. As opposed to the analytic approach, these methods are fully algorithmic. Nonetheless multi-loop multi-scale integrals generally present numerical instabilities that make their numerical integration challenging. A different numerical approach has been used in [52][53][54][55], where the solution of differential equations for Feynman integrals is obtained by using Runge-Kutta algorithms.
A third route of exploration has been methods based on series expansions. When it is difficult to obtain a closed form solution for a given integral, it is usually possible to obtain a (generalised) power series expansion of the solution. Series representations have a number of useful features. It is usually possible to compute several orders of the expansion and obtain an arbitrarily good approximation of the full solution. Moreover their numerical evaluation is virtually instant, since each term of the expansion is an elementary or a relatively simple function. All these features make series expansions a natural candidate to solve large classes of complicated Feynman integrals. Series expansion methods have been mostly applied to single scale problems in e.g. [41,[56][57][58][59][60][61]. On the other hand, for integrals depending on several scales, series expansions have been performed with respect to one variable, parametrising special kinematic configurations, while keeping the dependence on the remaining variables exact (see e.g. [62][63][64][65][66][67][68][69]), or to transport analytic boundary conditions in various regions [70]. Therefore, in the multivariate case, it is desirable to study a systematic approach to obtain results in all points of the kinematic regions.
In this paper we reduce the computation of a set of multivariate Feynman integrals [28] to a single scale problem, by defining differential equations along contours connecting two generic points of the kinematic regions. We then find generalised power series solutions by solving the (single-scale) differential equations with respect to the contour parameter (while replacing all the other variables with numbers). In this way the solution can be transported from a base point, where the integrals are assumed to be known, to a generic target point. We show that this approach is efficient, and can be repeated to compute the integrals in any point of the kinematic regions, with high numerical precision. More specifically, we apply this method to a family of planar (elliptic) Feynman integrals relevant for the two-loop QCD corrections to Higgs + jet production, below and above the heavyquark threshold. Previously [28] these integrals were computed in the Euclidean region by using integral representations. Our results are new, and provide at the same time the analytic continuation of these integrals to the physical region, and an efficient method for their numerical evaluation. Further applications of these methods to the non-planar integral topologies are presented in a companion paper [71].
The paper is organised as follows. In Section 2 we review general properties of the differential equations for dimensionally regulated scalar Feynman integrals, and their solution in terms of iterated integrals. In Section 3 we describe the series expansion strategy used in this paper. We show that after defining the (multi-scale) differential equations along a (one-dimensional) contour, the series solution can be obtained by series expanding the differential equations and iteratively integrating them up to the desired order of the dimensional regulator. In Section 4 we apply this strategy to a family of planar integrals relevant for Higgs + jet production in the physical region. We show how, once a series expansion is found, the analytic continuation is performed in a straightforward manner. We finally show high precision numerical results and timings, with comparisons to sector decomposition programs. In Section 5 we draw our conclusions.
2 Differential equations for dimensionally regulated Feynman integrals By using standard IBP reduction techniques, it is possible to identify a basis f ( x, ) for a set of dimensionally regulated scalar Feynman integrals, where f ( x, ) = {f 1 ( x, ), . . . , f n ( x, )} and x = {x 1 , . . . , x m } is the set of kinematic invariants. Given a basis, one is also able to define a system of first order linear differential equations satisfied by the basis, that in full generality takes the form where A x i is an n-by-n matrix that depends rationally on its variables. If the set of basis integrals is minimal the differential equations satisfy the integrability condition, where the last term is a commutator. Nonetheless the applicability of our method does not rely on this condition, and it can be applied also to an overcomplete set of integrals. The basis f is not unique. In [13] it was conjectured that, with a basis change, it is possible to cast differential equations for Feynman integrals in a canonical form, where the dependence on the dimensional regulator is factorised. In differential form the canonical equations have the following form In dimensional regularisation we are interested in a solution around = 0. By series expanding f ( x, ), the solution of Eq. (2.3) can be written in terms of iterated integrals [72]: where γ is a path with domain [0, 1] in the space of external invariants and f ( x 0 , ) is a vector of boundary conditions. If f (i) (x) admits a representation in terms of multiple polylogarithms [14], with G(, t) ≡ 1, the transformation matrix to the canonical basis is algebraic, and the matrixÃ(x) is a Q-linear combination of logarithms where d α is the number of linearly independent logarithms, C i are constant rational matrices, and α i ( x), i ∈ {1, . . . , d α } are called letters of the differential equations. The set of letters is referred to as the alphabet. If the letters admit a representation in terms of rational functions, then the iterated integrals of (2.5) can be directly expressed in terms of multiple polylogarithms. On the other hand, if a rational representation is not found, it is often possible to find a polylogarithmic representation by using the knowledge of the symbol of iterated integrals [14,73]. The symbol of the i-th basis element of (2.4) at order k is obtained by the following recursive formula Once the symbol of the solution is known, it is possible to find a corresponding polylogarithmic expression by using an ansatz for the set of polylogarithmic functions, and by imposing boundary conditions (see e.g, [28,74]). When considering integral sectors that do not admit a polylogarithmic representation, the properties of the transformation matrix to the canonical form are not yet well understood (but see e.g. [75] for recent progress in this direction). For this reason, when a canonical basis is not available, we will only assume that the differential equations are non-singular in the → 0 limit (note that it is always possible to find a basis of Feynman integrals satisfying such differential equations, see e.g. [76]).

Series expansion along a contour
Given a base point where the integrals are known, we show how the integrals are computed in a new point by series expanding the integrals along a contour connecting these points. We first define the differential equations along the contour, and then we show how a series expansion of the integrals is found by solving the differential equations around a set of points of the contour. The procedure can be repeated to reach any point of the kinematic regions. In this section we describe the general framework needed to perform series expansions along contours. In appendix B we discuss in detail a simple one-loop example.
Provided that a set of differential equations with respect to a complete set of kinematic invariants is available, we can define the differential equations along a contour γ(t) connecting two fixed points a = {a 1 , . . . , a m }, b = {b 1 , . . . , b m }. This is achieved by parametrising the contour with a parameter t: and by considering the differential equations with respect to t: where the new differential equations matrix can be readily obtained by using the chain rule, It is known that a set of master integrals, at a given order of the dimensional regulator, admit a solution in the vicinity of a singular point τ of the form (see e.g. [77]) where c (i,j 1 ,j 2 ,j 3 ) are vectors of dimension n (the number of master integrals), S i is a finite set of integers, w k is a complex constant (typically a rational number that accounts for the algebraic dependence of the matrix elements), and N i is the maximal power of the logarithms at order i . On the other hand, in the vicinity of a regular point τ , the integrals admit a standard Taylor series representation More simply, each Feynman integral, in the vicinity of a singular point τ , is expressed as a finite combination of terms of the form where ρ(t) is a Taylor series. On the other hand, in the vicinity of a regular point, Feynman integrals admit a Taylor series representation.
In the remainder of this section we show how the series solutions (3.4) and (3.5) can be found by series expanding the differential equations matrices, and by iteratively integrating them until the desired order of .
We know that, for many phenomenologically relevant processes, most integrals admit a polylogarithmic canonical basis. In the next subsection we discuss how we find a power series solution of a canonical set of differential equations. This form of the equations is usually much more compact then the differential equations for a generic basis of master integrals, and working with simplified differential equations renders their series expansion more efficient (for a discussion about timings see Section 4.4). In section 3.2 we discuss how we obtain series solutions for coupled sectors. We remark that the applicability of our method does not rely on (when it exists) a canonical basis of integrals, and a series solution for an arbitrary basis of integrals can be found by using the methods of section 3.2.

Canonical differential equations
When dealing with single scale problems, generalised power series solutions are usually obtained by defining generic power series, and by fixing the corresponding free coefficients by solving recurrence relations [58,59,78]. Here we proceed in a more direct way, which is particularly suited for differential equations in canonical form. Specifically, we consider a canonical system of differential equations of the form and we assume that the solution is known (analytically or numerically with very high precision, e.g. from a previous expansion) for some t = t 0 . From Eq. (2.5) it is easy to see that the solution is (3.8) Since we are interested in the solution in the vicinity of a point τ , we series expand the differential equations matrix around τ , where A (i) t are constant matrices 1 . For a canonical basis, w i , i ∈ N is expected to be the set of all half integer numbers with w i ≥ −1. However the following discussion holds for generic complex numbers w i . By plugging the previous expansion into (3.8), we get At each order k , we have to compute integrals of the form, (3.11) These integrals can be computed analytically in terms of integer powers of log(t − τ ) and (rational) powers of (t − τ ), by using recursively integration-by-parts identities 2 .
1 The series expansion for the matrix A (i) t can be obtained for example by using the built-in Mathematica function Series. 2 In practice, these integrals can be computed by using the built-in Mathematica function Integrate

Coupled sectors
The problem of finding a canonical basis for integrals that do not admit a polylogarithmic representation is still poorly studied in the literature. In practice, we often deal with differential equations where only a subset of the equations is in canonical form, while the other sectors admit a generic rational dependence on the dimensional regulator and an algebraic dependence on the kinematic invariants. In this case, by iteratively taking derivatives ∂/∂ t of equation (3.2), it is possible to obtain a k−th order differential equation for a single integral, say f where, for ease of notation, we denoted the generic master integral at order i as =j (t), and it is known at every iteration i. The homogeneous equation associated to (3.12) is and it admits k linearly independent solutions that we denote as h 1 (t), . . . , h k (t). We note that the homogeneous equation does not depend on the order under consideration. By defining the following fundamental matrix, the particular solution of (3.12) is where χ (i) j , j ∈ {1, . . . , k} is a set of boundary constants, the Wronskian is defined by W (Ĥ(s)) = det(Ĥ(t)), while W j (Ĥ(s)) is the determinant obtained by replacing the j-th column ofĤ(t) by (0, . . . , 0, 1). It can be shown that, for functions admitting a series solution of the form (3.4), the equation in the vicinity of a point τ can be written as, where the functions b i (t) are analytic in t = τ . In this case a series solution in the vicinity of t = τ can be found by applying the Frobenius method, which is discussed in detail in Appendix A (see also e.g. [79]). More specifically, the Frobenius method allows to find a complete set of k homogeneous series solutions in the vicinity of τ . These solutions have the form where c (j 1 ,j 2 ,j 3 ) i are complex constants, K i is the maximal power of the logarithm, S i ⊆ {1, . . . , k} and λ 1 , . . . , λ k are the roots of the indicial equation, defined as From the form of the homogeneous series solutions, it is clear that the full solution (3.2) will require, at each order i, the computation of integrals of the form (3.10) which, as explained in the previous section, can be done analytically in terms of integer powers of logarithms and complex powers of t. This shows that, also in the case of coupled differential equations, we can explicitly find a series representation for Feynman integrals in the vicinity of regular or singular points of the form of Eqs. (3.5) and (3.4) respectively. The Frobenius method is rather standard, and we review its general formulation for linear differential equations of generic order in Appendix A. Here we briefly show its application to second order equations, since this is the order encountered in this paper. Let us consider the case k = 2 of (3.13) and, without loss of generality, let us assume that t = 0 is a (regular) singular point, We have to distinguish two cases in order to proceed. Let us first assume that the two roots of the indicial equation λ 1 , λ 2 do not differ by an integer number, with λ 1 > λ 2 . In this case two linearly independent solutions are and the coefficients are fixed by requiring that the differential equation is satisfied orderby-order in t. If the two roots differ by an integer, the solution associated to λ 1 is obtained by (3.20) while the second solution is obtained by which can be expressed as a series by expanding all the integrands around t = 0 and performing the integrations term-by-term.

Matching
Given a base point where the integrals are assumed to be known, we want to transport the solution to a new point, by using series expansions along a contour connecting these points.
In the previous sections we have seen that, given a singular or regular point of the differential equations, we can find a series solution valid in the vicinity of these points. By construction, these solutions converge only in a region that does not contain any singularity other then the expansion point. In the following we will consider truncated series that, within a given accuracy, provide a good approximation of the full series. When considering a generic contour γ(t) for a range of values of t, we are interested in finding series expansions along the entire contour. The contour will in general contain multiple singular points of the differential equations, and it is necessary to find multiple series expansions and match them together in order to cover the entire contour. A good criterion for determining the domain of definition of a truncated series is that the series will converge fast enough only when considering it in a region such that the maximum distance from the expansion point is half the distance from the nearest singularity.
Let us assume that we have defined a contour γ(t), for t ∈ [0, 1] and we want to find truncated series expansions that approximate the full solution within a given accuracy. γ(0) is the known boundary point and γ(1) is the target point we want to compute. There might be real and complex singularities. Let us denote the real singularities as R = {τ i |i ∈ N s } and the complex singularities by Moreover we define the following set of real regular points We then find truncated power series around t 0 = 0 and t i for i ∈ {1, . . . , N e +1}. We denote these series as where t i is the expansion point, and r i is the radius of the series defined as the distance between t i and the closest point Since the t i are in general not equally spaced it can happen that some segments of the contour are not cover by any series. In this case we iteratively introduce new (regular) expansion points κ i , i ∈ N , and corresponding series , such that κ i is in the middle of an uncovered region and ρ i is the distance between κ i and the closest point among κ j =i , t j . The procedure is repeated until the entire contour is covered.
There are two special cases that need some care. If there are no points t i , we just expand around t 0 = 0 and the entire contour will be covered. If there are no points t Ne+1 , t Ne+2 we iteratively add new (regular) expansion points such that 0 < k i < t Ne . At the end of this procedure it can happen that the point t = 1 in not covered. We then set t Ne+1 = 1 and we add new (regular) expansion points as described above until the entire contour is covered.
We finally define the solution along the entire contour to be,

22)
r i ≤ r i , ρ i ≤ ρ i are such that there are no overlapping series, and S e ⊆ {1, . . . , N e }, S r ⊆ {1, . . . , N r } (N r is the total number of regular expansion points) are such that no series outside the unit interval t ∈ [0, 1] is considered in the sum. In Sections 3.1 and 3.2 we have seen that the series solutions depend on a set of integration constants that have to be fixed by imposing boundary conditions. When considering a set of series along a contour, the integration constants of one series are fixed by knowing the boundary conditions at a given point (e.g. when the boundary conditions are known analytically by other means, or they are known because the contour under consideration intersects another contour along which the series expansion is already known), while the other series are fixed by imposing that two consecutive series have the same value in the contact point, i.e. the point where they are both defined.

A planar elliptic family for Higgs+jet production
The planar two-loop QCD corrections to Higgs+jet productions are mediated by the four integral families depicted in Fig. 1. These integrals were computed in [28] in the non-physical region, where all the Mandelstam invariants are negative reals. For the polylogarithmic sectors, the solution was expressed in terms of logarithms and dilogarithms up to weight two, while the weight-three and four solutions were expressed in terms of one-fold integrals over the lower weight solutions. Family A contains two elliptic sectors. The first elliptic sector is I A 1,1,0,1,1,1,1,0,0 ; the homogeneous solutions are elliptic integrals, and the solution in dimensional regularisation requires iterated integrations over elliptic kernels. The second elliptic sector is I A 1,1,1,1,1,1,1,0,0 , and it admits a canonical form for the homogeneous differential equations, but it is coupled to the first elliptic sector via inhomogeneous terms. In [28] the solution for the elliptic sectors was expressed in terms of iterated integrals over elliptic kernels.
In this section we obtain generalised power series expansions for family A in the p 2 ↔ p 3 channel, which exhibits the most complicated (spurious, see section 4.3) singularity structure of the channels needed for the scattering amplitude. All the other families are simpler, as they do not involve elliptic structures and can be solved with the same methods.
The family under consideration is defined by the following set of integrals I Ap 2 ↔p 3 a 1 ,a 2 ,a 3 ,a 4 ,a 5 ,a 6 ,a 7 ,a 8 ,a 9 = propagators while d 8 and d 9 are numerators, therefore a 1 , . . . , a 7 are non-negative integers while a 9 , a 8 are non-positive integers. The kinematics is such that p 2 1 = p 2 2 = p 2 3 = 0 with where p 2 4 is the squared Higgs-mass, p 2 4 = m 2 H , and m 2 is the squared mass of the propagators. The relevant physical region is defined by We solve the kinematic constraints for s 12 , s 13 , p 2 4 and we define the following scaleless variables The family contains 73 master integrals, and our choice for the integral basis is provided in the ancillary files of the arXiv submission. The first 65 master integrals satisfy a set of differential equations in canonical form, The alphabet of the differential equations consists of 42 letters, depending on the following set of 6 square roots 3 , (4.7) We remark that our approach does not rely on the rational parametrisation of the set of square roots, and it works for general algebraic dependence of the differential equations. Integrals 66-73 are elliptic and satisfy coupled differential equations of the form where the vector g 66−73 ( x, ) depends on the polylogarithmic integrals f 1−65 ( x, ). Orderby-order in the polylogarithmic integrals f 1−65 ( x, ) satisfy completely decoupled differential equations. On the other hand, the elliptic sectors are coupled and, in general, one needs to solve higher order differential equations for single integrals. Specifically, the homogeneous matrix B (0) x i has the following schematic form where the lines separate sector I A 1,1,0,1,1,1,1,0,0 (first four rows) from sector I A 1,1,1,1,1,1,1,0,0 (last four rows). We see that integrals 67 and 69 are coupled. For each of these integrals we can define a second order differential equation. Once integral 67 or 69 is known, all the other integrals are decoupled and can be solved by considering first order equations only.

Series solution of the differential equations
In order to solve the differential equations we need a set of boundary conditions. We use the point x 1 = x 2 = x 3 = 0, and the boundary conditions are [28] f i ( 0, ) = 1 + π 2 We are interested in a solution in the relevant physical region (4.4). Therefore we transport the boundary condition to the point 2, 13 25 , −1 by series expanding the integrals along the contour 4 , Here we discuss in detail how we expand along the contour γ thr (t) defined as This contour is interesting because it allows to analytically continue the integrals above the physical threshold x 1 = 4. In order to achieve the decomposition of the contour described in Section 3.3, we need to know where the singularities lie. When considering differential equations with algebraic dependence (square roots in the present case), the singularities are all the non-analytic points of the differential equations, i.e. points where the differential equations diverge but also the zeros of the square roots (branching points). Once the path is fixed, the problem is one dimensional and we can find the position of the singularities by solving for the zeros of the denominators of the differential equations and for the zeros of the square roots (see Section 4.2 for a detailed discussion about the different classes of singular points of the differential equations). For the path γ thr (t), the relevant singularities are where τ 1 is the physical threshold and τ 2 is a spurious singularity outside the [0, 1] interval. We now partition the contour as described in Section 3.3. Specifically, we add the regular expansion points κ 1 = 0, κ 2 = 1 κ 3 = 1 8 which results in the following partitioning of the contour, where we replaced rational numbers with (exact) real numbers. All the boundary constants are fixed by imposing the following chain of boundary conditions,

Analytic continuation
In the previous sections we showed that in order to obtain converging power series expansions along a contour it is necessary to expand around the singular points of the differential equations (and regular points in order to ensure fast converging series, as described in Sec. 3.3). As already mentioned, by singular point we mean any non-analytic point of the differential equations and of the solution, i.e. power and logarithmic divergences, or branching points of the square roots. The singularities of the solution are a subset of the singularities of the differential equations. In this section we discuss the different classes of singularities encountered in the (series) solution of the differential equations and how to perform the analytic continuation across them. Moreover, for the sake of the following discussion, we remark that we solve the differential equations along contours entirely contained in the physical sheet. The first class of singularities are the so-called physical singularities, which correspond to the singularities predicted by unitarity cuts. In correspondence of physical singularities the solution develops branch cuts, and in order to analytically continue the solution across them we use Feynman prescription, i.e. we assign a small imaginary part to the contour parameter in such a way that the invariant crossing the singularity acquires a small positive imaginary part. In order to show this, let us write down explicitly the first few terms of the expansion around the heavy-quark threshold t = 1 2 , discussed in the previous section, for, e.g., the (elliptic) integral 68 at order 4 which reads 5 , f (4.16) When expanding along a path crossing a physical singularity, in this case x 1 = 4, the branch cuts are expressed by rational powers of the expansion parameter and logarithms. The analytic continuation is then performed by assigning a small imaginary part to the parameter t, consistently with Feynman prescription. In this case, since x 1 (t) = 2 + 4t, we have t → t + iδ, and, On the other hand, the expansion f The second class of singularities are the so-called non-physical singularities, which are not physical (they are not predicted by unitarity) but they appear in the solution of a given integral basis. In order to explain the origin of these singularities let us consider a generic basis element of the integral basis 6 , where I a 1ij ,...,a nij is a scalar Feynman integral of the form, If one of the coefficients of Eq. (4.18), say c ik ( x), is singular for x = x s , where x s is not a physical singularity and I a 1ik ,...,a nik is non-zero in x s , B( x) is singular for x = x s . This is what we call a non-physical singularity, since it originates only from the prefactors of the integral basis 7 . For integral bases with algebraic coefficients (as the ones considered in this paper) there are two kind of non-physical singularities: poles and branching points of square roots. Poles do not give rise to cuts, and no analytic continuation is needed across them. On the other hand, for square roots giving rise to non-physical branch cuts, we can choose an arbitrary analytic continuation across their branching points, since these cuts cancel at the level of the integrals I a 1ij ,...,a nij and at the level of the scattering amplitude. More specifically, for the integrals of family A the only physical singularities are for x 1 = 4 and x 2 = 4. According to Feynman prescription we define On the other hand, all the other square roots have no x 1 − 4, x 2 − 4 factors and we choose, for simplicity, to define them as √ a = i √ −a when a < 0, for any a and in any region. In Appendix B we discuss a simple one-loop example and we show how non-physical singularities appear in the series solution.
Finally we have the so-called spurious singularities. These are singularities of the differential equations that do not correspond to singularities of the solution. These singularities are not physical and the coefficients of the integral basis are regular in these points. Therefore, since we consider contours entirely contained in the physical sheet, the solution is regular in these points, and the singular terms of the series solution corresponding to spurious singularities cancel within the truncation error, and no analytic continuation is needed. In physical applications it is often necessary to consider integrals for very large or small values of the external invariants. When working with series expansions it is then convenient to map the relevant physical region to a finite region, so that it is possible to expand around limiting points in a straightforward manner. For two-to-two processes a convenient set of variables is the following, which, by keeping the Higgs mass fixed, maps the physical region (4.4) to the unit square in the l, z space. The l, z variables are also convenient when studying the singularities of the differential equations. The singularities are the zeros of the denominators of the differential equations and the zeros of the square roots. For our differential equations, the singular points are for l = 0, l = 1, z = 0, z = 1 and along the orbits defined by the following equations, 8 z = 13 100 , z = 12l 13 + 12l , z = −50 + 37l + 10 . (4.21) The orbit z = 0 corresponds to points where x 1 , x 3 are infinity and these are singular points for our integral basis (see below). l = 0 and z = 1 correspond to x 3 = 0 and are not physical singularities but the series expansions are not analytic there because some of the square roots vanish. l = 1 is a spurious singularity. The first singularity of Eq. (4.21) is the physical threshold. The second and last singular orbits of Eq. (4.21) are spurious. Indeed, the singular orbits do not correspond to physical singularities, and none of the square roots or the rational prefactors of the integral basis are singular along these orbits. The physical region, including the singular orbits, is represented in Fig. 3. The l, z variables are also the natural variables to perform expansions along contours reaching very large or small values of the invariants. Let us consider for example the contour, which corresponds to a contour from {x 1 = 24 In this case the only singular point is for t = 1, and it is sufficient to perform only one expansion around it to cover the full contour. t = 1 corresponds to a (singular) point at infinity. Nonetheless, we never cross the points at infinity and possible branch cut ambiguities are fixed imposing boundary conditions at a finite point and by treating the square roots as explained in Section 4.2. The plots of the elliptic integral sectors along γ ∞ are presented in Fig. 4.

Conclusion
We showed that by defining the differential equations for a set of multi-scale Feynman integrals along contours connecting two generic points of the kinematic regions, it is possible to systematically obtain one-dimensional series expansions for the integrals along the contours. Specifically, we applied this method to obtain new results for a planar family of integrals relevant for the two-loop QCD corrections to Higgs + jet production. When the expansion is performed along a contour such that only one invariant changes along it, the analytic continuation above the physical thresholds becomes straightforward. We demonstrated that performing an expansion along a contour is fast, and makes it possible to repeat the procedure to compute the integrals over the entire kinematic regions, with arbitrary numerical precision. Our approach is algorithmic, and it seems plausible to implement it in a computer code that can be applied to solve complicated integrals of phenomenological integrals.

A General formulation of the Frobenius method
The Frobenius method has general validity, and in this appendix we consider a generic order-k linear differential equation for an unknown function f (t) (the discussion closely follows [79]). The Frobenius method can be used to find a complete set of homogeneous series solutions to the equation in the vicinity of a singular (or regular) point that we assume being located at t = 0. The homogeneous equation can be written in general as: Since the b i (t) are analytic in t = 0 they admit a Taylor series representation of the form, We look for a solution of the form, where lambda is a (complex in general) parameter to be fixed. By substituting the formal series solution into the equation we obtain, Where the g(j) are linear in the c (1) , c (2) , · · · , c (j−1) with polynomial coefficients in λ, while the polynomial f (λ) is called the indicial polynomial and it reads, In order for f (t) to be a solution, the coefficients of the right-hand side of (A.4) need to satisfy which is a recursion relation that can be solved for c (j) with j ≥ 1 if σ(λ + j) = 0 for any j ≥ 1. The c (j) , j ≥ 1 are then uniquely determined functions of c (0) . In this case we get It is convenient to fix the value of c (0) , and we conventionally choose c (0) = 1. If λ = λ 1 where λ 1 is a root of the indicial polynomial, σ(λ 1 ) = 0, and σ(λ + j) = 0 for any j ≥ 1 then, f λ 1 ,1 (t) ≡ f (t)| λ=λ 1 , is a solution of the equation L[f (t)] = 0. If λ 1 has multiplicity 2, we need to find a second solution associated to it. If we differentiate Eq. (A.7) with respect to λ we get, If λ 1 is a root of multiplicity 2 then σ(λ 1 ) = ∂σ(t) ∂λ λ=λ 1 = 0, and ∂f (t) ∂λ is also a solution of Eq. (A.1). We have then that the second solution associated to the root λ 1 is If λ 1 has multiplicity m 1 , all the m 1 solutions associated to λ 1 are obtained by taking m 1 −1 derivatives.
Let us now suppose that λ 1 + k 2 = λ 2 , for a positive integer k 2 , i.e. the root λ 1 differs by an integer k 2 from another root λ 2 . Moreover we assume that λ 2 is the only root that differs from λ 1 by an integer. In this case the solution associated to λ 1 cannot be determined as explained above, since the recursion (A.6) becomes ill defined when j = k 2 . In order to find a solution we set c (0) = (λ − λ 1 ) m 2 , where m 2 is the multiplicity of the root λ 2 , (A.10) An explicit computation shows that Eq. (A.7) becomes and the recursion (A.6) can be written as, and where the functions G (i) are non-zero for λ = λ 1 . We note that the right hand side of (A.13) is now well defined when λ approaches λ 1 since σ(λ 1 + k 2 ) = σ(λ 2 ) = 0 and the denominator has a factor (λ − λ 1 ) m 2 that cancels against the numerator. This shows that The leading term of the series is t λ 2 , therefore this solution is linearly dependent from the solution one would obtain by considering the root λ 2 . In order to find a solution truly associated to λ 1 we take the m 2 -th derivative of f (t), which satisfies 16) and ∂ m 2 f (t) ∂λ m 2 | λ=λ 1 is indeed a solution. Moreover its leading term is m! t λ 1 so that ∂ m 2 f (t) ∂λ m 2 | λ=λ 1 is a genuine solution associated to λ 1 . If λ 1 has multiplicity m 1 > 1 the other solutions are obtained by taking m 1 − 1 derivatives of ∂ m 2 f (t) ∂λ m 2 with respect to λ evaluated at λ 1 . The last case we need to consider is when there are more than one root differing from λ 1 by an integer. More precisely, we consider r−1 roots λ 2 , . . . , λ r with Re(λ 2 ) < · · · < Re(λ r ) such that λ 1 + k i = λ i with k i a positive integer. Each root has multiplicity m i and we define M = r i=2 m i . We set c (0) = (λ − λ 1 ) M , so that By proceeding as in the case of only two roots differing by an integer, one obtains that a solution associated to λ 1 is given by In this appendix we apply the methods of Sec. 3 to a one-loop integral family, and we discuss in detail the derivation of the series expansion along a given contour. Specifically, we consider the integral family relevant for the leading order corrections to Higgs + one jet production in QCD depicted in Fig. 6, while the kinematics is p 2 1 = p 2 2 = p 2 3 = 0 with, where p 2 4 is the squared Higgs-mass, and m 2 is the squared mass of the propagators. We solve the kinematic constraints for s 12 , s 13 , p 2 4 . A canonical basis for the integral family is given by, The corresponding differential equations are, where the matrixÃ is, The letters are defined as l 1 = log (s 12 ) , l 2 = log (s 13 ) , l 3 = log p 2 4 , l 4 = log m 2 , l 5 = log 4m 2 − s 12 , l 6 = log 4m 2 − s 13 , l 7 = log 4m 2 − p 2 4 , l 8 = log s 13 − p 2 4 , l 9 = log s 12 − p 2 4 , l 10 = log −p 2 4 + s 12 + s 13 , l 11 = log s 12 s 13 − 4m 2 −p 2 4 + s 12 + s 13 , l 12 = log r 3 r 4 − p 2 4 , l 13 = log (r 5 r 6 − s 13 ) , l 14 = log (r 1 r 2 − s 12 ) , l 15 = log (r 1 r 5 r 7 + s 12 s 13 ) , l 16 = log 2m 2 −p 2 4 + 2s 12 + s 13 + r 1 r 6 r 7 − s 12 s 13 , l 17 = log 2m 2 −p 2 4 + s 12 + 2s 13 + r 2 r 5 r 7 − s 12 s 13 , l 18 = log −2m 2 p 2 4 s 13 − p 2 4 + s 12 p 2 4 + s 13 + p 2 4 s 12 s 13 + r 1 r 3 r 4 r 5 r 7 . (B.9) In order to discuss the solution we define, as in the main text, the following scaleless variables, We consider the contour, The differential equations along the contour are, It is easy to see that the differential equations along the contour are singular for t = 0 and t = 2/3, the latter being the physical threshold s 12 = 4 while the former is a non-physical singularity (see discussion below). Following the prescription of Sec. 3.3 we have , (B.14) We start with the computation of B [γ 1 ] (t, ) [0,0,1/3] . We series expand ∂Ã( x(t)) ∂t around t = 0.
By using the shorthand a i,j ≡ ∂a ij ∂t we have 9 , where C (i) is the boundary vector and the integral on the right-hand-side is an indefinite integral. We remark that by series expanding ∂Ã( x(t)) ∂t all the integrals in (B.16) are elementary and can be performed analytically. At order 1 we have [ The order 0 of the solution is given by the boundary conditions (B.11). By performing the integrals on the right-hand-side and by imposing the boundary condition B (1)  It is straightforward to iterate the procedure up to arbitrary order of . We see that the solution is singular in t = 0, since it develops a brunch cut for t < 0. This is what we call, in Section 4.2, a non-physical singularity, since the only physical singularity along the contour is s 12 = 4m 2 . The appearance of this singularity is due to the prefactors defining the canonical basis (B.4) and we see that, by inverting the basis, the integrals I a 1 ,a 2 ,a 3 ,a 4 admit a regular Taylor series expansion around t = 0. We remark that, according to (B.4), the differential equations depend originally on √ −t and, as discussed in Section 4.2, we consider the standard brunch for √ −t, i.e. √ −t > 0 for t < 0 and √ −t = i √ t for t > 0. The computation of B [γ 1 ] (t, ) [1/3,2/3,1] follows the same steps. One first expands ∂Ã( x(t)) ∂t around t = 2/3. The series solution is then obtained by, where t = 2/3 − t. We see that for t > 2/3 the solution develops a brunch cut, and the cut ambiguity is resolved by Feynman prescription which, by (B.12), is t → t + iδ with δ a small positive imaginary part.

C Plots
In this appendix we show plots of all the 73 master integrals at order 4 . The plots are obtained by series expanding along the contour γ thr defined in section 4.1 (s 13 = −1, p 2 4 = 13 25 , m 2 = 1, 2 ≤ s 12 ≤ 6). The solid dots represent numerical values computed with FIESTA 4.1. The real part of the integrals is in blue, the imaginary part is orange.