The two-mass contribution to the three-loop pure singlet operator matrix element

We present the two-mass QCD contributions to the pure singlet operator matrix element at three loop order in x-space. These terms are relevant for calculating the structure function $F_2(x,Q^2)$ at $O(\alpha_s^3)$ as well as for the matching relations in the variable flavor number scheme and the heavy quark distribution functions at the same order. The result for the operator matrix element is given in terms of generalized iterated integrals that include square root letters in the alphabet, depending also on the mass ratio through the main argument. Numerical results are presented.


Introduction
Starting at 2-loop order, massive operator matrix elements (OME) A ij , which are the transition matrix elements in the variable flavor number scheme (VFNS), receive two-mass contributions [1,2]. This also applies to the Wilson coefficients in deeply inelastic scattering. The single mass contributions to the pure singlet (PS) OME, A

PS,(3) Qq
, have been calculated in Ref. [3]. In the present paper we present the corresponding 2-mass contributions. They occur first at 3-loop order. Previously, we have calculated already the fixed moments of this OME for the Mellin variable N = 2, 4, 6, to O(η 3 ln 3 (η), with η < 1 the mass ratio of the heavy quarks squared in [2,4] using the package Q2E [5,6], as well as the two-mass contributions to the 3-loop OMEs in the flavor non-singlet cases A NS, (3) qq,Q and A NS−TR, (3) qq,Q and for A (3) gq,Q , cf. [2], both in N -and in x-space. The 2-mass contributions to A (3) gg,Q are in preparation [7]. Various other 3-loop single mass contributions have also been completed, cf. [8][9][10][11][12][13][14][15][16][17] and all the logarithmic contributions are known to this order [18]. Furthermore, for the OME A (3) Qg , all diagrams consisting of contributions that can be obtained in terms of first order factorizing differential or difference equations have been calculated [19][20][21]. For all OMEs a series of moments has been calculated in the single mass case in Ref. [22] using the code MATAD [23]. From the single pole terms of all the OMEs, we have derived the contributing 3-loop anomalous dimensions, cf. e.g. [3,24].
We perform the calculation of the 2-mass part of the OME A

PS,(3) Qq
mainly in x-space 1 , using only some elements of the N -space formalism. In the presnt case one obtains first order factorizable expressions in x-space, but not in N -space. This implies that the N -space solution cannot be given by sum and product expressions only.
The paper is organized as follows. In Section 2 we present the renormalized pure-singlet OME in the 2-mass case. Details of the calculation are given in Section 3. In Section 4 the result of the calculation is presented and numerical results are given in Section 5. Section 6 contains the conclusions. Details on new iterated integrals emerging in the representation, the G-functions, and a series of fixed moments as a function of the mass ratio of the two heavy quarks for the unrenormalized OME are given in the Appendix. 2 The renormalized 2-mass pure singlet OME The generic pole structure for the PS three-loop two-mass contribution is given by [2]  where we used the short hand notation 2 3) The tilde inÂ (3),PS Qq indicates that we are considering only the genuine two-mass contributions, and the double hat is used to denote a completely unrenormalized OME. Here the γ (l) ij 's are anomalous dimensions at l + 1 loops, β 0,Q = − 4 3 T F , and where m 1 and m 2 are the masses of the heavy quarks, and µ is the renormalization scale. Our goal is to compute the O(ε 0 ) termã (3),PS Qq (m 2 1 , m 2 2 , µ 2 ). In the MS-scheme, renormalizing the heavy masses on-shell, the renormalized expression is given byÃ The transition relations for the renormalization of the heavy quarks in the MS-scheme is given in [3], Eq. (5.100), but it only applies to the equal mass case since for the unequal mass case the first contributions emerge at 3-loop order. In Eqs. (2.1) and (2.5), a (2),PS Qq and a (2) gq represent the O(ε 0 ) terms of the two-loop OMEsÂ (2),PS Qq andÂ (2) gq , respectively, while a (2),PS Qq and a (2) gq represent the corresponding O(ε) terms, cf. Refs. [26][27][28][29][30]. Here and in what follows, ζ k = ζ(k), k ∈ N, k ≥ 2 denotes the Riemann ζ-function.
3 Details of the calculation 3

.1 The basic formalism
There are sixteen irreducible diagrams contributing toÃ , which are shown in Figure 1. The unrenormalized operator matrix element is obtained by adding all the diagrams and applying the quarkonic projector P q to the corresponding Green functionĜ ij Q , where p is the momentum of the on-shell external massless quark (p 2 = 0), ∆ is a light-like D-vector, with D = 4 + ε, i and j are the color indices of each external leg, and N c is the number of colors. The diagrams, D 1 , . . . , D 16 , are calculated directly within dimensional regularization in D dimensions. These diagrams have the following structure, . The dashed arrow line represents the external massless quarks, while the thick solid arrow line represents a quark of mass m 1 , and the thin arrow line a quark of mass m 2 . We assume m 1 > m 2 .
where N is the Mellin variable appearing in the Feynman rules for the operator insertions, cf. [2,22], m is the mass of the heavy quark where the operator insertion is not sitting, and with m 2 < m 1 , i.e. η < 1.
In previous publications where single-mass OMEs have been computed [3, 8-10, 12, 15, 16, 18, 20], and in our recently published two-mass calculations [2], we have always given the results both in N -and x-space. In the case of the two-mass pure singlet OME, finding a general N -space result at three loops, turns out to be rather cumbersome. We will, therefore, present our result only in x-space, see Eq. (3.3), which is anyway all we need in order to obtain the corresponding contribution to the structure function F 2 (x, Q 2 ) for large values of Q 2 , as well as the contribution to the variable flavor number scheme. In most of the applications one finally works in x-space.
The diagrams on the first row of Figure 1 all give the same result, and the diagrams on the second row are related to the diagrams on the first row by the exchange m 1 ↔ m 2 , i.e., (3.5) D i (m 1 , m 2 , N ) = D 1 (m 2 , m 1 , N ) for i = 5, 6, 7, 8. (3.6) All of these diagrams vanish for odd values of N .
A relation similar to (3.6) holds between the diagrams on the third row of Figure 1 and those on the fourth row. Furthermore, the diagrams on these two last rows which differ only in the direction of a fermion arrow are related by a factor of (−1) N relative to each other. Specifically one has (3.10) We can therefore write the whole unrenormalized pure singlet operator matrix element solely in terms of diagrams 1 and 9:  effectively massless by using a Mellin-Barnes integral [31][32][33][34][35] I µν,ab , (3.12) where µ and ν are the respective Lorentz indices of the external legs, a and b the color indices, k is the external momentum, m is the mass of the fermion, which can be either m 1 or m 2 , 4πα s is the strong coupling constant, and T F = 1/2 in SU (N c ), with N c the number of colors.
Diagram 1 can then be calculated using the following expression for the massive fermion bubble containing the vertex operator insertion (Figure 2b 1 ) : Likewise, diagram 9 can be calculated using the following expression for the bubble containing the operator insertion on the fermion line (Figure 2b 2 ): (3.14) After inserting these expressions, we introduced the corresponding projector for a quarkonic OME given in Eq. (3.1), and performed the Dirac matrix algebra and the trace arising in the numerator of the diagrams using the program FORM [36]. This leads to a linear combination of integrals, the denominators of which were then combined using Feynman parameters. The momentum integrals were then performed, and we obtained expressions where one of the Feynman parameter integrals is already in the form of a Mellin transform. Therefore, we left this Feynman parameter unintegrated, and integrated the remaining ones, after which we obtained the following expression for diagram 1, Here B 1 (ξ) is the following contour integral, For diagram 9, we get We want to compute the diagrams, and therefore the integrals (3.18, 3.28-3.30), as an expansion in ε up to O(ε 0 ). Notice that there is always a factor consisting on a ratio of Γ-functions depending on N in these integrals, and there are also additional factors of N in Eq. (3.19). After the aforementioned ε expansion is performed, some of the integrals will be left with a factor of the form 1 N + l , with l ∈ {−1, 0, 1}. (3.31) In order to get our results really in terms of a Mellin transform, these factors need to be absorbed into the integral in x, which can be accomplished using integration by parts, We will see later how this is manifested in the final result.

Computation of the contour integrals
We have seen that the diagrams where the operator insertion lies on the heaviest quark are related to the ones where the insertion lies on the lightest internal quark by the exchange m 1 ↔ m 2 , which means that in order to go from the former diagrams to the latter, we also need to do the change η → 1/η. Therefore, while in the case of diagrams D 1 , . . . , D 4 and D 9 , . . . , D 12 , we need to evaluate the contour integrals (3.18, 3.28-3.30) with , (3.34) in the case of diagrams D 5 , . . . , D 8 and D 13 , . . . , D 16 , we have instead As we will see, the contour integrals arising in the case of Eq. (3.34) will require a different treatment to those in the case of Eq. (3.35). The integrals (3.18, 3.28-3.30) can be computed with the help of the Mathematica package MB [37], together with the add-on package MBresolve [38]. These packages allow us to resolve the singularity structure in ε of these integrals by taking residues in σ, after which we can perform the expansion in ε. This leads to expressions of the form where the functions B i (ξ) are the sum of residues. For example, in the case of B 1 (ξ), we must take residues at σ = ε/2 and σ = 3ε/2, which leads to Once the residues are taken, it is possible to expand the original integrals in ε. This leads to the contour integrals that we denote by B In order to calculate the contour integrals (3.38-3.41), we need to close the contour either to the right or to the left, depending on the convergence of the resulting sum of residues. The integrand always contains a Γ-function in the denominator that can be rewritten using Legendre's duplication formula, which means that if we close the contour to the right, the sum of residues will be convergent if and only if ξ < 4, while if we close the contour to the left, the sum of residues will be convergent if and only if ξ > 4. Therefore, in the case where ξ = (ηx(1 − x)) −1 , we just need to close the contour to the left, since η < 1, and therefore ξ ≥ 4/η > 4. The calculation in the case where ξ = η/(x(1 − x)), on the other hand, is a bit more complicated, since in this case we need to split the integration range in x into the different regions where ξ can be bigger or smaller than 4, where Therefore, in this case, we will need to close the contour to the right when η − < x < η + , and to the left when 0 < x < η − or η + < x < 1.
Closing the contour to the right and summing residues, we obtain while closing to the left leads to (3.53) Here S a ≡ S a (N ) denotes the (nested) harmonic sums [39] S b, a (N ) = N k=1 (sign(b)) k k |b| S a (k), S ∅ = 1, b, a i ∈ Z\{0} . (3.54) In the above expressions ratios of Γ-functions are related to special binomial coefficients, like (3.55) All of the above sums can be performed using the Mathematica packages Sigma [40,41], Har-monicSums [42][43][44], EvaluateMultiSums and SumProduction [45]. The results are expressed in terms of generalized iterated integrals.
In principle, the letters in the alphabet of these iterated integrals (i.e., the functions f k (τ )) can be any function (or distribution), for which the iterated integral exists. In the particular case where the letters are restricted to 1 τ , 1 1−τ and 1 1+τ , these integrals correspond to the harmonic polylogarithms [46], which are defined by We will see that not only harmonic polylogarithms appear in our final result, but also iterated integrals with square roots in the letters. 4

Numerical Results
We compare the pure singlet 2-mass contributions to the complete O(T 2 F C A,F ) term as a function of x and µ 2 in Figure 3. Typical virtualities are µ 2 ∈ [30, 1000] GeV 2 . The ratio of the 2-mass contributions to the complete term of O(T 2 F C A,F ) grows in this region from slightly negative contributions to ∼ 0.36 for very large virtualities in most of the x-range. The behaviour of the ratio is widely flat in x, rising at very large x.

Conclusions
We have calculated the two-mass 3-loop contributions to the massive OME A PS,(3) Qq in analytic form in x-space for a general mass-ratio η. It contributes to the 3-loop VFNS and to the massive 3-loop corrections of the deep-inelastic structure function F 2 (x, Q 2 ) in the region m 2 Q 2 . As a function of x its relative contribution to the O(T 2 F C A,F ) terms of the whole matrix element A PS,(3) Qq behaves widely flat and grows with the scale µ 2 up to about ∼ 0.36. We used Mellin-Barnes techniques to obtain the x-space result by factoring out the Ndependence in terms of the kernel x N , and used integration by parts to absorb the N -dependent polynomial pre-factors. The result can be written as single limited integrals within the range x ∈ [0, 1] over iterated integrals containing also square-root valued letters. These integrals can be turned into polylogarithms of involved root-valued arguments, depending on η. The even Mellin moments of the OME exhibit a growing number of polynomial terms in η with growing values of N . Due to this structural property and the arbitrariness of η, which enters the ground field, the method of arbitrarily large moments [51] cannot be used to find the result in the present case.

B Fixed Moments
For fixed values of N = 2k, k ∈ N\{0}, we find the following moments: