On the next-to-next-to-leading order QCD corrections to heavy-quark production in deep-inelastic scattering

The contribution of quarks with masses m>>Lambda_QCD is the only part of the structure functions in deep-inelastic scattering (DIS) which is not yet known at the next-to-next-to-leading order (NNLO) of perturbative QCD. We present improved partial NNLO results for the most important structure function F_2(x,Q^2) near the partonic threshold, in the high-energy (small-x) limit and at high scales Q^2>>m^2; and employ these results to construct approximations for the gluon and quark coefficient functions which cover the full kinematic plane. The approximation uncertainties are carefully investigated, and found to be large only at very small values, x smaller about 10^-3, of the Bjorken variable.


Introduction
The production of heavy quarks in deep-inelastic lepton-hadron scattering (DIS) is an important process as it sheds light on a number of interesting issues in the theory of the strong interaction, Quantum Chromodynamics (QCD). First of all, the production mechanism can be described by standard renormalization-group improved perturbation theory for sufficiently large values of the momentum transfer Q 2 between the lepton and the hadronic states. For this purpose it is a great advantage that one can rely on the well-developed framework of the operator product expansion (OPE) which allows to understand essential features of heavy-quark DIS. In addition, the presence of the heavy-quark mass m introduces a second hard scale into the problem, m 2 ≫ Λ 2 QCD , where Λ QCD denotes the QCD scale. As a consequence, at any fixed order in perturbation theory (and neglecting bound-state effects), the heavy-quark mass defines the production threshold for the squared center-of-mass (CM) energy S, e.g., S ≥ 4m 2 for heavy-quark pair production. Moreover the heavy-quark mass acts as a regulator of collinear divergences, giving rise to large logarithms of the ratio Q 2 /m 2 at Q 2 ≫ m 2 , i.e., in the kinematic regime to be considered when matching QCD with n f light quarks and one heavy flavour to a theory with n f +1 light quarks.
Experimentally, heavy-quark production in DIS has been studied in fixed-target experiments and, in particular, at the electron-proton collider HERA. In neutral current reactions, measured with high accuracy at HERA, a considerable part of the inclusive DIS cross section at small Bjorken-x is due to the production of charm quarks, see, e.g., Refs. [1][2][3]. In this kinematic regime heavyquark DIS is dominated by the photon-gluon fusion process γ * g → cc X , and the respective highprecision measurements can provide invaluable information on non-perturbative parameters in the cross sections such as the gluon distribution of the proton and the strong coupling constant α s . Given sufficiently accurate data, they can even facilitate determinations of the heavy-quark masses in a reaction with space-like momentum transfer. With more data analyses from HERA run II still to be finalized, the constraints from heavy-quark DIS on global fits of parton distribution functions (PDFs) and its resulting consequences for collider phenomenology are, perhaps, the most important aspects in the era of the LHC.
It is thus of great importance to provide as accurate theoretical predictions for heavy-quark DIS as possible. Within the standard perturbative approach, this requires consideration of higher-order radiative corrections. At present the theory predictions for neutral current heavy-quark production rely on the complete next-to-leading order (NLO) QCD corrections [4] which are often used via the parametrizations of Ref. [5], see also Ref. [6] for a check and minor corrections. The complete contribution of the next-to-next-to-leading order (NNLO) is not known. However, partial information in various kinematic limits has been accumulated over the years. This includes in particular the logarithmically enhanced terms near threshold which are accessible by means of the soft-gluon exponentiation [7]. Also the high-energy limit of heavy-quark DIS has been known from small-x resummation to all orders at leading-logarithmic accuracy for a long time [8]. In the asymptotic regime of Q 2 ≫ m 2 fully analytic results have been obtained, and the NLO calculations [9,10] have been extended recently towards NNLO by the computation of a number of lowest even-integer Mellin moments [11][12][13], see also Ref. [14]. Finally, the dependence on the renormalization and mass-factorization scales is known exactly at NNLO from standard renormalization-group arguments [7,15] and can be used as an independent consistency check in all these limits.
In the present article, we provide approximate NNLO QCD corrections for the heavy-quark part of the dominant structure function F 2 in photon-exchange DIS. This is achieved by extending and, for the first time, combining the available NNLO information from the different kinematic regimes, i.e., from threshold (s ≃ 4m 2 , where s is the partonic CM energy), high energy (s ≫ 4m 2 ) and large scales Q 2 ≫ m 2 , in particular for the dominant gluon-initiated contribution. Specifically, we present results for all four soft-gluon enhanced logarithms near threshold (see Ref. [16] for a previous brief account), and we employ the all order-result in the high-s limit [8] to derive an analytic expression for the leading ln s term at NNLO. Finally, the known Mellin moments of the heavy-quark operator matrix elements (OMEs) [11][12][13] are employed together with the three-loop results for massless quarks [17] to construct approximate x-space expressions for the heavy-quark coefficient functions at Q 2 ≫ m 2 . By combining these individual constraints we construct NNLO coefficient functions for heavy-quark DIS which, while still being approximate, represent the most comprehensive results possible at this point. Below we will provide a detailed documentation of the required calculations as well as estimates of the accuracy of these approximate results. The resulting improvement of the predictions for heavy-quark DIS and the low-Q 2 small-x limitations of the present NNLO results are then illustrated in a brief phenomenological study.
The remainder of this article is organized as follows: In Section 2 we introduce our basic notations. Section 3 is devoted to the calculation of the new NNLO results for the gluon coefficient function near threshold (Section 3.1), at small-x (Section 3.2) and for scales Q 2 ≫ m 2 (Section 3.3). In the latter two cases we also address the light-quark coefficient function. In Section 4 we combine these results and provide an approximate NNLO expression for both coefficient functions for which we perform a number of quality checks. Section 5 addresses their impact on the charmand bottom-production structure functions F c 2 and F b 2 at NNLO and the limitations of the present results. Finally we conclude and briefly discuss possible future improvements in Section 6. The appendices contain, in Appendix A, the derivation of the analytical result for the gluon coefficient function at NLO in the limit s ≃ 4m 2 and, in Appendix B, the (up to two OMEs) exact, if somewhat lengthy results for the heavy-quark coefficient functions at asymptotic values Q 2 ≫ m 2 .

Setting the stage
We start by briefly reviewing the main results on heavy-quark production in deep-inelastic scattering mediated by neutral-current exchange. For the scattering of a charged lepton e off a proton P this reaction, e(l) + P(p) → e(l ′ ) + q h (p 1 ) +q h (p 2 ) + X , is dominated by virtual photon exchange if the momentum transfer Q 2 = −q 2 = −(l − l ′ ) 2 is much smaller than the Z-boson mass, Q 2 ≪ M 2 Z . The quark pair q hqh in the final state is heavy, m 2 ≫ Λ 2 QCD . The cross section, integrated over the momenta of the outgoing heavy (anti-) quarks, is written in terms of the heavy-flavour structure functions F k (x, Q 2 , m 2 ) with k = 2, L as where α is the electromagnetic coupling constant, m the mass of the heavy quark, and the wellknown DIS variables x and y are defined as x = Q 2 /(2p · q) and y = (p · q)/(p · l), respectively.
Disregarding power corrections, the heavy-quark contribution (2.2) to the DIS structure functions can be written in terms of a convolution of PDFs and coefficient functions as, see Ref. [5], where z max = 1/(1 + 4m 2 /Q 2 ) and e h is the heavy-quark charge. The strong coupling constant at the renormalization scale µ r is denoted by α s , and the PDFs for the parton of flavour i are f i (x, µ 2 f ) at the factorization scale µ f . The kinematic variables η and ξ in Eq. (2.3) are given by 4) and the partonic CM energy is s = Q 2 (1/z − 1). Instead of η, which measures the distance to the partonic threshold, one often uses the heavy-quark velocity β or the variable ρ, The coefficient functions of the hard partonic scattering process can be expanded in α s as where we have identified the renormalization and factorization scales, µ = µ f = µ r . This can easily be undone by expanding α s (µ f ) in terms of α s (µ r ) using the standard QCD beta-function.
In the above normalization, the coefficient functions at the leading order (LO) read [4,18,19]  L,g denote the contributions from transverse and longitudinal photon polarizations, and T F = 1 2 for the colour group SU (N c ). The radiative corrections to Eq. (2.6) at NLO have been known for a long time from Ref. [4]. Unlike the massless coefficient functions, the massive NLO coefficient functions c (1) k,i cannot be expressed in a simple analytic form. Instead, compact parametrizations for these functions were presented in Ref. [5], with minor corrections provided later in Ref. [6] (see Ref. [20] for parametrizations in Mellin-N space including complex values of N as required for a numerical Mellin inversion).
The NNLO coefficient functions c (2,0) k,i in Eq. (2.6) are not fully known, although, as outlined above, substantial partial information is available. On the other hand, all scale-dependent terms at this order, i.e., c (2,1) k,i and c (2,2) k,i , have been constructed by means of renormalization-group arguments [7,21] and are completely known in numerical form, using the parametrized NLO results and the well-known expressions for the NLO splitting functions, see Refs. [7,15]. In the following we will thus focus on improved predictions for the scale-independent parts of the NNLO coefficient functions. To be precise, we will confine ourselves to the gluon coefficient function c (2,0) 2,g , which is by far the most important contribution to the heavy-quark structure function F 2 after the convolution with the gluon distribution in Eq. (2.3), and to the pure-singlet light-quark coefficient function c (2,0) 2,q for the heavy-quark contribution proportional to e 2 h . For further reference we finally introduce the short-hand notations and we also note that throughout this article we will employ the on-shell scheme for the heavyquark mass m, i.e., the pole mass. Heavy-quark DIS with a running mass in the MS scheme has been considered in Ref. [15].

The NNLO corrections in three kinematic regions
In this section we provide improved predictions for the above two NNLO coefficient functions c (2,0) 2,g and c (2,0) 2,q . We address, in this order, the threshold region s ≃ 4m 2 , the high-energy regime s ≫ 4m 2 , and the high-scale region Q 2 ≫ m 2 . In all these limits the higher-order corrections exhibit well-known features and regularities which we briefly discuss and exploit. The results will be assembled into approximate coefficient functions for the whole kinematic plane in Section 4.

Threshold limit and soft-gluon resummation
Partonic cross sections generally receive large logarithmic corrections near threshold, which appear as the result of an incomplete cancellation of the soft (and collinear) contributions between the real and virtual corrections due to the reduced phase-space near threshold. In the present case of massive quarks, the corrections appear as a function of the heavy-quark velocity β at each order of the perturbation series in the form α n s ln m β with m ≤ 2n. The production threshold corresponds to β → 0, recall Eq. (2.5), and the highest logarithms can be resummed to all orders in perturbation theory, see, e.g., Refs. [22][23][24][25][26]. In addition to the logarithmic terms, heavy-quark pair production also receives corrections of the form α n s /β −m ln ℓ β with m ≤ n at all orders from the Coulombexchange of gluons between the heavy quarks. Also these terms can be resummed [27].
The raison d'être of phenomenology based on threshold resummation derives from the wellknown fact, see, e.g., Refs. [7,28,29], that, at not too large values of Q 2 , the convolution of the coefficient function for F 2 with the gluon PDF in Eq. (2.3) is often dominated by rather low partonic CM energies. Therefore the NNLO predictions of the threshold resummation can provide useful information on this numerically important contribution to F 2 . Previous research has already determined exactly the two highest [7] and approximately the third [30] threshold logarithms at this and all higher orders. In what follows, we derive exact expressions for all four logarithmically enhanced terms together with the complete Coulomb corrections for the NNLO gluon coefficient function c (2) 2,g in Eq. (2.6), a result that we have already reported in a numerical form in Ref. [16]. The (light-) quark coefficient function c 2,q for heavy-quark DIS is non-vanishing only from NLO and does not exhibit any enhancement near threshold.
According to the general formula for the threshold resummation [24][25][26], the gluon coefficient function in Mellin space (recall Eqs. (2.4) and (2.5) for the definitions of ξ and ρ) is of the form is the Mellin transform of the LO coefficient function of Eq. (2.7). In the large-N limit (which is dominated by the transverse component c T,g ) it is given by g 0 (α s , N) denotes the matching function to be discussed below, cf. Eq. (3.7). The exponential factor resums the threshold logarithms to all orders in α s . The exponent can be written in the standard form [16,26,31], Here the first term includes the gluonic cusp anomalous dimension A g , known to order α 3 s [32,33], which governs the process-independent soft-collinear gluon emission off the initial gluon. The second term, D γ * g→q hqh , collects process-dependent effects of soft-gluon emission from the initialand final-state particles and is built as from D g and D q hqh which are, respectively, the soft anomalous dimension for Higgs production in gluon-gluon fusion (known to order α 3 s [34,35]) and for the heavy-quark production in the coloroctet channel (known to order α 2 s [36][37][38], see also Refs. [39,40]). For reference we assemble D γ * g→q hqh , employing the usual notation (2.10) for the strong coupling, (3.5) Here and below ζ n stands for values of the Riemann zeta-function.
To next-to-next-to-leading logarithmic (NNLL) accuracy, as needed for the NNLO corrections we are considering, the exponent G N in Eq. (3.3) takes the form G N = ln N g 1 (λ) + g 2 (λ) + a s g 3 (λ) + . . . (3.6) with λ = β 0 a s ln N and, again, a s = α s /(4π). The functions g i (i = 1, 2, 3, . . .) are defined such that g i (0) = 0, with the first term ln N g 1 resumming the leading logarithms and so on. Explicit expressions for g i can be obtained from Ref. [31], even with the µ r and µ f dependence separated (for the computation see, e.g., Refs. [41,42]), with the obvious replacements A q → 1 2 A g , D q → 1 2 D g and D QQ → −D q hqh in Eq. (A.5) -(A.9) of Ref. [31] (note especially the unconventional sign for the soft anomalous dimension of color-octet heavy-quark production, there denoted by D QQ ). Below we will expand in terms of lnÑ = ln N + ln γ e , which absorbes the otherwise ubiquitous Euler-Mascheroni constant into the logarithms.
The process-specific new information needed for heavy-quark DIS resides entirely in the matching function g 0 (α s , N) in Eq. (3.1) which we need to discuss next. This function can be expressed as a product of a hard coefficient g h 0 (α s ) and a Coulomb coefficient g c 0 (α s , N), where the latter specifically induces a dependence on the Mellin variable N, Such a factorized form for the Coulomb corrections and the threshold logarithms has been discussed for top-quark pair production at hadron colliders before [43,44], see also Refs. [36,45] for studies in the framework of soft-collinear effective theory. For heavy-quark DIS a similar factorization applies to the gluon coefficient function, since the structure of radiative corrections to the DIS subprocess γ + g → q h +q h + X is similar to top-quark hadro-production in g + g → q h +q h + X , i.e., essentially the same up to replacements of the corresponding color factors.
The Coulomb coefficient g c 0 (α s , N) in Eq. (3.7) can be obtained to order α 2 s by making use of the result of the non-relativistic (NR) cross section in e + e − → q hqh calculated up to this order in Refs. [46,47]. The only difference in the present case is the color structure, i.e., the Coulomb corrections for the color-octet state require the colour factor replacement C F → (C F − C A /2) in the corresponding e + e − result [46,47]. Thus we identify the NR part of the NNLO coefficient function as where β 0 = 11/3C A − 2/3n f is the leading coefficient of the QCD beta function and we have removed the term corresponding to the product of the NLO hard and NLO Coulomb corrections to avoid double counting. Note that the single logarithmic term originates in the non-leading part of the NR QCD potential [48] and is not included in the soft-gluon exponent G N . Performing the Mellin transform, we obtain g The only missing information for the complete NNLL + NNLO resummation is the coefficient g h (2) 0 in Eq. (3.7), which requires the presently unknown NNLO constants in β, i.e., the two-loop analogues of c 0 (ξ) andc 0 (ξ) in Eqs. (3.10) and (3.11), which can be determined only from a full three-loop calculation of the heavy-quark contribution to DIS in the soft-gluon limit. Setting this constant to zero in β-space, we find g h (2) 0 = 135416 27 ln 2 − 35888 9 ln 2 2 + 11672 9 ζ 3 − 88856 27 + 1340 9 π 2 + 16 3 π 4 − 932 3 π 2 ln 2 − 1384 ζ 3 ln 2 + 488 3 π 2 ln 2 2 + 21584 9 ln 3 2 − 640 ln 4
Supplemented by an estimate of the missing β-space constant, e.g., from a Padé estimate as in Eq. (3.19) below, the results in Eqs. (3.6) -(3.14) can be used to provide an approximate NNLL + NNLO exponentiation of the gluon coefficient function in heavy-quark DIS, which is analogous to the case of hadronic top-quark pair production in Refs. [26,31].
In the present article, though, we are rather concerned with improved predictions for c 2,g to NNLO, essentially using Eq. (3.1) to generate results to fixed order in perturbation theory. To that end we now write down the β-space threshold behaviour for c 2,g up to NNLO. The LO and NLO results [4,18,19] are in terms of c 0 (ξ) andc 0 (ξ) in Eqs. (3.10) and (3.11). The convention for these constants derives from Refs. [4,5], where the logarithmic enhancement was expressed in terms of ln ℓ (8β 2 ), ℓ = 1, 2, and ln(4β 2 ) for the scale-independent and scale-dependent parts, respectively.
The corresponding NNLO result c (2) 2,g is obtained from our results above using the Mellin transform in Appendix A of Ref. [31], Several checks can be performed of this result which provides, for the first time, all logarithmically and Coulomb-enhanced contributions. First of all, the scale dependence at NNLO is fully known [7,15] from the renormalization group (see also Section 4), and we have verified that all terms proportional to powers of L µ in the threshold expansion in Eq. (3.18) agree. Another strong check is the agreement of Eq. (3.18) with the results presented in Eq. (A.1) of Ref. [38] for the hadronic heavy-quark production in the gluon-gluon-fusion channel, gg → q hqh , after replacing one gluon by a photon (together with the appropriate substitutions of the colour factors).
The NNLO constants in β multiplying the Born term are currently unknown. This includes the scale-dependent term c of which we will employ the scale-independent part at low ξ in Section 4 below. The main new NNLO result of this section, the scale-independent part of Eq. (3.18), is illustrated in Fig. 1 for two values of Q 2 . It is clear that a stable low-η result is established only by including our new results for the third and fourth logarithms. As a first estimate of the remaining uncertainly at η > ∼ 1 for low values ξ > ∼ 1 we have included the Padé estimate (3.19), assigning a 100% uncertainty (its small value at ξ = 20 is accidental and irrelevant in what follows).

High-energy expansion
Let us next discuss the high-energy limit of the coefficient functions c 2,g and c 2,q . To that end we build on the results of Ref. [8] from the small-x resummation, which had been derived even before NLO corrections were available and thus served as a check on Ref. [4]. At fixed Q 2 the high-energy limit s Therefore we consider the coefficient function c 2,g as a function of x throughout this section. In Mellin N-space, cf. Eq. (3.1) for the definition, the dominant contributions in this region behave as α n s /N n and can be resummed to all orders. One can define the partonic cross sectionσ 2,N in Mellin space corresponding to the gluon coefficient functions c 2,g (see Eq. (3.1)) as m 2σ 2,N = α α s e 2 q c 2,g (α s , N) , (3.20) which contributes to the structure function F 2 via Eqs. (2.2) and (2.3). In complete analogy, if multiplied with the Mellin moments of the gluon PDF f g,N , σ 2,N also defines the gluonic contribution to the total hadronic cross section. It is an interesting observation of Ref. [8] that the Q 2 -dependence in this quantity can be factorized if one considers the photo-production limit, i.e., the process γ g → q hqh X with a real photon. This leads to the ansatz The function K 2,N summarizes all ξ-dependence originating from the momentum transfer Q 2 of the off-shell photon in DIS electro-production of heavy quarks, and σ γ g,N denotes the Mellin moments of the gluonic contribution to the total hadronic cross section for photo-production of heavy quarks, see, e.g., Ref. [49,50].
The resummation of the logarithms at high energy is based on the framework of un-integrated PDFs in transverse momentum k t and the concept of k t -factorization, a procedure which involves two steps: first computing amplitudes with the initial particles off-shell in k t , and second performing the convolution with a gluon PDF which has the small-x corrections included. For heavy-quark photo-production, this allows to express σ γ g,N in Eq. (3.21) as the product [8] (cf. also Ref. [50]), that is, in terms of the gluon PDF and the impact factor h N (γ N ) depending on the anomalous dimension γ N . We are interested in the perturbative regime with γ N given by with a ≡ C A α s /π. This well-known result for γ N arises as the solution to the following condition on χ(γ), in a perturbative expansion for γ ≪ 1 and requiring γ(0) = 0. Here χ(γ) is the Mellin transform (with respect to the transverse momentum k 2 t ) of the lowest order BFKL kernel expressed through standard ψ-functions as (3.25) The anomalous dimension γ N governs the high-energy behavior of the gluon PDF in Eq. (3.22) as which up to terms relevant at NNLO simply gives In the perturbative regime the impact factor is needed only in the limit N ≪ γ −1 N , and can thus be approximated as (3.28) because h N (γ N ) is free of singularities and terms of order N are subleading. The calculation of the impact factor h is performed by considering off-shell amplitudes in a k t -factorization scheme. For photo-production of heavy quarks in γ g → q hqh , and with the incoming photon and gluon being off-shell by an amount k t , this results in [8,50,51] in terms of standard Euler Beta-functions. For γ N ≪ 1 this expression leads to h(γ N ) = π α α s e 2 Note that the precise definition of the k t -factorization scheme for the impact factor and the conversion to the MS scheme affects the result in the perturbative regime only at order γ 3 N , see Ref. [50]. This is beyond the NNLO accuracy in α s we are aiming at.
Finally the ratio K 2,N between the cross sections for DIS heavy quark electro-production and photo-production in Eq. (3.21) has been computed in closed analytical form in Ref. [8], where 2 F 1 denotes the hypergeometric function which can be conveniently expanded for small γ N in terms of harmonic polylogarithms H m (x) [52] (see e.g., Refs. [53][54][55][56] and Appendix B). To that end, we have used the HYPEXP2 package [57] in MATHEMATICA, which gives us  of which I(ξ) and J(ξ) have been introduced already in Ref. [8] where also an expression for I(ξ) in terms of normal (di-)logarithms can be found.
At this stage, all that remains is to perform the expansion of Eqs. (3.22) and (3.31) to second order for small γ N by assembling all formulae provided so far. Converting from N-space back to x-space we obtain for the gluon coefficient function, 2,g is a new result. For further reference (in Section 3.3, where we will also provide graphs of the above results), we have finally derived explicit results for the ξ ≫ 1 limits, i.e., Q 2 ≫ m 2 , of Eqs. (3.38) and (3.39). We obtain, at leading-logarithmic small-x accuracy, (3.41) Here the constant terms at NNLO, i.e., the terms independent of L Q and L µ (recall Eq. (2.10) above) in Eq. (3.41) are new. All other terms in this result agree with independent calculations. The scale-dependent L µ terms can, of course, be obtained by the standard renormalization-group behaviour, cf. Section 4. More importantly, we are also able to compare all L Q terms with the results obtained in the asymptotic limit Q 2 ≫ m 2 below, cf. Section 3.3. All powers of L Q and L µ match exactly with those results, which constitutes a strong check. Another test of Eq. (3.39) is possible by comparing with the small-x limit for c 2,g derived in Ref. [59] for the kinematic region Q 2 ≃ m 2 (see also Eqs. (10) and (11) in Ref. [30]), and we find agreement at the level of a few parts in a thousand for Q 2 /m 2 = 1.
In concluding this section, we briefly comment on the quark coefficient function for heavyquark DIS when the photon couples to the heavy quark, that is, c 2,q in Eq. (2.3) proportional to e 2 h and the quark-singlet PDFs. Its leading small-x term is related to that of c 2,g in Eqs. (3.38) -(3.41) by a simple colour factor replacement, c 2,q = C F /C A c 2,g , as pointed out already in Refs. [8,59].

High-scale limit
The results in the limit Q 2 ≫ m 2 build on an exact factorization of the heavy-quark coefficient functions into the respective coefficient functions with massless quarks and heavy-quark operator matrix elements (OMEs). The variable η in Eq. (2.4) factorizes as η → Q 2 /(4m 2 )(1/x − 1) + O (1), and, to denote this limit, we will be considering the coefficient functions as a function of x and the ratio Q 2 /m 2 throughout this section. At LO, one can quickly verify for c 2,g with the help of Eq. (2.7) that this limit reads where the term multiplying L Q is proportional to the lowest-order splitting function P (0) qg , while the L Q independent term (here written in terms of harmonic polylogarithms H m (x) [52], see Appendix B) is proportional to the one-loop gluon coefficient function with massless quarks.
Formally the underlying factorization may be expressed through massive partonic OMEs A i j (i.e., with massive quarks) and the DIS coefficient functions with massless quarks C light k, j together with their exact scale dependence from renormalization group invariance, see Ref. [21], i.e., at scales µ = Q and µ = m, Here ⊗ denotes the standard Mellin convolution, cf. Eq. (2.3). The dependence on the scale µ has been made explicit together with the separation of scales between the OMEs A i j and the coefficient functions C light k, j , which is realized through ln Q 2 /µ 2 = L Q + L µ . At NLO the formalism in Eq. (3.43) has first been applied in Ref. [9] to compute the corresponding two-loop OMEs for (unpolarized) heavy-quark DIS at asymptotic values Q 2 ≫ m 2 , see also the re-calculation in Ref. [10]. In the wider context, Eq. (3.43) implements the matching conditions for QCD with one massive and n f light quarks to QCD with n f +1 light quarks (note that n f denotes the number of massless quarks throughout this article). The matching of course also affects the strong coupling α s , for which the decoupling formulae are well known [60,61]. Unless indicated otherwise, we are always working in the scheme with α s (n f ), i.e., n f = 3 for charm quarks. Moreover, on the basis of Eq. (3.43), also the all-order resummation of the logarithms in Q 2 /m 2 has been discussed in the literature and has led to definitions of a so-called variable flavour number scheme (VFNS) [62], see Refs. [63,64] for implementations in phenomenological analyses.
To that end we rely on the extensions of the two-loop massive OMEs to higher orders in ε (the parameter of dimensional regularization) [11], on the two-loop gluonic OMEs [12] and, most importantly, on the even-integer Mellin moments of the three-loop heavy-quark OMEs [13] along with their complete n f -dependence [14]. The other important input consists of the results for the anomalous dimensions [32,33] and coefficient functions for DIS with massless quarks [17,65,66]  up to three loops. The necessary convolutions, Mellin transforms and inverse Mellin transforms are all performed with algorithms for harmonic sums [67] and harmonic polylogarithms [52] as implemented in the symbolic manipulation program FORM [68].
The coefficient functions C (being proportional to the colour factor d abc d abc ) affect the factorization for Q 2 ≫ m 2 in Eq. (3.43) in the following way. The available results for the massless coefficient functions in Ref. [17] (derived via the optical theorem) contain the sum of the contributions from all final state cuts which contribute to the inclusive massless structure function, say F 2 . However, these are not identical to the contributions for the heavy-quark structure function in Eq. (2.3), the latter being semi-inclusive and requiring a heavy-quark pair in the final state. This is because some admissible final state cuts go through the quark line while others do not, e.g., the cut of the two gluon lines in the right diagram of Fig. 2. Currently available information is, unfortunately, not sufficient to disentangle the different terms. For that reason, we have decided to omit the (usually small) threeloop contributions proportional to f l 11 and f l g 11 in C light k, j in the following, cf. also Appendix B. The general expansion of the heavy-quark OMEs in Eq. (3.43) in powers of α s reads where, at each order, the terms proportional to powers of L µ are determined by lower order OMEs and splitting functions (anomalous dimensions). The genuinely new ℓ-th order information resides in the expressions for a (ℓ,0) i j . The available results for those have been given in the literature for the heavy-quark mass in the on-shell scheme (pole mass) as well as in the MS scheme, see Ref. [13]. Note, that in the latter case, certain constants such as ζ 2 and ζ 2 ln 2 are absent in Mellin space, similar to DIS structure functions with massless quarks [17].
For the NNLO gluon coefficient function c (2) 2,g , we specifically need the heavy-quark OME A q h g (= A Qg in the notation of Ref. [13]); and for the NNLO light-quark coefficient function c (2) 2,q the pure-singlet heavy-quark OME A ps q h q (= A ps Qq in the notation of Ref. [13]). The three-loop expressions A (3) Qg and A (3),ps Qq are given in Eq. (4.39) and Eq. (4.26) of Ref. [13], respectively, with the µ-independent terms denoted by a (3) Qg and a (3),ps Qq . They can be decomposed in powers of n f as where the n f terms are now known exactly from Ref. [14]. For the n f = 0 terms, on the other hand, only a number of integer Mellin moments have been computed so far [13]. These are thus the only quantities missing for the construction of c (2) 2,g and c (2) 2,q via Eq. (3.43).
Therefore, leaving the functions a There are a number of checks on these results. First of all, up to NLO we agree with the previous results [9,10]. Next, in the convolution Eq. (3.43) for c (2) 2,g and c (2) 2,q all cubic powers L 3 µ from the individual renormalization of the OMEs A Qg and A ps Qq cancel, as they have to, against those from the massless coefficient functions C light 2,g and C light 2,q at three loops. Moreover, all remaining powers of L µ in c (2) 2,g and c (2) 2,q agree numerically with the exact results for the µ-dependence derived by renormalization-group methods in Refs. [7,15] which are valid for all ratios of Q 2 /m 2 , if the latter are evaluated in the present limit Q 2 ≫ m 2 , see also Section 4.
As a last step, in order to arrive at phenomenologically useful results, it thus remains to use the available information on the Mellin moments of a Qq, ps (x) for constructing (hopefully sufficiently accurate) approximate expressions together with estimates of their residual uncertainty. This proceeds along the lines of, e.g., Refs. [69][70][71], where Mellin moments of three-loop splitting and coefficient functions were successfully employed to derive useful approximation of these quantities prior to their exact computations [17,32,33,72,73].
A vital constraint on the functions a Qq, ps (x) is provided by the leading-logarithmic small-x behaviour of c (2) 2,g and c (2) 2,q determined above, see Eq. (3.41) and the last paragraph of Section 3.2. Since all other contributions to c (2) 2,g are known, that result can be used to deduce Taking into account also the double-logarithmic large-x terms, the OME a Qg can be written as is a complicated combination of harmonic polylogarithms which approaches a constant for x → 1 and includes powers of ln x at small x. For approximations based of the five known Mellin moments, N = 2, 4, . . . , 10, of Ref. [13] we select two of the large-x constants A m , the subleading small-x parameter B, and a two-parameter smooth function build from ln 1,2 x and various linear and quadratic functions in x.
The comparison of the Mellin transform of Eq. (3.48) with the known moments results in a system of linear equations for the chosen coefficients A m , B and the parameters of f (x). Varying the selected two A m and the two-parameter form of f (x), we arrive at a large number of approximations which indicate the remaining uncertainty of a Qg (x). As in Refs. [69][70][71], we discard a small number of functions for which the system of equations is almost singular, resulting in huge numerical coefficients and unrealistically large oscillations of the function. In this manner we have determined 50 to 100 acceptable approximations and finally selected two representatives, which reflect the error band for most of the x-range (and particularly at small x), for the final estimate of a Qg and its residual uncertainty.
The results of this process are shown in the left part of Fig. 3. The chosen two approximations (solid lines) are given by where the coefficient of x −1 ln x is a truncation of the exact result (3.47). The average of the two extremes can be used as the central result of our approximation procedure.
The same procedure has been applied to the function a Here the lowest six even-integer moments have been computed in Ref. [13]. Taking into account that this puresinglet quantity vanishes for x → 1, our approximations are built from a subset of the functions together with the small-x limit given by Qg of Eq. (3.47). The results are shown in the right part of Fig. 3, and the extremal representatives have been chosen as where, again, the coefficient of x −1 ln x arises from the truncated exact result, and the average of Eqs. (3.52) and (3.53) provides the central result.
Using these approximations, the high-scale coefficient functions c 2,g and c 2,q can now be assembled up to NNLO using Eqs. (B.4) -(B.10). The successive LO, NLO and NNLO results are shown in Fig. 4 at Q 2 = 100 GeV 2 , a scale that is low enough to be relevant to the HERA measurements and sufficiently high for the safe applicability of the Q 2 ≫ m 2 approximation (3.43). We observe a satisfactory convergence, especially of the dominant quantity c 2,g , over a wide range of x, together with the well-understood large-x region of soft-gluon enhancement for c 2,g , cf. Ref. [42]. For both coefficient functions the uncertainty band due to the heavy-quark OMEs A Qg and A ps Qq as given by Eqs. (3.49), (3.50) and (3.52), (3.53), respectively, is phenomenologically relevant only at small x, up to x ≈ 10 −2 . The uncertainty in this region is due to that of the (large) coefficient of the 1/x terms in Eq. (3.48) and its pure-singlet counterpart. The computation of more Mellin moments along lines of Ref. [13] would certainly help to improve the approximation also here; however the issue would be settled for all practical purposes by extending the results of Ref. [8] to the next-to-leading small-x terms at order α 3 s . c 2,g (2,1) c 2,g (2,2) : 2,g , ℓ = 0, 1, 2, to the NNLO gluon coefficient function defined in Eq. (2.6). Middle and bottom panels: the respective next-to-leading η 0 coefficients for c (2) 2,g and c (2) 2,q . The solid (black) lines are the exact all-Q 2 results, the dashed (blue) ones the high-scale asymptotic results; the dotted (red) low low-scale extrapolations are discussed in Section 4 below. Also illustrated, by the thin (green) lines in the bottom panels, is the small next-to-leading high-η deviation of c (2) 2,q from the 'Casimir-scaled' results for c (2) 2,g .
A crucial question is, obviously, down to which values of ξ = Q 2 /m 2 these high-scale asymptotic results are applicable. For the important high-energy (small-x) contributions, which one may expect to be least affected by m 2 /Q 2 corrections, this issue is addressed in Fig. 5 for the quantities c (2,ℓ) 2,i , i = q, g, in Eq. (2.6). Here the exact leading, cf. Eqs. (3.39) and (3.41), and next-to-leading high-η results are compared to the respective high-scale expressions. Except for the L 2 µ quantities c (2,2) 2,i in Eq. (2.6) fixed by LO information, the asymptotic results provide good approximations of the ln η parts of c (2,0) 2,i and the ln η and η 0 parts of c (2,1) 2,i down to ξ ≃ 4, but then deteriorate dramatically towards ξ = 1. This behaviour can be used to extrapolate the (unfortunately still rather uncertain) asymptotic coefficients of η 0 for c (2,0) 2,i to small ξ. This extrapolation will be discussed in the next section, but the corresponding curves have already been included in the figure.

Approximate coefficient functions at NNLO
We are now in a position to construct improved NNLO approximations for the genuinely new parts of the NNLO gluon and pure-singlet (e 2 h ) coefficient functions, c 2,q , for heavy-quark DIS. Specifically, we will discuss how to combine the above (approximate) expressions in the three kinematic limits, i.e., threshold, high-energy and large Q 2 , in order to arrive at functional forms which smoothly interpolate over the full relevant range in η and ξ (recall, again, Eq. (2.4)).
All logarithmic terms proportional to powers of L µ in Eq. (2.6) are known exactly from standard renormalization-group methods. For completeness, we give their explicit form in the MS scheme. Suppressing all arguments, these coefficients read up to NNLO [7,15,21] where, as in Eq. (3.43) above, ⊗ denotes the standard convolution. The splitting functions P (l) i j (x) can be taken from [32,33] and β 0 and β 1 are the standard coefficients of the QCD beta-function, normalized such that β 0 = 11/3 C A − 2/3 n f . The powers of 4π account for the normalization of Refs. [4,5] adopted in this article, cf. Eq. (2.6). Lacking more exact third-order information, in particular c (2,1) 2,g in Eq. (4.5) will provide important guidance below for the assembly of the approximate NNLO result for the most important quantity c Two issues need to be considered for the combination of the results from the various kinematic regions. Firstly, one has to merge the available information from threshold at low β (and η ) with the high-energy limit at large η. Since we use the complete Born coefficient function (2.7) which vanishes for η → ∞ as the prefactor in the threshold expansion (3.18), this can be done by multiplying the high-energy terms by sufficiently (but not excessively) high powers of β. Secondly, the asymptotic expressions of the large-ξ limit need to be joined with the low-ξ region. To that end we employ the function which provides a smooth transition between these two regimes. The parameters 2 and 4 in the exponent are chosen in such a way that the transition is sufficiently rapid at values close to ξ = 5. This choice is motivated by the above discussion of Fig. 5 and the finding of Ref. [9] that, e.g. for c (1,0) 2,g , the asymptotic limit (B.6) represents the exact result with a high accuracy already at ξ > ∼ 10.
These considerations lead the following ansatz for c (1) 2,g (cf. also the NLO parametrizations in Ref. [5]) which we test below, Here c (1)LLx 2,g and provides the high-η tail in the low-ξ region which is smoothly matched with the factor η γ /(C +η γ ). The values of γ, C and k in Eq. (4.8) will be specified below.
Correspondingly, for c (2) 2,g we use is the leading contribution (3.39) in the small-x limit, and we have divided out the factor ln x in order to be able to substitute ln x → − ln η, which is valid at high energies and determines the slope in η at all values of ξ. The next-to-leading large-η term, denoted by c . It is currently unknown in the low-ξ region. We will derive constraints on c (2)NLL 2,g below, although, even at large ξ, we still have to cope with the uncertainty of the heavy-quark OME A Qg estimated by Eqs. (3.49) and (3.50). Recall also that the unknown next term in the low-η threshold expansion (3.18) consists of terms proportional to the Born result as discussed in Section 3.1.
Therefore we adopt the following strategy for applying Eqs. (4.8) and (4.9). When merging the contributions from the various regions, we account for the residual uncertainties and construct two enveloping curves designed to span the uncertainty band in the entire kinematic plane. Thus at NLO we define Here c The resulting two approximations are compared, after adding the exact scale-dependent contribution, in the left parts of Fig. 6 to the parametrized [5] exact result of Ref. [4] at four typical values of ξ for the standard choice µ 2 = Q 2 + 4m 2 of the (identical) renormalization and factorization scale. It is worthwhile to note not only that the approximations (4.10) and (4.11) indeed span the exact results, but also that there is, at this scale, a considerable cancellation between the high-energy constants in c At NNLO we proceed in a similar manner. As indicated below Eq. (4.6), we first consider the known quantity c (2,1) 2,g and define with c and (3.39). For the next-to-leading term at small-x denoted by c The damping of this NLL high-energy tail can be approximated by (again using MINUIT) These approximations for c (2,1) 2,g are compared in the right parts of Fig. 6 with the exact result (4.5), for which these graphs supersede previous studies in Refs. [7,30]. Also these approximations are completely satisfactory, taking into account that a term corresponding to c 2,g , at four representative values of ξ = Q 2 /m 2 , with the approximations based on the threshold, high-η and high-Q 2 limits as respectively specified in Eqs. (4.10) and (4.11) and Eqs. (4.13) and (4.14). The first-order results are shown at a physically relevant scale µ after adding the exact scale-dependent contribution to all three curves.
Encouraged by these results, we now apply the same procedure to the main new NNLO coefficient function c (2,0) 2,g and write In Fig. 7 we display these approximate NNLO results for c 2,g are then added in Fig. 8 to these expressions for the standard scale µ 2 = Q 2 +4m 2 which we will also use in Section 5 below. The lack of a stronger constraint on the next-to-leading high-η coefficient leads to a large uncertainty in particular at low ξ, where it extends down to η ≈ 10 in the construction detailed above. We have investigated other approaches to this important mid-ξ region, but have been unable to derive a reliable stronger constraint with the available information on c (2,0) 2,g . The present NNLO results in Fig. 8 are consistent, albeit if with a large uncertainty, with the 'high-energy censorship' already seen for c (1) 2,g in Fig. 6, i.e., a small contribution of the high-energy tail to the total coefficient function despite large contributions to the individual components c For completeness, we apply the same procedure to the pure-singlet coefficient function c 2,q which exhibits a closely related high-energy behaviour, recall the end of Section 3.2 and the bottom parts of Fig. 5, but no low-η threshold enhancement. Suppressing the quantities corresponding to Fig. 6 for brevity, we directly turn to the new NNLO quantity c (2,0) 2,q which thus we approximate by The high-ξ asymptotic result c The resulting approximations for c (2,0) 2,q and the (small -compare the scales of the right parts to those of Fig. 8) complete coefficient function c (2) 2,q at µ 2 = Q 2 + 4m 2 are illustrated in Fig. 9. So far we have separately considered approximations to the NNLO gluon and pure-singlet quark coefficient functions. However, as mention above, the small-x/ large-η limits of these two functions are closely related; and while their C F /C A relation does not hold exactly beyond the leading logarithms, it still approximately holds for the NNLO η 0 -terms as shown already in the bottom panels of Fig. 5 in Section 3.3. In particular, the numerical breaking of this relation for these terms is much smaller than the current uncertainties of c (2,0) 2,g and c (2,0) 2,q in the low-ξ region. This opens up the possibility of one final improvement, i.e., using the C F /C A relation and the better constrained coefficient function to reduce the uncertainty of its more uncertain counterpart. For both the OME and the coefficient function the high-energy uncertainty is smaller in the quark case, as can be seen from Figs. 3 and 5, which is unsurprising given that one more high-ξ moment is known for the OME in this case [13], which is moreover an 'easier' function with a

Charm and bottom structure functions at NNLO
We finally illustrate the impact of the above approximate NNLO coefficient functions on the heavyquark contribution to the structure function F 2 in Eqs. (2.2) and (2.3). Since present-day parton distributions are rather well constrained for most of the kinematic range relevant to the charm and bottom structure functions, we confine ourselves to one set, that of Ref. [74]. See Refs. [63,[75][76][77][78][79] for other recent PDF fits including NNLO corrections. For the heavy-quark pole masses we use m 2 c = 2 GeV 2 and m b = 4.5 GeV , (5.1) and our default choice for the renormalization and mass-factorization scales is In Fig. 10 the resulting charm structure function F c 2 is shown for seven scales Q 2 which correspond to the ξ-values chosen in Figs. 7 -9 above. The dotted and short-dashed curves are the contributions of the LO and NLO coefficient functions to the NNLO results, i.e., these quantities have been calculated using the NNLO values of the strong coupling and the parton distributions as determined in Ref. [74]. The solid and dash-dotted curves show the NNLO corrections for the approximations 'A' and 'B' constructed in the previous section, including the improvement (4.25) of the latter. Where the approximations are sufficiently accurate, roughly speaking at x > ∼ 10 −3 , the overall effect of the NNLO coefficient functions is positive but considerably smaller than that of the NLO contribution, indicating a good convergence of the perturbation series. At x < ∼ 10 −3 the present results are inconclusive in this respect in particular, unfortunately, in the important lowand medium-Q 2 region. Also shown, by the long-dashed curves, are the physical NLO results involving the NLO coupling constants and parton distributions according to Ref. [74].
The corresponding results for the bottom-production contribution to F 2 are illustrated in Fig. 11 for the slightly shifted region 0.5 < ∼ ξ < ∼ 50. The pattern of the corrections is qualitatively similar to that for charm production in Fig. 10. Again the NNLO approximation uncertainty is largest at x < 10 −3 and 1 < ∼ ξ < ∼ 5, here corresponding to 20 GeV 2 < ∼ Q 2 < ∼ 100 GeV 2 , where a pure threshold approximation can be insufficient and the mid-η constraints of the high-scale expressions are not applicable yet. Note that this figure includes a bin at Q 2 < m 2 h , where the impact of the mediumand large-η parts of the coefficient functions is smaller than at higher scales. In general, the size and uncertainty of the NNLO corrections are smaller for F b 2 than for F c 2 at the same values of ξ due to the smaller values of α s and the steeper small-x PDFs at the corresponding higher scales.
Returning to the case of charm production which has a far larger impact on the determination of especially the gluon distribution from DIS data, Figs. 12 and 13 provide a more detailed look at the NNLO results and their remaining uncertainties for two typical values of Q 2 . The former figure shows the dependence on the scale µ for the range m c ≤ µ ≤ 2 µ std with µ std = µ of Eq. (5.2) which includes both hard scales entering the problem, 2m c and Q. Except at very low scales, µ < 2m c , (and possibly at very small x) the NNLO scale variation of F c 2 is smaller than its NLO counterpart given at µ = µ std by the long-dashed curves in Fig. 10. This is also obvious from Fig. 13  uncertainty which becomes important at x ≃ 10 −3 . The relatively small NNLO improvement at large x indicates the need for yet higher orders. Here a first step would be to include the (almost complete) NNLL threshold resummation discussed in Section 3.1.

Summary and outlook
The production of heavy quarks, in particular charm, contributes a sizeable fraction of the proton structure function F 2 measured at HERA, thus affecting fits of the parton distributions and hence the predictions for hard processes at other colliders such as the LHC. Unlike the case of massless quarks, the corresponding partonic cross sections (coefficient functions) are not fully known yet at the next-to-next-to-leading order (NNLO) of perturbative QCD. In this article, we have collected and extended the partial NNLO results for three kinematic limits, i.e., the threshold region with s > ∼ 4m 2 , the limit of high partonic CM energies, s ≫ m 2 , and the high-scale region Q 2 ≫ m 2 , and then combined these results into approximate expressions covering most of the kinematic plane.
For the NNLO gluon coefficient function c (2) 2,g we have determined all threshold contributions which are logarithmically enhanced in the heavy-quark velocity β, the so-called Sudakov logarithms, as well as the Coulomb corrections leading to powers of 1/β (Eq. (3.18) in Section. 3.1). In the high-energy limit, we have employed results of the corresponding 'small-x' resummation [8] to derive explicit expression for the (closely related) leading-logarithmic ln s contribution to the gluon and quark coefficient functions at order α 3 s (Eqs. (3.39) and (3.41) in Section 3.2). Finally, in the high-scale limit where m 2 /Q 2 power corrections can be disregarded, we have utilized the mass-factorization formula to derive exact expressions for the NNLO gluon and (pure-singlet) quark coefficient functions c (2) 2,g and c (2) 2,q (Section 3.3 and Eqs. (B.7) and (B.10) in Appendix B). The latter results are complete for the contributions proportional to n f , the number of light flavours, and rely on the three-loop operator-matrix elements (OMEs) of Ref. [14] for the Q 2 independent contribution. Their counterparts for the n f -independent contributions are not fully known yet, though, so we had to rely on the low even-N Mellin moments computed in Ref. [13] and our above small-x result for deriving the approximate expressions given in Eqs. (3.49) -(3.53) in Section 3.3. The largest uncertainty in these results is due to the next-to-leading small-x (at NNLO behaving as s 0 at large s) contributions which are not tightly constrained at this point.
We have then combined the results of the various regions to provide the best possible approximations for the coefficient functions for heavy-quark production in DIS at all values of s and Q 2 , carefully modeling and estimating the impact of missing information such as the NNLO threshold constant and of the uncertainty of the high-s constants and their dampening towards medium values of s. The corresponding results are given in Eq. (4.17) -(4.25) in Section 4. The implications of these uncertainties for the charm and bottom structure functions F c 2 (x, Q 2 ) and F b 2 (x, Q 2 ) have been illustrated in Section 5, where we have also addressed the dependence on the renormalization and factorization scale µ which can signal the importance of yet higher orders. The uncertainty due to the missing information on the NNLO coefficient functions is large at x < ∼ 10 −3 especially at low and medium scales Q 2 , but irrelevant at much larger values of x. The dominant uncertainty in that region (where the absolute values of F c 2 and F b 2 are much smaller) is due to higher orders, as shown by the rather modest NNLO improvement of the µ-dependence.
At high scales, Q 2 > ∼ 10 m 2 , the remaining small-x problem can be solved by extending the k t -factorization results of Ref. [8] to the next-to-leading high-energy contributions at order α 3 s or by extending the calculations of [14] to the n f -independent OMEs. Already the extension of the fixed Mellin-N computations of Ref. [13] by a couple of moments, however, would facilitate a considerable reduction of the uncertainties. These calculations could also be employed to considerably reduce, but presumably not remove, the large small-x uncertainties at low and medium scales, Q 2 < 10 m 2 . It is also worthwhile to note that there is a close relation between the nextto-leading high-energy contributions in the present case and in heavy-quark hadro-production, see Refs. [8,50]. This has been exploited in Ref. [80] to estimate the high-energy tail of the coefficient functions for top-pair production at the LHC, which complements numerically determined results on NNLO heavy-quark production in hadron-hadron collisions [81] valid in the region where the partonic CM energy s is not too large, s/4m 2 < ∼ 100.
The uncertainties at low Q 2 may be removed to a fully satisfactory extent only by an exact thirdorder calculation of heavy-quark production in DIS. This is a formidable task, including threeloop Feynman diagrams with two scales, i.e., an internal mass and an off-shell leg. A useful, but conceptually easier intermediate step would be a calculation in the threshold limit determining the missing threshold 'constant' (i.e., the contribution proportional to the Born cross section) at order α 3 s . This result would also provide the only piece missing for extending our results in Section 3.1 to a NNLO+NNLL (next-to-next-to-leading logarithmic) threshold resummation. This resummation, even without the NNLO constant, provides the natural framework for a first step beyond the third order, e.g., for addressing the still rather large scale uncertainty at large x.
At this point, though, our approximate third-order results represent the best and most complete predictions of heavy-quark DIS in the important context of NNLO determinations of parton distributions from inclusive and heavy-quark structure functions measured at HERA. Given the wider implications of those analyses on LHC physics, however, further improvements in particular at small x would be extremely welcome. For use until then, FORM files and FORTRAN subroutines with our main results can be obtained from the preprint server http://arXiv.org by downloading the source of this article. Furthermore they are available from the authors upon request.

A Gluon coefficient function at NLO for s → 4m 2
In this first appendix we derive the functions c 0 (ξ) andc 0 (ξ) in Eqs. (3.10) and (3.11) for the gluon coefficient function c 2,g at NLO. This derivation starts from the coefficient function in the differential kinematics of one-particle inclusive DIS [4], where the kinematic variables for the subprocess g(zp) The variable s 4 = s ′ + t 1 + u 1 = M 2 X − m 2 ≥ 0 measures the inelasticity of the subprocess, with s 4 → 0 corresponding to the limit of soft gluon emission. The upper limit of s 4 is given by The LO differential coefficient function is given by [4,18] The NLO differential coefficient function was calculated in Ref. [4]. In the threshold region s 4 ≃ 0, it can be simplified as [4,7,30] Here the factor K (1) (in the MOM scheme for the strong coupling α s ) is given by [4] and the +-distribution in s 4 defined as The conversion of K (1) to the standard MS scheme for α s proceeds via the well-known relation and affects the functionc 0 (ξ) for the scale-dependent matching constant, see Eq. (A.14) below. The scale-independent function R(s ′ ,t 1 , u 1 ) in Eq. (A.6) originates from the soft-and virtual-gluon contributions, for which the analytic expression is quite lengthy and available only via a FORTRAN code by the authors of Ref. [4].
In order to identify the leading contribution for s → 4m 2 , i.e., in the limit β → 0, we first consider Eq. (A.1) for the LO coefficient function with the integrand given by Eq. (A.3). After the trivial s 4 -integral, we have u 1 = −s ′ − t 1 , and the leading contribution in the t 1 -integral can be obtained by the simple replacement and multiplying with an overall factor s ′ β. This can be understood by observing where the first term on the right-hand side gives the leading contribution in the β unless the derivative of f (t 1 ) at t 1 = −s ′ /2 contains negative powers of β. Furthermore, we perform expansion in β such that s = 4m 2 (1 + O(β 2 )) and we recover from Eq. Using these formulae together with Eq. (A.10) and the expansions in β, we recover the threshold approximation of the gluon coefficient function c (1) 2,g in Eq. (3.16), now with explicit results for the the constant terms c 0 (ξ) andc 0 (ξ) at NLO.
The scale-dependent termc 0 (ξ) which is also fixed by renormalization-group constraints can be easily read off from Eq. (A.6), where the last term is due to the transformation of the strong coupling from the MOM scheme in [4] to α s in the standard MS scheme [12,62], cf. Eq. (A.9). The calculation of the scale-independent term c 0 (ξ) on the other hand requires some automated manipulations of R(s ′ ,t 1 , u 1 ), based on the FORTRAN code of [4], which finally lead to the new result (3.10) in Section 3.1.

B Exact results at asymptotic values Q 2 ≫ m 2
We close by presenting the exact expressions for the heavy-quark coefficient function H 2,g and H ps 2,q at high scales, Q 2 ≫ m 2 . In order to shorten the results, we have always separated the contribution of the coefficient functions with massless quarks in photon-exchange DIS, denoted as c (ℓ) 2,k below, as computed in Ref. [17] at ℓ ≤ 3. See also Refs. [65,66] for earlier results up to the second order.
For easier comparison with the literature, we will use the following established normalization of these asymptotic coefficient function, cf. Ref. [9,10], and correspondingly for H ps 2,q . As above we expand in terms of a s = α s /(4π) and, again following Refs. [9,10], this expansion is performed in the MS scheme with α s (n f + 1). All terms originating from subsequent matching α s (n f + 1) → α s (n f ) with the help of the decoupling formulae [60,61] will be presented separately as well. The relation of, e.g., the quantity (B.1) to the gluon coefficient function c 2,g , defined in Eq. (2.3) in the normalization of Refs. [4,5], is given by where we have suppressed the additional dependence on µ = µ f = µ r .
The explicit results are presented in terms of harmonic polylogarithms H m 1 ,...,m w (x) with m j = 0, ±1, where our notation follows Ref. [52] to which the reader is referred for a detailed discussion. For chains of indices zero we again employ the abbreviated notation in which also the argument of H m has been suppressed for brevity. The numerical evaluation of the the harmonic polylogarithms relies on the FORTRAN package of Ref. [82] and its weight-five extension provided by the authors.