Three loop massive operator matrix elements and asymptotic Wilson coefficients with two different masses

Starting at 3-loop order, the massive Wilson coefficients for deep-inelastic scattering and the massive operator matrix elements describing the variable flavor number scheme receive contributions of Feynman diagrams carrying quark lines with two different masses. In the case of the charm and bottom quarks, the usual decoupling of one heavy mass at a time no longer holds, since the ratio of the respective masses, $\eta = m_c^2/m_b^2 \sim 1/10$, is not small enough. Therefore, the usual variable flavor number scheme (VFNS) has to be generalized. The renormalization procedure in the two--mass case is different from the single mass case derived in \cite{Bierenbaum:2009mv}. We present the moments $N=2,4$ and $6$ for all contributing operator matrix elements, expanding in the ratio $\eta$. We calculate the analytic results for general values of the Mellin variable $N$ in the flavor non-singlet case, as well as for transversity and the matrix element $A_{gq}^{(3)}$. We also calculate the two-mass scalar integrals of all topologies contributing to the gluonic operator matrix element $A_{gg}$. As it turns out, the expansion in $\eta$ is usually inapplicable for general values of $N$. We therefore derive the result for general values of the mass ratio. From the single pole terms we derive, now in a two-mass calculation, the corresponding contributions to the 3-loop anomalous dimensions. We introduce a new general class of iterated integrals and study their relations and present special values. The corresponding functions are implemented in computer-algebraic form.


Introduction
The heavy flavor corrections to deep-inelastic scattering for pure photon exchange are known to leading [2] and next-to-leading order (NLO) [3] 2 . The present accuracy of the deep-inelastic world data requires next-to-next-to leading order (NNLO) QCD analyses in order to determine the strong coupling constant α s (M 2 Z ) [5][6][7] to ∼ 1% accuracy at NNLO, to obtain highly accurate values for the charm and bottom quark masses m c and m b , and to make precise determinations of the parton distribution functions. All of this is in turn needed to describe precision measurements at the LHC [8] and at facilities planned for the future [9,10].
In the region of large scales Q 2 m 2 , analytic expressions for the heavy flavor Wilson coefficients have been obtained at NLO [11,12]. A factorization relation valid in this asymptotic region was given in Refs. [11,13]. For the structure function F 2 (x, Q 2 ), the asymptotic corrections are sufficient at scales Q 2 /m 2 > ∼ 10, cf. [11]. The massless corrections at NNLO to the deepinelastic structure functions are available [14][15][16], while for the corresponding massive corrections in the asymptotic limit, a series of moments has been calculated in the single heavy mass case [1] for all contributing terms in neutral current deep-inelastic scattering. The calculation of the general expressions for the Wilson coefficients is still underway. The asymptotic Wilson coefficients for the structure function F L (x, Q 2 ) have been completed [17,18]. Here the first genuine two-mass contributions emerge at fourth order in the coupling constant. In the case of the structure function F 2 (x, Q 2 ), all corrections to the color factors O(N F T 2 F C A,F ) have been obtained in [19,20], which provides the complete results for two out of five contributing Wilson coefficients, cf. also [18]. The flavor non-singlet corrections have been calculated in Ref. [21] and the flavor pure singlet terms in Ref. [22]. The massive operator matrix elements (OMEs) calculated in [18,19,21,22] are also needed to describe the variable flavor number scheme (VFNS) in the case of a single heavy quark transition [13], for which also the gluonic contributions A gq,Q and A gg,Q are required and have been calculated at 3-loop order in [23] and in [20,24,25], respectively. 3 Technical aspects of these calculations have been described in [27][28][29]. Heavy quark corrections to charged current deep-inelastic processes have been dealt with in Refs. [30].
In the calculations mentioned above, besides internal massless fermion lines, only a single heavy mass is attached to massive fermion lines. However, starting at 3-loop order, there are also diagrams with two different masses attached to the massive lines. In the present paper, we consider corrections of this type. As before in the single heavy mass case [1], a series of finite moments for all massive OMEs and the Wilson coefficients in the asymptotic region Q 2 m 2 c,b is calculated. In some cases, we also compute the results at general values of the Mellin variable N and the momentum fraction z. Furthermore, we present the scalar two-mass integrals contributing to the OME A gg both in z-and N -space, in extension to the single mass case in Ref. [24]. In the present paper, we concentrate on the calculation of the two-mass effects in the case of massive OMEs, playing a central role in the variable flavor number scheme, and leave phenomenological studies of the contributions to various deep-inelastic structure functions for a separate publication.
The paper is organized as follows. In Section 2, the general formalism is outlined, describing the Wilson coefficients in the asymptotic region in the case of two massive quarks and the representation of the deep-inelastic structure functions. We also present the transition relations between a representation of three and five massless quarks to 3-loop order, which is governed by the massive OMEs and describes the matching conditions in the VFNS. In Section 3, the renormalization of the massive OMEs is described in the case of two massive flavors. Here we also derive the structure of the massive OMEs, which now receives logarithmic contributions depending on two masses. The fixed moments for N = 2, 4 and 6 are calculated for all massive OMEs in Section 4, for which we also present numerical illustrations. We have reported on a few results already briefly in [31][32][33]. In the flavor non-singlet and gq-cases, we have calculated the massive OMEs for general values of the Mellin variable N . These are presented in Section 5 and are numerically illustrated. In Section 6, we turn to the more involved case of the genuine two-mass contributions to the massive OME A (3) gg,Q , and outline the calculation strategy, which is significantly different from those of the easier cases being dealt with in Section 5. In the present paper, we limit the consideration to the calculation of all scalar 4 3-loop diagrams contributing to A (3) gg,Q , both in N -and z-space, leading to new functional structures. Unlike the case for the moments, cf. Section 4, where we can expand in the mass ratio of the heavy quarks, this is in general not possible in the case of the diagrams contributing to A (3) gg,Q for general values of N . Therefore, as in Section 5, we derive the analytic solution for general values of the mass ratio. Section 7 contains the conclusions. The z-space results of a series of OMEs are given in the Appendix A, and a collection of new root-valued iterated integrals is presented in Appendix B.

Massive OMEs and Wilson Coefficients with two masses
Starting at 3-loop order, Feynman diagrams carrying internal fermion lines of different mass contribute to the OMEs. The relevant masses are those of the charm and bottom quark, m c and m b . In the following, we will work in the on-shell scheme. Here the masses are given by [34,35]  with m 2 = m c , m 1 = m b , amounts to η ∼ 1/10. Later we will also use the symbol η 1 = √ η. The two masses do not form a strong hierarchy and charm cannot be assumed to be massless at µ 2 = m 2 b . The asymptotic decoupling thus rather proceeds under the condition with Q 2 the virtuality of the exchanged gauge boson in the deep-inelastic process and µ 2 the factorization scale, which we will set equal to the renormalization scale in the following. The transition relation to the MS-scheme for the mass renormalization will be given in Section 3.7. We refer to the on-shell scheme in the following for computational reasons, rather than giving preference to this scheme. In any data analysis, the mass effects shall be expressed in the MS-scheme, which provides perturbative stability. In view of this, the associated variable flavor number scheme (VFNS) differs from the one in which only a single heavy quark is decoupled at the time [1,13], which also works up to 2-loop order since there no diagrams containing fermion lines of different mass contribute.
In the following, we will mainly work in Mellin space to take advantage of the simplicity of the emerging convolution formulae, which are given by ordinary products. The Mellin transform of a function f (x) is defined by (2.5) The convolution formula of two functions reads In what follows, we will use the Mellin transform to map between the z-and the Mellin N -spaces.
Let us now derive the massive Wilson coefficients for deep-inelastic scattering in the kinematic range of large virtualities Q 2 , cf. (2.4). We generalize the considerations in the case of a single heavy quark mass in Refs. [1,13] and obtain the following factorization relation in the non-singlet case: (2,L) N, N F , (2,L) N, N F + 2, Here N F denotes the number of massless flavors (with N F = 3 in QCD). C j i and A kl are the massless Wilson coefficients, cf. [14,36,37] and massive operator matrix elements (OMEs), respectively.
Above and in what follows we use the notatioñ (2.14) The double tilde inH PS q,(2,L) andH g, (2,L) should not be interpreted as applying Eq. (2.13) twice. Instead, it is used to differentiate these Wilson coefficients from those of the single mass case, indicating now the required sum over charges as made explicit later in Eqs. (2.23) and (2.24).
The massive operator matrix elements are the expectation values of the local twist τ = 2 operators O j , obtained in the light cone expansion [38] of the products of electromagnetic currents, The partonic states |i(p) , with i = q (quark) or i = g (gluon), are on-shell with p 2 = 0. In Eqs. (2. 16-2.18) Sp is the color-trace, and S denotes the symmetrization operator 5 of the Lorentz indices µ 1 , . . . , µ N . D µ is the covariant derivative, ψ and ψ are the quark and anti-quark fields, and F a µν the gluonic field strength tensor, with a the color index in the adjoint representation. Furthermore, λ r is the flavor matrix of SU (N F ). The labels q, g on the left-hand side of Eqs. (2.16-2.18) distinguish quarkonic and gluonic operators.
For convenience we will express the strong coupling constant by a s = α s /(4π) ≡ g 2 s /(4π) 2 in the following. Expanding the expressions (2.8-2.12) up to O (a 3 s ) we obtain: g, (2,L) qg,Q (N F + 2) δ 2 + A (1) gg,Q (N F + 2)N FC (2) g,(2,L) (N F + 2) +A (2) gg,Q (N F + 2)N FC (1) g,(2,L) (N F + 2) + A (1) Qg (N F + 2)N FC gg,Q (N F + 2)C (1) g,(2,L) (N F + 2) +C (2) g, (2,L) Qg (N F + 2) C (1),NS q,(2,L) (N F + 2) + A (2) gg,Q (N F + 2)C (1) g,(2,L) (N F + 2) + A (1) Qg (N F + 2) C (2),NS q,(2,L) (N F + 2) +C (2),PS q,(2,L) (N F + 2) + A (1) gg,Q (N F + 2)C (2) g,(2,L) (N F + 2) +C (3) g,(2,L) (N F + 2) . (2.24) Here the symbol δ 2 takes the values δ 2 = 1 for F 2 0 for F L . (m 1 , m 2 ) denotes the part for which the current couples to the fermion-loop of the heavy quark of mass m 1 . This line is carrying the respective local operator. In general, the following representation holds (2.27) The charge-weighted OME is thus given bỹ In the case of the structure function F L (x, Q 2 ), the asymptotic massive 3-loop corrections are obtained by the massive OMEs up to 2-loop order only and therefore do not contain genuine two-mass contributions, cf. [17,18]. The inclusive deep inelastic structure functions F i (x, Q 2 ), i = 2, L can be represented in the fixed flavor number scheme in terms of their purely massless contributions and the remaining terms consisting of the real and virtual heavy quark contributions, (2.29) Since the parton distribution functions are related to massless partons only, F massless i (x, Q 2 ) is obtained in a completely massless calculation. One finds with Σ and ∆ k the flavor singlet and non-singlet distributions given by and G denoting the gluon density. The heavy quark part is given by (2,L) x, N F + 2, (2,L) x, N F + 2, (2,L) x, N F + 2, (2,L) x, N F + 2, The presence of diagrams with c-and b-quarks at 3-loop order also yields power corrections in η to the massive operator matrix elements 6 . One obtains the following transition relations decoupling both the charm and bottom contributions at high scales µ 2 : The flavor singlet, non-singlet and gluon densities for (N F + 2) flavors are given by Here f k(k) (N F ), Σ(N F ) and G(N F ) denote the massless quarkonic parton densities. Note that the above process independent leading twist OMEs A i,j for fixed moments N contain besides logarithmic corrections in η also power corrections. For general values of N the η-dependence is more involved and requests at least generalized harmonic sums [39,40] and binomially weighted generalized harmonic sums [41] as will be shown below. 7 We would like to mention, that although ∆ k is the genuine flavor non-singlet distribution, sometimes the combination f k + f k may be considered to take its role, [11,21].
The presence of 2-mass terms in Eqs. (2.34-2.38) only allows to define the new parton densities at (N F +2) out of those at N F at sufficiently high decoupling scales µ 2 m 2 1 , m 2 2 at 3-loop order, while up to 2-loop order, flavors can technically be decoupled one by one, if m 2 2 m 2 1 (which is not the case, however for b-and c-quarks). The picture of an individual charm and bottom quark density does therfore not hold from 3-loop order onwards. The quantities f k + fk, Σ, ∆ k and G are not affected, as they depend on all heavy quark masses in a symmetric way. The two-mass generalization (2.35) of the single mass case [1,13], is a formal relation as it stands. It can be rewritten expressing the charm and bottom quark densities in the variable flavor scheme, still requesting We turn now to the calculation of the massive two-mass OMEs and discuss first their renormalization in the case of two heavy quark masses.

Renormalization of the Massive Operator Matrix Elements
The Feynman integrals contributing to the various operator matrix elements contain mass, coupling, ultraviolet operator singularities, and collinear divergences, due to massless sub-graphs. They are regularized by applying dimensional regularization [44] in D = 4 + ε dimensions. The singularities appear as poles in the Laurent series in ε, with the highest pole corresponding to the loop order. At one and two loop order the two-mass massive operator matrix elements A ij are given in terms of the known single mass contributions since they do not contain more than one internal massive fermion line [11-13, 17, 18, 45-47]. The first single particle irreducible diagrams with two masses emerge at O(α 3 s ). In the following, we consider the renormalization of the two mass contributions in individual terms together with the genuine two-mass contributions. The latter terms will then be obtained subtracting the former ones, cf. Ref. [1]. The unrenormalized OMEs are given bŷ are the single-mass OMEs [1] andÂ (l) ij are the new two-mass contributions. The last term in Eq. (3.1) for l = 3 contains a factor (m 1 m 2 /µ 2 ) 3/(2ε) . Furthermore, a change in the renormalization scheme as in Eqs. (3.41, 3.42) generally introduces a mixing between the different components of Eq. (3.1).
In the main steps we follow the renormalization procedure outlined in Ref. [1], incorporating the necessary modifications for the two-mass case. We consider the case of N F massless and two massive quark flavors as this covers the physical case of contributions e.g. due to the charm and bottom quarks.
We first consider mass and coupling constant renormalization, followed by the renormalization of the ultraviolet singularity of the local operators, and the factorization of the collinear singularities.

Mass Renormalization
The schemes most frequently used for the mass renormalization are the MSand the on-mass shell scheme (OMS). In the following, we renormalize the mass in the OMS and provide the finite renormalization to switch to the MS-mass at a later stage, cf. Eq. (3.140). We perform the mass renormalization first, i.e. the respective expressions are still containing the bare couplinĝ a s =ĝ 2 s /(4π) 2 . 8 The bare massesm i , i ∈ {1, 2} are expressed by the renormalized on-shell masses m i viâ and Here δm 0 2 is the single mass-contribution, whereasδm 2 i denotes the additional contribution emerging in the case of two massive flavors. Note that from order O(â 2 s ) onward the Z-factor renormalizingm 1 depends on m 2 and vice versa. For the massive operator matrix elements this can be observed at 3-loop order for the first time. The coefficients δm 1 and δm 2 have been derived in [52,53] up to O(ε 0 ) and O(ε −1 ), respectively. The constant part of δm 2 was given in [48,54,55] and the O(ε)-term of δm 1 in [1]. One obtains has been dropped as they are independent of the renormalized mass m i . Furthermore, H a (ζ) are the harmonic polylogarithms (HPLs) [56] H 0 (ζ) = ln(ζ) (3.11) H −1,0 (ζ) = Li 2 (−ζ) + ln(ζ) ln(1 + ζ) (3.12) (3.13) Eq. (3.9) states the complete analytic form of the contribution of the respective other massive flavor to the renormalization of the bare masses. In the present analysis we will focus on m 1 , m 2 being the masses of the bottom and charm quarks, respectively. Due to the size of the ratio it is enough to do the expansion up to O (η 4 ln(η)), as we will do in general for the fixed Mellin moments of the OMEs. The mixed-mass terms are given bỹ

Renormalization of the Coupling
When renormalizing the coupling constant, it is important to note that the factorization relation (2.8-2.12) strictly requires the external massless partonic legs of the operator matrix elements to be on-shell, i.e.
with p the external momentum of the OME. This condition would be violated by naively applying massive loop corrections to the gluon propagator. We follow [1] and absorb these corrections uniquely into the coupling constant by using the background field method [58][59][60] to maintain the Slavnov-Taylor identities of QCD. In this way, one first obtains the coupling constant in a MOM-scheme. A finite renormalization to transform to the MS-scheme is applied subsequently.
The light flavor contributions to the unrenormalized coupling constant in terms of the renormalized coupling constant in the MS-scheme read Here the coefficients δa MS s,i (N F ) are given by with β k (N F ) the expansion coefficients of the QCD β-function [61][62][63][64][65][66] We split the renormalized gluon self-energy Π into the purely light and the heavy flavor contributions, Π L and Π H , The heavy quarks are required to decouple from the running coupling constant and the renormalized OMEs for µ 2 < m 2 1 , m 2 2 which implies [11] Π H (0, m 2 1 , m 2 2 ) = 0 . We apply the background field method, which has the advantage of producing gauge-invariant results also for off-shell Green's functions, to compute the heavy flavor contributions to the unrenormalized gluon polarization function [58,67]. Applying the respective Feynman rules [68] one obtainŝ where the masses m 1 and m 2 have been renormalized in the on-shell scheme (3.2). In order to write (3.27) in a more compact form we use the notation (3.28) and keep this factor unexpanded in the dimensional regularization parameter ε for the moment. Furthermore, we denote the contributions to the QCD β-function coefficients by β Eq. (3.27) differs from the sum of the two individual single-mass contributions [1] by the last term only, which is due to additional reducible Feynman diagrams in the cases of two heavy quark flavors of different mass. The background field is renormalized using the Z-factor Z A which is split into light and heavy quark contributions, Z A,L and Z A,H . It is related to the Z-factor renormalizing the coupling constant g via Concerning the light flavors, we require the renormalization to correspond to the MS-scheme with N F light flavors The heavy flavor contributions are fixed by condition (3.25) which implies The Z-factor in the MOM-scheme is read off by combining Eqs. (3.39) Finally, we express our results in the MS-scheme. For this transition we assume the decoupling of the heavy quark flavors. The transformation to the MS scheme is then implied by   (3.17) we obtain the OME after mass and coupling renormal-izationÂ where we have suppressed the dependence on the masses, ε and N in the arguments of the OMEs.

Operator Renormalization
Next we remove the ultraviolet divergence of the different local operators defined in Eqs. (2.16-2.18) by introducing the respective Z-factors In the singlet case, the operator renormalization introduces a mixing between the different operators as they carry the same quantum numbers. Analogously to the OMEs, here the Z-factors are split into the flavor pure singlet (PS) and non-singlet (NS) contributions Each Z-factor is associated with an anomalous dimension γ ij via Here both the anomalous dimensions and the operator Z-factors obey perturbative series expansions in the coupling constant In order to renormalize the respective operators, we first consider operator matrix elements with off-shell external legs as a sum of massive and massless contributions: Here the massless contribution depends on a MS s since the MOM-scheme, cf. Section 3.2, has been constructed in such a way that it corresponds to the MS-scheme concerning the renormalization of the light quark flavor and gluon contributions.Â Q ij denotes any massive OME we consider. The term δ ij , which appears in the expansion of the OMEs (see Eqs. (3.17) and (3.43)), does not have any mass-dependence and is considered a part of the light flavor partÂ ij −p 2 µ 2 , a MS s , N F . We first consider the renormalization of the purely massless contribution in the MS-scheme In the non-singlet and pure singlet cases one has respectively. The Z-factors describing the ultraviolet renormalization of the complete operator matrix elementsÂ ij p 2 , m 2 1 , m 2 2 , µ 2 , a MOM s , N F + 2 are obtained by inverting (3.55-3.57) and replacing N F → N F + 2. Finally, the transformation (3.42) is applied. The resulting operator Z-factors read: (3. 60) Here and in the Eqs. (3.55-3.57) we have dropped the N F -dependence of the anomalous dimensions γ ij and β i for brevity. The inverse Z-factors for the purely light-parton case correspond to (3.58-3.60) after substituting N F + 2 → N F and δa MOM s,i → δa MS s,i . We are only interested in performing the ultraviolet renormalization for the massive contributions to the operator matrix element in (3.52) and thus subtract the contributions stemming from purely light parts agaiñ Finally, the limit p 2 → 0 is performed. Since scale-less diagrams vanish if computed in dimensional regularization, only the Born piece of the massless OME contributeŝ One obtains the UV-renormalization prescriptioñ Here Z-factors at N F + 2 flavors describe the massive case (3.58-3.60) while those with argument N F denote the Z-factors for the massless case.

Collinear Factorization
At this point only collinear singularities remain. They arise from massless subgraphs only and are therefore independent of the additional heavy quark flavor considered in these analyses. We thus follow [1] directly and remove the collinear singularities via mass factorization Note that in a fully massless scenario the transition functions Γ ij would be related to the light flavor renormalization constant via cf. [11]. However, in the presence of one or more heavy quark flavors the transition functions stem from the corresponding massless subgraphs only. Due to this and the subtraction of the δ ij -term in the OMEs after ultraviolet renormalizationÃ Q ij the transition functions contribute up to O(α 2 s ) only. The renormalized OME is then obtained by Eq. (3.66) differs from the corresponding renormalization and factorization prescription for one heavy quark flavor [1] only by the definition of the renormalization constants Z . Now the term δ ij is added back to the massive OME. In a final step, the coupling constant is transformed to that in the MS-scheme via Eq. (3.41).

One-particle reducible contributions
We will perform the renormalization of the massive operator matrix elements starting from the set of Feynman diagrams which also include the one-particle reducible contributions. These terms contribute from O(α 2 s ) onward and are obtained by quark and gluon self-energy contributions to the external legs of lower order one-particle irreducible diagrams. From 3-loop order onward the reducible contributions to the OMEs A Qg and A gg,Q may contain three different heavy flavors, while this is not the case for the irreducible contributions. Note that the inclusion of the top quark in a loop of the irreducible terms for A (3) ij would demand to consider the energy range Q 2 m 2 t . At a scale µ 2 m 2 t , both charm and bottom can be dealt with as effectively massless. The emergence of massive top loops in the reducible contributions is accounted for by renormalization. In the following we will strictly consider the case of two heavy flavors only.

Self-energy contributions
The scalar self-energies are obtained by projecting out the Lorentz-structurê We decompose the irreducible two-mass self-energies into contributions which depend on one mass only and an additional part stemming from diagrams containing both heavy quark flavorŝ Up to two-loop order no diagrams with two heavy flavors contributễ The single-mass contributions for the gluon are known from [1,70,71] and for the quark self-energy, Similarly to other massive processes [1,72] the constant irreducible diagrams with two different masses contribute for the first time. For the gluonic case, we compute the respective diagrams up to O(η 3 ) using the codes Q2E/Exp [73,74], The quarkonic self-energy contributions have been computed analytically in η,

The reducible operator matrix elements
As in Eqs. (3.71-3.72) we define the two-mass OMEs at one-loop order and the irreducible OMEs at O(α 2 s ) byÂ

The General Structure of the Massive Operator Matrix Elements
In the following, we present the structure of the different unrenormalized and renormalized OMEs for the genuine two-mass contributions.
In the case of only one heavy quark flavor with mass m [1], the mass dependence of the unrenormalized massive operator matrix element at order α l s is given bŷ Here the OMEÂ (l) ij ε, N does not depend on the mass explicitely anymore. It exhibits poles in the dimensional parameter ε up to ε −l (3.95) We adopt the notation of Ref. [1] and denote The unrenormalized operator matrix elements with two massive fermion flavors with masses m 1 = m 2 are split into the respective single-mass contributions (3.94, 3.95) and a part The two-flavor contributionsÂ to the massive OMEs do not obey a factorization relation as (3.94) and the mass dependence is pulled into the coefficients of the Laurent expansionÂ In the following, a (l,k) , a (l) , a (l) without argument will denote the single mass-quantities corresponding to the definitions in   It is technically advantageous to perform the renormalization on the complete two-flavor OMEsÂ µ 2 , ε, N only, which is obtained after subtracting the respective single-mass contributions [1,75].
The analytic expressions for the respective single mass contributions and renormalization constants to two-loop order, which appear in subsequent relations, have been given in Refs. [1,12,15,16,46] and references therein.
The lowest non-trivial flavor non-singlet (NS) contribution is of O(a 2 s ), Starting from O(a 3 s ) it exhibits a non-trivial two-mass contributioñ The renormalized two-mass OME in the MOM-scheme is obtained from the bare quantities combining Eqs. (3.43, 3.66). It is given by After a finite renormalization to the MS-scheme and the subtraction of the single-mass contributions one obtains the pole-structure of the two-flavor piece bŷ The renormalized expression in the MS-scheme is given bỹ For N = 1 the OME vanishes due to fermion number conservation; this applies both for the anomalous dimensions γ

Qq
Depending on whether the operator couples to a heavy or a light fermion, there are two puresinglet contributions [1] Up to O(a 3 s ) only the OME A Qq contains a generic two-mass contribution, since A PS qq,Q emerges only at O(a 3 s ) and contains one internal massless fermion line. One has The combined renormalization relation at third order is given by This yields the generic pole structure for the PS two-mass contributioñ In the MS-scheme one obtains the renormalized expression bỹ

A Qg
Like in the PS case, there are two different contributions to the OME A Qg Of these OMEs only A Qg contains two-flavor contributions starting from O(a 2 s ) In Eq. (3.114) the O(a 2 s ) contribution consists of one-particle reducible diagrams only, see Eq. (3.86). As a consequence the flavor dependence factorizes in the O(a 2 s ) terms. The renormalized MOM-scheme two-loop contribution is obtained by (3.115) The unrenormalized terms are given bŷ The coefficientsã The renormalized expression at 2 loops then reads The renormalized 3-loop OMEs in the MOM-scheme are obtained from the charge-and massrenormalized OMEs by The structure of the unrenormalized OME is more complex than in the NSor PS case. It is given bŷ For the renormalized operator matrix element in the MS scheme we finally obtain, The matrix element A gq,Q contains contributions starting at O(a 2 s ), Diagrams with two different masses, however, contribute only from O(a 3 s ) The renormalization in the MOM-scheme is performed using Applying Eq. (3.126) yields the unrenormalized expressioñ and the renormalized operator matrix element reads (3.128)

A gg,Q
Finally, the matrix element A gg,Q obeys the expansion with two-mass contributions starting at O(a 2 s ), The renormalization formulae in the MOM-scheme read After subtracting all single-mass contributions we obtain the unrenormalized two-flavor contribution at 2 loopŝ and the renormalized expressioñ The O(a 2 s ) contribution consists of one particle reducible contributions only and the coefficients follow from Eq. (3.87) The unrenormalized 3-loop contribution from two masses readŝ The renormalized result in the MS-scheme is given bỹ

Mass renormalization schemes
The heavy quark masses in the MS and on-shell renormalization schemes are related bŷ where m denotes the mass in the MS scheme and m in the OMS scheme. The ratio of these two masses for a quark of mass m 2 in the presence of a second heavy quark of mass m 1 is given by with [48][49][50][51] 10 , with L µ = ln(µ 2 /m 2 ) and x = m 1 /m 2 .
In data analyses one usually fits the MS-mass m, which is free of infrared renormalon ambiguities, unlike the on-shell mass, which grows significantly order-by-order in perturbation theory [51].

Fixed Moments of the Massive Operator Matrix Elements
In Ref. [1] a series of fixed Mellin moments of all massive operator matrix elements at 3-loop order have been calculated in the single mass case by projecting the corresponding integrals onto massive tadpoles and evaluating them using the code MATAD [76]. These moments serve as important reference points for the general N solution. In the following, we will calculate the Mellin moments N = 2, 4, 6 in the case of unequal masses. The number of moments is less than in the equal mass case, where values of N = 10...14 could be reached, which is due to the presence of the second variable η and the performance of the codes Q2e/Exp [73,74], which we are going to use. The full calculation took about one CPU year. We still obtain very useful reference points by this. The Feynman diagrams are generated using the code QGRAF [77]. In order to take into account the local operator insertions, we introduce new additional propagators which either carry an operator insertion or which generate an operator on an attached vertex. In the case of operator insertions on a gauge boson line, this method leads to a double counting of some vertex diagrams which has to be removed. For the calculation of the color algebra of the expressions we used the code Color [78].
After inserting the Feynman rules, cf. Section 8.1 [1], and the projection operators, the momentum integrals take the form Here p denotes the external momentum, p 2 = 0, ∆ is an arbitrary light-like vector ∆ 2 = 0 and q i are linear combinations of the loop momenta k j and the external momentum p. The exponents n i are integer-valued and obey n i = N , while the function f (k 1 . . . k l , p, m 1 , m 2 ) contains the remaining numerator structure and denominators. In Eq. (4.1), we have omitted possible summations over indices on which the exponents n i might depend.
We may represent (4.1) as Since N j=1 ∆ µ j constitutes a completely symmetric tensor only the purely symmetric part of I (l) µ 1 ,...,µ N contributes. We thus symmetrize by shuffling the indices, [79], and normalize it by dividing by the number of terms. For the general integral (4.1) the symmetrized tensor is given by where S is the symmetrization operator given in Eq. (2.19). The result of the original integral (4.1) may then be obtained again by applying the projection operator [1] The pre-factors F (N ) and the combinatorial factors C(i, N ) for odd values of N are given by , (4.6) and read for even values of N . The pre-factors F odd (N ), F even (N ) are chosen such that the projector (4.4) is normalized The integrals with a local operator insertion for fixed values of N are thus represented in terms of tadpole diagrams with a modified numerator structure. The projection operators (4.4) become sizeable for large values of N , which leads to an exponential increase in the computation time.
In the calculation, the projected Feynman integrals are first expanded in the mass ratio η by an expansion in subgraphs [80][81][82][83] using the codes Q2e/Exp [73,74], which also rely on MATAD to evaluate the single-mass tadpole diagrams, using Form and TForm [84].
The pole structure of the unrenormalized OMEs corresponds to the one which was deduced from the renormalization prescription given in Section 3. As a by-product of the present calculation, also the terms in these 3-loop anomalous dimensions for the moments N = 2, 4, 6, ∝ T F are obtained, cf. [85], here in a two-mass calculation.
The moments of the OMEs calculated in the following depend on the logarithms We expand up to remaining terms of The pole terms in the dimensional parameter ε do not contain any power corrections in η.
In the following, we present the moments N = 2, 4 and 6 for the two-flavor contributions to the constant parts of the various operator matrix elements as defined in Eq. (3.1). 11 The flavor non-singlet contribution two-mass contribution is obtained bỹ 11 We have presented a few of these results before in [31,32]. (4.14) The constant two-mass contribution to the OME A PS,(3) Qq is given bỹ   In Table 1 we illustrate the ratio of the constant parts of the unrenormalized 3-loop twomass OMEsã gq for general values of N . Therefore, these ratios can be used as a first estimate for these OMEs in case the two-mass contribution is only known for some moments.  In order to obtain the results shown above, we have expanded the constant parts of the 3-loop unrenormalized OMEs for fixed even integer values of N . This is a valid representation for some but not for all of the OMEs also at general values of N , as is shown in Sections 5 and 6. In case the expansion in η exists, one might try to reconstruct the η-expanded solution from the moments using guessing methods [86], which have been successfully applied in other cases [28,87]. However, many more moments are needed in this case. They cannot be provided using Q2e/Exp [73,74], and usually require at least the analytic solution of part of the integrals and possibly generating function methods [28,88].
gg,Q andã gq,Q toã NS (3) qq,Q as a function of Q 2 and N .

The Non-Singlet and gq-Contributions at general Values of N
All non-singlet diagrams at 3-loop order contain two massive fermion bubbles. One of these may be rendered effectively massless by using the Mellin-Barnes representation [90][91][92][93][94], see Figure 1. This yields similar integrals as in the case with one massive and one massless fermionic line [19].
One may now introduce a Feynman parameter representation, integrate the momenta and perform the Feynman parameter integrals in terms of Euler Beta-functions The remaining contour integral is then of the general form dξ Γ g 1 (ε) + ξ, g 2 (ε) + ξ, g 3 (ε) + ξ, g 4 (ε) − ξ, g 5 (ε) − ξ g 6 (ε) + ξ, g 7 (ε) − ξ η ξ , where the f j and the g j are linear functions. Furthermore, the notation is applied. After closing the contour in (5.2) and collecting the residues a linear combination of generalized hypergeometric 4 F 3 -functions [95] is obtained For the flavor non-singlet (NS) contributions, and for A gq , the arguments of the hypergeometric P F Q -function are completely independent of the Mellin variable N , and each term factorizes into contributions that describe the operator insertions and the generalized hypergeometric functions covering the mass structure of the diagrams. Due to the fact that the parameters of the hypergeometric functions depend on the dimensional regularization parameter ε only, their respective expansion may be performed with the code HypExp 2 [96]. The results of these expansions are then given in terms of the following (poly)logarithmic functions [43,[97][98][99], The pre-factor C j (ε, N ) may contain a sum stemming from the operator insertion on the vertex, see Section 8.1 [1]. This sum is easily evaluated in terms of single harmonic sums using the summation package Sigma [100,101] Applying these methods we calculate the two-mass contributions in the flavor non-singlet cases and for the OME A gq . In the following, we denote byã ij the remainder constant part of the genuine two-mass term of the unrenormalized matrix element, omitting the terms coming from the expansion of the factor (m 1 m 2 /µ 2 ) 3ε/2 in the OMEs for brevity, cf. (3.1), i.e. the terms ∝ L 1 , L 2 , which are given in the remainder part of the OMEsÃ ij instead together with other the terms of this kind. The expressions forã ij depend on η and are symmetric under the interchange Also the OMEsÃ ij are symmetric in interchanging m 1 ↔ m 2 . One may furthermore check, calculatingã ij (N ) for N = 2, 4, 6 and expanding in η = m / c m 2 b < 1 up to O(η 3 ) that the values for the fixed moments in the corresponding parts ofã ij (N ) are obtained. The latter ones do not obey the symmetry interchanging the masses anymore, since we have chosen to expand for η < 1. To obtain the representations in Section 4 L η is given by − ln(η) expanding the expressions in the remainder part of this section or the z-space expressions given in Appendix A. Since the present expressions obey the symmetry (5.6) a choice has to be made.
For the single mass contributions the different OMEs receive a 3-loop correction changing from the on shell mass m to the MS mass expanding the OME in a MS s . For the two-mass contributions at 3-loop order this is not the case for A gg,Q (see also Eqs. (3.117, 3.118, 3.135, 3.136). They are not dealt with in the present paper.

The flavor non-singlet contribution
The general pole structure for the unrenormalized two-mass contribution to the OME A NS qq,Q is given in Eq. (3.103). The only contribution which is not determined by the renormalization prescription is the constant part, for which we obtaiñ Here S a ≡ S a (N ) denote the (nested) harmonic sums [102] S b, a (N ) = N k=1 (sign(b)) k k |b| S a (k), S ∅ = 1, b, a i ∈ Z\{0} . (5.8) The polynomials R i read The two-mass part of the renormalized OME A

(3),NS
qq,Q is given bỹ Both the constant part of the unrenormalized two-mass OME (5.7) and the OME (5.11) vanish for N = 1 due to fermion number conservation for any value of the heavy quark masses. In the Appendix we present the the corresponding z-space expressions for Eqs. (5.7) and (5.11). The analytic continuation of the N -space result may also be obtained by expressing the contributing sums in the asymptotic region |N | → ∞ and using their recurrence relations, cf. [103]. One may derive semi-numeric representations, cf. [104]. The inversion to z-space is then done by a contour integral around the singularities of the problem.  In Figure 2 we show the ratio of the genuine 3-loop 2-mass contributions to A

The transversity contribution
The pole structure of the unrenormalized transversity OME corresponds to the one in Eq. (3.103) after substituting the anomalous dimensions γ NS qq → γ NS,trans qq . The constant contribution is given by !!ã with The two-mass part of the renormalized OME A The corresponding expressions for (5.7,5.11, 5.14, 5.17) in z-space are given in Appendix A.
As before in the equal mass case [21] and for the O(N F T 2 F ) contributions [19], we obtain the O(T 2 F C A,F ) terms of the 3-loop flavor non-singlet contributions to the anomalous dimensions in the vector and transversity case from the single pole terms of the unrenormalized non-singlet OMEs, confirming once again the result in [15], see also [105]. In Figure 3 we show the ratio of the genuine 2-mass contributions to the complete T 2 F 3-loop term for transversity as a function of x and Q 2 . The spikes are due to a zero in the denominator of this ratio. Except for a small region of x around these spikes, the ratio takes values between 1.5 and -0.6. For Q 2 not too low, mostly values between 0 and 0.6 are obtained.

The gq-contribution
The genuine two-mass contributions to the OME A ln 3 (η) and the polynomials The two-mass contribution to the OME is then given bỹ The corresponding z-space expressions are given in Appendix A. Also in this case we obtain as before in Ref. [23] the corresponding contributions to the 3-loop anomalous dimension [16]. In Figure 4 we show the ratio of the genuine 2-mass contribution to the complete T 2 F 3-loop result for A (3) gq,Q for typical values of Q 2 and x. The ratio varies between 0 and 0.5. At higher values of Q 2 , an almost flat behaviour is observed.
6 Scalar A gg,Q diagrams with m 1 = m 2 The factorization into parts depending purely on the Mellin variable N and contributions depending only on the mass ratio η, which has been observed for the non-singlet diagrams, constitutes a very special case. Normally both variables appear in a more intertwined form and more advanced methods are required to perform the calculation. Since the complexity of the mathematical structures contributing to a Feynman diagram depends on the denominator functions and on the form of the operator insertion, we will first study the scalar topologies contributing to the OME A gg,Q in this paper. Due to the nesting between the Mellin variable and the mass ratio, novel η-dependent sums and integrals will emerge. In particular, it turns out that the expansion in η is not possible in general, unlike the case for fixed integer moments. Therefore, the integrals have to be calculated for general values of η.

The Calculation Strategy
As we expect new functions to appear in the results and since the construction of the inverse Mellin transforms for these functions turns out to be a non-trivial task, we opt for an approach where we derive the z-space representation of the respective diagrams first. The N -space representation 13 is then obtained in a final step by using the generating function method, constructing a difference equation and solving it using the package Sigma [100,101]. These representations can be then evaluated at fixed integer moments in N , be expanded in the parameter η and compared to the fixed moments having been calculated using the code Q2e/Exp [73,74].
First we introduce Feynman parameters and perform the momentum integration for one of the closed fermion lines. This leads to an effective propagator, the mass of which we can detach using the Mellin-Barnes representation [90][91][92][93][94] 1 Then we perform the remaining momentum integrals, which leads to an expression where the Feynman parameter integrals are now of the generalized hypergeometric type [95] and the appropriate application of techniques used earlier in Refs. [12,19,46] allow to integrate all Feynman parameter integrals as Beta-functions, of which only one depends on both the Mellin variable N and the Mellin-Barnes variable ξ and is kept in unintegrated form. We obtain a representation of the general form where a k , d k , β and δ are integers, b k , e k , α and γ are integers or half-integers, c k ∈ {−1, 1} and f k ∈ {−2, −1, 1, 2}, with i k=1 c k = j k=1 f k . The dependence on N of the function C(N, m 1 , m 2 , ε) arises from gamma functions that depend on N (and possibly on ε) but not on ξ.
Mellin-Barnes integrals of the form are usually solved by closing the contour either to the left or to the right and by applying Cauchy's theorem If we close the integration contour in (6.3) to the left(right) the residue sum only converges for Z > 1 (Z < 1), respectively. In (6.2) we have which covers both ranges for possible values of η < 1 and η > 1. In the calculations we applied the code MB [108]. We follow the method applied in the equal mass case in Ref. [24,109], split the integration range and remap the individual parts to the domain [0, 1]: A further advantage of this procedure is that the contour integration decouples the η-dependence which now only enters through the T -integration. We follow the well known procedure of deforming the contour integral in order to separate the ascending from the descending poles 14 applying Cauchy's theorem. At this point we are left with only one integral and no overlapping singularities anymore. If necessary, we map T → 1−T in order to have singularities regulated by ε only at T = 0. They appear as ε-poles after applying the following integration-by-parts relation: We may then perform the Laurent series expansion around ε = 0.
In the next step we rewrite the sums obtained using the package Sigma [100,101] 15 . The sums are then expressed in terms of generalized harmonic sums [39,40] at infinity, (6.8) which have to be rewritten in terms of generalized harmonic polylogarithms (GHPLs) [40] at argument x = 1 using HarmonicSums [40,106,107]. These functions are iterated integrals over the following alphabet : In order to process them, we want the remaining integration variable T to only appear in the argument of the HPLs. Because of the emergence of letters with non-linear denominators, we cannot apply the methods used in Ref. [28,88] directly, although extensions of it, as it is described below, should suffice to transform these HPLs. However, due to the relatively simple structure of the letters in Eq. (6.9), there is a way based on applying the shuffle relations, cf. [79], and rescaling the internal integration variables, to rewrite the corresponding iterated integrals in the desired form.
Instead of computing the remaining integrals, we rather aim at transforming them into a Mellin transform from which one can then read off the z-space representation. Next, we absorb rational N -dependent factors into the integral, which appear both in the numerator and denominator. These factors stem from the integration of the Feynman parameters, and are now pulled into the T -integration by performing a partial fraction decomposition and then applying the following partial integration identities repeatedly, , (6.10) Relation (6.10) has to be especially handled with care, as its application may introduce new divergences in each term. This issue is solved by regularizing the remaining integral in (6.10) by a +-type distribution which cancels these additional singularities, as e.g.
We now rewrite the remaining integral as the following Mellin transform: , for g(x, η) < 0, 0 < x < 1, η < 0 . (6.13) Note that the function g is monotonous (cf. Eq. (6.6)) and thus the inverse function g −1 exists. The class of harmonic polylogarithms is not sufficient to perform this step and generalizations are required to allow for quadratic forms in the denominator. One such generalization is given by the cyclotomic harmonic polylogarithms [106]. We use the label {4, i} to denote the following letter where Φ 4 (τ ) = τ 2 + 1 is the fourth cyclotomic polynomial, and dτ indicates that the iteration proceeds over τ . For example, H 0,{4,1} (x) represents the iterated integral More generally, we write In the calculation of some of the diagrams shown in the next subsection, we thus performed simplifications such as where in the last step we removed the η-dependence of the argument by again applying a rescaling of the inner integration variables. At this point, it is desirable to remove the square roots in the arguments of the HPLs and to obtain iterated integrals with the argument x only. In order to obtain this representation, we once again exploit the property that taking the derivative reduces the transcendental weight of a hyperlogarithm and use a method similar to the one given in [28,88]. For example, However, not all the occurring HPLs can be expressed in terms of generalized HPLs of the previous kind and new, root-valued letters have to be introduced. To perform this in a systematic way, we introduce a more general class of iterated integrals as follows: G ({f 1 (τ ), f 2 (τ ), · · · , f n (τ )} , z) = z 0 dτ 1 f 1 (τ 1 )G ({f 2 (τ ), · · · , f n (τ )} , τ 1 ) , (6.20) with the special cases and Here f i (τ ) are real functions, with τ ∈ [0, 1]. At the moment we do not discuss matters like algebraic or structural independence of these quantities, cf. [40,79,106], but rather consider (6.20) as a placeholder. Algebraic and other relations are applied later in the concrete cases appearing. These functions are given in explicit form in Appendix B.
Using these generalized iterated integrals we rewrite the HPLs with root-valued functions in the argument. For example one has In the present computation, similar HPLs up to weight w = 3 had to be transformed. Due to the size of the expressions and the necessity to cancel spurious terms, all relations obeyed by these quantities have to be used. These are • shuffle relations • integration by parts relations, such that only factors with exponents ∈ {1/2, −1} contribute to the different letters • shuffling of single square root terms to the end and performing the integrals e.g.: dx or 1 f (η) dx due to the splitting of the X-integration in Eq. (6.6). We therefore substitute (6.25) As it would have been expected, the integrals f (η) 0 dx G(x) completely cancel up to trivial integrals of the form .

(6.26)
We now use HarmonicSums [40,106,107] to perform the inverse Mellin transform for terms that do not contain any x-integration. They usually stem from integration-by-parts applied in steps (6.7), (6.10) or (6.11). We are left with a z-space representation for our diagram. This representation usually also includes a part proportional to a δ-distribution and one term proportional to a +-distribution. As a last step, we want to generate a N -space representation for our result, for which the last remaining integration has to be performed. This is done with the help of a generating function representation mapping the integral into generalized HPLs and then generating a recurrence relation for the N th coefficient of this result. This procedure is automated within the package HarmonicSums [40,106,107]. The resulting recurrences were solved using the package Sigma [100,101]. The result contains many generalized HPLs at argument x = 1, which stem from the upper integration limit. In case their letters are free of the mass ratio η, they can be evaluated in terms of special constants like π, ln(2), the Catalan number C, ζ 2 and ζ 3 by using standard integration methods or applying the internal integration algorithms of computer algebra packages like Mathematica or Maple. In case these generalized HPLs are not entirely free of η, it is desirable to rewrite them as iterated integrals with argument η in order to obtain algebraic independence and an easier access to series representations. Rewriting these generalized HPLs cannot be done by rescaling integration variables or by just applying the methods of [28,88], since due to the root valued letters the derivative with respect to an inner variable in general does not lead to a weight reduction in this case. There is, however, an extension to the ideas in [28,88]: Taking the derivative with respect to inner variables we observe, that only GHPLs of a lower weight, GHPLs independent of this variable and the original GHPL itself contribute, as e.g.: Therefore, the linear first order differential operator does lead to a weight reduced expression when applied to the GHPL The weight reduced expression can be rewritten with the same method and we have to undo the effect of the differential operator by using the general solution for the linear first order differential Applying this method to the GHPL considered above we obtain For all the GHPLs considered in this section, it is always possible to construct a linear first order differential operator 16 which does yield a weight reduced expression when applied to the corresponding generalized HPL, and all the GHPLs could thus be rewritten in terms of GHPLs with argument η. See Appendix B for a list of relations for the GHPLs.

The results for individual diagrams
In the following, we present the results for all scalar two-mass topologies contributing to A gg,Q both in z-and in N -space. Up to a global pre-factor, all results are expressed as functions of the mass ratio η. We consider only the cases where the operator insertion is located on a line, and not on a vertex, since the latter case can be easily derived from the former. The powers of the propagators are taken to be the highest ones appearing in the corresponding physical diagrams (this means that in all of the diagrams the sum of powers of propagators equals 9). We define the following functions which appear often in the z-space expressions of the diagrams, 16 First order linear differential operators could be used instead of pure differentiation in order to extend the parametric integration method. However, remapping parameters might be a more suitable method to integrate Feynman parameter integrals which are not a priori reducible. Both methods become inapplicable when noniterative integrals appear, as e.g. genuine elliptic integrals and others.
with the polynomials P i (η, N ) Here the factor 1 2 (1 + (−1) N ) comes from the operator insertion Feynman rule. This factor is removed from the z-space results in all diagrams, due to the analytic continuation from the even moments. Although topologically very similar to diagram D 1 , diagrams D 2a and D 2b exhibit a much more involved mathematical structure. As we restrict ourselves to a representation within the class of iterated integrals of argument z, additional root-valued integration kernels had to be introduced. Furthermore, iterated integrals depending on both variables η and z contribute.
The summands of many of these sums diverge for η → 1 due to factors as (1 − η) −j , where j is a summation index which assumes positive integer values. Furthermore, also contributions ∝ (1 − η) −N emerge. Physically the limit η → 1 represents the equal mass case m 1 = m 2 [24] and thus the diagrams are expected to be convergent in this limit. Due to the many individually divergent terms this is highly non-trivial to prove for general values of N . However, evaluating a series of Mellin moments N = 2 . . . 30, yields convergent results for η = 1, which agree with the results given in Ref. [24] previously. This indicates that these apparent divergences are just a relic of this specific representation which has been applied. By induction one may prove that the result is valid at general values of N . The diagrams (D 2a , D 2b ), (D 4a , D 4b ), (D 5a , D 5b ), (D 6a , D 6b ) and (D 8a , D 8b ) have all been computed independently. One notes that as expected the respective z-and Mellin-space results can be translated into each other by interchanging the masses m 1 ↔ m 2 , η → 1/η. Furthermore, the results for the mass-symmetric diagrams D 1 , D 3 and D 7 turn out to be invariant under this interchange, which constitutes a further check of these results.
For all scalar A gg,Q -topologies, series expansions up to O(η 3 ln 3 (η)) for a series of fixed Mellin moments (N = 2, 4, 6) have been computed using the code Q2e/Exp [73,74]. All the general N and general-η results agree with these expansions.

Conclusions
Genuine two-mass contributions to the Wilson coefficients and the transition matrix elements in the VFNS occur at 3-loop order in QCD. We derived the renormalization of these contributions, which extends the single mass case considered earlier in Ref. [1]. Although the new contributions manifest themselves as two-mass contributions in single diagrams carrying local operators, it is possible to assign a diagram to either of the heavy flavor distributions in the VFNS by the quark species carrying the operator. The diagrams arise from separating off the massless Wilson coefficient in the light-cone expansion. Through this, one knows the charge assignment for the corresponding diagram. In this way an asymmetric separation of the otherwise symmetric OMEs under m 1 ↔ m 2 occurs. This only applies to the OMEs A Qg . All other OMEs enter the VFNS in a mass-symmetric way.
In a first step we have calculated a series of moments (N = 2, 4, 6) for all contributing massive OMEs and presented the constant part of the unrenormalized genuine two-mass OME. With current technologies [73,74], the 6th moment required one CPU year of computational time. For a series of OMEs, the solution for general values of the mass ratio η, and at general values of the Mellin variable N , could be derived along with its z-space representation. This is the case for the OMEs A qg,Q . The corresponding expressions depend on harmonic sums, weighted with a (poly)logarithmic dependence on the mass ratio. In these cases we presented also numerical results studying their relative contribution to the complete O(T 2 F )-term of the OMEs A ij in a wide range of x and Q 2 , in order to illustrate the two mass effects compared to the single mass contributions. In all cases these ratios vary between 0 and ∼ 0.5 in part of the kinematic region, exhibiting scaling violations.
We have also calculated all the scalar topologies appearing in the more involved case of the OME A (3) gg,Q . Here, more advanced computation methods were required. The corresponding integrals do not allow an expansion in the mass ratio at general values of N , so we calculated these integrals exactly. In z-space the corresponding integrals could be represented in terms of iterated two-variate and partly root-valued integrals, the G-functions, see also Appendix B. Associated to it, one obtains in Mellin-N space, sum representations containing functions of η in denominators, with a formally divergent behaviour as η → 1. However, since N ∈ N, one obtains convergent representations for each individual integer N in this limit. Also because of this behaviour, the inverse Mellin transform to z-space requires a series of special steps, which we have outlined. It is expected that the corresponding representation in the case of the two-mass contributions to the OME A Qg is even more involved, since already in the equal mass case elliptic integrals and iteration of other letters over them contribute.

A Massive Operator Matrix Elements in z-Space
In the following, we present a series of genuine two-mass contributions in z-space. These are distribution-valued and consist of the three parts A δ ij , A + ij (z) and A reg ij (z). The Mellin convolution of the OMEs with a function f (z) is defined by, cf. e.g. [110],