A semi-analytic method to compute Feynman integrals applied to four-loop corrections to the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document}-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 MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} 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 fiveloop 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.

JHEP09(2021)152
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 "normalized 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,19]. 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. [20,21] and implemented in the code DiffExp [21]. 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 [22][23][24][25][26][27][28]. 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 JHEP09(2021)152 type can be found in figure 1. The case 0 = m 2 = m 1 was considered in refs. [23,29] at two-loop and in refs. [30,31] at three-loop order (see also refs. [32,33]).
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 : +8 further colour structures involving n m +23 further colour structures without n m , (1.4) 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. (1.4) are given in section 4. For Z OS m we assume the same decomposition as in eq. (1.4).
For the class of integrals considered in this paper we can obtain analytic results expressed in terms of Goncharov polylogarithms [34]. 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 < ∞. One has to be able to compute the master integrals for (at least) one limiting behaviour of x. For our application the limit x 1 leads to a straightforward asymptotic expansion and to integrals which are available

JHEP09(2021)152
in the literature. Other special points in our calculation are the equal-mass case x = 1 and the limit of vanishing mass m 2 , x = 0. In general there are also other exceptional points depending on the physical application under consideration. 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,35] or apply the methods developed in refs. [12,36,37]. 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 extent, 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.
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 [38]. 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.

JHEP09(2021)152
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 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. (1.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 ) [31] and presented in ref. [37]. 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 [39,40] with FireFly [41,42] 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 JHEP09(2021)152 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 [44] 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 into 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 block-triangular form, i.e. the diagonal elements are square matrices with possible nonvanishing 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 [45], which is based on Sigma [46,47], 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 [13,[48][49][50][51][52][53][54][55][56][57][58][59]. 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.
As an illustrative and representative example we consider a coupled 2 × 2 sub-system of our differential equations which is given by The inhomogeneities R 1 (z, ) R 2 (z, ) depend on previously calculated master integrals and are known as an expansion in . They are given by where H a ≡ H a (z) are harmonic polylogarithms. For simplicity we only show the quartic and cubic poles. Using OreSys we can decouple the system in eq. (3.1) and obtain a differential equation of second order for F 1 (z, ), and an equation which can be used to compute F 2 (z, ), The new inhomogeneities are given bỹ To solve the differential equation in eq. (3.3) we can expand in . In each order in the homogenous differential equation is the same and can be factorized as The next task is to solve the factorized differential in eq. (3.6). The solution of a first order differential equation is straightforward: the homogenous solution is given by 8) and the particular solution can be given by variation of constants as The general solution is then given by with an arbitrary constant c 1 . For our application the lower integration limits in eqs. (3.8) and (3.9) can be chosen to be 0. In general one might have to choose a different value which

JHEP09(2021)152
modifies the constant c 1 in eq. (3.10). This procedure can be generalized to a differential equation with n factors of the form Denoting the homogenous solutions of the individual factors by h k , k = 1, . . . , n the homogenous solutions of the differential equation in eq. (3.11) are given by , (3.12) and the particular solution can be written in the form This algorithmic way of solving the differential equation is implemented in the package HarmonicSums and naturally leads to iterated integrals. The letters can be read off from the denominators of eqs. (3.12) and (3.13). In our application we find a factorization as in eq. (3.11) for all differential equations and the results for the master integrals contributing to the amplitude can be written in terms of iterated integrals (3.14) with letters drawn from the set The first three letters define the harmonic polylogarithms, the following two letters are also needed for the calculation at O(α 3 s ) [31] while the last two cyclotomic letters are needed for the n m · n l colour factors at O(α 4 s ). The application of this algorithm to the example in eq. (3.6) leads to the two homogenous solutions The higher order terms in quickly become more complicated. The general solution for F 1 can be cast into the form where the constants c 1,(k) and c 2,(k) have to be fixed by comparing to boundary conditions which in general are obtained from an explicit calculation of the Feynman integrals in a certain limit. In the following paragraph we describe our approach for the computation of the boundary conditions. For the present calculation 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 [60,61] which generates for each master integral all relevant sub-and co-subgraphs according to the rules of the hard mass procedure [38]. 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. [62][63][64][65][66][67][68][69][70].
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.
For the integral in eq. (3.17) we find for the boundary condition

JHEP09(2021)152
Using eq. (3.4) we can immediately obtain results for the integral F 2 (z, ). Note that we only need the boundary conditions of F 1 (z, ) but obtain solutions for both integrals. The boundary conditions calculated for the integral F 2 (z, ) can be used to check the consistency of the solutions. Iterating this procedure to the next orders in we find (3.20) and This provides results for F 1 and F 2 up to O( −2 ):

JHEP09(2021)152
Note that the homogeneous and particular solutions individually depend on the square root valued letter f w 2 . However, this dependence drops out in the pole terms of F 1 (z) and F 2 (z) and only starts to contribute to the solutions from O( 0 ) onwards. We fix the boundary conditions for m 2 m 1 . 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 x → 1/x, so we do not need to introduce new letters for the physical region. In the following we show an example of how to obtain the necessary analytic continuation.
Let us consider the integral I({f w 2 , f w 1 }, z). From the definition of iterated integrals we have Now we change the variable to x = 1/z with 0 < x < 1 and obtain We see that if we are able to find a solution of I({f w 1 }, z = 1/x) from eq. (3.25) we can insert it into eq. (3.24) and obtain the desired solution of I({f w 2 , f w 1 }, z = 1/x) also in terms of iterated integrals. Integrating eq. (3.25) we find with an arbitrary constant C 1 . We choose to fix the constant at x = 1/z = 1. Here we have where the last term comes from the limit x → 1 of eq. (3.26).
To fix the boundary values it is necessary to compute the iterated integrals for special values, in our case for z = 1. For harmonic polylogarithms this is straightforward. For the iterated integrals containing the letters f w 1 and f w 2 we perform the variable change z = 2t/(1 + t 2 ). In this new variable t the iterated integrals turn into cyclotomic harmonic polylogarithms with the letters f {4,0} and f {4,1} for which the constants at argument 1 have been calculated [71]. The variable change can be done using the same method of differential equations shown above. However, since no analytic continuation is needed we

JHEP09(2021)152
can fix the boundary constants by fixing the limit z = 0 = t where no new constants are generated.
In our example we obtain Inserting eq. (3.28) into eq. (3.24) and performing the integration over x we now get (3.29) In the limit x = 1/z = 1 we have On the other hand, the limit x → 1 on the r.h.s. of eq. (3.29) leads to Thus, the final solution for the analytic continuation of I({f w 2 , f w 1 }, z) is given by setting C 2 = −iπ 2 /8 in eq. (3.29).
To obtain the final results for our calculation we had to analytically continue iterated integrals up to weight 5, which can be done iteratively starting from lower weights as demonstrated above.
As an illustration for an analytic result we show the colour structure C F n m n 2 h of z m which is given by

JHEP09(2021)152
(3.32) In this expression only harmonic polylogarithms are present. It is interesting to note that the unrenormalized expressions of all contributions with three closed fermion loops, except z F M LL m , receive contributions from iterated integrals over the letters f w 1 and f w 2 . However, these contributions are canceled against terms coming from renormalization so that the final results are given by the usual harmonic polylogarithms. This is not the case for z F F M L m and z F AM L m . The corresponding full analytic results are too long to be printed here but are given in the ancillary file to our paper [72].
In the following we present the first three expansion coefficients around x → 0 and x → 1 for the eight color factors under consideration. They have been obtained from the analytic results. In the ancillary file one also finds the first 51 terms of the respective expansions. For the expansions x → 0 we obtain with l 2 = log (2), l x = log(x), a i = Li i (1/2) and Riemann's ζ-function ζ i = has been chosen.
In the limits x = 0 and x = 1 the results in eqs. (3.33) and (3.34) determine certain colour factors in the one-mass limit computed in refs. [5,6]. Using the notation of ref. [5] we have  [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 figure 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 threeloop 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. [64] 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 figure 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,73]. 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   figure 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 figure 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 JHEP09(2021)152 from these expressions. In the expansions for the terms proportional to the color factor n m n l also Catalans constant contributes. For some of the remaining 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. [74,75]). 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 [72].