The non-first-order-factorizable contributions to the three-loop single-mass operator matrix elements $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$

The non-first-order-factorizable contributions (The terms 'first-order-factorizable contributions' and 'non-first-order-factorizable contributions' have been introduced and discussed in Refs. \cite{Behring:2023rlq,Ablinger:2023ahe}. They describe the factorization behaviour of the difference- or differential equations for a subset of master integrals of a given problem.) to the unpolarized and polarized massive operator matrix elements to three-loop order, $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$, are calculated in the single-mass case. For the $_2F_1$-related master integrals of the problem, we use a semi-analytic method based on series expansions and utilize the first-order differential equations for the master integrals which does not need a special basis of the master integrals. Due to the singularity structure of this basis a part of the integrals has to be computed to $O(\varepsilon^5)$ in the dimensional parameter. The solutions have to be matched at a series of thresholds and pseudo-thresholds in the region of the Bjorken variable $x \in ]0,\infty[$ using highly precise series expansions to obtain the imaginary part of the physical amplitude for $x \in ]0,1]$ at a high relative accuracy. We compare the present results both with previous analytic results, the results for fixed Mellin moments, and a prediction in the small-$x$ region. We also derive expansions in the region of small and large values of $x$. With this paper, all three-loop single-mass unpolarized and polarized operator matrix elements are calculated.


Introduction
The heavy-flavor contributions both to the unpolarized and polarized deep-inelastic structure functions form an essential part of these quantities.Their scaling violations are different from those of the massless contributions.Since the experimental precision reached the 1% level in the unpolarized case with HERA [3], which will also be the case for polarized deep-inelastic scattering at EIC [4] and the proposed LHeC [5], the three-loop heavy-flavor corrections are needed in the QCD analysis of these data.
In a previous paper [2], we have calculated the first-order factorizable contributions to the constant parts of the unrenormalized three-loop massive operator matrix elements (OMEs) A Qg ), which are based on 1009 of a total of 1233 Feynman diagrams.For all contributions at least 1000 non-vanishing Mellin moments are known.Furthermore, 15 of the 25 color-ζ contributions of the OMEs have been calculated by using the method of arbitrarily high Mellin moments [6].Here ζ n , n ∈ N, n ≥ 2 denote the values of Riemann's ζ-function at integer values of n.Their associated recurrences were computed by using the guessing method [7,8] and solved using the package Sigma [9,10] in all first-order-factorizing cases.
The major new aspect of the present calculation concerns the contribution of higher transcendental letters in the iterated integrals forming the master integrals given by 2 F 1 -solutions [11] 2 of Heun differential equations [21,22].The corresponding basic integrals have been computed to O(ε 0 ) in Ref. [1].Here ε = D − 4 denotes the dimensional parameter.The other master integrals can be obtained by iterating  and square-root valued letters [2,33] on the former solutions, by which all master integrals can be obtained.Future work will be devoted to this method of iterating non-iterative integrals, including higher transcendental letters.
In the present paper, we are following a different avenue.Here we do not compute the analytic results for the previously mentioned OMEs, but we obtain highly precise semi-analytic results.To derive these results, we start with the coupled first-order differential equation system of master integrals, which is obtained from the integration by parts (IBP) reduction [34,35].The solution of these equations is performed in t-space [36,37] via series expansions around various points which are numerically matched in overlapping regions of convergence.Here t denotes the resummation parameter.The initial conditions are given at t = 0 as the Mellin moments to the required order in ε.In our set of master integrals one has to calculate up to O(ε 5 ) in individual cases.The initial values have been computed already before for determining the corresponding recurrences in Mellin-N space, cf.[2,38].After analytic continuation from t to x-space, cf.[1], one obtains the final expression.The analytic continuation has to pass a series of pseudo-thresholds and thresholds from x → ∞ (t = 0) to x = 0 (t → ∞) and matching conditions have to be evaluated.The present approach uses large mantissa rational matching in this process.In this way, we finally obtain the constant parts of the unrenormalized massive OMEs A Qg and ∆a (3) Qg .This formalism is only applied to the part of the amplitude which is affected by 2 F 1 -related letters.
The paper is organized as follows.In Section 2 we describe the basic computation method.In Section 3 we compute a  Qg in x-space, compare to previous partial results in the literature, and present numerical results.The results for small and large values of Bjorken x are presented in Section 4 and Section 5 contains the conclusions.

The main steps of the calculation
The calculation of the contributing Feynman diagrams from their generation to the reduction to the master integrals has been described in Ref. [2].Here we use the packages QGRAF, Form, Color and Reduze 2 [34,35,[39][40][41][42], and apply the Feynman rules given in Refs.[43,44].In the polarized case we compute the OME in the Larin scheme [45].The OMEs are calculated using the method described in Ref. [1].This means the operators defined for discrete integer values of N are resummed into a generating function which depends on the continuous real variable t.The master integrals are computed in this variable by solving linear systems of coupled differential equations, see also Refs.[2,[46][47][48][49].The initial values are provided by the Mellin moments [38], which are the expansion coefficients at t = 0.The system of differential equations is solved at a series of thresholds and pseudo-thresholds in the region t ∈ [0, ∞[.The variables t and the Bjorken variable x are related by The set of thresholds and pseudo-thresholds in the differential equations of all master integrals, resp.the necessary expansion points, given the convergence radius of the respective local series, are for x ∈ [0, 1] in the present problem.Imaginary parts of the amplitude develop only at the transition point x = 1, and for no other points in the regions x ∈]0, 1[ and x ∈]1, ∞[.The expansion points for x > 1 were 3) The analysis starts at x = ∞ using the previously computed moments as initial values.The expansion order in ε depends on the individual master integral.In the present case one needs to expand up to O(ε 5 ) for some integrals.One performs series expansions around the (pseudo-)thresholds by inserting a suitable ansatz into the differential equation.Comparing coefficients in ε, the expansion parameter t − t 0 (and possible powers of logarithms) one obtains a large system of equations for the symbolic expansions.This system of linear equations is solved with FireFly [50,51] using modular methods in terms of a small number of boundary constants.These are determined by matching two neighboring expansions in the middle with 250 digits accuracy.In order to solve this linear system, we rationalize the arising floating point numbers, which allows for a stable solution.Except for the expansions around x = 0 and x = 1, we compute 100 expansion coefficients, while at the latter points, which contain in addition powers of the logarithms ln(x) and ln(1 − x), respectively, 50 expansion terms are used.
The initial values at x → ∞ are real, as are also the coefficients of the linear differential equation systems.For the expansion points for x > 1, the series expansions are real-valued Taylor series in x and their contribution to the massive OMEs vanish, see [1].The analytic continuation at x = 1 implies logarithmic-modulated series containing powers of ln k (1 − x) and thus an imaginary part after analytic continuation.The imaginary part is proportional to the x-space representation of the massive OMEs, see Ref. [1].
In Ref. [2], we computed all color-ζ contributions which can be obtained by solving difference equations that factorize into first-order factors.Furthermore, we calculated all remaining irreducible Feynman diagrams with contributions from master integrals, whose associated differential equations factorize into first-order factors.The remaining 224 Feynman diagrams are related to 2 F 1 -solutions [1,11] and are calculated in the present paper by solving the first order differential equations obtained from the IBP-relations directly in a highly precise numerical approach, adding the previous results to the complete solution.
On top of the irreducible Feynman diagrams mentioned before, there are also reducible Feynman diagrams and ghost contributions to the amplitudes A  Qg (x), contributing to the final result, which we would like to characterize briefly.In the unpolarized case, in Mellin N -space, they are spanned by the harmonic and generalized harmonic sums [32,52,53] with rational and 2 N -prefactors.In x-space, these terms convert to harmonic polylogarithms [54] at argument x or 1−2x, as in Ref. [55].The generalized sums stem from the ghost contributions, which are absent in the polarized case.
Expanding their contribution around x = 0 terms of order are present, which are not expected in the complete result.Here the color factors are ) and N c = 3 for Quantum Chromodynamics (QCD).Indeed, the calculation shows that these terms are canceled at a relative accuracy of in the complete result numerically.Contributions of this kind do not emerge in the polarized case.
Our present results can be tested also in various other ways.Next we compare the result in x-space with the moments computed in Ref. [44] by a totally different method, using MATAD [56].For the moments N = 2, 4, 6, 8, 10 we obtain agreement up to relative accuracies of Qg turns out to be zero, we compare here the relative deviation of the moments N = 3, 5, 7, 9, 11, for which we obtain Qg (x) as a function of x, rescaled by the factor x(1 − x).Left panel: smaller x region.Full line (red): a (3) Qg (x); dashed line (blue): leading small-x term ∝ ln(x)/x [57]; dotted line (green): ln(x)/x and 1/x term; dash-dotted line (black): all small-x terms, including also ln k (x), k ∈ {1, . . ., 5}.Right panel: larger x region.Full line (red): a A further test of accuracy consists in the comparison of the present differential equation method with the analytic results obtained by N -space techniques for the N F terms in a           For smaller values of x, the deviations are even smaller.A comparable accuracy has been obtained for the pole terms of the unrenormalized amplitudes A Qg (x).In Figures 1 and 2 we illustrate the analytic result for a (3) Qg (x) for QCD and by setting N F = 3 in the smaller and larger x regions.Here we add different lines for the small-x and large-x contributions to show their validity.The so-called leading small-x result [57] turns out not to describe the physical quantity a (3) Qg (x) quantitatively; see, however, the discussion in Section 4. This is, as in all known other cases, see e.g.[55,[58][59][60][61], due to sub-leading terms which cancel the leading behaviour.
Here the inclusion of the 1/x term, not predicted by small-x methods, leads to a description up to x ∼ 10 −4 .To describe the region to x ∼ 2 • 10 −2 one needs also all contributing ln k (x)terms, with k ∈ {1, . . ., 5}.In the large-x region no expansion terms have been predicted.Here one obtains a description down to x ∼ 0.9 by considering all ln k (1 − x) terms for k ∈ {1, . . ., 4}, including the (1 − x) 0 and (1 − x) contributions.
In Ref. [62] estimates on the size of the charm quark contributions in F 2 (x, Q 2 ) were made based on five moments for A Qg and six moments for A (3),PS Qq calculated in Ref. [44], the two-loop contributions of Refs.[63,64], and the N F -terms from our calculation in Ref. [65].Furthermore, the small-x behaviour from [57] for Qg and a corresponding color-rescaled leading small-x term for A (3),PS Qq were assumed.The latter has only later been proven in Ref. [55] by calculating A (3),PS Qq in complete form analytically.In [62] the three other contributing OMEs A ,PS qq,Q , as well as the two-mass corrections, were not taken into account.In Figure 2 we illustrate the former estimate on a Qg (x) in the region of smaller and larger values of x (gray band) and compare it to the exact result (red lines), which lie close to the upper end of the former estimate in the region of small values of x.
Let us now turn to the polarized case.The quantity ∆a Qg (x) is shown in Figure 3, where we also indicate the small-and large-x terms, cf.Section 4. We note that the latter approximations match the exact result only in the extreme regions.∆a (3) Qg (x) shows an oscillatory behaviour as also known for the polarized structure functions g 1,2 (x, Q 2 ).One reason for this behaviour is that the first moment of ∆a (3) Qg (x) vanishes.Also at the first and second order in the strong coupling constant the first moment vanishes, cf.Ref. [66]. 4 The small-and large-x limits In this section, we discuss in more detail the small-and large-x limits which have already been illustrated in Section 3. The leading small-x contribution to the unpolarized quantity a (3) Qg (x) has been predicted in Ref. [57], within a leading order calculation based on k ⊥ -factorization.The result is given by As we saw in Section 3, this result is interesting for the theoretical comparison to the corresponding term in the complete calculation, but cannot be used for phenomenology due to the destructive sub-leading corrections.We obtained the term ∝ ζ 2 in Ref. [2] since it results from first order factorizing contributions only.From the small-x expansion of the present result, we obtain an agreement on the purely rational and ζ 3 term of Eq. ( 4.1) at a relative accuracy of This is the first independent recalculation of the result of Ref. [57] using a different method, and it also establishes the rescaling to the corresponding analytic result in the pure-singlet case [55].
Since we knew the coefficient of the ζ 2 -term, cf.[2], in (4.1), the question is whether integer relations allow to determine the other two terms in (4.1) at the level of 15 known digits.At least there could be a conditional answer in assuming a certain rational prime factor pattern, starting out with 2, 3, 5, 7 to reasonably small powers watching out for matches by using the LLL algorithm [67] lindep in the package Pari [68].David Broadhurst has been so kind to do this for us.He obtained very quickly, by which (4.1) can be considered to be confirmed using methods of experimental mathematics. 3he small-x terms for a Qg are given numerically by The alternating sign of the first two coefficients is the main reason why the leading small-x contributions are not sufficient to describe the final result even at very small values of x.For a precise description also the sub-leading terms are needed.In the large-x limit a Qg is given by a Here the coefficients of the terms ln(x)/x and 1/x have been shown to be zero in [2].The N F term ∝ ln 5 (x) has been derived in Ref. [2] as C F T 2 F N F (4/3) ln 5 (x).The coefficients in (4.6) are alternating, except for the last term.(4.7) The large-x expansions, Eqs.(4.5, 4.7) are the same in the unpolarized and polarized case.The same behaviour has already been seen for the factorizing contributions before, cf.[2].The results in the unpolarized and polarized cases were obtained by separate calculations.

Conclusions
We have calculated the non-first-order-factorizable contributions to the three-loop massive operator matrix elements A Qg and ∆A (3) Qg in the single-mass case.This completes the computation of these matrix elements and thereby of all of the three-loop single-mass unpolarized and polarized OMEs [2,55,65,[70][71][72][73][74][75][76].Also the two-mass three-loop corrections [77][78][79][80][81][82], except those for (∆)A (3) Qg , have already been computed.The solution of the first-order differential equation system of master integrals in different sub-intervals of x ∈]0, ∞[ at very high numerical precision and high precision matching using the methods [1,46] allowed us to derive the three-loop corrections tied up to iterated non-iterative integrals containing 2 F 1 -letters in terms of local series expansions.The latter are logarithmic-modulated with powers of ln(x) around x = 0, and ln(1 − x) around x = 1.
We confirm the leading small-x prediction for the O(ln(x)/x) term in the unpolarized case in an independent calculation using a different method for the first time.We compared our results with the moments of Ref. [44] and other terms, which were calculated by us using different methods and found agreement.
The present results are important for future measurements of the strong coupling constant α s (M 2 Z ) [83][84][85][86], the charm quark mass, m c , [87], and the parton distribution functions, see e.g.[88,89].All the three loop single-and two-mass corrections to deep-inelastic scattering will be released in form of a numerical code in a forthcoming publication.
Qg (x) before, where N F denotes the number of massless flavors.

( 3 )
Qg (x) as a function of x, rescaled by the factor x(1−x).Left panel: smaller x region.Full line (red): a

( 3 )
Qg (x) as a function of x, rescaled by the factor x(1 − x).Left panel: full line (red): ∆a