The Two-mass Contribution to the Three-Loop Gluonic Operator Matrix Element $A_{gg,Q}^{(3)}$

We calculate the two-mass QCD contributions to the massive operator matrix element $A_{gg,Q}$ at $\mathcal{O} (\alpha_s^3)$ in analytic form in Mellin $N$- and $z$-space, maintaining the complete dependence on the heavy quark mass ratio. These terms are important ingredients for the matching relations of the variable flavor number scheme in the presence of two heavy quark flavors, such as charm and bottom. In Mellin $N$-space the result is given in the form of nested harmonic, generalized harmonic, cyclotomic and binomial sums, with arguments depending on the mass ratio. The Mellin inversion of these quantities to $z$-space gives rise to generalized iterated integrals with square root valued letters in the alphabet, depending on the mass ratio as well. Numerical results are presented.


Introduction
The massive operator matrix elements (OMEs) A ij for partonic transitions rule the matching conditions in the variable flavor number scheme. They start to receive contributions from two different massive quarks starting at two-loop order through one-particle reducible graphs and at three-loop order due to irreducible contributions. Since the mass ratio between charm and bottom is not sufficiently small (m 2 c /m 2 b ∼ 0.1 1 ), there is no single heavy quark dominance in the mass region of charm and bottom and one has to account for both mass effects at the same time, cf. also [3].
In the present paper we calculate these two-mass contributions to the OME A gg,Q up to O(α 3 s ). The O(T 2 F ) single mass contributions have been computed in [4] and the O(T 2 F N F ) terms in [5]. The complete single mass OME is given in [6]. Previously, fixed moments of the OME for the Mellin variable N = 2, 4, 6 as a series expansion up to O(η 3 ln 2 (η)) with η = m 2 c /m 2 b < 1, the ratio of the heavy quark masses squared, were calculated in Ref. [3,7] using the package Q2E/EXP [8,9]. There, also all contributing scalar prototype graphs were calculated in z-space, and the results were converted to N -space by an automated Mellin transform. Furthermore, the analytic results for A NS, (3) qq,Q , A NS-TR, (3) qq,Q and A (3) gq,Q both in N -and z-space have been obtained. Recently, also the z-space result in the pure-singlet case A PS, 3 Qq has been calculated in Ref. [10]. Additionally, various three-loop single mass OMEs have been completed, cf. [5,[11][12][13][14][15][16][17][18]. The logarithmic contributions to all OMEs have been computed to 3-loop order in [19]. For A (3) Qg , all contributions which can be expressed via first order factorizable differential or difference equations in N -or z-space have been obtained in Refs. [20][21][22].
We perform the calculation of the two-mass contribution to A (3) gg,Q first in N -space and then use an automated inverse Mellin-transform to arrive at the z-space result. This is a change of paradigm from earlier work on scalar prototype diagrams [3], where first a representation in zspace was derived and the N -space solution was found by a Mellin-transform. We have checked the present method for the scalar diagrams and found agreement with the previous results. Since both the corresponding difference equations in Mellin N -space and the differential equations in z-space for the contributing diagrams are first order factorizable, one can choose to perform the calculation either in N -or z-space without further difficulties, cf. [21]. There is even no need to refer to special bases.
The paper is organized as follows. In Section 2 we recall essentials for the representation of the renormalized OME A gg,Q in the two-mass case, cf. [3]. In Section 3 the fixed moments for the constant part of the two-mass contributions are given in complete form for later comparison 2 . We outline the general steps of the computation in Section 4 and illustrate them in detail for a particular diagram in Section 5. In Section 6 the result of the two-mass contributions to A gg,Q are given both in N -and z-space. In the latter case, we use in part single-valued integral representations in order to still allow for root-valued iterated integrals (appearing as integrands) that have representations in terms of harmonic polylogarithms (HPLs) of more involved arguments. These integrals can be performed in principle, but lead to voluminous functional expressions for parts of which new numerical implementations would have to be developed as two-variable functions, which we tried to avoid. Numerical results are presented comparing the two-mass contribution to the whole T 2 F term for the heavy quark contributions. The conclusions are given in Section 7. Special integrals and functions of the momentum fraction z and the mass ratio squared η appearing in intermediate and the final result of the calculation are listed in the Appendix, also for further use in other two-mass projects. Here we have expressed the appearing iterated integrals in terms of harmonic polylogarithms at more complicated arguments, which is thoroughly possible in these cases. This allows particular fast numerical implementations. We also present the renormalized OMEÃ gg,Q in N and z-space. 2 The renormalized two-mass OMEÃ gg,Q The complete renormalization and factorization procedure for all OMEs up to O(α 3 s ) in the case of a single massive quark has been presented in Ref. [7], while that for the two-mass case has been derived in Ref. [3]. In order to perform a separated treatment for the two-mass effects, one splits the OMEs into parts stemming from graphs with only one massive flavor and one part containing the contributions from both viâ Here we briefly summarize the main results for A gg,Q . In the following we use the mass ratio η given by We abbreviate a series of logarithms by Here µ 2 denotes the renormalization and factorization scales, which we set equal. The physical masses are denoted by m i , while the unrenormalized quantities are denoted bym i . 3 The operator matrix element A gg,Q receives two-mass contributions beginning at two-loop order. At O(α 2 s ) these contributions stem from one-particle reducible contributions only, while from 3-loop order onwards, genuine two-mass effects appear. The unrenormalized OME is given byÂ with a s = α s /(4π) = g 2 s /(4π) 2 and g s the strong coupling constant. The generic pole structure of A gg,Q up to three-loop order is given by [3] A (2) gg,Q Here we used the notation 4γ ij are anomalous dimensions at (l + 1)-loops. The quantities a (2) gg,Q andā (2) gg,Q denote the O(ε 0 ) and O(ε) terms of the two-loop OMEÂ (2) gg,Q , cf. Refs. [23][24][25][26][27], and β 0,Q = − 4 3 T F , (2.8) 1 + O(ε 2 ), (2.10) 12) 4 In Eqs. (3.137, 3.138) of Ref. [3] unfortunately only the shift N F + 1 → N F has been used, which we correct here. 5 In Eqs. (3.137) of Ref. [3] the notationγ ij was used, which does not reproduce the N 2 F term if the anomalous dimension starts at O N 0 F .
Ref. [28], with r 1 = √ η and r 2 = 1 √ η . (2.13) We work in D = 4 + ε space-time dimensions. Furthermore, ζ i = ∞ j=1 1 j i , i ∈ N, i ≥ 2 denotes the Riemann ζ-function, N F denotes the number of massless flavors and C A , C F and T F are the color factors which for a general SU (N ) take the values T F = 1 2 , C A = N and C F = N 2 −1 2N . The H i (x) are simple harmonic polylogarithms [29]. (2.14) In the following the relationδ gg,Q m 2 1 , m 2 2 , µ 2 . (2.17) These expressions already include contributions from one-particle reducible graphs. They can be written by lower order terms of A gg,Q and the gluon vacuum polarization, defined viâ Π ab µν (p 2 ,m 2 1 ,m 2 2 , µ 2 ,â s ) = iδ ab −g µν p 2 + p µ p ν Π (p 2 ,m 2 1 ,m 2 2 , µ 2 ,â s ) , We split the gluon self-energy into parts depending only on one mass and one part stemming from graphs containing both masses, the same way as we did in the case of the massive OMEŝ Up to two-loop order the gluon self-energy does not obtain contributions from graphs with two massesΠ (k) (p 2 ,m 2 1 ,m 2 2 , µ 2 ) = 0 for k ∈ {1, 2} . (2.21) The single-mass on-shell vacuum polarization of the gluon is given by 6 is a constant frequently encountered in massive calculations [7,30] and Li n (z) denotes the polylogarithm [31]. The two-mass contributions toΠ (3) have been given before as power series up to O(η 3 ln 2 (η)) in [3]. The full dependence on η can be implicitly found in [32] and has been independently calculated in Ref. [33]. The quantity is given bŷ With these ingredients we can split the reducible contributions as followŝ in the initial and final state only. When using the projector, which will be introduced in Eq. (4.1), also diagrams with a ghost in the initial and final state contribute toÂ (2) gg . These, however, must not be included in the reducible contributions for the three-loop OME. This statement does also directly apply to the one loop OMEÂ (1) gg , but since no ghost contributions are present here, we can identifyÂ 3 Fixed moments ofÂ (3) gg,Q In Ref. [3] the fixed moments N = 2, 4, 6 of all two-mass OMEs at 3-loop order were presented as series expansions up to O(η 3 L 2 η ). For the constant part ofÂ (3) gg,Q , the irreducible contributions were given. To allow for a direct comparison with the general N results presented later, we list in the following these terms, including the reducible parts. They are given bỹ

Details of the calculation
There are 76 irreducible diagrams contributing to the OMEÃ (3) gg,Q , out of which 6 contain external ghost lines. Since the value of a diagram is not changed by moving the operator insertion to a different gluon line with the same momentum, we are left with 12 topologically different diagrams. We checked these identities for fixed moments N = 2, 4, 6 with the help of Q2E/EXP [8,9]. Half of these diagrams are symmetric under the exchange m 1 ↔ m 2 , while the other half has to be evaluated for both possible mass assignments. One representative of each of the twelve topologies is shown in Figure 1.
The unrenormalized OME is obtained by applying the gluonic or ghost projector to the Green's functions with external gluon or ghost lines, respectively, adding all contributions up, including the one-particle reducible contributions from Eq. (2.28) as well. For the Feynman rules we follow Ref. [34]. Here p denotes the momentum of the external on-shell gluon, ∆ is a light-like D-vector (i.e., p 2 = ∆ 2 = 0) and a, b are the color indices of the external gluons (ghosts). For the ghost diagrams an additional factor of 2 has to be included. Furthermore, special care has to be taken when including the reducible contributions. Here the irreducible Figure 1: The twelve different topologies forÃ (3) gg,Q . Curly lines: gluons; dotted lines: ghosts; thin arrow lines: lighter massive quark; thick arrow lines: heavier massive quark; the symbol ⊗ represents the corresponding local operator insertion, cf. [35].

two-loop contributionÂ
(2),irr gg,Q in Eq. (2.28) has also to be calculated using the projector of Eq. (4.1), excluding the ghost contributions which enter the complete two-loop result.

Computation Strategy
Since the number of diagrams we have to calculate is small and the reduction to master integrals can introduce spurious terms which only cancel in the final result, we aim at computing the diagrams without reduction to master integrals. Moreover, the reduction to master integrals in the two-mass case with local operator insertions requires a substantial computational time. The direct calculation of the Feynman diagrams in N -space will require the treatment of a large amount of terms as well, due to large numerator structures in the gluonic case. The result is obtained directly, without having to calculate two-parameter master integrals, e.g. by solving differential equations. As we will see later, in this way we would obtain very involved expressions, which can be avoided by introducing efficient one-dimensional integral representations, see also Ref. [10]. They can be found most easily working in Mellin N space.
Furthermore, our first goal, in contrast to the treatment in Ref. [10], is to derive first the Nspace solution and to obtain the z-space result via an analytic Mellin-inversion thereafter. This is possible since all occurring difference equations turn out to be first order factorizable, so closed form solutions of these sums can be found using established difference field techniques using the package Sigma [36,37]. In the following paragraph, we outline the basic computational strategy to calculate the diagrams. After that, we give a more detailed description of the calculation of a particular diagram as an example.
The 76 contributing irreducible diagrams have been generated using QGRAF [38], in the version given in Ref. [35] which includes local operator insertions. After identifying the 12 different topologies, we set up dedicated FORM [39] routines to perform the Dirac algebra and traces. The color algebra is done using the FORM program COLOR [40]. For fermionic bubble insertions we use the identity Next, the Feynman parameterization was performed on the full numerator and denominator structure, i.e., we do not cancel structures appearing in the numerator against the denominator. This provides us with a uniform Feynman parameterization for the whole diagram. At last the momentum integrations were performed one after another, starting from loops without the operator insertion. The resulting tensor integrals were reduced to scalar ones according to the rules stated in Appendix A and thus mapped to the basic one-loop integral [34] It is important to perform the integration of the momentum with the operator insertion as the last one. In this way only the additional scalar product p.k can appear, which simplifies the reduction to scalar integrals drastically, since only a single term of the binomial decomposition of (k.∆ + R 0 p.∆) N can contribute to the integral. After these steps we are left with a linear combination of up to 7-fold Feynman parameter integrals, with the general structure Here R 1 and R 2 are simple rational functions of x i and 1 − x i and R 0 is a polynomial in x i stemming from the local operator insertion. In the next step we split the rightmost factor by means of a Mellin-Barnes integral [41-45] where the real part of the integration contour has to be chosen such that the ascending poles are separated from the descending ones. Our next aim is to compute the Feynman parameter integrals. To do this, the operator polynomial R 0 can be decomposed with the help of the binomial theorem This splitting has to be performed as often as necessary to obtain hyperexponential terms in x i and 1 − x i only. In the present case, we had to split the polynomial up to three times.
Attempts to combine the expression into a linear combination of higher transcendental functions in order to keep the additional summations as few as possible have failed, because overlapping divergencies of the Γ-functions appeared, preventing to choose a proper path for the Mellin-Barnes integral. This indicates that these transformations cannot be performed naively after the Mellin-Barnes representation has been applied. Applying these transformations, all Feynman parameter integrals can be expressed by Euler's Beta-functions For example, we encountered the integral After applying the Mellin-Barnes integral and integrating the Feynman parameters we find . (4.10) Note that the summands arising from the binomial decomposition in Eq. (4.6) appear naturally in nested form. We have not yet specified m a or m b to the physical masses, since there are diagrams with both possibilities. In the following we choose to exploit the symmetry of the Mellin-Barnes integral to arrive at two different representations either proportional to (m 2 a /m 2 b ) σ or to (m 2 b /m 2 a ) σ . In this way we can choose m 2 a /m 2 b = η or m 2 b /m 2 a = η and close the contour to the right in both cases. At this point we could have followed earlier approaches by applying the packages MB [46] and MBresolve [47] to resolve the singularity structure of the integrals and expand the final integral in ε. However, the additional dependence on N and up to four summation quantifiers renders the automated finding of a suitable integration contour non-trivial. Therefore, we calculated these integrals by summing up the residues of the ascending poles of the integrand keeping the ε-dependence and are expanding afterwards. In general, residues had to be taken at σ = k, σ = k + ε/2 and σ = k + ε, where k is an integer larger than an integral specific minimum. In the end, each integral is represented by a linear combination of three infinite sums, over which additional binomial sums have to be performed. Nevertheless, we used the packages MB and MBresolve to check our sum representations for fixed values of the Mellin variable N .
The final multi-sum can now be handled by the packages Sigma [36,37], EvaluateMultiSums and SumProduction [48]. Here additionally HarmonicSums [49][50][51] was used for limiting procedures and operations on special functions and numbers. The sum representation of each integral, which can take up to O(4MB), was crushed to a optimal representation using SumProduction. This representation contains constants from taking out points from summation boundaries and multi-sums with large summand structures. These multi-sums were then handled by EvaluateMultiSums, which uses Sigma and HarmonicSums. The results were expressed in terms of nested harmonic-, generalized harmonic-, cyclotomic-and binomial-sums. Furthermore, generalized harmonic-and cyclotomic-sums at infinity contribute. These can be expressed in terms of HPLs depending on η in the argument with the help of HarmonicSums.
Prior to the solution for general values of N , our sum representations also allow to calculate fixed even moments, without expanding in the parameter η, cf. Section 3. They also serve as input values for the general N -solution.  for some nonnegative integers l, and the remaining 1140 sums consist of double sums of similar type. One of these typical triple sums is

Summation techniques
. (4.12) In the following we will present our summation toolbox that enables one to compute the ε-expansions of the arising sums and thus of the desired integrals by summing up all these εexpansions. More precisely, we will utilize summation algorithms that succeed in representing the coefficients of the ε-expansion of sums like (4.12) in terms of hypergeometric products and indefinite nested sums defined over such products that can be defined as follows.
Definition. Let f (N ) be an expression that evaluates at non-negative integers (from a certain point on) to elements of a field 7 K of characteristic 0. Then f (N ) is called an indefinite nested sum (over hypergeometric products) w.r.t. N if it is composed by elements from the rational function field K(N ), the three operations (+, −, ·), hypergeometric products of the form N k=l h(k) with l ∈ N and h(k) being a rational function in k and being free of N , and sums of the form N k=l h(k) with l ∈ N and with h(k) being an indefinite nested (over hypergeometric products) w.r.t. k and being free of N .
Note that this class covers as special cases the harmonic sums [52], the generalized harmonic sums [51,53] cyclotomic harmonic sums [50] or finite nested binomial sums [54]. 8 Here the variable d may depend on √ η. Furthermore, generalized harmonic-and cyclotomic-sums at infinity contribute.
These can be expressed in terms of HPLs depending on η in the argument with the help of HarmonicSums.
In Subsection 4.2.1 the basic summation mechanism of simplification for such definite sums to indefinite nested sums is presented using the packages Sigma [36,37] and EvaluateMultiSums. As it turns out, this general tactic is not sufficient for the explicitly given expressions: the expressions are scattered into too many pieces of sums and in the intermediate calculation steps for the individual sums, clumsy sums arise that cannot be handled properly with our summation tools. Therefore, we will utilize in Subsection 4.2.2 in addition the package SumProduction [48], which merges the input sums to appropriate forms that can be handled with our summation techniques.

Definite summation tools
In the following we present a survey of the crucial summation tools that assist in the calculation of an ε-expansion for the triple sum (4.12). First, we compute the first coefficients of the ε-expansion of the summand and f 0 in terms of such factorials, the harmonic sums S 1 (i − k), S 1 (k), S 1 (−i + j + k), S 1 (N ) and the cyclotomic harmonic sum Next, we apply the summations over each coefficient, and get the ε-expansion of the triple sum: One is now faced with the task of simplifying T −1 and T 0 , both being free of ε. The simpler coefficient T −1 (N ) can now be simplified by the summation machinery of Sigma [36,37] based on difference ring theory [57]. Namely, one transforms from inside to outside the arising objects in (4.16) to the desired indefinite nested sum form. E.g., for the innermost sum h(η, N, j, i) = i k=0 f −1 (η, N, j, i, k) of T −1 (η, N ) we proceed as follows: (1) We compute a linear recurrence of h(η, N, j, i) in i of order 2: with polynomial coefficients a 0 , a 2 , a 2 in η, N, j, i. The right hand side r is given in terms of a linear combination of hypergeometric products depending on η, N, j, i. This machinery is based on the creative telescoping paradigm [58] in the setting of difference rings [36,37,57].
(2) Next, we solve the found recurrence (4.17) in terms of indefinite nested sums [36,37,59]: we find 2 linearly independent solutions of the homogeneous version of the recurrence and one particular solution of the recurrence itself. One of the most complicated indefinite nested sums w.r.t. i is: where (x) k = x(x + 1) . . . (x + k − 1) denotes the Pochhammer symbol.
(3) Finally, we compute the 2 initial values h(η, N, j, l) with l = 0, 1 and combine the solutions of the recurrence such that they match with the given initial values. This yields an alternative expression of the sum h(η, N, j, i) where by construction the occurring objects are indefinite nested w.r.t. i.
For all three steps, 195 seconds are used and we obtain an expression where 5 indefinite nested sums w.r.t. i appear; one of them is (4.18). Now we apply the next summation quantifier j i=0 to this expression and repeat the same machinery: compute a linear recurrence of this new sum w.r.t. the j (which is the summation variable of the final sum), solve the recurrence in terms of indefinite nested sums w.r.t. j and combine the solutions to find an alternative representation of the double sum which now is indefinite nested w.r.t. j. For this calculation step, 1210 seconds are needed and 9 indefinite nested sums w.r.t. j arose in the found representation. Finally, we repeat this once more in 295 seconds and obtain an expression of the single pole term T −1 (η, N ) in terms of 19 indefinite nested sums w.r.t. N . Summarizing, we needed about 1700 seconds to transform the triple sum (4.12) to an expression in terms of indefinite nested sums which turn out to be harmonic sums, generalized harmonic sums and generalized binomial sums, like, e.g., In general, this summation mechanism of recurrence finding and recurrence solving (see the summation steps (1)-(3) from above) has to be applied slightly more carefully: • If there are poles at the summation bounds (coming from the internal summation representations), the summation range has to be be refined and extra terms have to be treated by another round of our definite summation tools.
• If the summand is too large (e.g., more than 100 MB of size or composed by more than 300 indefinite nested sums), it is split into appropriate smaller parts. Then the summation mechanism is applied separately to them, and the results are combined properly before the next summation quantifier is applied.
• In the case of infinite summation quantifiers (see, e.g., the second sum in (4.11)), also limits have to be handled. Here asymptotic expansions are computed using the summation package HarmonicSums [49].
All these steps are skillfully combined within the package EvaluateMultiSums [48] using the difference ring machinery of Sigma and the special function algorithms of HarmonicSums. E.g., executing the command the single pole term T −1 (η, N ) of (4.15) is simplified to an expression in terms of indefinite nested sums as sketched above. Performing all 9122 triple sums in this way (and ignoring double sums) indicates that already the single pole terms for all these sums require more than 180 days of calculation time. Calculating the constant term in this way seems hopeless. Besides, the following intrinsic problem arose: when we tried to calculate the constant term for our triple sum (4.12), we encountered internally definite sums that are not expressible in terms of indefinite nested sums 9 . However, such alien sums can be avoided by merging all these scattered sums as much as possible and treating them in one stroke. This observation will be utilized in the next subsection.

The full tactic
In order to cure the problems mentioned at the end of the last subsection, we proceed by performing the following three steps.
Step I: The arising sums are crushed to an optimal representation using SumProduction. In this way, one only obtains very few master sums that have to be treated.
Step II: These remaining multi-sums are then handled by EvaluateMultiSums, which uses Sigma and HarmonicSums.
Step III: The results of the master sums are combined to get the final result. Since the calculations of the master sums are carried out independently, the found indefinite nested sums between different master sums are not synchronized, i.e., many relations among them exist. Thus all relations among the arising sums are computed with Sigma and the final result is given in terms of indefinite nested sums that are all algebraically independent among each other. As a consequence, most of the arising sums vanish and a rather compact expression remains.
In the following we work out further details for the two most challenging diagrams: diagram 7 and diagram 11b, i.e. with the bubble-mass for the heaviest quark m = m 2 .

Details for Diagram 7
We proceed with the input expression of diagram 7 introduced already in Subsection 4.2.
Step I: Using the package SumProduction [48] the 10262 sums are merged to the following six sums plus extra terms that only depend on hypergeometric products. More precisely, if the Mathematica variable D7 contains the input expression of D7, the SumProduction-command synchronizes the summation ranges, and maps the arising hypergeometric products (Γ-functions, factorials, Pochhammer symbols, powers) to products that are algebraically independent among each other. This computation required 37892 seconds (10.5 hours) and produced an alternative expression of diagram 7 with 156 GB size. Next, we apply in Step II our summation technologies to these 6 sums and compute the εexpansions up to the constant term in terms of indefinite nested sums, and combine in Step III the found coefficients to obtain the complete ε-expansion of diagram 7. Here the arising indefinite nested sums are algebraically independent among each other. To perform these last two steps, we needed 1 hour for the triple pole term, 2 hours for the double pole term, 2 days for the single pole term, and 20 days for the constant term. Together with step (I), this amounts to 23 days of computation time to obtain the desired sum representation of diagram 7.
In the following we give some further details of Steps II and III for the computation of the constant term of the ε-expansion.
Step II: For instance, consider the sum F 1 whose sum requires 35.6 GB memory; after expanding the summand in ε, the constant term uses 47 MB of memory. Then activating the summation machinery from above to the given triple sum, one needs 605563 seconds (7 days) to obtain the constant term of the ε-expansion of F 1 . The result can be given in terms of 280 indefinite nested sums. This information of F 1 and of the other sums F 2 , . . . , F 6 can be also found in Table 1. In summary, the total computation time for the simplification of all 6 sums required 19 days.
Step III: Combining the constant coefficients of all 6 master sums yields an expression using 23 MB memory consisting of 788 sums and products. Finally, we eliminate all algebraic relations In total, we needed 26 hours to rewrite the found expression in terms of only 46 basis sums that are all algebraically independent from each other. The algebraic independence follows by difference ring theory [57]; for a connection to the underlying quasi-shuffle algebras of the arising sums we refer also to [60]. As a consequence, the expression of 788 sums collapsed in the last step to an expression in terms of 46 sums that requires in total only 0.365 MB. Here one of the most complicated sums is

Details for Diagram 11b
After carrying out the transformations of Subsection 4.1, diagram 11b can be represented by an expression in terms of 14865 sums that requires in total 95 MB of memory. More precisely, the expression consists of 150 single sums, 1000 double sums, 12160 triple sums and 1555 quadruple sums.
Step I: We utilize first the package SumProduction and crunch in 8640 seconds the expression to an expression of 377 MB size consisting only of 8 sums where the summation ranges are given in the first column in Table 2. Next, we calculate the ε-expansions for each of the 8 sums and combine the results to get the ε-expansion of diagram 11b itself. For the triple pole term, this amounts to 89 minutes, for the double pole term to 19 hours, for the single pole term to 6.9 days, and for the constant term to 77.7 days. In the following some extra information is given for the calculation of the constant term.
Step II: In Table 2 further details are given for the treatment of the 8 multi-sums. E.g., for the quadruple sum involving one infinite sum (first row) the input summand uses 17.7 MB of memory. After its expansion the constant term requires 266 MB of memory and the total computation time to transform this sum in terms of indefinite nested sums requires 2.1 days (using the package EvaluateMultiSums utilizing Sigma and HarmonicSums). Let us be even a bit more specific: carrying out the infinite sum requires 32.1 hours and leads to an expression of 8.6 MB size in terms of 41 indefinite nested sums. Carrying out the next sum quantifier sum size of sum summand size of time of number of (with ε) constant term calculation indefinite sums  Dealing with the summation sign i 4 −2 i 3 =0 produces in 18.9 hours an expression of size 4.1 MB in terms of 222 indefinite nested sums, and finally, processing the last summation quantifier N −3 produces in 18.6 hours the final result of 15.2 MB size in terms of 1188 indefinite nested sums 10 .
Step III: Combining the constant coefficients of all 8 master sums yields an expression using 154 MB of memory consisting of 4110 sums and products. Finally, the elimination of all algebraic relations among the arising sums needed 32.5 days and yield a compact expression in terms of 74 sums/products that requires in total only 8.3 MB of memory.
We remark that during the calculations of diagram 7 and in particular of diagram 11b the hardest calculations arose that the summation packages Sigma and EvaluateMultiSum have faced so far. Various sub-algorithms and sub-routines had to be improved and optimized in order to compute recurrences for such gigantic summands and to solve the found recurrences efficiently in terms of indefinite nested sums and products. In particular, the elimination of algebraic relations among the arising sums where pushed to the limit: the underling tower of difference rings contained up to 1500 extension variables (so-called RΠΣ-extensions [57]) and the underlying algorithms were heavily optimized to work with them efficiently.

Deriving the z-space solution
The general method to go from the N -space to the z-space is elaborated in [61,62]. The main idea is to find a recurrence in N for the quantity under consideration and from that to derive a differential equation for the solution in z-space which can be afterwards solved. In the frame of the current project, we used an improved version of the method presented in [61,62], which we will sketch in the following. A detailed description of this enhanced method will be given in [63].
For a given nested sum of the form we look for a representation in the form Note that for the sums under consideration the v j from (4.20) can be read off from F (N ), and that it is straightforward to find recurrences After deriving such recurrences R j , we can use the algorithm from [61] to derive a differential equations D j (with differentiation in z) for the inverse Mellin transforms of v −N jF j (N ). The f i,j (z) are precisely the solutions of the differential equations D j . Hence, after solving the differential equations it remains to fix the d i,j by checking a sufficient amount of initial values.
This method is implemented in HarmonicSums and with the help of the HarmonicSumscommand we find for instance: (4.23) Here G denotes an iterated integral defined in (5.20). In total we computed the inverse Mellin transforms for about 50 sums using this method, which took around 2 hours on a standard desktop PC. The following sum is one of the most complicated sums we had to consider: Note that similar to the final representation of Ref.
[10], we do not include all polynomial prefactors in N but leave these to be included by a Mellin convolution. In this way the inner generalized iterated integrals can be evaluated as HPLs with involved arguments.

An Explanatory Example
In this section we want to illustrate the computational steps in more detail considering diagram 2 of Figure 1. The small numerator structure of this diagram allows to present the calculation in detail. Since here the η and N structures do not factorize, they give rise to more involved structures compared to the single mass case. After inserting the Feynman rules, applying the gluonic projector, performing the Diracalgebra and combining the denominators via Feynman parameters, one obtains where A(B) represent different mass assignments. We normalize the functions J i according to In the following we use the short hand notation

The N -space solution
The functions J 1 to J 3 are given by the following functions The contour integrals are evaluated by taking residues at the ascending poles, which are added up. One obtains where T i,1 follows from the residue at σ = k, T i,2 from the residue at σ = ε + k and T i,3 from the residue at σ = ε/2 + 2 + k. The explicit expressions read Here we applied Legendre's duplication and Euler's reflection formula to the Γ-ratios. The expressions for J B i look similar. In the following we concentrate on the calculation of D A 2 . It is worth mentioning, however, that care is needed at taking the residues for the other mass assignment. Here structures like develop residues at isolated boundary points, i.e., in this example the residues at σ = ε, 1 + ε have to be treated differently than the ones at σ = 2 + n + ε + k with k ∈ N. Therefore, the final representation for D B 2 does not only contain sums but also terms from separately taken residues. In Mellin N -space we use harmonic sums [52] and generalized harmonic sums [51,53] to represent the result. In z-space the corresponding functions are harmonic polylogarithms H a (z) [29] and generalized iterated integrals, G[{ b}, z] over alphabets of the kind discussed in [54], which we find algorithmically [51,54], and special values thereof. The sum representation, moreover, also contains harmonic polylogarithms of the mass ratio √ η.
Here the functions g i , h are arbitrary functions for which the respective integral (5.20) exists. The full expression for D A 2 can now be handled with SumProduction, EvaluateMultiSums, Sigma and HarmonicSums. For the complete diagram we obtain with the polynomials The diagram explicitly fulfills the symmetry We calculated all diagrams which differ for the different mass assignments separately and checked that the symmetry relation holds. For mass symmetric diagrams, we checked the independence of the mass assignment explicitly.

The z-space solution
The result in z-space for diagram 2, split into the usual contributions, reads: Here we use the notations n = N − 1, Furthermore, we will set µ = m 1 for brevity to shorten the expression. The logarithmic dependence on the mass can be easily restored by using the full N -space result and will be entirely given in terms of HPLs. One obtains with the polynomials The functions F k are given by H 0 (η) + H 0 (y) + H 1 (y) H 0 (y) + H 1 (y) The functions G i and K i are given in the appendix. The additional polynomials read We note that although single sums have a different support other than z ∈ 0, 1 , for example

Results
The renormalized 2-and 3-loop OMEsÃ gg,Q (2.16, 2.17) can be obtained from the different contributions to the renormalized masses, the expansion coefficients of the β-function and anomalous dimensions, together with the constant part of the unrenormalized 3-loop OMEã (3) gg,Q in the twomass case. In the following we present this function both in N -and z-space and will give the corresponding results forÃ 2(3) gg,Q in Appendix B.

Transformation to the MS scheme
Since there is a finite two-mass contributionÃ (2) gg,Q , which depends on both heavy quark masses, at 3-loop order the OMEÃ (3) gg,Q differs if calculated in the on-mass shell scheme (OMS) or the MS scheme. The change from the OMSto the MS-mass is best performed on the complete, renormalized OME, which we present in the following.
The relation between the OMS-mass m i and the MS-mass m i for, i = 1, 2 reads with the coefficients, cf. Section 2,

Numerical results
In Figure 2  The contribution of the two-mass term to the whole T 2 F -contribution is significant. At lower values of µ 2 the ratio in Figure 2 shows a profile varying with the momentum fraction z. It flattens at large µ 2 due to the dominating logarithms and reaches values of O(0.4) at µ 2 1000 GeV 2 .

Conclusions
We have calculated the two-mass 3-loop contributions to the massive OME A gg,Q in analytic form both in Mellin N -and z-space for a general mass ratio η. The OME contributes to the two-mass variable flavor number scheme. The close values of the charm and bottom quark masses make it necessary to use this extended scheme. The relative contribution of the two-mass contributions to the whole massive T 2 F contributions of A gg,Q are significant and they amount to values of O(0.4) for µ 2 1000 GeV 2 .
The OME has been first calculated in N -space by direct integration of the contributing Feynman integrals, which made one Mellin-Barnes representation necessary. The problem was thus turned into a nested summation problem, in which the mass ratio η appeared as fixed parameter in the ground field. The corresponding sums could be calculated using the packages Sigma, EvaluateMultiSums and SumProduction, being the largest and most demanding computation we have ever performed as a summation project. For the infinite sums the limit N → ∞ was performed using procedures of the package HarmonicSums. The overall computational time in the summation part amounted to four to five months, including runs needed for code optimization. The N -space result contains harmonic sums, generalized harmonic sums due to the η dependence, and (inverse) binomial extensions thereof. The Mellin variable N also occurs as exponent in η-ratios. We proved analytically that the evanescent poles at N = 1/2 and N = 3/2 vanish. The package HarmonicSums provides algorithms to calculate the inverse Mellin transform of the N -space expressions, which are needed for a series of phenomenological and experimental applications. This is the case because not all parameterizations of parton densities have a simple Mellin space representation, even not at the starting scale Q 2 0 , cf. [64]. The z-space representation can finally be given in terms of general iterated integrals over root-valued letters, also containing the parameter η. These can be reduced to (poly)logarithms of involved arguments, up to one integral in some cases. We were choosing this representation to obtain a fast numerical implementation. The corresponding integrals can in principle be performed within the G-iterated integrals. However, corresponding fast numerical implementations would have still to be worked out for part of these functions. We have checked that our general N -results and those in z-space are in accordance with the moments we have calculated before for N = 2, 4, 6 using different techniques.
With this contribution only one further two-mass OME, A Qg , has to be calculated to complete all two-mass quantities of the VFNS to 3-loop order. At 2-loop order the study of the two-mass VFNS has already been performed in Ref. [65].
During the calculation we obtained a series of analytic integrals, which are listed in the appendix. They are of use in further 3-loop two-mass calculations. One more result of the present analytic calculation is that special numbers, appearing in intermediary steps, and which are not multiple zeta values, cancel in the final result. This is as well the case for one singular Mellin transform due to the behaviour ∝ N 2 , which cancels between different Feynman diagrams.
(A.5) Furthermore, integrals in which the local operator insertion contributes are calculated using Other terms vanish, since they turn out to be ∝ ∆.∆ = 0. gg,Q are given in N -and z-space by :

B The OMEs A
with the polynomials The corresponding contributions in z-space, again separated into δ-, + and regular components, read Also special constants contribute, e.g.: The hypergeometric 4 F 3 -constant can be expressed by (C.8) Similar expressions of this kind often appear in iterated integrals over root-valued alphabets, cf. also [66]. There are many other, partly lengthy, relations more, which we all had to reduce to a suitable level to prove that the evanescent poles at N = 1/2 and N = 3/2 vanish. D Representation of the functions G l and K l Before the absorption of a few rational pre-factors in N , all emerging integrals first written in G-functions can be expressed in terms of polylogarithms at algebraic arguments in z and η.
In cases it leads to simplifications, we also use arcus-and area-functions instead of logarithms, which belong to the harmonic (poly)logarithms of complex-valued argument. The different functions G l ≡ G l (z, η) and constants K l = K l (η) are given by Furthermore, the functions K l (η) ≡ K l contribute. For the more complicated among them we first obtained a longer representation, which finally could be reduced. In these cases we present both representations, since they contain relations between polylogarithms. Structures like this are particularly obtained by integrating using Mathematica. The comparison of both these cases my be helpful in other calculations to obtain more compact results.