Exact results for ZmOS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Z}_m^{\mathrm{OS}} $$\end{document} and Z2OS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Z}_2^{\mathrm{OS}} $$\end{document} with two mass scales and up to three loops

We consider the on-shell mass and wave function renormalization constants ZmOS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Z}_m^{\mathrm{OS}} $$\end{document} and Z2OS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Z}_2^{\mathrm{OS}} $$\end{document} up to three-loop order allowing for a second non-zero quark mass. We obtain analytic results in terms of harmonic polylogarithms and iterated integrals with the additional letters 1−τ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{1-{\tau}^2} $$\end{document} and 1−τ2/τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{1-{\tau}^2}/\tau $$\end{document} which extends the findings from ref. [1] where only numerical expressions are presented. Furthermore, we provide terms of order O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(ϵ2) and O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(ϵ) at two- and three-loop order which are crucial ingredients for a future four-loop calculation. Compact results for the expansions around the zero-mass, equal-mass and large-mass cases allow for a fast high-precision numerical evaluation.


Introduction and notation
Once quantum corrections to quantities, which involve heavy quarks, are computed to higher orders in perturbation theory the renormalization of the mass and wave function has to be performed. The corresponding renormalization constants are usually denoted by Z OS m and Z OS 2 , respectively. They are defined through where m 0 and ψ 0 stand for the bare quark mass and wave function. The superscript "OS" refers to the on-shell scheme, which for QCD corrections is used synonymous to the pole scheme. Within QCD, analytic results for both renormalization constants are available up to three loops [2][3][4][5][6][7][8][9]. At four-loop order [10][11][12][13] semi-analytic methods were used. 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 figure 1. The case m i = 0 (i = 1, 2) and m 1 = m 2 was considered in refs. [3,4] and [1] at two-and three-loop orders (see also refs. [14,15]). In this work we re-consider these contributions to Z OS m and Z OS 2 up to threeloop order and provide analytic results including O( ) terms. In ref. [1] only expansions for m 2 /m 1 → 0 and numerical results have been provided up to the constant term in .
In the Standard Model the largest quark mass ratio comes from charm and bottom quarks. In fact, with increasing experimental precision in B physics observables it becomes more and more important to include charm mass effects even at higher orders in perturbation theory. A prominent example is the semileptonic decay b → c ν which is currently known to O(α 2 s ) [16][17][18]. The extension to O(α 3 s ), necessary to determine |V cb | with a relative uncertainty of about 1%, requires the two-mass renormalization constants considered in this paper. They must be know analytically in all three limits m 2 /m 1 → 0, m 2 /m 1 → 1 and m 2 /m 1 → ∞. Indeed, to renormalize the mass and wave function of the JHEP10(2020)087 In addition, the future analysis of the data from the MUonE experiment [19,20] will likely require the knowledge of the dominant α 3 em corrections to muon-electron scattering [21][22][23]. As it often happens for QED processes, it is not possible to set the electron mass m e = 0 in the evaluation of virtual corrections, since charged leptons are experimentally distinguishable from their collinear photon radiation (contrary to what happens in QCD to quarks in jets). Therefore in QED, keeping finite fermion masses is often a necessity in order to regularize collinear singularities. Even in the case when finite electron mass effects are restored via massification [24] of virtual amplitudes computed for m e = 0, the procedure employs the lepton's wave-function renormalization constant with finite m e effects to the relevant order in α em . The results presented in this paper thus apply also to such processes, with the proper translation of QCD color factors to QED.
A further, at first sight not obvious, application of the three-loop two-mass result for Z OS m is the renormalon analysis of the relation between the top quark pole and MS mass [25,26] with the aim to determine the ultimate uncertainty of the top-quark pole mass. Note that the typical loop momentum at order α n+1 s scales as m t e −n and thus light quark mass effects are important at higher loop orders although at first sight m b /m t and m c /m t seem to be negligibly small.
It is convenient to introduce the variable We furthermore adopt the notation from [1] and write (i = m, 2) with the SU(N c ) colour factors C F = (N 2 c − 1)/(2N c ), C A = N c and T F = 1/2. We have introduced the quantities n l , n h and n m to label closed quark loops with mass zero, m 1 and m 2 , respectively, see also table 1. We have n f = n l + n m + n h = n l + 1 + 1 active quark flavours. Note that only the terms proportional to n m and n 2 m have a non-trivial dependence on x. This is the main subject of the present paper.
For the quark mass renormalization constant we also introduce the ratio which is finite since both Z OS m and the MS renormalization constant Z MS m only contain ultra-violet poles, which cancel in the ratio. Note that Z OS 2 contains both ultra-violet and infra-red poles. z m has an analogous perturbative expansion as Z OS m and Z OS 2 in eq. (1.3). In ref. [1] the n m -dependent terms of Z OS m and Z OS 2 were computed up to three loops. At two-loop order analytic results were obtained. However, at three-loop order, for the complicated master integrals only an expansion for x → 0 could be obtained. For larger values of x a numerical evaluation was necessary. For most practical purposes this is sufficient. However, in some cases analytic expressions or expansions are useful. In this work we extend the result of [1] in the following aspects: • We extend the expansion by one order both at two and three loops, which is necessary input for a future four-loop calculation of the n m terms of the on-shell renormalization constants.
• We provide analytic results in terms of iterated integrals with the letters τ,
In the next section we briefly describe the approach which we use to obtain the analytic results and the expansions in the various limits. In section 3 we discuss our results for the renormalization constants and give our conclusions. In the appendix we provide details to the supplementary material attached to this paper.

Technicalities
We base our three-loop calculation on intermediate expressions obtained in ref. [1]. In particular, we use the results where Z OS m and Z OS 2 are expressed in terms of the 28 master integrals shown in figure 2. They are obtained from three integral families introduced in refs. [1,27]. We have implemented the integral families in LiteRed [28,29] and redone the redution to the 28 master integrals shown in figure 2. After using the additional relation [27] one ends up with 27 master integrals, which have to be computed. We utilize LiteRed [28,29] to derive a closed system of differential equations. Out of the 27 master integrals 23 can be solved in terms of harmonic polylogarithms (HPLs). Their analytic results are given in [27]. There it was also noted that four master integrals (M 20 , M 21 , M 22 , M 23 ) cannot be expressed in terms of HPLs at higher orders in = (4 − d)/2. For these integrals expansions around x = 0 were obtained and numerical values for larger values of x were calculated. In this work we obtain analytic results for the missing master integrals and, furthermore, extend all master integrals by one order in such that we obtain the renormalization constants to O( 2 ) at two and to O( ) at three loops, respectively. This will be a crucial input for a future four-loop calculation of the renormalization constants.
We solve the coupled system of differential equations using the algorithmic approach presented in ref. [30]. For the convenience of the reader we outline in the following the main steps of this approach. The differential equation can be written in the form where M (x, ) is the vector of our 27 master integrals. It can be chosen such that the matrix A is in upper-block diagonal form, i.e. the diagonal elements are square matrices with possible non-vanishing entries to the left. The square matrices on the diagonal represent coupled sets of master integrals, which only depend on themselves and integrals from lower sectors. One can therefore solve the system successively starting from simpler systems and insert the solutions as inhomogeneities into the more involved ones. In total we find one 3 × 3 and seven 2 × 2 systems. The x dependence of seven of the remaining ten master integrals factorizes. JHEP10(2020)087

JHEP10(2020)087
We decouple the systems of differential equations with the package OreSys [31], which is based on Sigma [32,33], to obtain a single differential equation of higher order for one of the master integrals in the system. Furthemore, OreSys provides rules to construct the other master integrals from the solution of the differential equation. The higher order differential equations are then expanded in and iteratively solved order by order. To solve the differential equations we make use of the solver implemented in HarmonicSums [34][35][36][37][38][39][40][41][42][43][44][45][46], which is particularly well suited to find solutions in terms of iterated integrals. In a first step we consider the homogeneous part of the differential equation and try to write it in factorized form. If it fully factorizes into first order factors the solution in terms of iterated integrals can be obtained in a straightforward way. If second order factors remain Kovacic's algorithm [47] is used to find all solutions of the differential equation, which can be written in terms of iterated integrals. In our case the homogenous solutions in terms of iterated integrals exist and thus also the particular solutions can be expressed in terms of iterated integrals. The construction and simplification of the homogeneous and particular solution is automated in HarmonicSums. To fully solve the differential equations we still need to fix the integration constants multiplying the homogeneous solutions. Boundary values for all integrals at x = 0 and x = 1 can be extracted from the on-shell integrals given in ref. [48]. To fix all integration constants we need both limits since for some master integrals the homogeneous solutions vanish at x = 0 or x = 1.
Four master integrals cannot be expressed in terms of usual HPLs. We want to illustrate this for the system of differential equations of the integrals M 22 and M 23 . After decoupling the homogenous differential equations for the master integral M 22 we obtain where d = 4 has been used in the coefficients since the -dependent terms enter the inhomogeneous part. Equation (2.3) has the two solutions where I denotes a generalized iterated integral over the specified integration kernels indicated in the curly brackets with the help of the variable τ , i.e.
Note that a regularization is needed for letters which lead to divergent expressions for t → 0. This is in complete analogy to HPLs [49]. Equation (2.4) illustrates that one has to introduce the new letter √ 1 − τ 2 /τ in order to solve the differential equation. Analogously the system of master integrals M 20 and M 21 introduces the letter √ 1 − τ 2 . After fixing the boundary conditions it turns out that for all master integrals the generalized letters are JHEP10(2020)087 only needed from O( ) onwards. Note that the O( ) terms enter the finite contribution of Z OS m and Z OS 2 . Since the additional letters only introduce one square root it is possible to rationalize the letters with a suitable variable transformation. One possibility is the so-called trigonometric substitution which introduces the letters 1 1 + τ 2 , τ 1 + τ 2 . (2.7) Iterated integrals over these kinds of letters have been studied in ref. [37]; the corresponding iterated integrals are called cyclotomic HPLs. Alternatively one can factor the polynomial over the complex numbers and introduce Goncharov polylogarithms [50] with letters taken from the set of the 4th root of unity. For example one has which shows that the variable transformation in eq. (2.6) converts HPLs with argument x into cyclotomic HPLs with argument y. Note, however, that the transformation in eq. (2.6) significantly increases the complexity of the rational functions in the differential equations. Thus, we have chosen to solve them in the variable x. However, eq. (2.6) is needed to fix the boundary conditions at x = 1, since this requires the evaluation of the iterated integrals at this point. The corresponding results up to weight 5 are conveniently obtained by transforming the iterated integrals to cyclotomic HPLs for which the values at x = 1 are known up to weight 6 [37]. This leads to relations like For the expansion around x = 1 we first map the argument of the iterated integrals to 1 − x. This can be achieved iteratively with the formula which can be easily proven from the integral representation. In our case this step does not introduce new letters, but introduces the iterated integrals at argument x = 1. The same constants were already needed to fix boundary conditions for the differential equations. Afterwards we can expand easily around 1 − x.
The expansion for x → ∞ is more involved since the letters involving square roots develop a branch cut for x > 1. Thus, in a first step we have to construct the analytic continuation for the iterated integrals, i.e., the relations for the corresponding functions with argument x < 1. We use differential equations to do this. Let us for illustration consider an iterated integral of weight one. Then we have Now we change the variable to z = 1/x and find where f (z) is the analytic continuation of the iterated integral in eq. (2.15). We assume 0 < z < 1 in accordance with x > 1. Note that in our case the change of variables again does not introduce new letters. The differential equation can be easily solved by integrating JHEP10(2020)087 the right hand side over z and fixing the integration constant for x = z = 1. This again only requires the knowledge of the iterated integrals at argument x = 1. For our example, we obtain (2.17) For higher weights one can proceed iteratively, since the derivative of an iterated integral of weight w with respect to its argument only depends on iterated integrals of weight w − 1.
Note that the analytic continuation of the individual iterated integrals introduces imaginary parts (cf. eq. (2.17)). However, after inserting the analytic continuations for all iterated integrals into the expressions for Z OS m and Z OS 2 all imaginary parts cancel analytically and the expansion around 1/x = z = 0 can be obtained in a straightforward way.

Results and conclusions
In this section we briefly discuss our results for Z OS m , z m and Z OS 2 . After inserting the exact master integrals into the corresponding amplitudes we renormalize the quark masses m 1 and m 2 in the on-shell scheme, the strong coupling constant in the MS scheme and expand in such that we obtain results up to 2 at two-loop and 1 at three-loop order. Whereas the two-loop results are still quite compact (see, e.g., eqs. (15) and (28) of ref. [1]), the threeloop expressions are too big to be printed. Instead we provide the analytic expressions in the supplementary material attached to this paper. We also provide transformation rules which map the iterated integrals introduced in the previous section to Goncharov polylogarithms which can be evaluated numerically with the help of GiNaC [51]. Note that our final threeloop result contains iterated integrals up to weight five and six in the 0 and 1 term.
More compact expressions are obtained after expanding for x → 0, x → 1 or x → ∞. For illustration we show for the n m dependent terms of z m , which we define via the first three expansion terms at two and three loops. To keep the expressions compact we specify the colour factors to QCD, i.e. C A = 3, C F = 4/3 and T F = 1/2. Furthermore we set µ = m 1 , n h = 1 and restrict ourselves to the 0 term. For x → 0 we obtain  be found in ref. [1]. In the supplementary material attached to this paper we provide for z m and Z OS 2 different variants for the expansions in the three limits x → 0, x → 1 and x → ∞.
We update the Mathematica routines provided in ref. [1] for the numerical evaluation of z m and Z OS 2 . In the supplementary material attached to this paper one finds the functions zmnum[x,m1,mu1,mu2[,scheme]] and Z2OSnum[x,m1,mu1,mu2[,scheme]] (see appendix) which implement the expansion for x → 0, x → 1 and x → ∞. We switch between the first two expansions at x = 1/2 and between the latter two at x = 3/2. The justification for this choice is illustrated in figure 3, where we show the expansions for z (3),M m for = 0, n l = 3, n h = n m = 1 and µ = m 1 . In the regions where the expansions converge (x < 1/2, 1/2 < x < 3/2, 3/2 < x for x → 0, x → 1 and x → ∞, respectively) we plot solid lines and outside these region we switch to dotted lines. One observes that both around x = 1/2 and x = 3/2 there is a large overlap among at least two expansion, which justifies that we use the expansion results to contruct the functions zmnum[x,m1,mu1,mu2[,scheme]] and Z2OSnum[x,m1,mu1,mu2[,scheme]]. Let us also mention that we observe an agreement with the exact result to at least 8 digits over the whole range in x. Similar results are obtained for the O( ) term of z m and for Z OS 2 .

JHEP10(2020)087
To summarize, we have obtained analytic results of all 27 master integrals which are needed to obtain the three-loop contributions for the on-shell renormalization constants Z OS m and Z OS 2 with dependence on two different quark masses, m 1 and m 2 . Our final result includes terms of O( ), which are relevant for a future four-loop calculation. Furthermore, we have obtained 26 expansion terms for three cases m 1 m 2 , m 1 ≈ m 2 and m 1 m 2 .
JHEP10(2020)087 api lmm1 cf ca tr nl nm nh α s (µ)/π log(µ 2 /(m OS 1 ) 2 ) C F C A T F n l n m n h x xOSOS xMSOS xOSMS xMSMS mu1 mu2 x Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.