A semi-analytic method to compute Feynman integrals applied to four-loop corrections to the $\overline{\rm MS}$-pole quark mass relation

We describe a method to numerically compute multi-loop integrals, depending on one dimensionless parameter $x$ and the dimension $d$, in the whole kinematic range of $x$. The method is based on differential equations, which, however, do not require any special form, and series expansions around singular and regular points. This method provides results well suited for fast numerical evaluation and sufficiently precise for phenomenological applications. We apply the approach to four-loop on-shell integrals and compute the coefficient function of eight colour structures in the relation between the mass of a heavy quark defined in the $\overline{\rm MS}$ and the on-shell scheme allowing for a second non-zero quark mass. We also obtain analytic results for these eight coefficient functions in terms of harmonic polylogarithms and iterated integrals. This allows for a validation of the numerical accuracy.


Introduction
The techniques used for the computation of single-scale integrals are quite advanced. It has become standard to compute massless propagator and massive vacuum integrals up to four loops (see, e.g., Refs. [1,2]). Recently even the master integrals for massless five-loop propagator integrals have been computed [3]. On-shell integrals up to four-loop order have been considered in Refs. [4][5][6], where, however, most integrals are only available in numerical form.
In this paper we discuss an approach to compute multi-loop integrals which involve two dimensionful scales, m 1 and m 2 , and thus depend on the dimensionless quantity x = m 2 /m 1 . Among other things (which are described below) it requires single-scale integrals as input and thus a routine application at four-loop order is possible. Our method allows to obtain numerical results with the help of differential equations in the whole kinematic region of x. Details are presented in Section 2.
There are a number of algorithms in the literature which can be used to obtain analytic and/or numeric results of multi-loop integrals. Some of them make heavy use of integration-by-parts relations and many approaches exploit the power of differential or difference equations [7][8][9][10]. For example, in Ref. [11] difference equations are constructed by raising one of the denominators of a given single-scale master integral to an arbitrary power x. The solution of the difference equations leads to high-precision numerical results of the original single-scale integral.
The main idea of Ref. [12] is to construct difference equations for the coefficients of the Laurent expansion in x of the master integrals by plugging in a suitable ansatz into the system of differential equations. Once the difference equations are established, a large number of expansion terms can usually be obtained rather efficiently so that one can try to find a closed form solution using other methods (see e.g. Ref. [13]). However, the method presented in Ref. [12] cannot be applied for points where the master integrals obey power-log expansions.
An interesting approach to obtain numerical results of loop integrals has been presented in Ref. [14] where an imaginary mass is added to all propagators. The differential equations with respect to this mass are solved numerically using the infinite-mass limit as boundary.
In the approach of Ref. [15], the differential equations are solved numerically after obtaining boundary conditions by performing expansions in a suitable limit. Afterwards, the numerical results can be used to determine the coefficients of expansions around arbitrary points [16].
Our approach is close to the method presented in Ref. [17]. This approach uses expansions around singular points together with the system of differential equations to transfer the information of the integrals from a starting point, x 0 , where boundary conditions are available, to another point, x 1 . Note, however, that the program which comes together with Ref. [17], DESS.m, requires the system of differential equations in the so-called "nor-malized global Fuchsian form", which is not required in our algorithm. Furthermore, the authors of Ref. [17] aim for high-precision results with several hundred significant digits, necessary to reconstruct analytic expressions with the help of the PSLQ algorithm [18]. On the other hand, the goal of our method is the construction of an approximation for all values of x.
Another similar approach is discussed in Refs. [19,20] and implemented in the code DiffExp [20]. It is more general than our approach in the sense that the differential equations are solved without the need of an appropriate ansatz in a respective kinematic point. It aims at integrating multi-scale integrals along line segments. Our approach is taylored to problems which depend on only one dimensionless parameter with known analytic properties and is optimized for multi-loop problems with large coupled systems of differential equations in mind.
As an application of our method we consider the relation between a heavy-quark mass defined in the pole (or on-shell) scheme (m OS 1 ) and the MS scheme (m 1 ), which is given by where z m is finite and depends on the renormalization scale µ. For convenience we also introduce the on-shell renormalization constant via where m 0 1 is the bare quark mass. For the pertubative expansion of z m we write with an analogous definition of Z OS m . Within QCD, analytic results for z m are available up to three loops [21][22][23][24][25][26][27]. At fourloop order semi-analytic methods were used [4][5][6]. Starting from two loops there are contributions with closed quark loops, which can either be massless, have the mass of the external quark (m 1 ), or have a different mass (m 2 ). Sample Feynman diagrams of this type can be found in Fig. 1. The case 0 = m 2 = m 1 was considered in Refs. [22,28] at two-loop and in Refs. [29,30] at three-loop order (see also Refs. [31,32]).
In this work we concentrate on the four-loop contributions which involve at least one additional closed quark loop with mass m 2 . We introduce the symbol n m to count the number of such loops. In analogy n h counts the closed loops of fermions with mass m 1 and n l the massless ones. To demonstrate our method we consider the following eight (out of 16 in total) four-loop colour structures which involve n m : with C F = (N 2 c − 1)/(2N c ), C A = N c and T F = 1/2 for an SU(N c ) gauge group. Note that in practical applications we have n f = n l + n m + n h = n l + 1 + 1 active quark flavours. Numerical result for the coefficients of Eq. (4) are given in Section 4. For Z OS m we assume the same decomposition as in Eq. (4).
For the class of integrals considered in this paper we can obtain analytic results expressed in terms of Goncharov polylogarithms [33]. We describe our calculation in Section 3.

Method
In this section we describe a method to obtain numerical results of Feynman integrals. To be concrete we consider a set of master integrals which depend on the variable x (usually a ratio of two kinematic invariants) and the dimension d = 4 − 2 . Let us assume that we are interested in the results for the integrals for 0 ≤ x < ∞ and that an analytic calculation of the master integrals for x 1 is possible.
In such cases one often proceeds as follows: One establishes the differential equations with respect to the variable x for the master integrals [7,8]. If they are sufficiently simple, a direct integration is possible and with the help of the boundary conditions for large x an exact solution can be constructed. Often it is helpful to transform the differential equations into -(or canonical) form [10,34] or apply the methods developed in Refs. [12,35,36]. There are a number of non-trivial examples where these approaches have provided analytic results in terms of harmonic polylogarithms or even more general functions. However, there are severe limitations. For example, it is not always possible to construct an -form and thus a simple solution of the differential equations is harder to obtain. The method, which is described in the following, does not have such limitations. In fact, it is, to a large extend, insensitive to the complexity of the differential equations since only expansions around certain kinematical points are considered. In particular, it can be applied in case the differential equations contain elliptic sub-systems.
To apply our method one has to be able to reduce the integrals to a set of master integrals and to establish a system of differential equations for the latter. Furthermore, it must be possible to compute the master integrals in a given limit of the variable x. In this limit it is allowed that the master integrals obey a power-log expansion. It is not necessary to bring the system of differential equations into a particular form (e.g. Fuchsian form) or require that the x-and -dependence factorizes in the denominators. Furthermore, it is not necessary that the set of master integrals is minimal.
Our algorithm consists of the following steps: 1. Reduce all contributing Feynman integrals to master integrals.
2. Establish the system of differential equations for the master integrals.
3. Compute boundary conditions for the system of differential equations, i.e., evaluate the master integrals for x approaching some limit. For clarity, let us consider x 1.
In our case, since m 2 is an internal mass scale, one can apply the hard mass procedure [37]. This leads to vacuum integrals of the considered loop order and products of lower-loop integrals.
4. Expand the differential equations in this limit and insert an ansatz for the master integrals. In general the ansatz is a power-log expansion with even and odd powers of x. One can use the boundary conditions to fix the constants in the leading term(s) of the ansatz. Afterwards the expanded differential equations are used to obtain a deep expansion in 1/x.
In our application only even powers are present in the limit x 1. Typically we compute 50 expansion terms. 5. Expand the differential equations for x ≈ 1 and insert an ansatz for the master integrals in this limit.
In our application the expansion around x = 1 is a simple Taylor expansion.
6. Choose a value x 1 where both the expansions in 1/x and around x = 1 converge. We evaluate the (known) 1/x-expanded master integrals for x = x 1 and use the results as boundary conditions for the expansions around x = 1. A typical value for our application is x 1 ≈ 1.5.
Proceeding this way one, in principle, ends up with an over-determined system of linear equations which in general has no solution due to the numerical errors introduced by truncating the expansions. To circumvent this issue we proceed as follows: We start with the simplest master integrals and fix the constant for the leading pole. We then apply this relation to the whole expression and proceed with the next term in the -expansion and repeat this for all master integrals. If more than one unknown constant appears, we solve for one of them. At some point we encounter equations which are linearly dependent on equations solved before. In this case we get a numerical value which would be equal to 0 if we had the exact boundary conditions at x = x 1 . We store such values and use their absolute value to examine the validity of our procedure. One observes that they get smaller and smaller the more terms in the expansion are used. With the help of the differential equations we can again compute 50 terms in the expansion around x = 1.
7. In a next step we repeat the same procedure to perform a matching between the expansions around x = 0 and x = 1.
Note that the expansion of the master integrals around x = 0 contains again both monomials in x and log(x). In contrast to the 1/x expansion both even and odd powers in x are present in general. For our application a proper matching point is It might happen that the differential equations have further singularities at x = x s with x s = 0, 1 or ∞, even in case the physical amplitude does not have thresholds for this value of x. In that case we perform a similar matching at the intermediate value x = x s . For the integral families considered in this paper we did not encounter such additional singularities.

MS-OS relation at four loops: analytic results
For the colour factors which we consider (see Eq. 4) it is possible to obtain analytic results in terms of iterated integrals. We use the same techniques already discussed for the calculation at O(α 3 s ) [30] and presented in Ref. [36]. Let us briefly discuss our procedure.
After generating the amplitudes we map each diagram to an integral family and express it as a linear combination of scalar functions with 14 arguments, the powers of the scalar propagators. We use Kira [38,39] with FireFly [40,41] for the reduction to 339 master integrals out of 18 contributing integral families. 1 At that point we could attempt for a symmetrization of the master integrals over all different integral families. However, we prefer to consider each integral family separately. Of course, there are several master integrals which are present in different families. We use the comparison as a cross check.
To obtain the differential equations for this set of master integrals we use LiteRed [43] for the differentiation with respect to x and reduce the result again with the help of Kira.
To obtain a closed set of differential equations we have to consider a larger set of integrals than the one present in the amplitude. In total we compute the analytic solution of 520 master integrals. Since our boundary constants are computed in the limit x → ∞ we cast the differential system in the form with z = 1/x and M (z, ) the vector of master integrals. Working with the variable z instead of x allows to fix the boundary at z → 0 and thus no analytic continuation has to be performed when solving the differential equations. In the end, however, an analytic continuation to the region z > 1 (i.e. x < 1) is necessary.
The vector M (z, ) can be chosen in such a way that the matrix A(z, ) is in lower blocktriangular form, i.e. the diagonal elements are square matrices with possible non-vanishing entries to the left. The square matrices represent coupled systems of master integrals. We find at most 5 × 5 systems in our calculation. To solve the coupled systems of differential equations we utilize OreSys [44], which is based on Sigma [45], to decouple the systems of equations and obtain a higher-order differential equation for one of the master integrals in the respective system. Furthermore, we obtain rules which allow to construct the other master integrals of the coupled system from its solution. The higher order differential equation is solved with the help of HarmonicSums [46]. Internally the solver factorizes the differential equation and, if successful, finds the solution in terms of iterated integrals without the need of specifying an alphabet.
The boundary constants are fixed using the expansions in the limit x 1, i.e. the limit in which m 2 is much larger than m 1 . For the computation of the boundary conditions for m 2 m 1 we use the program exp [47,48] which generates for each master integral all relevant sub-and co-subgraphs according to the rules of the hard mass procedure [37]. In some cases up to eleven subdiagrams are generated. By construction the subdiagrams are one-to four-loop vacuum integrals where the relevant scale is given by m 2 . On the other hand the co-subgraphs are propagator-type on-shell integrals up to three loops. They only depend on the mass scale m 1 . All relevant integral families are well studied in the literature and the master integrals can be found in Refs. [49][50][51][52][53][54][55][56][57].
We compute the expansion for each master integral up to order 1/x 4 . Note that only a subset of this information is needed in order to fix the constants in the x → ∞ ansatz. Since we decouple the differential equations, we can choose to fix the boundary constants by considering the leading term in the limit x → ∞ of every master integral in the system or by considering higher orders in the expansion for the master integral which remains after decoupling. We chose the latter approach. Still we need at most expansions up to z 2 = 1/x 2 . All remaining expansion coefficients of the master integrals that we do not need to fix the boundary constants are used to cross-check the consistency of our results.
The results for the master integrals contributing to the amplitude can be written in terms of iterated integrals with letters drawn from the set The first three letters define the harmonic polylogarithms, the following two were also needed for the calculation at O(α 3 s ) [30] while the last two cyclotomic letters are needed for the n m · n l colour factors at O(α 4 s ). To arrive at a result in the physical region z = 1/x > 1 we need to analytically continue the iterated integrals. For the iterated integrals involving the square roots, i.e. f w 1 and f w 2 , we use differential equations to obtain the analytic continuation while the iterated integrals involving the other letters can be analytically continued using HarmonicSums. Note that the above letters are closed under the analytic continuation, so we do not need to introduce new letters for the physical region. An example how to obtain the analytic continuation can be found in Ref. [30].
In the following we present results for x = 0 and x = 1, which correspond to particular colour factors in the one-mass limit computed in Refs. [5,6]. The analytic results in the variable x are too lengthy to be printed here and can be found in the ancillary material. We find with l 2 = log (2) [5] which shows that the uncertainty estimate made there is correct.

MS-OS relation at four loops: numerical results
In the following we apply the algorithm described in Section 2 to the on-shell-MS relation. We start from the same systems of differential equations and boundary conditions already needed for the analytic calculation discussed in the previous section. Note that also in this approach only a subset of the boundary conditions are needed in order to fix the constants in the x → ∞ ansatz. In fact, in general we need at most expansions up to z 2 = 1/x 2 . Furthermore, to fix the boundary conditions of the differential equations not all master integrals are needed. All remaining coefficients are used to cross-check the consistency of our results. For example for family d4L456 (see Fig. 2), the family which introduces the cyclotomic letters in the analytical result, we find 33 master integrals. Only the boundary constants of 18 master integrals, expanded at most up to the constant contribution in the large-m 2 limit, are needed.
For some families we observe spurious poles from the reduction up to 1/ 4 , which requires an expansion of the boundary integrals up to order 4 . For most of the three-loop on-shell and four-loop tadpole integrals, which appear in the boundary conditions, an expansion to such high order is not available. We parameterize the unknown coefficients and check that they drop out in the physical result. Alternatively, we could have used the algorithm in Ref. [51] in order to construct an -finite basis. However, for the colour structures considered in this paper this was not necessary.
We are now in the position to discuss our results. For the renormalization scale we again use µ = m OS 1 . We start with the bare expressions and show in Fig. 3 the results of the 0 term for the colour structures C 2 F n m n l and C F C A n m n l for 0 ≤ x ≤ 1. The light and dark blue curves show the expansion results around x = 0 and x = 1, respectively. The red dots denote the known results for m 2 = 0 and m 2 = m 1 [5,58]. Note that in the case of C 2 F n m n l the limit x → 0 does not exist for the bare expression. In fact, there are logarithmic divergences which arise from diagrams containing a fermion self energy (see Fig. 1(g) for an example). The corresponding analytic expression is given by   Fig. 3 shows that there is a substantial cancellation of more than an order of magnitude between the bare expressions and the counterterm contributions. Nevertheless, we manage to reproduce at least five digits of the exact result.
The renormalized results for the six n 3 f colour factors are shown in Fig. 6. Also here we reproduce the exact result with a precision between 5 and 10 digits.
The results discussed above are based on expansions around x = 0, 1 and ∞ involving 50 terms. We have repeated the analysis also with fewer expansion terms and observe a fourto six-significant digit agreement with the exact expressions in case we use 45 terms for the matching. One quickly loses precision in case even less terms are used which is mainly due to the inability to match the x → ∞ and x → 1 expansions properly. This problem could probably be cured by introducing further matching points. We have checked that, in case we increase the number of expansion terms in the matching step from 50 to 60, the agreement with the exact result increases by about one significant digit.

Conclusions
In this paper we present a numerical method to compute multiloop integrals with two mass scales, i.e., one dimensionless parameter. The prerequisites necessary to apply our algorithm are quite simple: It is necessary that the reduction problem can be solved, that the differential equations can be established and that boundary conditions can be computed for some limit of the dimensionless parameter. We have shown that our approach works for systems of differential equations which involve a few hundred master integrals.
As an application we have considered eight colour factors for the four-loop relation between a heavy-quark mass defined in the MS and on-shell scheme where a second quark mass, m 2 , is present in a closed loop. Analytic boundary conditions are obtained in the large-m 2 limit.
For the considered colour structures we were able to obtain analytic results which allowed us to quantify the numerical precision to about ten digits in the whole kinematic range. For some of the other colour factors we have observed that not all results can be expressed in terms of iterated integrals and thus an analytic calculation is much more involved or even impossible with current techniques. However, our numerical approach can be applied.
We have demonstrated that our approach can reproduce logarithmic divergences with high accuracy. It could thus also be applied to compute corrections to the electron contribution of the anomalous magnetic moment of the muon which develops log(m e /m µ ) terms in the limit m e m µ (see, e.g., Refs. [59,60]). Furthermore, one can consider non-fermionic on-shell integrals by introducing an artificial mass in such a way that analytic boundary conditions can be obtained. We defer further applications to future work.
Our analytic and most of the numeric results can be found in the ancillary file to our paper [61].
after the Higgs Discovery". The Feynman diagrams were drawn with the help of Axodraw [62] and JaxoDraw [63].