The unpolarized two-loop massive pure singlet Wilson coefficients for deep-inelastic scattering

We calculate the massive two--loop pure singlet Wilson coefficients for heavy quark production in the unpolarized case analytically in the whole kinematic region and derive the threshold and asymptotic expansions. We also recalculate the corresponding massless two--loop Wilson coefficients. The complete expressions contain iterated integrals with elliptic letters. The contributing alphabets enlarge the Kummer-Poincar\'e letters by a series of square-root valued letters. A new class of iterated integrals, the Kummer-elliptic integrals, are introduced. For the structure functions $F_2$ and $F_L$ we also derive improved asymptotic representations adding power corrections. Numerical results are presented.


Introduction
The complete massive two-loop Wilson coefficients for deep-inelastic scattering corresponding to the structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ) were only available in numerical form [1][2][3] 1 for a long time. Later the flavor non-singlet Wilson coefficients have been calculated analytically in [5] in the tagged-flavor case and recalculated for the inclusive case [6] to obtain a representation consistent with the associated sum rules.
In the present paper we calculate the massive pure singlet two-loop Wilson coefficients analytically. Due to the corresponding graphs, the formulae are structurally the same for the charm and the bottom contributions. In the numerical illustrations we will concentrate on the charm contributions, considering the first three quarks as massless. The knowledge of the complete analytic expressions allows to derive important limiting cases such as the limit of large virtualities Q 2 m 2 , m being the heavy quark mass, or the threshold expansion in a direct way. In the former case it is possible to derive systematic expansions in m 2 /Q 2 with coefficients represented in terms of harmonic polylogarithms, while the complete result depends on much more general functions. Harmonic polylogarithms can be easily calculated numerically [7][8][9]. Furthermore, they can be directly transformed to Mellin space [10,11]. It has been observed numerically in Ref. [5] that the limit of large virtualities is approached beyond some process-dependent scale Q 2 0 . The Wilson coefficient in this limit can be calculated with the help of massive operator matrix elements (OMEs) and massless Wilson coefficients, cf. [5]. It is important to prove this analytically. At three-loop order the massive Wilson coefficients are only known in the asymptotic region [12][13][14][15][16][17][18][19][20][21][22][23]. We also recalculate the corresponding massless two-loop Wilson coefficients given in [24][25][26][27][28][29][30][31] before and compare to these results.
The analytic calculation of the massive pure singlet Wilson coefficient has been envisaged by W.L. van Neerven and one of the authors (J.B.) 20 years ago, after the non-singlet contribution had been obtained in [5]. In retrospect, however, adequate mathematical techniques to perform this task have only become available very recently. This includes the elimination of all functional relations in the final result and techniques to obtain a compact representation. The massive Wilson coefficient is given by a four-fold non-trivial phase space integral. Three of the integrals can be carried out using standard techniques. The integrand of the last integral is obtained as a polynomial of rational terms, logarithms and polylogarithms [32,33] with an involved argument structure. Therefore, the last integral is performed after determining the contributing irreducible structure of letters of the contributing iterated integrals, using the techniques described in [34,35]. The Wilson coefficient can finally be obtained as a d'Alembertian integral over a finite alphabet. The analytic results allow to perform expansions in m 2 /Q 2 including power corrections, which is of particular importance for the structure function F L (x, Q 2 ). Here the corresponding expansion coefficients are then harmonic polylogarithms. Such a representation is easily envisaged for the two-loop non-singlet Wilson coefficients given in [5,6], since there the whole Wilson coefficient depends at most on classical polylogarithms.
The paper is organized as follows. In Section 2 we first illustrate the asymptotic factorization using the example of the O(α s ) calculation. The corresponding scattering cross sections will be used in the two-loop massless and massive calculation later. In Section 3 the massless two-loop pure singlet Wilson coefficients are calculated. The mathematical method used to prepare for the last analytic integral in the massive case is described in Section 4 and in Section 5 we present the analytic results for the massive Wilson coefficients. The asymptotic and threshold expansions are derived in Section 6 and numerical results are presented in Section 7. Section 8 contains the conclusions. Some technical aspects of the calculation are given in the Appendix.

Asymptotic cross section factorization
The massive Wilson coefficients are calculated by factorizing the massless initial states (quarks and gluons). In the unpolarized case and for longitudinal polarization the factorization is longitudinal, i.e. by setting p = zP, z ∈ [0, 1]. Here P denotes the incoming hadron momentum and p the quark momentum. In the transversal polarized case one has to use the covariant parton model [48], see [49][50][51][52]. As an illustrative example we consider the unpolarized one-loop heavy flavor contribution to deep-inelastic scattering [53][54][55][56][57]. As for all the massive Wilson coefficients, it can be written in three parts: the massive operator matrix element, the massless Wilson coefficient and a remainder part. The last one vanishes in the limit Q 2 /m 2 → ∞ in the case of asymptotic factorization. A simple prediction on the structure of this term is not easily possible, but usually requires the calculation of the whole process followed by the expansion in m 2 /Q 2 . This term depends on the structure of the phase space and it is a process-dependent quantity. In Figure 1 the contributing Feynman diagrams are shown. The massive Wilson coefficients have the following series representation where i denotes the incoming parton and 2(L) refer to the associated structure functions and a s ≡ a s (µ R ) = g 2 s /(4π) 2 denotes the strong coupling constant at the renormalization scale µ R . We work in d = 4 + ε space-time dimensions. Since we also need the O(ε) term of the LO result later on, we further define where we dropped the arguments of the coefficient functions for brevity.
Let us consider the leading order contribution for the process γ * + g → QQ as an example, cf. [53][54][55][56][57]. In the following we use the variable The Wilson coefficients H L,g z, h 2,g z, with θ(x) the Heaviside function, a = 1/(1+4m 2 /Q 2 ) and Here we refer to the harmonic polylogarithms [58] defined by and the letters f c are using the definitionf Note that Eqs. (12,13) hold for z ∈ [0, 1]. Here C (1) g,2(L) denote the massless two-loop Wilson coefficients and A (1) Qg the massive one-loop operator matrix element (OME) with external gluons [5,19,36] A (1) The massless one-loop Wilson coefficients read [59][60][61] whereP is a one-loop splitting function [62,63]. 2 It can now be seen that the massive Wilson coefficients can be decomposed in terms of the part obtained at large virtualities Q 2 m 2 , Eqs. (12,13), consisting of massive OMEs and massless Wilson coefficients, and a remainder part vanishing in the limit Q 2 /m 2 → ∞. Whenever this is the case one calls the respective process asymptotically factorizing. The factorization scale µ cancels in the cross sections (12,13) since they are free of collinear singularities. As a peculiarity in this case, the massive OME only contributes to the pure logarithmic term. This, however, is due to its vanishing constant part and is generally not the case. Numerically it is interesting to see from which value of Q 2 0 /m 2 onward the asymptotic representation holds, say at the accuracy of O(2%) or better, cf. [5,6] and Section 7.

The massless Wilson coefficients
The massless pure singlet Wilson coefficients obey the expansion with δ 2 = 1 for C 2 and δ 2 = 0 for C L . Throughout this paper we will identify the factorization scale µ F and the renormalization scale µ R .
The unrenormalized Wilson coefficients F L(2),q are related to the hadronic tensor of deeply inelastic scattering in the partonic sub-system,Ŵ µν , by Here p denotes the incoming parton momentum and q the space-like momentum of the virtual photon with q 2 = −Q 2 .
In the massive case we will also consider the Wilson coefficient as a subsidiary function in order to avoid redundancies in the calculation. Note that this Wilson coefficient does not correspond to the structure function F 1 , cf. [64].
The following expressions will be given in Mellin-N space. They are obtained from the momentum fraction z-space by a Mellin transform The unrenormalized Wilson coefficients F L,g + c (2),PS L,q withâ s the unrenormalized coupling constant, the spherical factor and γ E the Euler-Mascheroni constant. We work in the MS-scheme and set S ε = 1 at the end of the calculation. Here the factors of 1/2 in Eq. (25) emerge since for the splitting into the upper quark-antiquark pair, the quarks are produced correlated. Since the pure singlet contributions start at O(a 2 s ) only, the renormalized Wilson coefficients C (2),PS L, (2) are obtained after mass factorization with In z-space the functions in Eqs. (24,25) read see as well Eqs. (16,17) for µ 2 = Q 2 . The splitting functions are The massless Wilson coefficients C PS,(2) L and C PS,(2) 2 are thus given by We agree with the results given in [30,31] and note a typo in [27], Eq. (13), where the next-tolast term should read (448/27)x 2 . In Appendix A.1 we present details of the calculation in the massless case. The massless two-loop pure singlet contribution to the structure functions F 2(L) for pure virtual photon exchange is given by where µ denotes the factorization scale, Q H = 2/3 for charm and Q H = −1/3 for bottom, and denotes the quark singlet distribution for three light quarks.

Systematic integration in the massive case
We will express the scattering cross sections in terms of a minimal number of special functions.
In the case of single scale quantities, various methods have been worked out in the past to achieve this; for a recent survey see [65]. In the present case, we deal with a two-scale process, since the cross sections depend on z and m 2 /Q 2 in a non-factorizing way. The complete massive Wilson coefficients are represented in terms of four non-trivial integrals. The first three integrations are evaluated in terms of logarithms and polylogarithms at various complex arguments involving square-roots and trigonometric functions. What remains is a one-fold integral with respect to an angular variable ϕ of a function that also depends on the parameters z and β. The overall aim is to write this integral in terms of nested integrals. To this end, we first write its integrand in terms of nested integrals. First, we apply the change of integration variables t = sin(ϕ).
As a result, we get rid of the trigonometric functions in the integrand. In addition, we introduce the quantity which satisfies √ z < k < 1. We use it to express β as Altogether, the integrand is then an expression in terms of z, k, and t as well as logarithms and dilogarithms with arguments expressed in terms of square-roots involving these quantities.
Next, we eliminate redundancies among square-root expressions to express the integrand using only the roots In order to facilitate the conversion of the logarithms and dilogarithms appearing in the integrand to nested integrals, we exploit the argument relations to avoid arguments on branch cuts. After these pre-processing steps, all the following steps for computing the integral are done by our code [66] in Mathematica, which also uses the routine DSolveRational of the package HolonomicFunctions [67]; see [34,68] for the general theory underlying [66]. We also refer to [69] for the simpler case when no singularities are present at the endpoints of integration, which, however, does not apply here.
First, the logarithms and dilogarithms are converted to nested integrals, which is based on repeated differentiation followed by expressing the integrands of these nested integrals in the form developed in (3.16)-(3.19) of [35]. In fact, a generalized version of those forms is used to avoid the necessity of introducing new square-roots in terms of z and k in addition to √ 1 − k 2 above. Then, a normal form of the integrand is computed. This affects all parts of the representation, also those that do not depend on t. For the nested integrals we use the shuffle relations and also for their coefficients we compute normal forms in terms of the logarithms and square-roots.
As a result, we obtain a representation of the integrand as a linear combination of nested integrals evaluated at t whose integrands also depend on z and k. Their coefficients only contain z, k, t, as well as all other logarithms and dilogarithms depending on z and k, do not appear in this representation anymore. Moreover, since both the integrand as a whole and all integrands of the nested integrals in its representation are real, all complex expressions drop out of the coefficients as well and we have a completely real representation. This is ensured since the integrands in (3.16)-(3.19) of [35], and also their generalization used here, were designed so that the corresponding nested integrals all are linearly independent.
Finally, the integral over t from 0 to β is computed as a linear combination of nested integrals evaluated at β, again in normal form. Like before, their integrands also depend on z and k and their coefficients only contain z, k, t, , ln(k +z), and ln(k −z). The following letters contribute in the present case: The set of letters span the Kummer-Poincaré iterated integrals [70] defined as The letter f w 9 can be rewritten into Kummer-Poincaré letters [70], which we, however, avoid here. Some of the above letters contain the elliptic letter as a factor. Therefore, one expects that in iterated integrals the incomplete elliptic integrals of the 1st, 2nd, and 3rd kind cf. [71], are emerging, over which further Kummer-Poincaré letters are iterated. We call iterated integrals of this type Kummer-elliptic integrals. Their alphabet is Note that integrals of depth 1 over the letters f w 1 to f w 12 are (poly)logarithmic, since one may change variables t → √ t, cf. Eqs. (52)(53)(54)(55)). Yet Kummer-elliptic integrals appear in the iterated case. Therefore, iterated integrals of depth 2 formed out of some of these letters will form results containing incomplete elliptic integrals in part. These iterative integrals cannot be reduced to the Kummer-Poincaré iterated integrals for general values of k. As also the incomplete elliptic integrals, they belong to the d'Alembert class, unlike the complete elliptic integrals [71], which also appear in various higher order calculations, cf. e.g. [72], as letters in other iterated integrals.

The massive Wilson coefficients
The unrenormalized two-loop massive pure singlet Wilson coefficients H i,q with i = 1, 2, L, see also Eq. (22), are given in Mellin space by The functions h 1,g are given by Since the two heavy quarks do not induce collinear divergences the mass factorization in the massive case reads Therefore, we find Identifying the renormalization and factorization scale, µ = µ F , we finally obtain Note that in the pure singlet case the coupling constant is not renormalized at two-loop order.

The asymptotic and threshold expansions
The complete expressions calculated in Section 5 allow now to perform the asymptotic expansion for Q 2 m 2 and the threshold expansion for β 1. In the asymptotic limit Q 2 m 2 the massive pure singlet Wilson coefficient have the following representations [5,36] Here the massless Wilson coefficientsC The constant part of the unrenormalized OME a in z-space.
Expanding the fully massive result given in Section 5 in the asymptotic limit Q 2 m 2 and setting µ 2 = Q 2 we find with the polynomials We note that the asymptotic terms are exactly reproduced, cf. [5,12,36], proving the asymptotic factorization in this process. The additional power suppressed terms can be used to obtain fast numerical implementations for the heavy quark Wilson coefficients which are valid for lower values of Q 2 . The reach of this approximations is discussed in Section 7.
The threshold expansion of the Wilson coefficients for β 1 is given by compare to Ref. [16] for H 2,PS 2,q . Next we study the ratios cf. also [5], comparing the full (69,127)   to χ = 10 does not introduce an error larger than 5% in this region. At larger z (here z = 1/2) the asymptotic representation begins to deviate significantly from the full calculation beginning at χ ∼ 1000. However, the Wilson coefficients are very small in this region. As it was already noted earlier [5] the asymptotic representation for H 2,PS L,q is only valid for much higher values of χ. Demanding an agreement of ≤ 2% requires χ > 900 for the small values of z and even higher values for larger z. Similar to the ratio of the full and asymptotic Wilson coefficient we define the ratio whereF (2),PS i,q is the structure function obtained by using the expansion of the respective Wilson coefficient up the desired level. The corresponding results are depicted in Figure 4. We use the parameterization of the parton distribution [73] at NNLO to better compare previous numerical results [16]. We used the LHAPDF interface [74]. Demanding an agreement within ±2% for In Figures 5 we show the complete results for the two-loop pure singlet contributions to F 2 and F L as a function of x for a series of Q 2 -values. At large values of Q 2 the corrections are negative and turn to positive values around Q 2 ∼ 10 GeV 2 . In the small x region the corrections are large and grow with Q 2 . The absolute corrections to F L are smaller in size than those to F 2 .
In Figure 6 we illustrate the ratios Eq. (148) as a function of x for different values of Q 2 for F 2 and F L comparing the asymptotic result to the full result. The corrections behave widely flat in x, turning to lower values in the large x region. For F 2 the ratios are larger than 0.96 for Q 2 ≥ 500 GeV 2 . At Q 2 = 100 GeV 2 , values of ∼ 0.85 are obtained. For lower values of Q 2 the ratio is even smaller.
For F L the corrections are generally larger. At Q 2 = 10 4 GeV 2 one obtains a ratio of 0.96, for Q 2 = 10 3 GeV 2 0.85, and for Q 2 = 500 GeV 2 ∼ 0.75, with even larger deviations from one for lower values of Q 2 .
In Figure 7 we depict the ratio of the full result over the O((m 2 /Q 2 ) 2 ) improved asymptotic results for F 2 and F L as a function of x for a series of Q 2 -values. In the region x < 0.1 the ratios for F 2 are larger than 0.98 for Q 2 > 50 GeV 2 and grow for larger values of x. Stronger deviations are observed for lower Q 2 values. For F L the corrections are larger. In the region x < 0.3 and Q 2 > 100 GeV 2 the ratio is larger than 0.97, while for lower scales Q 2 the deviations are larger. We limited the expansion to terms of ∼ O((m 2 /Q 2 ) 2 ) in the present paper, but higher order terms can be given straigtforwardly. The expanded expressions do also allow direct Mellin transforms and provide a suitable analytic basis for Mellin-space programmes.   [73] .

Conclusions
We have calculated the massless and massive two-loop unpolarized pure singlet Wilson coefficients of deep-inelastic scattering for the structure functions F 2 and F L . In the massless case, we confirmed earlier analytic results in the literature, which can be expressed by harmonic polylogarithms. In the massive case, the Wilson coefficients are calculated analytically for the first time. They are also given in terms of iterative integrals, including now, however, Kummerelliptic integrals. The corresponding alphabets contain also elliptic letters. All integrals can be represented by classical (poly)logarithms with involed arguments with partly one more (elliptic) letter iterated upon. This representation is very well suited to obtain numerical results.
We have studied systematic expansions in the ratio m 2 /Q 2 in the asymptotic region and the velocity parameter β in the threshold region. In the former case the leading asymptotic result has been recovered, known form calculations based on massive OMEs and massless Wilson coefficients, proving asymptotic factorization in the present case. We have obtained a series of power corrections. Here the expansion coefficients are also spanned by harmonic polylogarithms. Retaining these terms extends the validity of the cross sections to lower scales of Q 2 , which is relevant for experimental analyses. In particular, the predictions for the structure function F L (x, Q 2 ) are significantly improved. In general, the Kummer-elliptic integarals, also obeying shuffling relations, span a wide class of iterative integrals which play a role as well in other multi-scale calculations.

A Details of the calculation
Our calculation closely follows classical calculations in the literature, cf. e.g. [61,[76][77][78]. Although these calculations are typically well documented, we encountered subtleties at several points of our calculation. Therefore, we provide a more detailed discussion of our calculation in the massless and massive case in this Appendix. First we will give the parametrization of the phase space we used in the massless and massive case, then we will proceed by explaining the angular integration and give explicit results for the angular integrals in d dimensions. In the end, we will comment on our resolution of the poles in ε and subtleties encountered in the massless case.

A.1 Phase Space Parametrization
The 2 → 2 Process In the 2 → 2 case in Figure 1 we refer to the invariants We will also use the notation β = 1 − 4m 2 /s. In the cms of the outgoing particles, k 1 + k 2 = 0, the scattering angle θ is defined by with The phase space integral is given by The limit m → 0 is easily obtained by setting m = 0 and β = 1. (166) In the limit m → 0, we recover the parameterization given in [61].
In a next step we want to introduce dimensionless variables with support over the unit cube. Here it is advantageous to distinguish between the massless and the massive case. In the massless case, we follow [61] and introduce the new variables The massless three-particle phase space then reads In the massive case the change to the following variables is useful The new parameterization then reads The limit m → 0 is not easily recovered, because of the mass dependent transformation.

A.2 Angular Integrals
The massless case There are four angle dependent denominator structures appearing for the pure singlet process: with Using partial fractioning we can express all angular integrals via We only encounter integrals with k ≤ 0, however, it is possible to find closed form solutions for k ≤ 0 and l ≤ 0 in the massless case. In the following we will list the result for these angular integrals in d-dimensions. l negative: For l = 0 this reduces to k negative: For k = 0 this reduces to Expanding these results around ε = d − 4 dimensions we recover the integrals given in [76].

The massive case
In the massive case the four denominator structures read Therefore, we have to consider the more general angular integral in the following. For l ≥ 0 and arbitrary k (the only case we encounter), we find:

A.3 Regularization
In order to perform the ε-expansion of the functions we use a simple subtraction term for y = 0. However, there is a subtlety hiding in this limit. The hypergeometric functions of interest are all of the argument Inserting the coefficients from Eqs. (176), we see that which means that there is a potential logarithmic singularity for y → 0 in the massless case. This divergence can be made explicit by transforming the 2 F 1 's from argument The new hypergeometric functions have Taylor expansions around y = 0. The only singular behavior can now occur for y → 0. This means that we can resolve the divergences via where we used the notation In the massive case we have which means that this divergence is regulated by the quark mass. The subtraction term (B) can be trivially integrated over y, which will lead to poles in ε. In the massless case the expansion in ε can be performed afterwards and the last integration over z can be carried out. In the massive case there can be additional singularities hiding in the z → 1 limit. Therefore, term (B) has to be regularized accordingly. Term (A) is not singular in the limit y → 0 and can be expanded in ε and then integrated over y and z.

B Contributing expressions due to renormalization
In the following we list some Mellin-convolutions, which occurred in Eqs. (69,70). These are convolutions with leading order splitting functions, using the parameter κ = m 2 /Q 2 .