Massive 2- and 3-loop corrections to

hard scattering processes in QCD

zur Erlangung des Doktorgrades an der Fakultät

für Mathematik, Informatik und Naturwissenschaften

Fachbereich Physik

der Universität Hamburg

vorgelegt von

Marco Saragnese

Gutachter der Dissertation: 					    Prof. Dr. habil. Johannes Blümlein

Prof. Dr. Sven-Olaf Moch

Zusammensetzung der Prüfungskommission:				    Prof. Dr. habil. Johannes Blümlein

Prof. Dr. Sven-Olaf Moch

Prof. Dr. Bernd Kniehl

Prof. Dr. Katerina Lipka

Prof. Dr. Gleb Arutyunov

Vorsitzender der Prüfungskommission:				    Prof. Dr. Bernd Kniehl

Datum der Disputation:								    14 Juli 2022

Vorsitzender Fach-Promotionsausschusses Physik:		    Prof. Dr. Wolfgang Parak

Leiter des Fachbereichs Physik:						    Prof. Dr. Günter H. W. Sigl

Dekan der Fakultät MIN:								    Prof. Dr. Heinrich Graener

§ ABSTRACT

This thesis deals with calculations of higher-order corrections in perturbative quantum chromodynamics (QCD). The two-mass contributions to the 3-loop, polarized twist-two operator matrix elements (OMEs) $A_{Qq}^{(3),PS}$ and $A_{gg,Q}^{(3)}$ are calculated. The $N$-space result for $A_{gg,Q}^{(3)}$ is obtained analytically as a function of the quark mass ratio, which for $A_{Qq}^{(3),PS}$ is not yet possible. In the $z$-space representation, one obtains for both matrix elements semi-analytical representations in terms of iterated integrals, whereby for reasons of efficiency an additional integral is necessary for some terms.

These universal (process-independent) massive OMEs govern the asymptotic behaviour of the Wilson coefficients in deep-inelastic scattering at large virtualities $Q^2\gg m_{c,b}^2$, with $m_{c,b}$ the charm and bottom quark masses. These corrections are also required to define the variable flavour number scheme. This scheme describes the transition from massive quark corrections to the massless ones for very high momentum scales, which is relevant to the description of collider data.

In the single-mass, polarized case, we derive the logarithmic corrections for the Wilson coefficients of the structure function $g_1$ in the asymptotic region $Q^2\gg m_{c,b}^2$. This is done using the known OMEs and massless Wilson coefficients, using the renormalization group equations.

For the non-singlet structure functions $F_2^{NS}$ and $g_1^{NS}$ we revisit the scheme-invariant evolution operator known for massless quarks and extend it to the massive case with single- and two-mass corrections. In this case, the evolution can effectively be described up to $\mathcal O(a_s^3)$ in the Wilson coefficients, where $a_s=\alpha_s/(4\pi)$ denotes the strong coupling constant. The influence of the hitherto not fully known 4-loop non-singlet anomalous dimension can be described effectively. It turns out that the effect of the theory error in question can be completely controlled. A representation by a Padé approximant proves to be sufficient.

We consider the class of functions of multivariate hypergeometric series and study systems of differential equations obeyed by them. We describe an algorithmic method to solve some classes of such differential systems which delivers a hypergeometric series solution having nested hypergeometric products as summand; we discuss the relationship between these products and Pochhammer symbols. For a number of classical hypergeometric series we derive differential systems and their associated difference equations. We present some examples of series expansions of such functions and of the mathematical objects which arise therein. We also present a Mathematica package which implements algorithms related to the solution of partial linear difference equations, focusing in particular on bounding the degree of the denominator of solutions which are rational functions. These methods are of particular importance when solving multi-leg calculations for Feynman diagrams, but also come into play when hypergeometric methods for multi-loop integrals are used.

We describe a numerical implementation of an $N$-space library for the calculation of scaling violations for  structure functions, which can perform the evolution of parton distribution functions up to NNLO from a parametrization chosen by the user, and encodes massless and massive Wilson coefficients for the structure functions $F_2$ and $g_1$ in the case of photon exchange, and for the structure functions $F_3^{W^+\pm W^-}$ in the case of charged-current exchange. The library contains analytic continuation of the relevant harmonic sums in Mellin-space up to weight 5 and many weight-6 harmonic sums.
The numerical representation in $x$ space is performed by contour integration around the singularities of the solution of the evolution equations in $N$ space.

§ ZUSAMMENFASSUNG

Die vorliegende  Arbeit beschäftigt sich mit Berechnungen von Korrekturen höherer Ordnung in
der perturbativen Quanten Chromodynamik (QCD). Es werden die
zweimassigen Beiträge zu den 3-loop massiven, polarisierten twist–2 Operatormatrixelementen (OMEs),
$A_{Qq}^{(3),PS}$ und $A_{gg,Q}^{(3) }$, berechnet. Das $N$-Raum-Ergebnis für $A_{gg,Q}^{(3)}$
erhält man analytisch als Funktion des Massenverhältnisses der schweren Quarks, was für
$A_{Qq}^{(3),PS}$ nicht durch iterierte Inegrale möglich ist. In der sog. $z$-Raum Darstellung
erhält man für beide Matrixelemente
analytische Darstellungen durch iterierte Integrale, wobei aus Effizienzgründen hierüber für
manche Terme ein weiteres Integral notwendig ist.

Diese universellen (prozeßunabhängigen) massiven OMEs bestimmen das asymptotische
Verhalten der Wilson-Koeffizienten bei tief–inelastischer Streuung für große Virtualitäten
$Q^2\gg m_{c,b}^2$. Hier bezeichnen  $m_{c,b}$ die Charm- und
Bottom-Quark-Masse. Diese Korrekturen sind auch erforderlich, um das Variable Flavor Number
Scheme zu definieren. Dieses Schema beschreibt für sehr hohe Impulsskalen den Übergang massiver
Quark–Korrekturen in den masslosen Fall, was für die Beschreibung von Kollider–Daten von
Bedeutung ist.

Im einmassigen Fall leiten wir die logarithmischen Korrekturen für die
Wilson-Koeffizienten der Strukturfunktion $g_1$ in der asymptotischen Region $Q^2\gg m_{c,b}^2$ ab.
Dies geschieht unter Verwendung der bekannten OMEs und der masselosen Wilson-Koeffizienten, unter
Verwendung der Renormierungsgruppengleichungen.

Für die nicht-singulett Strukturfunktionen $F_2^{NS}$ und $g_1^{NS}$ berechnen wir den
schema-invarianten Evolutionsoperator, der für masslose Quarks bekannt war und erweitern ihn für
den massiven Fall mit ein- und zwei-massigen Korrekturen. Hierdurch kann die Evolution in diesem
Fall effektiv bis zur $O(a_s^3)$ in den Wilsonkoeffizienten beschrieben werden. Hierbei bezeichnet
$a_s = \alpha_s/(4\pi)$ die starke Kopplungskonstante.
Der Einfluss der bisher nicht vollständig bekannten 4–loop nicht-singulett anomalen Dimension
kann effektiv beschrieben werden. Es stellt sich heraus, daß der Effekt des betreffenden
Theoriefehlers vollständig kontrolliert werden kann. Eine Darstellung durch eine
Padé-Approximation zeigt sich als ausreichend.

Wir betrachten die Klasse der Funktionen multivariater hypergeometrischer Reihen und untersuchen
Systeme von Differentialgleichungen und Differenzengleichungen, welche diese beschreiben.
Wir beschreiben ein algorithmisches Verfahren zur Lösung einiger Klassen solcher
Differentialgleichungssysteme, welche eine hypergeometrische Reihenlösung mit verschachtelten
hypergeometrischen Produkten als Summanden liefert und  diskutieren die Beziehung zwischen
den Strukturen rationaler Monome aus Pochhammer-Symbolen. Für eine Reihe
generalisierter klassischer hypergeometrischer Reihen leiten wir Differentialgleichungssysteme
und die zugehörigen Differenzengleichungen her. Wir stellen einige Beispiele für Reihenentwicklungen
solcher Funktionen und der darin auftretenden mathematischen Objekte vor. Es wird ein
Mathematica-Paket beschrieben, welches Algorithmen implementiert, die sich auf die Lösung partieller
linearer Differenzgleichungen beziehen, wobei der Schwerpunkt
insbesondere auf der Begrenzung des Grades der Nenner von Lösungen liegt, die rationale Funktionen
sind. Diese Methoden haben besondere Bedeutung bei der Lösung von sog. multi-leg Berechnungen
bei Feynman Diagrammen, kommen jedoch auch bei Anwendung der hypergeometrischen Methoden für
multi-loop Diagramme zum Einsatz.

Wir beschreiben eine numerische Implementierung einer analytischen $N$-Raum Bibliothek zur
Berechnung von Skalenverletzungen von Strukturfunktionen, welche die Evolution von
Parton-Verteilungsfunktionen bis zu NNLO aus einer vom Benutzer gewählten Parametrisierung
durchführen kann und masselose und massive Wilson-Koeffizienten
für das Photon für die Strukturfunktionen $F_2$ und $g_1$, im Falle des Photonaustausches, und für
die Strukturfunktion $F_3^{W^+\pm W^-}$ im Falle geladener Ströme beschreibt. Die Bibliothek
enthält im Mellin-Raum analytische Fortsetzungen der relevanten harmonischen Summen bis zum
Gewicht $w=5$. Die numerische Darstellung im $x$-Raum erfolgt durch eine Kontur-Integration um die
Singularitäten der vollkommen analytischen Lösung der Evolutionsgleichungen im $N$ Raum.

§ LIST OF PUBLICATIONS

Chapters of this thesis have been published in part in:

§.§  Journal articles

J. Blümlein and M. Saragnese,
The $\mathrm N^3\mathrm{LO}$ scheme-invariant QCD evolution of the non-singlet structure functions $F_2^{NS}(x,Q^2)$ and $g_1^{NS}(x,Q^2)$,
Phys. Lett. B 820 (2021) 136589
[arXiv:2107.01293 [hep-ph]].

J. Blümlein, A. De Freitas, M. Saragnese, C. Schneider and K. Schönwald,
Logarithmic contributions to the polarized O($\alpha_s^3$) asymptotic massive Wilson coefficients and operator matrix elements in deeply inelastic scattering,
Phys. Rev. D 104 (2021) no.3, 034030
[arXiv:2105.09572 [hep-ph]].

J. Ablinger, J. Blümlein, A. De Freitas, A. Goedicke, M. Saragnese, C. Schneider and K. Schönwald,
The two-mass contribution to the three-loop polarized gluonic operator matrix element $A_{gg,Q}^{(3)}$,
Nucl. Phys. B 955 (2020) 115059
[arXiv:2004.08916 [hep-ph]].

J. Ablinger, J. Blümlein, A. De Freitas, M. Saragnese, C. Schneider and K. Schönwald,
The three-loop polarized pure singlet operator matrix element with two different masses,
Nucl. Phys. B 952 (2020) 114916
[arXiv:1911.11630 [hep-ph]].

§.§  Proceedings

J. Ablinger, J. Blümlein, A. De Freitas, M. Saragnese, C. Schneider and K. Schönwald,
New 2- and 3-loop heavy flavor corrections to unpolarized and polarized deep-inelastic scattering,
[arXiv:2107.09350 [hep-ph]].

§.§  Preprints

J. Blümlein, M. Saragnese and C. Schneider,
Hypergeometric Structures in Feynman Integrals,
[arXiv:2111.15501 [math-ph]].

§ INTRODUCTION

The first experiment in subatomic particle physics was arguably performed by J. Thomson in 1897 [1], see also [2] for historical accounts of the period.  Thomson was the first to isolate an electron beam (a cathodic ray in the language of the time), speculating that it was composed of charged particles. By measuring their charge to mass ratio he could conclude that those particles were distinct from any known ion, therefore discovering the electron.
Further insights on the structure of matter came from Rutherford's experiment culminating in the establishment of the presence of a positively charged nucleus in 1911 [4, 5, 3] and therefore discovering the proton. The famous gold-foil experiment is arguably the first scattering experiment in the modern sense. In 1932, Chadwick [6, 7] discovered the neutron.

The notion of a particle-like interpretation of the photon was introduced in the same era, by Planck's theory of blackbody radiation and Einstein's theory of the photoelectric effect, gaining widespread acceptance following Compton's experiments on photon-electron scattering in 1923 [8].

Given the known components of the nucleus, it was necessary to deduce the existence of a strong force liable to hold the positively charged particles together. The first theory of the strong force is due to Yukawa in 1934 [9]. Yukawa postulated the existence of a new particle acting as a mediator of the strong force and formulated a prediction of its mass. Yukawa's theory proved inadequate to settle the question of the nature of the strong force completely, as in the following two decades, a large number of new particles and antiparticles, both weakly and strongly interacting, were discovered mainly in the study of cosmic rays.

In the 1960s, the spectra and quantum numbers of the known strongly-interacting particles were classified by Gell-Mann and Ne'eman into multiplets based on a $SU(3)$ symmetry (the “Eightfold Way” [10]). Based on this model, Gell-Mann successfully predicted the existence of an undiscovered resonance, the $\Omega^-$, its charge, mass and decay rate. It was found as predicted in 1964 [11].
The success of this model led Gell-Mann [12] and Zweig [13, 14] to the interpretation of hadrons as composite objects, mesons and baryons, composed of two or three quarks. At the time, the model contained three flavours of quarks, $u$, $d$, and $s$, explaining the $SU(3)$ symmetry of the  Eightfold Way. A new quantum number, color, was introduced by Greenberg in 1964 [15] to explain the existence of the $\Delta^{++}(uuu)$, $\Delta^-(ddd)$ and $\Omega^-(sss)$ resonances. Without color as a new quantum number, the wavefunction of these resonances would have been totally symmetric, and hence prohibited for Fermions by Pauli's exclusion principle.
Color was included formally in a Yang-Mills theory by Han and Nambu [16].

Early experimental evidence for the presence of substructure in the hadrons was obtained by measurements of the anomalous magnetic moment of the proton by Frisch and Stern in 1933 [17] and of the neutron by Alvarez and Bloch in 1939 [18, 19]. Hofstadter's team measured the charge distribution inside the nucleons in the 1950s [20] and measured their size to be of the order of $10^{-15}\text{~m}$.

Further evidence for nucleon substructure came from the deep-inelastic scattering experiments performed at SLAC in the 1960s at higher resolution [21, 22, 23, 24, 25, 26, 27]. In these experiments, based on electron-proton scattering, the measured cross-section was found to be incompatible with the proton being point-like, but was compatible with an internal structure formed by three point-like constituents. The hypothesis of the proton having a uniform charge distribution was definitely excluded.
On the theoretical front, before the deep-inelastic scattering experiments of the late 1960s, the existence of quarks was not universally accepted, because the hypotheses of fractional charge and of confinement were considered too arbitrary; the fact that no quark could be observed directly was for many physicists a reason to accept them only as a book-keeping method. Theoretical research was more focused on the deduction of general properties of the $S$-matrix than on perturbation theory or even field theory. The SLAC experiments were crucial for establishing a baseline of experimental observations that a successful theory of the strong interactions would need to explain. Among these were Bjorken scaling [28], the observation that the structure functions are approximately independent of the exchanged momentum (it is intimately connected to asymptotic freedom) and the Callan-Gross relation [29], a relation between the structure functions $F_L$ and $F_2$ of the proton. Both are successfully explained by Feynman's parton model [30, 31].

Kinematically, deep-inelastic scattering is described by the exchanged momentum $q^2=-Q^2$ and by $x=Q^2/2p.q$, with $p$ the nucleon momentum (see Section <ref>). Bjorken scaling refers to the observation that the structure functions, to first approximation, are independent on $Q^2$. The Callan-Gross relation corresponds to the vanishing of the longitudinal structure function, $F_L\ll F_2$. In QCD, these properties are recovered to leading order, and higher-order corrections are calculable in perturbation theory.

Bjorken scaling and the Callan-Gross relation had been predicted in 1969. At the time of Feynman's formulation of the parton model, which provided an explanation of both properties, the existence of quarks started to be commonly accepted in the scientific community.

Further theoretical developments were the proof of renormalizability of Yang-Mills theory by 't Hooft in 1971 [32], the formulation of QCD as a $SU(3)$ Yang-Mills theory by Gell-Mann and Fritzsch in 1972 [33, 34] and the computation of the $\beta$-function by Gross, Wilczek and Politzer in 1973 [35, 36] establishing the firm theoretical footing for asymptotic freedom.

The predictions of QCD are, to first order in the coupling constant, the same as those of the parton model, but become different at higher order. In particular, QCD predicts a pattern of scaling violations: a specific dependence of the experimental observables on the scale (typically the exchanged virtual momentum).
For example, Bjorken scaling is not valid beyond lowest order, and the dependence of the structure functions on $Q^2$ is calculable. Similarly, the Callan-Gross relation is not valid beyond lowest order. This feature of QCD provides a path to very precise tests of the theory.

The calculation of the theoretical ingredients for this type of analysis started when Gross and Wilczek [37, 38] and Georgi and Politzer [39] first computed the anomalous dimensions of the twist-2 operators to leading order. In the framework of QCD, they account for the anomalous scaling of the parton densities.
The Wilson coefficients, which account for the partonic scattering amplitudes, were first computed in 1978 [40], and completed in [41] enabling the calculation of the leading order QCD corrections to the parton model.

The scaling of the parton densities was formulated in $x$-space as a set of coupled integro-differential equations in [42, 43] and, in a fermion-pseudoscalar theory [45, 46, 47, 48], in [44]. The parton densities and the splitting functions are given in this context an intuitive interpretation related to the partonic content of the hadron and to the amplitude for collinear splitting of the partons. The splitting functions $P_{ij}$ are just Mellin transforms of the already known anomalous dimensions $\gamma_{ij}$,

\begin{equation}
\gamma_{ij}(N) = -\int_0^1 dx \ x^{N-1} P_{ij}(x) ~.
\end{equation}

Advances in the formulation of theoretical predictions are also due to the factorization technique whereby the infrared and collinear singularities due to initial state partons are absorbed into the bare parton densities to obtain a finite renormalized parton density. The first instance of this idea is due to Politzer in 1977 [49]; systematic study on factorization theorems has been performed by a number of authors, among whom Amati, Petronzio and Veneziano [50, 51], Libby and Sterman, [52, 53], in several studies [54, 55, 56, 57, 58, 59]; see also the reviews [60, 61].

Higher order results have been obtained for the deeply inelastic process since the pioneering works of the 1960s.
The Wilson coefficients for the structure functions $F_2$ and $F_L$ and the anomalous dimensions have been computed up to three-loop order in [62, 63, 64, 65, 66, 67, 68, 69, 70, 76, 71, 72, 75, 74, 73] in massless QCD. For the polarized case, the anomalous dimensions were calculated at LO in [42, 79, 80], at NLO in [81, 82] and at NNLO in [83, 84, 85]; a calculation and definition of the M-scheme commonly used for polarized studies appears in [86, 85].
Polarized Wilson coefficients for the structure function $g_1$ were computed in [87, 88]. The Wandzura-Wilczek sum rule was presented in [89], see also [90, 91] for a modern perspective on polarized sum rules.
Higher-order results are also available for structure functions in charged-current exchange in massless QCD [92].

The coefficients of the $\beta$-function in QCD are known to NLO [93, 94], NNLO [95, 96], N$^3$LO [97, 98] and to five-loop order [99, 100, 101, 102].

The experimental basis for the Standard Model gained another building block in 1974 with the discovery of the $J/\Psi$ resonance [103, 104]. The new meson was convincingly interpreted as a $c\bar{c}$ resonance in the quark model, where $c$ denotes the charm quark. The existence of a fourth quark had been conjectured earlier in [105, 106] in order to preserve a symmetry between leptons and quarks: this fact contributed, together with the discovery of many other charmed particles, to the general acceptance of the quark model. The bottom quark $b$ was identified in 1975 [107] as a constituent of the $\Upsilon$ meson, and the top quark $t$ was discovered in 1995 [108, 109] at the Tevatron.

In Table <ref> the masses of the six quarks in the $\MS$ scheme are summarized. The $u,d$ and $s$ quarks have a mass which lies in the non-perturbative regime typically characterized by $\Lambda_{QCD}\sim\mathcal O(200\mathrm{~MeV})$. For this reason, in perturbative calculations, they are typically treated as massless. The heavier $c$, $b$ quarks have a mass which is not negligible at the energies probed by scattering experiments and cannot be neglected, while the top quark has a mass much higher still, and is treated as decoupling in most perturbative calculations.

$u$     $2.2_{-0.4}^{+0.5}$ MeV

$d$     $4.7_{-0.3}^{+0.5}$ MeV

$s$     $95_{-3}^{+9}$ MeV

$c$     $1.275_{-0.035}^{+0.025}$ GeV

$b$     $4.18_{-0.03}^{+0.04}$ GeV

$t$     $160.0_{-4.3}^{+4.8}$ GeV

Quark masses quoted in [77]. The $u$, $d$, $s$ masses are $\MS$ masses at $\mu\approx2$ GeV; the $c$, $b$ masses are the running $\MS$ masses. The $\MS$ $t$ mass is quoted in [78].

At intermediate energies, the heavy $c$, $b$ quarks are often treated in a variable flavour number scheme, where, depending on the value of $Q^2$, a definite number of quarks is considered light and attributed a parton density. As $Q^2$ increases, when a threshold for quark production is reached, one more quark is treated as light and attributed a parton density, while the other parton densities are redefined according to a set of matching conditions. This type of prescription enables to formulate predictions for a wide range of values of $Q^2$ across the thresholds for the heavy quarks and to limit the occurrence of large logarithms in the perturbative series, which are instead resummed by the renormalization group equation. Various such variable flavour number schemes (VFNS) have been defined in the literature [110, 111, 112, 114, 118, 116, 115, 119, 113, 117, 120]. The matching conditions in the VFNS have been computed in [121, 122] for the two-mass case.

Experimentally, the largest kinematic range has been achieved by the HERA experiment at DESY [123, 124, 125, 126], for which the extreme values of $x$ and $Q^2$ covered were $10^{-6}<x<0.65$ and $0.45\mathrm{~GeV^2}<Q^2<50000\mathrm{~GeV^2}$. At $Q^2=10\mathrm{~GeV^2}$, the region up to around $x=10^{-4}$ could be probed [126].

Reviews of QCD and deep-inelastic scattering can be found in [127, 128, 129, 130, 131, 132, 133, 61, 134].
Deep-inelastic scattering experiments have been performed with polarized leptons and hadrons since the early experiments at SLAC in the 1970s [135, 136], and more recently at CERN with the EMC experiment[137, 138], at SLAC [139, 140, 141, 142, 143, 144, 145], at CERN [151, 146, 147, 148, 149, 150], CLAS [152, 153] and at the HERMES experiment at DESY [154, 155].
In such experiments, the lepton is longitudinally polarized and the hadron can have longitudinal or transverse polarization. For polarized scattering, a larger number of structure functions than in the unpolarized case appears: in the most general electroweak case, which was studied in [156], five unpolarized structure functions $F_i$ and nine polarized ones $g_i$ contribute to the hadronic tensor; in the case of pure photon exchange, however, only $F_{1,2}$ and $g_{1,2}$ appear, see also [156, 90, 157].

The contribution due to $g_2$ is suppressed, compared to the contribution due to $g_1$, by a factor of $M^2/s$, with $M$ the mass of the hadron and $s$ the energy in the center-of-mass frame, making it harder to measure experimentally, especially in the kinematic range where higher-twist contributions are less pronounced. Because the theoretical treatment of the polarized structure functions in terms of the operator product expansion closely follows that of the unpolarized case, an opportunity exists for the extension to the structure function $g_1$ of many results already known in the literature.

In particular, the operator matrix elements of twist-two operators are important ingredients in the VFNS, where they enter in the redefinition of the parton densities across a quark mass threshold. Matching relations in the VFNS have been studied at $\mathcal{O}(a_s^2)$ in [112] and at $\mathcal{O}(a_s^3)$ in [158] for the single-mass case, where heavy quarks are decoupled one by one under the assumption $Q^2\gg m_b^2,m_c^2$.

However, because the mass ratio between the charm and the bottom quark is $\eta=m_c^2/m_b^2\sim0.1$, it is desirable for precision applications to build a VFNS with the purpose of matching the high-energy region where both quarks are decoupled and the low-energy region where they are active. To this end, it is necessary to compute operator matrix elements involving both quarks in the loops.
In QCD, operator matrix elements will contain contributions from diagrams with two different massive quarks in the loops starting at $\mathcal{O}(a_s^3)$. This causes the emergence, in their calculation, of classes of special functions and iterated integrals which depend on the ratio of the two masses.
The extension of the VFNS to the decoupling of two quarks has been studied in [122] and the renormalization of the OMEs in [159] to $\mathcal{O}(a_s^2)$ and in [158, 121] to $\mathcal{O}(a_s^3)$ for the single-mass and the two-mass case respectively.

The calculation of the matrix elements of twist-two operators has historically been one of the methods by which their anomalous dimensions have been obtained, starting with [38, 39].
The calculation of the unpolarized single-mass OMEs has been performed to $\mathcal{O}(a_s^2)$ in [159, 112, 163, 160, 161, 162] and to $\mathcal{O}(a_s^3)$ in [164, 166, 167, 169, 165, 73, 168, 170], and of the two-mass OMEs in [173, 174, 121, 171, 172]. The polarized OMEs have been computed in [176, 177, 175] and to $\mathcal{O}(a_s^3)$ in [179, 178, 180]. Currently, up to $\mathcal O(a_s^3)$, only the OME $A_{Qg}^{(3)}$ is not fully known, and is likely to fall in the function space of elliptic integrals, at least for some color structures.

From the factorization theorems of QCD, it is possible to obtain the asymptotic form of the Wilson coefficients involving a massive quark, i.e. the coefficients of the logarithmic factors $\ln(Q^2/\mu^2)$ and $\ln(Q^2/m^2)$, from the factorization into massive OMEs and massless Wilson coefficients. This method was pioneered in [112] where the asymptotic two-loop charm contributions to the Wilson coefficients were calculated. Those asymptotic results were found to be in agreement with an analytic calculation in [181] and with the small-$x$ limit derived in [182]. To $\mathcal O(a_s^3)$ the logarithmic terms of massive Wilson coefficients have been calculated in [183, 165, 73, 184].
For charged-current interactions, the massive Wilson coefficients have been computed with the method of massive OMEs in [185, 186, 187], correcting results of [188].

The factorization of mass singularities has also found applications in QED: initial-state radiative corrections to the process $e^+ e^-\rightarrow\gamma^*/Z^*$ have been calculated with the method of massive OMEs [189], which has historically been important to cross-check and to correct a direct computation performed in [190] and to obtain radiator functions to high orders [191, 192, 193]. Nowadays, these QED results allow a very precise determination of the $Z$-boson width and of the asymmetry in $e^+ e^-$ annihilation [194, 195].

This thesis reports on computations which fit into the long-running project of calculating OMEs and heavy-quark corrections to the Wilson coefficients. It is organized as follows:
in Chapter <ref> we review the theoretical basis of deep-inelastic scattering, of the renormalization of the OMEs of twist-two operators, and the definition of the VFNS. We also briefly review the mathematics of nested sums, iterated integrals and of the Mellin transform.
In Chapter <ref>, we describe the calculation of two of the polarized OMEs, namely the contributions to $A_{Qq}^{(3),\text{PS}}$ and $A_{gg,Q}^{(3)}$ with two different quark masses.
For $A_{Qq}^{(3),\text{PS}}$ the result, given in Chapter <ref>, is obtained by Feynman parametrization and Mellin-Barnes integration. The Mellin-Barnes integrals are calculated by the residue theorem, which turns the integral into a sum, which is then treated using techniques in summation theory. The calculation closely follow that of the unpolarized OME [196]. Working in the Larin scheme [197] for the treatment of $\gamma_5$ in dimensional regularization, we computed the OME and compared its poles in the dimensional parameter $\ep$ to the known prediction which is obtainable from the knowledge of OMEs to lower perturbative order and of the renormalization structure of the theory. We could confirm the pole structure of the OME. The constant part in $\ep$ is new and is given in $x$-space in semi-analytic form as iterated integrals over an alphabet which contains root-valued expressions. The $N$-space result is not given, since the respective recurrences are not first-order factorizable and hence the solution falls outside of the function space under consideration. The calculation was published in [196].

In Chapter <ref> we calculate the polarized OME $A_{gg,Q}^{(3)}$ in the Larin scheme. The method is the same as for $A_{Qq}^{(3),\text{PS}}$, and closely follows the unpolarized calculation [172]. Here, we could compute the analytic $N$-space result as well, as a function of the ratio $\eta$ of the squared quark  masses. The OME turns out to be expressible in nested harmonic and binomial sums. We also derive the $x$-space result in semi-analytic form using iterated integrals, in a form suitable for numerical evaluations. We also derive a number of identities which re-express many of these iterated integrals into the more familiar harmonic polylogarithms and multiple polylogarithms, which are less cumbersome for numerical evaluation. We review the mathematical objects in which these results are expressed, namely nested sums containing binomial coefficients involving $\eta$ as summands, and iterated integrals whose alphabet contains square roots, at arguments containing $\eta$. These results were published in [198].

In Chapter <ref> we present the asymptotic form, for $Q^2\gg m^2$, of the polarized single-mass Wilson coefficients for the structure function $g_1$ in the Larin scheme. These single-mass Wilson coefficients were obtained by the factorization theorems which state that the massive corrections can be obtained asymptotically as the product of the massive OMEs and the massless Wilson coefficients. The known OMEs and Wilson coefficients allow us to write the logarithmic terms at $\mathcal O(a_s^3)$ with the exclusion of the constant part. These results were published in [199].
In Chapter <ref> we describe the the scheme-invariant evolution of the structure functions $F_2^{NS}$ and $g_1^{NS}$ in the asymptotic region accounting for the effects of the $c$ and $b$ quarks to N$^3$LO. These results were published in [200].

In Chapter <ref> we describe methods to classify some differential systems of hypergeometric type. These systems are obeyed by multivariate hypergeometric series and are in correspondence with difference equations with shifts, obeyed by the summand, involving one variable only. We discuss such a case and describe an algorithm to solve the system of differential equations and recover the summand as a nested hypergeometric product. We discuss how this algorithm works in the classical cases of functions studied by Appell, Horn,  Lauricella, and Exton, as these functions have been used in the physics literature in the context of calculating Feynman integrals. In [201] we provide a computer algebra package which implements the algorithm and a computer-readable list of these classical functions and of the systems obeyed by them. We also give some examples of how the series expansion of some hypergeometric series can be obtained using the package Sigma by Schneider [202, 203, 204] and of the types of functions arising in those examples.

The calculation of OMEs to high loop order requires the solution and classification of nested sums, arising for instance from the Mellin-Barnes integration method. These sums, which in physics applications can be very numerous and very complicated, need to be solved, or, in other words, classified in a minimal set of simpler objects, the simplest of which are the harmonic sums. This classification is often done in summation theory by deriving and solving difference equations, through a number of techniques falling under the name of telescoping (see [205] for a survey and [206, 208, 207, 160] for applications to OMEs). The package Sigma encodes such techniques for solving univariate difference equations in this context. In principle, the Laporta algorithm [331] also gives rise to difference equations. Motivated by the application to physics of univariate difference equations in deep-inelastic scattering, in Chapter <ref> we review the problem of partial linear difference equations in several variables. To date, only a limited number of algorithms are available towards a solution of such equations. We describe a computer algebra implementation of one known approach which targets the solution space of rational functions, possibly containing harmonic sums or Pochhammer symbols in the numerator. The package has been released in [201].

In Chapter <ref> we describe a numerical library in Fortran which encodes the splitting functions up to NNLO and massless and asymptotic massive Wilson coefficients for the structure functions $F_2$ and $g_1$ for photon exchange, and $F_3^{W^+\pm W^-}$. It also includes Fortran routines for the Wilson coefficients of the polarized and unpolarized Drell-Yan process, as well as for scalar and pseudoscalar Higgs boson production. The library works in $N$-space by encoding the analytic continuation of the Mellin transforms of harmonic polylogarithms, sufficient to evaluate the analytic continuation, through the even or odd moments, of the harmonic sums up to weight 5 and in several cases to weight 6. The library is suitable for the numerical evolution of singlet and non-singlet PDFs given a programmable initial parametrization. It is usable in principle for experimental fits of these structure functions. We study the numerical precision attained over the computation of low integer moments and show the evolution of a test input set of unpolarized and polarized PDFs.

§ BASIC FORMALISM

§.§ Deep-inelastic scattering

Deep-inelastic scattering (DIS) is the scattering process between a lepton, with momentum $k$, and a hadron, with momentum $p$, in a specific kinematic regime (Figure <ref>). The lepton interacts with a quark by exchanging an electro-weak boson and a final state is produced. We denote the momentum of the outgoing lepton by $k'$ and the momentum of the final hadronic state by $p_X$.
In many inclusive experiments, only $k'$, but not $p_X$, is measured. A measurement of also $p_X$ has been first possible in the collider experiments H1 and ZEUS at HERA.

Kinematics in DIS

We define $q=k-k'$ the momentum of the vector boson. In the reference
frame of the hadron, $p=(M,0)$, $k=(E_{k},\vec{k}).$ Neglecting
the electron mass, we have:

\begin{equation}
\end{equation}

It is customary to use the variables

\begin{eqnarray}
Q^2   &=& -q^2 > 0, \\
\nu   &=& \frac{p. q}{M}, \\
x     &=& \frac{Q^2}{2p. q}, \\
y     &=& \frac{p.q}{p.k}, \\
W^2   &=& p_X^2
\end{eqnarray}

to define the kinematics of the process.
The variables $x$ and $y$ are called Bjorken variables.
One refers to deep-inelastic scattering if the process occurs in a kinematic region where $Q^2 \gtrsim 4\text{~GeV}^2$ and $W^2 \gtrsim 4\text{~GeV}^2$ [134], a region where the picture offered by perturbative QCD becomes applicable. Further cuts are typically applied in order to limit the size of power corrections $\mathcal O(M^2/Q^2)$ which would otherwise be visible in the experimental data, but subleading with respect to the logarithmic terms most readily obtainable in QCD.
Target mass corrections, which are one class of such subleading terms, have been discussed in [209, 210] and in the context of experimental fits in [211, 212]. In this thesis, they will be neglected.

Computing the amplitude for the process outlined above, one obtains:

\begin{equation}
i\mathcal{M}=L^{\mu}\frac{-i}{q^{2}}\int dxe^{iq.x}\langle X|J_{\mu}(x)|P\rangle=L^{\mu}\frac{-i}{q^{2}}W_{\mu},
\end{equation}

where $L^{\mu}$ represents the leptonic contribution to the amplitude,
$J^{\mu}$ is the electromagnetic current $J^{\mu}=e_q\bar{q}\gamma^{\mu}q$ in case of pure photon exchange,
and $e_q$ is the quark charge. Squaring the amplitude requires us to
examine the quantities $L^{\mu\nu}$ and $W_{\mu\nu}$. The leptonic
tensor has the form $L^{\mu\nu}\propto\mathrm{Tr}[/\!\!\!k\Gamma^{\mu}/\!\!\!k'\Gamma^{\nu}]$,
$\Gamma^{\mu}$ being the lepton-boson vertex coupling. Because the
hadron is a composite object, such an explicit formula cannot be written
about the hadronic tensor

\begin{align}
W_{\mu\nu}(p,q) & =\sum_{X}(2\pi)^{4}\delta^{4}(p_{X}-p-q)\langle P|J_{\mu}(x)|X\rangle\langle X|J_{\nu}(x)|P\rangle\\
& =\int d^{4}xe^{iq.x}\langle P|J_{\mu}(x)J_{\nu}(0)|P\rangle.\label{eq:w-tensor}
\end{align}

Nevertheless, its form can be constrained by considering its Lorentz
structure, because, due to the Ward identity and the conservation of the electromagnetic current, $q^{\mu}W_{\mu\nu}=q^{\nu}W_{\mu\nu}=0$.
This restriction allows us to write in general

\begin{eqnarray}
W_{\mu\nu}(p,q) & = & \Bigl(-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{q^{2}}\Bigr)F_1(x,Q^2)\nonumber \\
&& + \frac{2x}{Q^2} \Bigl(p_{\mu}-\frac{p.q}{q^2}q_{\mu}\Bigr)\Bigl(p_{\nu}-\frac{p.q}{q^2}q_{\nu}\Bigr)F_2(x,Q^{2}) \nonumber \\
&& + i \varepsilon_{\mu\nu\alpha\beta} \frac{p^\alpha q^\beta}{p.q} F_3(x,Q^2).\label{eq:w-tensor-1}
\end{eqnarray}

The scalar functions $F_i(x,Q^{2})$ are known as structure functions and

\begin{equation}
\ep_{\mu\nu\alpha\beta} = 	\begin{cases}
\text{sign}(\sigma) 	& \text{if~} (\mu\nu\alpha\beta) = \sigma(0123), \\
0						& \text{otherwise},
\end{cases}
\end{equation}

for $\sigma$ a permutation, is the Levi-Civita tensor.

Formula (<ref>) is valid if we restrict the analysis to electromagnetic currents. In the case of unpolarized nucleon targets, the term with the Levi-Civita tensor does not contribute, unless weak interactions are considered. In the case of polarized nucleon targets, the hadronic tensor also acquires an antisymmetric part, and the structure functions $g_1$ and $g_2$ are defined by [91, 156]

\begin{equation}
W_{\mu\nu}^A = i \ep_{\mu\nu\lambda\sigma} \Bigl[ \frac{q^\lambda S^\sigma}{p.q} g_1(x,Q^2) + \frac{q^\lambda }{(p.q)^2} (p.q \ S^\sigma - S.q \ p^\sigma) g_2(x,Q^2) \Bigr],
\end{equation}

with $S$ the nucleon spin 4-vector normalized as $S^2=-M^2$.

In the case of pure photon exchange on unpolarized targets, the structure functions can be mapped to the differential cross section of deep-inelastic scattering by [134]

\begin{equation}
\frac{d^2\sigma^{\gamma,\text{unpol.}}}{dx dy} = \frac{2\pi\alpha^2}{xy Q^2} \Big\{\big[1+(1-y)^2\big] F_2(x,Q^2) -y^2 F_L(x,Q^2) \Big\}
\end{equation}

and therefore can be measured experimentally; they are observables. In the literature, different definitions and normalizations for the structure functions are used; here we follow [134].
For completeness we repeat, for the case of pure photon exchange on unpolarized targets, the relation between the structure functions and the differential cross-section, [91]

\begin{align}
\frac{d^2\sigma^{\gamma,\text{pol.}}(\lambda,\pm S_L)}{dx\ dy} &= \pm 2 \pi s \frac{\alpha^2}{Q^4}
\left [
-2\lambda y \left( 2-y-\frac{2 x y M^2}{s} \right) x g_1(x,Q^2) +
8 \lambda \frac{y x^2 M^2}{s} g_2(x,Q^2)
\right],
\label{eq:dsigpolL} \\
\frac{d^3\sigma^{\gamma,\text{pol.}}(\lambda,\pm S_T)}{dx\ dy\ d\phi} &= \pm \frac{2\alpha^2 s}{Q^4} \sqrt{\frac{M^2}{s}} \sqrt{xy\Bigl(1-y-\frac{xyM^2}{s}\Bigr)} \cos(\theta-\phi)
\nonumber\\&
\times \bigl[-2\lambda xy g_1(x,Q^2) -4\lambda x g_2(x,Q^2) \bigr],
\label{eq:dsigpolT}
\end{align}

with $\alpha$ the fine structure constant, $s$ the energy in the center of mass frame, $\lambda$ the helicity of the incoming lepton, $S_{L,T}$ the spin vector of the longitudinally or transversally polarized nucleon, which are

\begin{align}
S_L &= (0,0,0,M), \\
S_T &= M(0,\cos\theta,\sin\theta,0)
\end{align}

in the nucleon rest frame, $\phi$ is the azimuthal angle.

Instead of $F_1$, in the literature it is common to study the longitudinal
structure function

\begin{equation}
F_L(x,Q^2) = F_2(x,Q^2) - 2x F_1(x,Q^2) .
\label{eq:fl}
\end{equation}

By projecting the hadronic tensor (<ref>) with $g_{\mu\nu}$ and with $p^\mu p^\nu$ and setting $p^2=0$ one can write [175]

\begin{eqnarray}
g^{\mu\nu}  W_{\mu\nu}(p,q) &=& \frac{2-D}{2x}   F_2(x,Q^2) + \frac{D-1}{2x} F_L(x,Q^2) , \\
p^\mu p^\nu W_{\mu\nu}(p,q) &=& \frac{Q^2}{8x^3} F_L(x,Q^2) ,
\end{eqnarray}

\begin{equation}
\end{equation}

the dimensions of spacetime. These are inverted as

\begin{eqnarray}
F_2(x,Q^2) &=& \frac{2x}{D-2} \Big[(D-1)\frac{4x^2}{Q^2} p^\mu p^\nu W_{\mu\nu}(p,q) - g^{\mu\nu}W_{\mu\nu}(p,q) \Big]~, \\
F_L(x,Q^2) &=& \frac{8x^3}{Q^2} p^\mu p^\nu W_{\mu\nu}(p,q)~.
\end{eqnarray}

The process known as deep-inelastic scattering refers to the kinematic
region of $Q^{2}\rightarrow\infty$ as $x$ is kept finite, the Bjorken
limit [28]. In this region it is possible to apply the
methods of perturbative QCD, due to the asymptotic freedom of the

§.§ Light-cone dominance

Consider [133] the quantity

\begin{equation}
\int d^{4}xe^{iq.x}\langle P|J_{\nu}(0)J_{\mu}(x)|P\rangle\label{eq:current-product-is-zero}
\end{equation}

where the two currents have been interchanged with respect to their order
in $W_{\mu\nu}$. In the physically allowed region, $q^{0}=E-E'>0$.
It can be shown that in this region the quantity (<ref>)
is zero. Inserting a complete set of states, one obtains

\begin{equation}
\int d^{4}xe^{iq.x}\langle P|J_{\nu}(0)J_{\mu}(x)|P\rangle=\sum_{X}(2\pi)^{4}\delta^{4}(q-p+p_{X})\langle P|J_{\mu}(x)|X\rangle\langle X|J_{\nu}(x)|P\rangle.
\end{equation}

One can prove that if $q^{0}>0$ then $q-p+p_{X}\neq0$: in the rest
frame of the proton, assuming the equality holds,

\begin{equation}
(q-p)^{2}=p_{X}^{2}\rightarrow q^{2}-2q^{0} M+M^{2}=p_{X}^{2}\rightarrow q^{0}=\frac{1}{2M}(q^{2}+M^{2}-p_{X}^{2}).
\end{equation}

In the physically allowed region, $q^{2}<0$ and $p_{X}^{2}>M^{2}$,
so it is impossible to have $q^{0}>0$.

As a consequence we can rewrite (<ref>) as

\begin{equation}
W_{\mu\nu}=\int d^{4}xe^{iqx}\langle P|[J_{\mu}(x)J_{\nu}(0)]|P\rangle.
\end{equation}

Because of causality, the commutator must vanish for $x^{2}<0$. Additionally,
it can be shown that in the deep-inelastic scattering limit where
$q^{2}\rightarrow-\infty$ and $Q^{2}/2p.q\rightarrow\textrm{constant}$,
the dominant contribution to $W_{\mu\nu}$ is due to the region $0\leq x^{2}<1/Q^{2}$.

This statement is known as light-cone dominance: the hadronic
tensor $W_{\mu\nu}$ receives contributions from the product of currents
$J_{\mu}(x)J_{\nu}(0)$ which are dominated by the region $x^{2}\sim0$.

§.§ The operator product expansion

The product of two composite operators $O(x_{1})O(x_{2})$ can become
singular in certain limits. Wilson [213] considered the
limit $x_{1}\rightarrow x_{2}$ and postulated that such a product, when singular, can be expanded as a linear combination of all other operators $O_{i}$
appearing in the theory which are finite in the limit, with the singular
behaviour encoded in singular coefficient functions $C_{i}(x)$:

\begin{equation}
O(x_{1})O(x_{2})\rightarrow\sum_{i}C_{i}(x)O_{i}(x)\qquad\textrm{as }x_{1}\rightarrow x_{2}.
\end{equation}

For the application to deep inelastic scattering [214, 215, 216],
such an expansion is needed for the light-cone region of the product
of two currents:

\begin{align}
J_{\mu}(x)J_{\nu}(0) & \rightarrow g_{\mu\nu}\left(\frac{\partial}{\partial x}\right)^{2}\sum_{i,n}C_{i,1}^{(n)}(x^{2})x^{\mu_{1}}\cdots x^{\mu_{n}}O_{\mu_{1}\cdots\mu_{n}}^{(i)}(0)\nonumber \\
& +\frac{1}{x^{2}}\sum_{i,n}C_{i,2}^{(n)}(x^{2})x^{\mu_{1}}\cdots x^{\mu_{n}}O_{\mu\nu\mu_{1}\cdots\mu_{n}}^{(i)}(0)+\cdots\label{eq:ope}
\end{align}

In a massless, free-field theory, it is possible to find the behaviour
of the singular functions $C_{i}(x^{2})$ around $x^{2}=0$ by a power-counting
argument: call $d_{J}$ the mass dimension of the current $J(x)$
and $d_{O}^{(i)}(n)$ that of $O_{\mu_{1}\cdots\mu_{n}}^{(i)}$. Then,
from (<ref>) it follows that the mass dimension of $C_{i}(x^{2})$
is $[C_{i}(x^{2})]=2d_{J}+n-d_{O}^{(i)}(n)$ and

\begin{equation}
\end{equation}

Thus, the operators corresponding to the minimum value of

\begin{equation}
\tau_{n}^{(i)}=d_{O}^{(i)}(n)-n,
\end{equation}

a quantity called twist, will be dominant in the light-cone expansion. The dominant operators $O_{\mu_{1}\cdots\mu_{n}}^{(i)}$ therefore have twist 2. They are traceless and symmetric, and have definite spin $n$ and dimension $n+2$. The operators can be explicitly written:

\begin{align}
O_{q;\mu_1\cdots\mu_n}^{S}(x) & =i^{n-1} \mathbf S \bigl[ \overline{\psi}(x)\gamma_{\mu_1}D_{\mu_2}\cdots D_{\mu_n}\psi(x) \bigr] - \text{trace terms} \nonumber \\
O_{q;\mu_1\cdots\mu_n}^{NS,(i)}(x) & =i^{n-1} \mathbf S \bigl[ \overline{\psi}(x)\gamma_{\mu_1}D_{\mu_2}\cdots D_{\mu_n}\lambda^{(i)}\psi(x) \bigr] - \text{trace terms} \label{eq:operators}\\
O_{g;\mu_1\cdots\mu_n}^{S}(x) & =i^{n-1} \mathbf {S\,S_p} \bigl[ F_{\mu_1\nu}^a(x) D_{\mu_2}\cdots D_{\mu_{n-1}}F_{\mu_n}^a{}^{\nu}(x) \bigr] - \text{trace terms} \nonumber
\end{align}

where $\mathbf S$ stands for a symmetrization in the Lorentz indices

\begin{equation}
\mathbf S f_{\mu_{1}\cdots\mu_{n}}=\frac{1}{n!}(f_{\mu_{1}\cdots\mu_{n}}+\textrm{permutations}),
\end{equation}

$\mathbf{S_p}$ is  the trace over $SU(N_c)$ and $\lambda^{(i)}$ are the $SU(N_f)$ generator matrices if we assume the theory to have such a flavour symmetry. In equations (<ref>), the subtraction of trace terms contains factors of $g_{\mu_i\mu_j}$; it is needed in order to make the operators have definite spin. Here, $\psi(x)$ and $F_{\mu\nu}$ denote respectively the quark field and the electromagnetic field strength, and $D_\mu$ is the covariant derivative,

\begin{eqnarray}
D_\mu &=& \partial_\mu-igt^aA_\mu^a,\\
F_{\mu\nu}^a &=&\partial_\mu A_\nu^a-\partial_\nu A_\mu^a +gf^{abc}A_\mu^b A_\nu^c,
\end{eqnarray}

where $t^a$ are the generators of $SU(3)$ in the fundamental representation and $f^{abc}$ are the structure constants.

In the case of polarized scattering, the contributing twist-two operators are [81, 79]

\begin{align}
O_{q;\mu_1\cdots\mu_n}^{S,5}(x) & =i^n \mathbf S \bigl[ \overline{\psi}(x)\gamma_5 \gamma_{\mu_1}D_{\mu_2}\cdots D_{\mu_n}\psi(x) \bigr] - \text{trace terms} \nonumber \\
O_{q;\mu_1\cdots\mu_n}^{NS,(i),5}(x) & =i^n \mathbf S \Bigl[ \overline{\psi}(x) \gamma_5 \gamma_{\mu_1}D_{\mu_2}\cdots D_{\mu_n} \frac{\lambda^{(i)}}{2} \psi(x) \Bigr] - \text{trace terms} \label{eq:pol-operators}\\
O_{g;\mu_1\cdots\mu_n}^{S,5}(x) & =i^n \mathbf {S\,S_p} \Bigl[ \frac{1}{2} \varepsilon^{\mu_1\alpha\beta\gamma} F_{\beta\gamma}^a(x) D_{\mu_2}\cdots D_{\mu_{n-1}}F_{\alpha\mu_n}^a(x) \Bigr] - \text{trace terms} .\nonumber
\end{align}

§.§ The forward Compton amplitude

The optical theorem states that the imaginary part of an amplitude can be related to the amplitude for scattering into all possible final states. The theorem can be applied [215] in the case of the hadronic tensor to relate it to forward virtual Compton scattering, whose amplitude is determined by

\begin{equation}
T_{\mu\nu}=i\int d^{4}x\,e^{iq.x}\langle P|T[j_{\mu}(x)j_{\nu}(0)]|P\rangle.\label{eq:t-mu-nu-definition}
\end{equation}

The theorem states that

\begin{equation}
\frac{1}{2\pi}\textrm{Im}(T_{\mu\nu})=W_{\mu\nu},
\end{equation}

since $\textrm{Im}(T_{\mu\nu})$ is equal to the discontinuity of the tensor in the $q^{0}$ plane. The forward Compton scattering can be more readily computed in a perturbative expansion, and can be related to the hadronic tensor by a dispersion integral: by calling

\begin{equation}
\omega=\frac{1}{x}=\frac{2p.q}{Q^{2}},
\end{equation}

we have

\begin{equation}
2\int_{0}^{1}dx\,x^{n-2}W_{\mu\nu}=\frac{2}{\pi}\int_{1}^{\infty}\frac{d\omega}{\omega^{n}}\text{Im}(T_{\mu\nu})=\frac{1}{2\pi i}\int_{C}d\omega\frac{T_{\mu\nu}}{\omega^{n}},\label{eq:mellin-w}
\end{equation}

and $C$ is a contour that circles around the branch cuts of $T_{\mu\nu}$ in the $\omega$ plane. Noting the formula

\begin{equation}
\frac{1}{2\pi i}\int_{C}d\omega\,\omega^{m-n}=\delta_{m,n-1}
\end{equation}

it follows that Eq. (<ref>) will pick one power in a series expansion of $T_{\mu\nu}$.

To the forward Compton scattering tensor $T_{\mu\nu}$ one can apply arguments related to an operator product expansion in a way analogous to those applicable to $W_{\mu\nu}$. We report here the results of this operator product expansion as presented in the textbook by Muta [133], see also [217]:

\begin{align}
T[j_{\mu}(x)j_{\nu}(x')]= & (\partial_{\mu}\partial_{\nu}'-g_{\mu\nu}\partial.\partial')O_{L}(x,x')\nonumber \\
& +(g_{\mu\lambda}\partial_{\rho}\partial_{\nu}'+g_{\rho\nu}\partial_{\mu}\partial_{\lambda}'-g_{\mu\lambda}g_{\rho\nu}\partial.\partial'-g_{\mu\nu}\partial_{\lambda}\partial_{\rho}')O_{2}^{\lambda\rho}(x,x')\\
& +\textrm{ antisymmetric terms.}\nonumber
\end{align}

The antisymmetric terms only contribute to polarized scattering processes.
In analogy with (<ref>), the operators have the light-cone expansion

\begin{align}
O_L(x,x') & =\sum_{i,n}C_{L,n}^{(i)}(y^2)y^{\mu_1}\cdots y^{\mu_n}O_{L\mu_1\cdots\mu_n}^{(i)}\Big(\frac{x+x'}{2}\Big)\\
O_2^{\lambda\rho}(x,x') & =\sum_{i,n}C_{2,n}^{(i)}(y^2)y^{\mu_1}\cdots y^{\mu_n}O_{2\mu_1\cdots\mu_n}^{(i)\lambda\rho}\Big(\frac{x+x'}{2}\Big)
\end{align}

where $y=x-x'$, and $i$ runs over the operators in (<ref>). Applying the expansion to (<ref>),
one obtains

\begin{equation}
T_{\mu\nu}=2\sum_{i,n}\omega^n\left[\left(g_{\mu\nu}-\frac{q_\mu q_\nu}{q^2}\right)A_{L,n}^{(i)}C_{L,n}^{(i)}(Q^2)+d_{\mu\nu}A_{2,n}^{(i)}C_{2,n}^{(i)}(Q^2)\right],\label{eq:t-tensor-a-c}
\end{equation}

\begin{equation}
d_{\mu\nu} = - g_{\mu\nu} - \frac{p_\mu p_\nu}{(p.q)^2} q^2 + \frac{p_\mu q_\nu+p_\nu q_\mu}{p.q} ,
\end{equation}

\begin{equation}
C_{(2,L),n}^{(i)}(Q^{2})q^{\mu_{1}}\cdots q^{\mu_{n}}/(q^{2})^{n}=i\int d^{4}x\,e^{iq.x}C_{(2,L),n}^{(i)}(x^{2})
\end{equation}

\begin{equation}
A_{L,n}^{(i)}p_{\mu_1}\cdots p_{\mu_n}+\textrm{terms with }g_{\mu_i\mu_j}=\langle P|O_{L\mu_1\cdots\mu_n}^{(i)}(0)|P\rangle,
\end{equation}

with a similar formula for $A_{2,n}^{(i)}$.

From equations (<ref>), (<ref>),
(<ref>), (<ref>), and applying the replacement
$n\rightarrow N+1$ to harmonize different conventions present in
the literature, one finally deduces the master formula

\begin{equation}
\int_0^1 dx\,x^{N-1}F_{(2,L)}(x,Q^2) = \sum_i A_{(2,L),N}^{(i)} \Bigl( \frac{p^2}{\mu^2} \Bigr) \, C_{(2,L),N} ^{(i)}\Bigl(\frac{Q^2}{\mu^2}\Bigr) ~.\\
\label{eq:fact-struct-func}
\end{equation}

The operation on the left-hand side is called the Mellin transform,
and can be inverted to give an explicit formula for the structure

\begin{equation}
F_{(2,L)}(x,Q^{2})=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}dN\,x^{-N}\sum_iA_{(2,L),N}^{(i)}\Bigl( \frac{p^2}{\mu^2} \Bigr) C_{(2,L),N}^{(i)}\Bigl(\frac{Q^2}{\mu^2}\Bigr).
\label{eq:f-inv-mellin}
\end{equation}

An interesting property with practical consequences for actual computations [91, 218] is the symmetry
In the unpolarized case, this has the consequence that $C_{(2,L),N}(Q^{2})=0$ for odd $N$. In the case of scattering over a polarized target, instead, $C_{g_1,N}(Q^{2})=0$ for even $N$.

The importance of formula (<ref>) lies in the fact that the functions $C_{N}(Q^{2})$, called Wilson coefficients, depend only on short-distance physics, and, for this reason, can be computed in perturbative QCD because of the property of asymptotic freedom. Asymptotic freedom refers to the fact that the running coupling constant goes to zero at high energy: this makes the computation of physical quantities possible as the expansion parameter $a_s = \alpha_s/(4\pi)$ is small. By contrast, the quantities $A_{N}$, whose Mellin inverse is related to the parton density through a factorization procedure, depend on long-distance physics: they are the matrix elements between the target hadron state of the operators appearing in the operator product expansion.

We will discuss later in greater detail how parton densities are defined and also treat possible collinear and mass singularities present in the operator matrix elements. Broadly speaking, in calculations these singularities arise from the initial-state partons, over which the deep-inelastic scattering process is not inclusive. In reality, the initial state partons are confined inside the hadron, and are not themselves asymptotic states. It is therefore reasonable to assume that the initial state infrared and collinear singularities that arise when we compute with initial-state partons are manifestations of the inapplicability of perturbation theory to the long-range nonperturbative regime, and are actually shielded by non-perturbative physics. Therefore, it is sensible to reabsorb them into the unobservable, bare parton density, which encodes the long-distance physics. After this step, the structure function is postulated to be a Mellin convolution of a finite parton density and a finite Wilson coefficient along the lines of (<ref>), and the renormalized parton density becomes scale dependent. The procedure, first applied in [49], is analogous to what happens when renormalizing the coupling constant.

The parton densities cannot be computed perturbatively and must be extracted from experimental data. This can be accomplished by measuring the structure functions experimentally and computing the Wilson coefficients, then applying (<ref>) in a suitable form.

Still, there is predictive value in factorizing the structure function (or its Mellin moments), as one can assume that the $A_{N}$ are universal (see for a discussion [127, 131]): they will not depend on which scattering process is being studied, and, once they are deduced for a given hadronic target and leptonic probe, they can be used in predictions for processes which involve a different leptonic probe, provided that the virtuality $\mu^2$ does not change sign.

One important remark is the existence of Carlson's theorem [220, 219], which implies that the analytic continuation of the Mellin transform in the complex plane can be derived given the knowledge of the even or odd moments. Extensions and applications of Carlson's theorem to this physical case are derived in [221].

§.§ Scaling violation and renormalization group

Let us describe Eq. (<ref>) in greater detail. To the operators (<ref>) correspond three parton densities, classifiable according to their symmetry properties under the flavour symmetry. The gluon density $G(N,\mu^2)$ corresponds to the gluon operator. Assuming $N_F$ massless quarks, the singlet combination $\Sigma(N,\mu^2)$ and the non-singlet combination $\Delta_k(N,\mu^2)$ are given by [64]

\begin{eqnarray}
\Sigma(N,\mu^2)   &=& \sum_{k=1}^{N_F} [ f_k(N,\mu^2) + f_{\bar k}(N,\mu^2) ] , \label{eq:singletPDF}\\
\Delta_k(N,\mu^2) &=& f_k(N,\mu^2) + f_{\bar k}(N,\mu^2) - \frac{1}{N_F} \Sigma(N,\mu^2) , \label{eq:NSPDF}
\end{eqnarray}

with $f_k$ and $f_{\bar k}$ the densities of quarks and anti-quarks. The scale $\mu^2$ is the factorization scale, which separates the high-energy and the non-perturbative contributions. In the following, it will be set equal to the renormalization scale, although in principle the two scales can be treated separately. The structure functions then can be written as [64]

\begin{eqnarray}
F_{(2,L)}(N_F,N-1,Q^2) &=& \frac{1}{N_F} \sum_{k=1}^{N_F} e_k^2 \Bigl[ \Sigma(N_F,N,\mu^2) C_{q,(2,L)}^S \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr)
\nonumber\\
&& +G(N_F,N,\mu^2) C_{g,(2,L)}^S \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr)
\nonumber\\
&& + N_F \Delta_k(N_F,N,\mu^2) C_{q,(2,L)}^{NS} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr)
\Bigr]~.
\end{eqnarray}

The quarkonic singlet contribution is split into a non-singlet part and a pure-singlet part,

\begin{equation}
C_{q,i}^S = C_{q,i}^{NS} + C_{q,i}^{PS}, \quad i=2,L.
\end{equation}

The perturbative expansion of the Wilson coefficients is as follows:

\begin{eqnarray}
C_{g,i}^S \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) &=& \sum_{k=1}^\infty a_s^k C_{g,i}^{(k),S} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) ~, \\
C_{q,i}^{PS} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) &=& \sum_{k=2}^\infty a_s^k C_{q,i}^{(k),PS} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) ~, \\
C_{q,i}^{NS} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) &=& \delta_{i,2} + \sum_{k=1}^\infty a_s^k C_{q,i}^{(k),NS} \Bigl(N_F,N,\frac{Q^2}{\mu^2}\Bigr) ~,
\end{eqnarray}

for $i=2,L$. We call

\begin{equation}
a_s = \frac{\alpha_s}{4\pi} = \Bigl(\frac{g_s}{4\pi}\Bigr)^2.
\end{equation}

the renormalized strong coupling constant.

In order to investigate the behaviour of the Wilson coefficients in the high-$Q^2$ regime, and to formulate predictions for the dependence of the structure functions on $Q^2$, it is necessary to consider the renormalization of the operator matrix elements.
The operators (<ref>) are renormalized by [38]

\begin{eqnarray}
O_{q;\mu_1\ldots\mu_n}^{NS,(i),\text{bare}} &=& Z^{NS}(\mu^2)   O_{q;\mu_1\ldots\mu_N}^{NS,(i),\text{ren}} , \\[1mm]
O_{i;\mu_1,\ldots,\mu_N}^{S,\text{bare}}    &=& Z_{ij}^S(\mu^2) O_{j;\mu_1,\ldots,\mu_N}^{S,\text{ren}}, \qquad i,j=q,g .
\end{eqnarray}

The renormalization produces a scale dependence in the renormalized operators, whose anomalous dimension is

\begin{eqnarray}
\gamma_{qq}^{NS} &=& \mu \bigl(Z^{NS}(\mu^2)\bigr)^{-1} \frac{\partial}{\partial\mu} Z^{NS}(\mu^2) ~,\\
\gamma_{ij}^S    &=& \mu \bigl(Z_{il}^S(\mu^2)\bigr)^{-1} \frac{\partial}{\partial\mu} Z_{lj}^S(\mu^2) ~.
\end{eqnarray}

Because the structure functions are observables, they must be independent on the renormalization scale $\mu$; therefore their total derivative with respect to $\mu^2$ must vanish,

\begin{equation}
\mathcal D(\mu^2) F_{(2,L)}(N,Q^2) = 0,
\label{eq:RGE}
\end{equation}

\begin{eqnarray}
\mathcal D(\mu^2) &=& 	\mu^2 \frac{\partial}{\partial\mu^2} + \beta\bigl(a_s(\mu^2)\bigr) \frac{\partial}{\partial a_s(\mu^2)} -\gamma_m\bigl(a_s(\mu^2)\bigr) m(\mu^2) \frac{\partial}{\partial m(\mu^2)}, \\
\beta(a_s(\mu^2)) &=& \mu^2 \frac{\partial a_s(\mu^2)}{\partial\mu^2}, \\
\gamma_m(a_s(\mu^2)) &=& -\frac{\mu^2}{m(\mu^2)} \frac{\partial m(\mu^2)}{\partial\mu^2}.
\end{eqnarray}

Eq. (<ref>) is the renormalization group equation [214, 216, 222].
In a massless theory, it follows that [127, 223]

\begin{eqnarray}
\sum_{i=q,g} \Bigl[\Bigl\{ \mu^2 \frac{\partial}{\partial\mu^2} + \beta(a_s(\mu^2)) \frac{\partial}{\partial a_s(\mu^2)} \Bigr\} \delta_{ij} -\frac{1}{2}\gamma_{ij}^{S,NS} \Bigr] C_{i,(2,L)} \Bigl(\frac{Q^2}{\mu^2},a_s(\mu^2)\Bigr) = 0, \\
\sum_{j=q,g} \Bigl[\Bigl\{ \mu^2 \frac{\partial}{\partial\mu^2} + \beta(a_s(\mu^2)) \frac{\partial}{\partial a_s(\mu^2)} \Bigr\} \delta_{ij} +\frac{1}{2}\gamma_{ij}^{S,NS} \Bigr] \langle l|\mathcal O_j(\mu^2)|l \rangle = 0.
\end{eqnarray}

(In the non-singlet case, the indices $i,j$ only take the value $q$.) From these equations one can obtain the evolution equations for the parton densities [37, 38, 39],

\begin{eqnarray}
\frac{\partial}{\partial\ln\mu^2}\begin{pmatrix}
\Sigma(N_F,N,\mu^2) \\
\end{pmatrix}
\begin{pmatrix}
\gamma_{qq} & \gamma_{qg} \\
\gamma_{gq} & \gamma_{gg}
\end{pmatrix}
\begin{pmatrix}
\Sigma(N_F,N,\mu^2) \\
\end{pmatrix}
\\
\frac{\partial}{\partial\ln\mu^2} \Delta_k(N_F,N,\mu^2) &=& -\frac{1}{2} \gamma_{qq}^{NS} \Delta_k(N_F,N,\mu^2)
\end{eqnarray}

and correspondingly

\begin{eqnarray}
\frac{\partial}{\partial\ln\mu^2}\begin{pmatrix}
C_{q,i}^S (N_F,N,Q^2/\mu^2) \\
C_{g,i}   (N_F,N,Q^2/\mu^2)
\end{pmatrix}
\frac{1}{2}
\begin{pmatrix}
\gamma_{qq} & \gamma_{gq} \\
\gamma_{qg} & \gamma_{gg}
\end{pmatrix}
\begin{pmatrix}
C_{q,i}^S (N_F,N,Q^2/\mu^2) \\
C_{g,i}   (N_F,N,Q^2/\mu^2)
\end{pmatrix} ,
\\
\frac{\partial}{\partial\ln\mu^2} C_{q,i}^{NS}(N_F,N,Q^2/\mu^2) &=& \frac{1}{2} \gamma_{qq}^{NS} C_{q,i}^{NS}(N_F,N,Q^2/\mu^2) .
\end{eqnarray}

§.§ Renormalization in the presence of heavy quarks

Because of the wide difference in the quark masses, it is natural to divide the quarks into “light" and “heavy". Typically, the $u$, $d$, and $s$ quarks are considered light, because they have a mass comparable or smaller than $\Lambda_{QCD}$, which separates perturbative and non-perturbative physics. In perturbative calculations, their mass is neglected.
The masses of the $c$, $b$ and $t$ quarks are outside of the non-perturbative regime, and their treatment in calculations will depend on the scale of the process being considered.
In general, the $t$ quark, because of its extremely high mass, is not considered to produce any effects at the lower energies available for DIS experiments.

The possibility to disregard particles much heavier than the scale under consideration is a consequence of the Appelquist-Carazzone theorem [224]. Its physical meaning is that physics at very high energies should not affect the physics at low energies and should not be discernible only from low-energy phenomena. More precisely, it states that if a low-energy effective theory is constructed by removing the heavy-particle fields, then the effect on the Green's functions that involve only light particles is equivalent to a finite renormalization of the couplings, up to terms suppressed by the heavy mass, $\mathcal{O}(1/M)$. In other words, removing the heavy particles does not produce observable effects. (An exception to this rule applies when the low-energy theory has different symmetries than the high-energy one; in such cases the effect of removing the heavy particles will be visible).

In practice, whether it is appropriate to include the effects of the $c$, $b$ quark, or both, will depend on the observable under consideration and the experimental precision available. However, if calculations are performed in the $\MS$ scheme, the decoupling of heavy particles is not manifest order by order. This is because the $\MS$ renormalization prescription is mass-independent. In dimensional regularization, it prescribes the removal of the $\varepsilon$ poles and of the universal spherical factor

\begin{equation}
\end{equation}

for each perturbative order, in $D$ dimensions, with $\gamma_E$ the Euler-Mascheroni constant,

\begin{equation}
\gamma_E = \lim_{n\rightarrow\infty} \Bigl( \sum_{k=1}^n \frac{1}{k} - \ln(n) \Bigr) ~.
\end{equation}

In the $\MS$ scheme, the decoupling of heavy particles is only manifest after all perturbative orders are summed. However, it is highly desirable in practice to adopt a renormalization scheme which exhibits the decoupling order by order. An example of such a scheme is the CWZ prescription [110].

The CWZ prescription corresponds to a set of subschemes related to each other by matching conditions. If the masses of $n_f$ quarks are

\begin{equation}
m_1, \ldots\, m_{n_\ell}, m_{n_\ell+1},\ldots, m_{n_f}
\end{equation}

and the virtuality of the process is $Q\sim\mu\sim m_{n_\ell}$,
the CWZ prescription is to divide the quarks into $n_\ell$ light or “active" and $n_h=n_f-n_\ell$ heavy quarks. The light quarks are considered massless and are renormalized in the $\MS$ scheme. Graphs containing heavy quark lines are renormalized using zero-momentum subtraction, also known as BPHZ. Explictly, it is demanded that $\Pi(p^2, m^2)|_{p=0}=0$ for the heavy quarks, where $\Pi(p^2,m^2)$ refers to the contribution from the heavy quarks to the gluon self-energy.

In the CWZ scheme, the decoupling of the heavy quarks is manifest, and the numerical value of the $\beta$-function is the same as in the effective theory with $n_\ell$ quarks renormalized in the $\MS$ scheme.
When implementing a variable flavour number scheme, as a way of resumming the large logarithms that occur, the choice of which subscheme to use is determined by the virtuality $Q$, which should be of the order of $m_{n_\ell}$.

A one-loop calculation in the CWZ scheme can be found in [225]. Multi-loop calculations can be found in [176, 159, 175].

Let us review the renormalization procedure for the massive OMEs developed in [159, 158, 175] for the single-mass OMEs; the two-mass case was discussed in [121, 178]. This whole section follows the exposition in those References.

The unrenormalized OMEs are first obtained by 2-point functions including self-energies in the external legs, amputated of the external fields. The external legs are kept on-shell, thus avoiding the potential problem of the mixing of non-gauge-invariant operators [228, 229, 230, 231, 75, 227, 226]. We call the momentum flowing through them $p$, with $p^2=0$.

The trace terms present in the twist-two operators are projected out by multiplying by an external source

\begin{equation}
J_{\mu_1\ldots\mu_N} = \Delta_{\mu_1}\cdots\Delta_{\mu_N}~,
\end{equation}

with $\Delta_\mu$ an auxiliary light-like vector, $\Delta^2=0$.

Calling the unrenormalized OMEs, which are denoted by two hats,

\begin{equation}
\hat{\hat A}_{ij} \Bigl(\frac{m^2}{\mu^2},\varepsilon,N\Bigr) = \langle j | O_i | j \rangle~,
\end{equation}

the general structure of these Green's functions is then [159, 112]

\begin{eqnarray}
G_{\mu\nu,l,Q}^{ab} &=& \hat{\hat A}_{lg} \ \delta^{ab} (\Delta.p)^N \Bigl( -g_{\mu\nu}+\frac{\Delta_\mu p_\nu+\Delta_\nu p_\mu}{\Delta.p} \Bigr) ~,\\
G_{l,Q}^{ij} &=& \hat{\hat A}_{lq} \ \delta^{ij} (\Delta.p)^N /\!\!\!\!\Delta ~,
\end{eqnarray}

depending on whether the external states are gluons or quarks. Here $l$ can indicate a light parton or the heavy quark. One projects out the OMEs through

\begin{eqnarray}
\hat{\hat A}_{lq} &=& P_q \ G_{l,Q}^{ij} = \frac{\delta^{ij}}{N_c} (\Delta.p)^{-N} \text{Tr}(/\!\!\!p G_{l,Q}^{ij} ) ~,
\end{eqnarray}

for the case of external quarks, with $N_c$ the number of colors. In the case of external gluons, one should distinguish between the two projectors

\begin{eqnarray}
P_g^{(1)} &=& -\frac{\delta^{ab}}{N_c^2-1}\frac{1}{D-2} (\Delta.p)^{-N} g_{\mu\nu} ~,\\
P_g^{(2)} &=& \frac{\delta^{ab}}{N_c^2-1}\frac{1}{D-2} (\Delta.p)^{-N} \Bigl(-g_{\mu\nu} + \frac{\Delta_\mu p_\nu + \Delta_\nu p_\mu}{\Delta.p} \Bigr) \ .
\end{eqnarray}

The projector $P_g^{(2)}$ enforces the transversity of the gluon polarizations and does not require the inclusion of diagrams with external ghosts. Instead, they must be considered if one chooses to perform the calculation using $P_g^{(1)}$. In either case, one has

\begin{equation}
\hat{\hat A}_{lg} = P_g^{(1,2)} \ G_{\mu\nu,l,Q}^{ab} \ .
\end{equation}

In the case of polarized OMEs, different projectors are used.

The renormalization of the OMEs occurs in four steps: mass renormalization, coupling constant renormalization, operator renormalization and mass factorization. In the first step, the unrenormalized mass $\hat m$ is replaced with the renormalized mass. Considering for illustration the case of one massive quark, the relation reads in the on-mass shell (OMS) scheme [232, 233, 234, 235, 236]

\begin{equation}
\hat m = Z_m m = m \bigg[ 1 + \hat a_s \ \Big(\frac{m^2}{\mu^2}\Big)^{\varepsilon/2} \delta m_1 + \hat a_s^2 \ \Big(\frac{m^2}{\mu^2}\Big)^{\varepsilon} \delta m_2 + \mathcal O(\hat a_s^3) \bigg] ~,
\end{equation}

\begin{eqnarray}
\delta m_1 &=&C_F
\left[\frac{6}{\ep}-4+\left(4+\frac{3}{4}\zeta_2\right)\ep
\right] \label{delm1}  \\
&\equiv&  \frac{\delta m_1^{(-1)}}{\ep}
+\delta m_1^{(0)}
+\delta m_1^{(1)}\ep~, \label{delm1exp} \\
\delta m_2 &=& C_F
\Biggl\{\frac{1}{\ep^2}\Big[18 C_F-22 C_A+8T_F(n_f+N_h)
\Big]
\nonumber\\&&
\nonumber\\&&
\N\\ &&
\left(\frac{69}{2}-14\zeta_2\right)\right]\Biggr\}
\label{delm2}  \\
&\equiv&  \frac{\delta m_2^{(-2)}}{\ep^2}
+\frac{\delta m_2^{(-1)}}{\ep}
+\delta m_2^{(0)}~, \label{delm2exp}
\end{eqnarray}

where $N_F$ is the number of light flavours and $N_H$ that of heavy flavours.
The constants $\zeta_k$ refer to the Riemann $\zeta$-function at integer values,

\begin{equation}
\zeta_k = \sum_{n=1}^\infty \frac{1}{n^k}, \qquad k\in\mathbb N, k\ge 2.
\label{eq:zetak}
\end{equation}

After this replacement one can write the mass-renormalized OMEs as [158]

\begin{eqnarray}
\Ahathat_{ij}\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
&=& \delta_{ij}
\Ahathat_{ij}^{(1)}\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
\N\\ &&
+ \hat{a}_s^2 \left[
\Ahathat^{(2)}_{ij}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
+ {\delta m_1}
\Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
m \frac{d}{dm}
\Ahathat_{ij}^{(1)}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
\right]
\N\\ &&
+ \hat{a}_s^3 \Biggl[
\Ahathat^{(3)}_{ij}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
+{\delta m_1}
\Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
m \frac{d}{dm}
\Ahathat_{ij}^{(2)}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
\N\\ &&
+ {\delta m_2}
\Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep}
m \frac{d}{dm}
\Ahathat_{ij}^{(1)}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
+ \frac{\delta m_1^2}{2}
\Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep}
m^2  \frac{d^2}{{dm}~^2}
\Ahathat_{ij}^{(1)}
\Bigl(\frac{m^2}{\mu^2},\ep,N\Bigr)
\Biggr]~.
\nonumber\\
\label{maren}
\end{eqnarray}

The coupling is renormalized in the $\MS$ scheme as follows:

\begin{eqnarray}
\hat a_s &=& \Big( Z_g^{\MS}(\ep,N_F) \Big)^2 a_s^{\MS}(\mu^2) \\
&=& a_s^{\MS}(\mu^2) \Bigl[ 1 + \delta a^{\MS}_{s, 1}(N_F)
+ \delta a^{\MS}_{s, 2}(N_F)
\big( {a^{\MS}_s}(\mu)\big)^2
+ \mathcal O\Big(\big({a_s^{\MS}}\big)^3\Big) \Bigr].
\label{eq:asMS}
\end{eqnarray}

The coefficients in this expansion are related to the $\beta$-function as follows: in dimensional regularization, the renormalization scale $\mu$ is defined through

\begin{eqnarray}
\hat g_{s,(D)} &=& \mu^{-\ep/2} \hat g_s ~,\\
\hat g_s^2 &=& (4\pi)^2 \hat a_s ~.
\end{eqnarray}

From the independence of the bare coupling on $\mu$, one derives

\begin{eqnarray}
0 &=& \frac{d}{d\ln\mu^2} \hat a_{s,(D)} = \frac{d}{d\ln\mu^2} (\mu^{-\ep} \hat a_s) = \frac{d}{d\ln\mu^2} \big[\mu^{-\ep} Z_g(\ep,N_F,\mu^2) a_s(\mu^2) \big]
\end{eqnarray}

from which it follows that

\begin{equation}
\beta = \frac{\ep}{2} a_s(\mu^2) -2a_s(\mu^2) \frac{d}{d\ln\mu^2} \ln Z_g(\ep,N_F,\mu^2) ~.
\end{equation}

Specializing to the $\MS$ scheme for the coupling constant,

\begin{equation}
\beta^{\MS}(N_F) = -\beta_0(N_F) \big(a_s^{\MS}\big)^2 -\beta_1(N_F) \big(a_s^{\MS}\big)^3 +\mathcal O\Big(\big(a_s^{\MS}\big)^4\Big)
\end{equation}

and one obtains [36, 35, 237, 93, 94]

\begin{eqnarray}
\delta a^{\MS}_{s, 1}(N_F) &=& \frac{2}{\ep} \beta_0(N_F) \\
\delta a^{\MS}_{s, 2}(N_F) &=& \frac{4}{\ep^2} \beta_0^2(N_F) + \frac{1}{\ep} \beta_1(N_F)
\end{eqnarray}

\begin{eqnarray}
\beta_0(N_F)
&=& \frac{11}{3} C_A - \frac{4}{3} T_F N_F \label{beta0}~, \\
\beta_1(N_F)
&=& \frac{34}{3} C_A^2
- 4 \left(\frac{5}{3} C_A + C_F\right) T_F N_F \label{beta1}~.
\end{eqnarray}

In order to preserve the condition of having on-shell massless external particles and the decoupling of the massive quarks in the running of the coupling constant, one demands that the gluon self-energy receives no contribution from the heavy quark at zero momentum,

\begin{equation}
\Pi_{H,BF}(p^2=0,m^2) = 0 ~.
\end{equation}

This condition is enforced in the background field method [238, 239, 240, 158]. In this renormalization scheme the renormalization factor of the coupling constant is defined through

\begin{eqnarray}
Z_g^{\MOM} (\ep,N_F+1,\mu^2,m^2) &=& \frac{1}{(Z_{A,l}+Z_{A,H})^{1/2}} ~,\\
Z_{A,l} &=& \big(Z_g^{\MS} (\ep,N_F) \big)^{-2} ~,
\end{eqnarray}

with $Z_{A,(l,H)}$ the contributions to the renormalization factor of the background field due to the light quarks and the heavy quark.

In this $\MOM$ scheme, a formula analogous to (<ref>) holds, but reads [158]

\begin{eqnarray}
\hat a_s &=& \Big( Z_g^{\MOM}(\ep,N_F+1,\mu^2,m^2) \Big)^2 \ a_s^{\MOM}(\mu^2,m^2) \\
&=& a_s^{\MOM}(\mu^2,m^2) \Bigl[ 1 + \delta a^{\MOM}_{s, 1}
+ \delta a^{\MOM}_{s, 2}
\big( a_s^{\MOM}(\mu)\big)^2
+ \mathcal O\Big(\big({a_s^{\MOM}}\big)^3\Big) \Bigr] ~,
\nonumber\\
\end{eqnarray}

with [158]

\begin{eqnarray}
\delta a_{s,1}^{\MOM}&=&\Bigl[\frac{2\beta_0(N_F)}{\ep}
\Bigr]~,\label{dela1} \\
\delta a_{s,2}^{\MOM}&=&\Bigl[\frac{\beta_1(N_F)}{\ep}+
\Bigl\{\frac{2\beta_0(N_F)}{\ep}
\Bigl(\beta_{1,Q}+\ep\beta_{1,Q}^{(1)}
\Bigr)\Bigr] + O(\ep^2)~.\label{dela2}
\nonumber\\
\end{eqnarray}

\begin{eqnarray}
\beta_{0,Q} &=&-\frac{4}{3}T_F~, \label{b0Q} \\
\beta_{1,Q} &=&- 4 \left(\frac{5}{3} C_A + C_F \right) T_F~, \label{b1Q} \\
\beta_{1,Q}^{(1)}&=&
+15T_FC_F~, \label{b1Q1} \\
\beta_{1,Q}^{(2)}&=&
+T_FC_F\right) ~,\label{b1Q2}
\end{eqnarray}

\begin{eqnarray}
\Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
\exp \Bigl(\sum_{i=2}^{\infty}\frac{\zeta_i}{i}
\Bigl(\frac{\ep}{2}\Bigr)^{i}\Bigr)~. \label{fep}
\end{eqnarray}

From the invariance of the unrenormalized coupling,

\begin{equation}
\big( Z_g^{\MS}(\ep,N_F+1) \big)^2 \ a_s^{\MS}(\mu^2) = \big( Z_g^{\MOM}(\ep,N_F+1,\mu^2,m^2) \big)^2 \ a_s^{\MOM}(\mu^2)
\end{equation}

one can obtain the relation between $a_s^{\MS}$ and $a_s^{\MOM}$:

\begin{eqnarray}
\delta a^{\MOM}_{s, 1}
-\delta a^{\MS}_{s, 1}(N_F+1)
\biggr]
\delta a^{\MOM}_{s, 2}
-\delta a^{\MS}_{s, 2}(N_F+1)
\N\\ &&
-2\delta a^{\MS}_{s, 1}(N_F+1)\Bigl[
\delta a^{\MOM}_{s, 1}
-\delta a^{\MS}_{s, 1}(n_f+1)
\Bigr]
\biggr]+O({a_s^{\MOM}}^4)~, \label{asmsa}
\end{eqnarray}

Here, $a_s^{\MS}=a_s^{\MS}(N_F+1)$. Using the MOM scheme for the coupling, one can write the mass- and coupling-renormalized OME, indicated by ${\hat{A}}_{ij}$ [158]

\begin{eqnarray}
{\hat{A}}_{ij} &=&  \delta_{ij}
+ a^{\MOM}_s \Ahathat_{ij}^{(1)}
+ {a^{\MOM}_s}^2 \left[\Ahathat^{(2)}_{ij}
+ \delta m_1 \Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
m  \frac{d}{dm} \Ahathat_{ij}^{(1)}
+ \delta a^{\MOM}_{s,1} \Ahathat_{ij}^{(1)}\right]
\N\\ &&
+ {a^{\MOM}_s}^3 \Biggl[ \Ahathat^{(3)}_{ij}
+ \delta a^{\MOM}_{s,2} \Ahathat_{ij}^{(1)}
+ 2 \delta a^{\MOM}_{s,1} \left( \Ahathat^{(2)}_{ij}
+  \delta m_1 \Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
m \frac{d}{dm} \Ahathat_{ij}^{(1)} \right)
\N\\ &&
+ \delta m_1 \Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep/2}
m  \frac{d}{dm} \Ahathat_{ij}^{(2)}
+ \delta m_2 \Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep}
m  \frac{d}{dm} \Ahathat_{ij}^{(1)}
+ \frac{\delta m_1^2}{2} \Bigl(\frac{m^2}{\mu^2}\Bigr)^{\ep}
m^2  \frac{d^2}{{dm}~^2} \Ahathat_{ij}^{(1)}
\Biggr]
\N\\&&
\label{macoren}
\end{eqnarray}

At this stage, ultraviolet singularities are still present due to the composite operators. In the massless case, they are removed by imposing

\begin{equation}
A_{ij} \Bigl(\frac{-p^2}{\mu^2},N_F,N \Bigr) = Z_{ij} (\ep, a^{\MS},N_F,N) \ \hat A_{ij} \Bigl(\frac{-p^2}{\mu^2},\ep,N_F,N \Bigr) ~.
\label{eq:opZ}
\end{equation}

The pole structure in the operator $Z$-factors $Z_{ij}$ can be determined by the requirement that

\begin{eqnarray}
\gamma_{ij} = \sum_{k=0}^\infty \gamma_{ij}^{(k)} \big(a_s^{\MS}\big)^{k+1} = \mu Z_{il}^{-1}(\mu^2) \frac{\partial}{\partial \mu} Z_{lj}(\mu^2) ~.
\label{eq:opGammas}
\end{eqnarray}

Equations (<ref>) and (<ref>) have to be specialized to the singlet, non-singlet and pure-singlet cases appropriately, i.e.

\begin{eqnarray}
Z_{qq}^{-1}&=&Z_{qq}^{-1, {\sf PS}}+Z_{qq}^{-1, {\sf NS}}
\label{ZPSNS1}~,\\
A_{qq}     &=&A_{qq}^{\sf PS}+A_{qq}^{\sf NS} \label{ZPSNS2}~.
\end{eqnarray}

From Eq. (<ref>) one determines the pole terms in the $Z$-factors, which can be found written in explicit form in [158] and are not repeated here. One obtains the $Z$-factors in the case of $(N_F+1)$ flavours by taking them at $N_F+1$ flavours and applying the scheme transformation between $a_s^{\MS}$ and $a_s^{\MOM}$, i.e. the inverse of (<ref>). The OMEs are split into a light and a heavy flavour part,

\begin{equation}
\hat A_{ij}(p^2,m^2,\mu^2,a_s^{\MOM},N_F+1) = \hat A_{ij} \Big(\frac{-p^2}{\mu^2},a_s^{\MS},N_F\Big) + A_{ij}^Q (p^2,m^2,\mu^2,a_s^{\MOM},N_F+1)
\end{equation}

and the operator renormalization of the heavy contribution is [158]

\begin{eqnarray}
\bar A^Q_{ij}(p^2,m^2,\mu^2,a_s^{\MOM},N_F+1)&=&
\hat{A}^Q_{ij}(p^2,m^2,\mu^2,a_s^{\MOM},N_F+1)
\N\\ &&
\hat{A}_{ij}\Bigl(\frac{-p^2}{\mu^2},a_s^{\MS},N_F\Bigr)
\N\\ &&
\hat{A}_{ij}\Bigl(\frac{-p^2}{\mu^2},a_s^{\MS},N_F\Bigr)~.
\label{eqXX}
\end{eqnarray}

Taking now $p^2=0$, the contribution of the unrenormalized massless OMEs reduces to their tree level value because scaleless loop integrals vanish in dimensional regularization, so they reduce to

\begin{equation}
\hat A_{ij}(0,a_s^{\MS},N_F) = \delta_{ij}~.
\end{equation}

The remaining collinear singularities in $\bar A^Q_{ij}$ are removed by mass factorization,

\begin{equation}
A_{ij}^Q = \bar A_{il}^Q \Big(\frac{m^2}{\mu^2},a_s^{\MOM},N_F+1 \Big) \Gamma_{lj}^{-1}~.
\end{equation}

The transition functions $\Gamma_{ij}$ correspond to the inverse of the $Z$-factors in the massless case. We repeat here from [158] the formula for the renormalized OMEs which one obtains after these steps:

\begin{eqnarray}
&& A^Q_{ij}\Bigl(\frac{m^2}{\mu^2},a_s^{\MOM},n_f+1\Bigr)=
\N\\&&\phantom{+}
\hat{A}_{ij}^{(1),Q}\Bigl(\frac{m^2}{\mu^2}\Bigr)
\Biggr)
\N\\&&
\hat{A}_{ij}^{(2),Q}\Bigl(\frac{m^2}{\mu^2}\Bigr)
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\N\\
+\Bigl[ \hat{A}_{il}^{(1),Q}
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\Bigr]
\Gamma^{-1,(1)}_{lj}(n_f)
\Biggr)
\N\\
\hat{A}_{ij}^{(3),Q}\Bigl(\frac{m^2}{\mu^2}\Bigr)
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\nonumber\end{eqnarray}

\begin{eqnarray}
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\hat{A}_{il}^{(1),Q}
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\N\\ &&\phantom{+{a^{\MOM}_s}^3\Biggl(}
\Bigr]
\Gamma^{-1,(2)}_{lj}(n_f)
\hat{A}_{il}^{(2),Q}
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\N\\
\Bigl(\frac{m^2}{\mu^2}\Bigr)
\Bigr]
\Gamma^{-1,(1)}_{lj}(n_f)
\Biggr)~. \label{GenRen3}
\end{eqnarray}

By applying (<ref>) and (<ref>) as well as the pole expansions of the $Z$ and $\Gamma$ factors, the coefficients of the $\ep$-expansion of the OMEs have been predicted in terms of the renormalization constants up to $\mathcal O(a_s^3)$. These predictions are derived by demanding that (<ref>) be free of poles in $\ep$.

In analogy with the process described above, the renormalization of the OMEs in the presence of two heavy quarks has been preformed in [121, 178]. The main differences to the single-mass case summarized above are due to the fact that the renormalization of one mass will depend on the other; also, two quarks have to be decoupled in the running of the coupling. The procedure is not repeated here and only the relevant results are printed in the next sections for the OMEs respectively under consideration, where we focus on the genuine two-mass contribution to the OMEs, i.e. on graphs which contain two different masses.

§.§ Variable flavour number scheme

In the so-called fixed-flavour number scheme, heavy quarks are treated as if they were radiatively generated, and are not assigned a parton distribution. Such calculations give rise to logarithms of the type $\ln(Q^2/m^2)$, with $m$ the mass of the heavy quark. In principle, for very large virtualities, these logarithms could become large and ruin the convergence of the perturbative series. In the presently accessible kinematic range at HERA, however, the fixed-flavour number scheme has been shown to be numerically stable [241, 242].

The variable flavour number scheme (VFNS), by contrast, exploits the possibility to factorize the structure functions for $Q^2\gg m^2$ using the massive operator matrix elements and the massless Wilson coefficients, and devises a transition between, e.g., $N_F$ and $N_F+1$ flavours, however, only at very high scales $\mu^2=Q^2$. At the transition scale $\mu^2$, a new PDF is introduced for the heavy quark, which from then on is treated as massless. The transition is designed such that the structure functions are unchanged asymptotically after the new PDF is introduced and the light quark PDFs are adjusted appropriately. The VFNS was discussed in [110, 111, 112, 114, 118]; the matching conditions were given at $\mathcal O(a_s^2)$ in [112], at $\mathcal O(a_s^3)$ in [158] and the VFNS where two quarks are decoupled in [121]. This VFNS prescription is also known as the “zero-mass VFNS".

For concreteness we reproduce here the single-mass VFNS relations given in [158]:

\begin{eqnarray}
\label{HPDF1}
f_k(N_F+1,\mu^2,m^2,N) + f_{\overline{k}}(N_F+1,\mu^2,m^2,N)
&=& A_{qq,Q}^{\rm NS}\left(N_F,\frac{\mu^2}{m^2},N\right)
\cdot \bigl[f_k(N_F,\mu^2,N)
\nonumber\\ && \hspace*{3.8cm}
+ f_{\overline{k}}(N_F,\mu^2,N)\bigr]
\nonumber\\
& & \hspace*{-4mm} + \tilde{A}_{qq,Q}^{\rm
\cdot \Sigma(N_F,\mu^2,N)
\nonumber\\
& & \hspace*{-4mm} + \tilde{A}_{qg,Q}\left(N_F,\frac{\mu^2}{m^2},N\right)
\cdot G(N_F,\mu^2,N),
\nonumber\\
\\
\label{fQQB}
f_Q(N_F+1,\mu^2,m^2,N) + f_{\overline{Q}}(N_F+1,\mu^2,m^2,N)
{A}_{Qq}^{\rm PS}\left(N_F,\frac{\mu^2}{m^2},N\right)
\cdot \Sigma(N_F,\mu^2,N)
\nonumber\\ && \hspace*{-4mm}
+ {A}_{Qg}\left(N_F,\frac{\mu^2}{m^2},N\right)
\cdot G(N_F,\mu^2,N)~.
\nonumber\\
\end{eqnarray}

\begin{eqnarray}
\Sigma(N_F+1,\mu^2,m^2,N)
&=& \sum_{k+1}^{N_F+1} \big[f_k(N_F+1,\mu^2) + f_{\bar k}(N_F+1,\mu^2) \big]
\nonumber\\
&=& \Biggl[A_{qq,Q}^{\rm NS}\left(N_F, \frac{\mu^2}{m^2},N\right) +
N_F \tilde{A}_{qq,Q}^{\rm PS}\left(N_F, \frac{\mu^2}{m^2},N\right)
\nonumber \\ &&
+ {A}_{Qq}^{\rm PS}\left(N_F, \frac{\mu^2}{m^2},N\right)
\Biggr]
\cdot \Sigma(N_F,\mu^2,N) \nonumber
\\
%\end{eqnarray} \begin{eqnarray}
& & \hspace*{-3mm} + \left[N_F \tilde{A}_{qg,Q}\left(N_F,
\frac{\mu^2}{m^2},N\right) +
{A}_{Qg}\left(N_F, \frac{\mu^2}{m^2},N\right)
\right]
\cdot G(N_F,\mu^2,N)
\nonumber\\
\\
\Delta_k(N_F+1,\mu^2,m^2,N) &=&
+ f_{\overline{k}}(N_F+1,\mu^2,m^2,N)
\nonumber\\ &&
- \frac{1}{N_F+1}
\Sigma(N_F+1,\mu^2,m^2,N)\\
\label{HPDF2}
G(N_F+1,\mu^2,m^2,N) &=& A_{gq,Q}\left(N_F,
\frac{\mu^2}{m^2},N\right)
\cdot \Sigma(N_F,\mu^2,N)
\nonumber\\
+ A_{gg,Q}\left(N_F, \frac{\mu^2}{m^2},N\right)
\cdot G(N_F,\mu^2,N)~.
\end{eqnarray}

where $f_{k,\bar k}$ are the light quark and antiquark PDFs, $f_{Q,\bar Q}$ refer to the new PDF of the heavy quark, $G$ is the gluon distribution, $\Sigma$ is the singlet distribution defined in Eq. (<ref>)
and $\Delta_k$ the non-singlet distribution, Eq. (<ref>).

The two-mass VFNS derived in [121] is also reproduced below, for completeness:

\begin{eqnarray}
\label{eq:VFNS1}
\hspace*{-2.4cm}
f_k(N_F+2,N,\mu^2,m_1^2,m_2^2) + f_{\overline{k}}(N_F+2,N,\mu^2,m_1^2,m_2^2)=
\NN\\ &&  \hspace{20mm}
A_{qq,Q}^{\sf NS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
\cdot\bigl[f_k(N_F,N,\mu^2)+f_{\overline{k}}(N_F,N,\mu^2)\bigr]
\NN\\ &&  \hspace{20mm}
+\frac{1}{N_F}A_{qq,Q}^{\sf PS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
\cdot\Sigma(N_F,N,\mu^2)
\NN\\ &&  \hspace{20mm}
\cdot G(N_F,N,\mu^2), \label{HPDF12M} \\
\hspace*{-2.4cm}\label{fQQB2M}
f_Q(N_F+2,N,\mu^2,m^2_1,m_2^2) + f_{\overline{Q}}(N_F+2,N,\mu^2,m^2_1,m_2^2)=
\NN\\ &&  \hspace{20mm}
A_{Qq}^{\sf PS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},
\frac{m_2^2}{\mu^2},
\right)
\cdot \Sigma(N_F,N,\mu^2)
\NN\\ &&  \hspace{20mm}
\cdot G(N_F,N,\mu^2)~.
\end{eqnarray}

In this case, the flavor singlet, non-singlet and gluon densities for $(N_F+2)$ flavors are given by

\begin{eqnarray}
\Sigma(N_F+2,N,\mu^2,m^2_1,m_2^2)
&=& \Biggl[
A_{qq,Q}^{\sf NS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
+A_{qq,Q}^{\sf PS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
\NN \\ && \hspace*{-20mm}
+A_{Qq}^{\sf PS}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
\Biggr]  \cdot \Sigma(N_F,N,\mu^2)
\NN \\ && \hspace*{-23mm}
\right]   \cdot G(N_F,N,\mu^2)~,
\NN\\ \\
\Delta_k(N_F+2,N,\mu^2,m_1^2,m_2^2)
&=& f_k(N_F+2,N,\mu^2,m_1^2,m_2^2)+f_{\overline{k}}(N_F+2,N,\mu^2,m_1^2,m_2^2)
\NN\\ &&
-\frac{1}{N_F+2}\Sigma(N_F+2,N,\mu^2,m_1^2,m_2^2)~, \\
\label{HPDF22M}
G(N_F+2,N,\mu^2,m^2_1, m^2_2)
&=& A_{gq,Q}\left(N,N_F+2,\frac{m_1^2}{\mu^2},\frac{m_2^2}{\mu^2}\right)
\cdot \Sigma(N_F,N,\mu^2)
\NN\\ &&
\cdot G(N_F,N,\mu^2)~.
\end{eqnarray}

One can observe how the two-mass contributions to the OME enter the definitions of the parton densities. Due to the presence of diagrams with two different quark masses, decoupling the charm and bottom quarks one at a time becomes theoretically an ill-defined procedure at higher order. This is one main motivation for the computation of the two-mass contributions to OMEs.

§.§ Mathematical methods

§.§.§ The Mellin transform

The Mellin transform is a central mathematical operation in the study of deep-inelastic scattering, because the Bjorken variable $x$ and the variable $N$ are conjugate to each other with respect to it. In Eq. (<ref>) we saw an example of this fact.

The Mellin transform [243, 244, 245] of a function $f(x)$ is defined as

\begin{equation}
\mathbf{M}\left[f\right](N)=F(N)=\int_{0}^1 x^{N-1}f(x)dx ,
\end{equation}

whenever the integral exists.
In our physical applications, poles of $\mathbf{M}\left[f\right](N)$ arise along the negative real axis, and the rightmost pole will be located at $N=1$.

In Table <ref>, a few common cases of Mellin transforms are summarized.

$f(x)$     $\mathbf M [f(x)](N)=F(N)$$x^a$     $\frac{1}{N+a}$ $(1+x)^{-a}$     $\frac{\Gamma(s)\Gamma(a-s)}{\Gamma(a)}$$\delta(1-x)$     $1$$\ln^a x$     $(-1)^a N^{-a-1} a!$ $\ln^a f(x)$     $F^{(a)}(N)$ $f^{(a)}(x)$     $\frac{\Gamma(a+1-N)}{\Gamma(1-N)}F(N-a)$

Some elementary properties of the Mellin transform.

From the Mellin transform of a function it is possible to recover the original function using the inverse Mellin transform

\begin{equation}
f(x)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}x^{-N}F(N)dN .
\end{equation}

The Mellin convolution of two functions is defined as

\begin{equation}
\left[f\otimes g\right](x)=\int_0^\infty dx_1 \int_0^\infty dx_2 \, \delta(x-x_1 x_2)f(x_1)g(x_2).
\end{equation}

A property of the Mellin convolution is:

\begin{equation}
\mathbf{M}\left[f\otimes g\right](s)=\left(\mathbf{M}\left[f\right](s)\right)\left(\mathbf{M}\left[g\right](s)\right),
\end{equation}

or, in other words, the Mellin convolution becomes an ordinary product in Mellin space.

The functions of relevance to physics which we are concerned with, such as parton distribution functions, are defined in the interval $[0,1]$, and vanish outside of this interval.

§.§.§ Nested sums

From the earliest computations of anomalous dimensions in QCD, sums have appeared in physical quantities [38, 39], even though a more systematic study of these objects in the context of high-energy physics appeared only much later [246, 247].

The first class of sums to be studied systematically [246, 247], see also [208] is that of the harmonic sums, defined recursively as

\begin{eqnarray}
S_{n_1,\ldots,n_k}(N) &=& \sum_{i=1}^{N} \frac{(\text{sign}(n_1))^i}{i^{|n_1|}} S_{n_2,\ldots,n_k}(i), \qquad n_i \in \mathbb Z \backslash\{0\}, \label{eq:defHSums1}\\
S_\emptyset &=& 1. \label{eq:defHSums2}
\end{eqnarray}

The quantity $w = \sum|n_i|$ is called the weight of the sum, and $k$ is called the depth. For any given weight there are $2\cdot 3w-1$ different harmonic sums.

Generalizations of these objects which have been encountered in QCD computations have the form [221, 208, 248]

\begin{equation}
S_{n_1,\ldots,n_k}(x_1,\ldots,x_k;N) = \sum_{i=1}^N \frac{x_1^i}{i^{|n_1|}} S_{n_2,\ldots,n_k}(x_2,\ldots,x_k;i), \qquad x_i\neq 0 .
\label{eq:genSsums}
\end{equation}

A further generalization is the class of cyclotomic sums [249], defined as

\begin{equation}
S_{\{a_1,b_1,c_1\},\ldots,\{a_k,b_k,c_k\}}(x_1,\ldots,x_k;N) = \sum_{i=1}^N \frac{x_1^i}{(a_1 i+b_1)^{c_1}} S_{\{a_2,b_2,c_2\},\ldots,\{a_k,b_k,c_k\}}(x_2,\ldots,x_k;i) .
\end{equation}

In three-loop massive calculations, in addition, more complex sums, both finite and infinite, involving binomial summands have been found to contribute, [250, 253, 251, 252, 254]. Examples of such summands are

\begin{eqnarray}
&&\binom{2i}{i}^r \frac{x^i}{i}, \qquad r=\pm 1,  \\
&& \frac{4^i}{i \binom{2i}{i}} \Bigl(\frac{\eta}{\eta-1}\Bigr)^i, \qquad 0<\eta<1,
\end{eqnarray}

where in our applications $\eta$ is the ratio of the squares of two quark masses.
No specific symbol has been used for this type of sums in the literature and they are typically written explicitly.

Nested sums obey algebraic and structural properties [257, 256, 255], most importantly the algebraic quasi-shuffle algebra [259]. This algebra can be derived from relations of the type

\begin{equation}
\Bigl(\sum_{i=1}^N a_i\Bigr) \Bigl(\sum_{i=1}^N b_i\Bigr) = \sum_{i=1}^N a_i \sum_{j=1}^i b_j + \sum_{i=1}^N b_i \sum_{j=1}^i a_j - \sum_{i=1}^N a_i b_i .
\end{equation}

Such relations allow to reduce the sums to a smaller set of “basis” sums which are algebraically independent.

The classes of sums described above have been an object of interest for mathematicians, and Mathematica packages exist to perform algebraic reductions and inverse Mellin transforms. We have extensively used HarmonicSums [260] in this thesis.

§.§.§ Iterated integrals

Iterated integrals are a powerful way to represent many classes of functions which arise naturally in the computation of Feynman integrals. They naturally appear, for example, when the method of differential equations is used to compute the Feynman integrals [264, 265], as well as in other areas [221, 208]. In this thesis, they occur in the $x$-space representation of OMEs, as the result of inverse Mellin transforms.

Iterated integrals are functions defined as

\begin{equation}
=\int_0^z  d\tau_1~f_1(\tau_1)
\label{eq:Gfunctions}
\end{equation}

\begin{equation}
\cdots,\frac{1}{\tau}}_{\text{n times}}\Biggr\},z\Biggr)
\equiv
\frac{1}{n!} \ln^n(z)~.
\label{eq:Gfunctions2}
\end{equation}

The set $\{f_i\}$ of functions appearing in a given physical problem is called the alphabet.

One example of iterated integrals studied in the context of high-energy physics is that of the harmonic polylogarithms or HPLs [266, 267, 268, 269], which can be considered a special case of iterated integrals. They are defined as

\begin{eqnarray}
\HA_{b,\vec{a}}(x) &=& \int_0^x dy f_b(y) \HA_{\vec{a}}(y),~~~\HA_\emptyset = 1,~a_i, b~\in~\{0,1,-1\}~,
\label{eq:HPLdef1}
\end{eqnarray}

\begin{eqnarray}
f_0(x) = \frac{1}{x},~~~f_1(x) = \frac{1}{1-x},~~~f_{-1}(x) = \frac{1}{1+x}~,
\label{eq:HPLdef2}
\end{eqnarray}

\begin{equation}
\HA_{\underbrace{\text{\scriptsize 0},\ldots,\text{\scriptsize 0}}_{\text{\scriptsize n times}}}(x)
\frac{1}{n!} \ln^n(x)~.
\label{eq:HPLdef3}
\end{equation}

The number of indices in $H_{\vec a}(x)$ is called the weight of the HPL. A numerical library for the evaluation of HPLs up to weight 5 has been published in [270] and to weight 8 in [271].

Harmonic polylogarithms are closely related to the Mellin transform of harmonic sums. For example, defining the “+”-distribution as

\begin{equation}
\int_0^1 dx \, f(x) \bigl(g(x)\bigr)_+ = \int_0^1 dx \, \bigl(f(x)-f(1)\bigr) g(x)
\end{equation}

one has

\begin{equation}
S_1(N) = \mathbf M\Big[\Bigl(\frac{1}{x-1}\Bigr)_+\Big] .
\end{equation}

By repeated integration-by-parts, one can obtain the Mellin transform of HPLs in terms of harmonic sums (possibly evaluated at infinity) provided that appropriate regularizations such as the +-distribution are used in the Mellin transform.

Another special class is that of cyclotomic multiple polylogarithms [249], whose alphabet is

\begin{equation}
\Bigl\{\frac{1}{x}\Bigr\}\cup\Big\{ \frac{x^b}{\Phi_n(x)} \mid n \in \mathbb N_+,0\leq b\leq\varphi(k) \Big\}
\end{equation}

\begin{equation}
\Phi_n(x) = \prod_{\substack{1\le k\le n\\\text{gcd}(k,n)=1}} \Bigl(x-e^{2i\pi \frac{k}{n}} \Bigr)
\end{equation}

where $\varphi(n)$ is Euler's totient function. The polynomials $\Phi_n(x)$ are called cyclotomic polynomials. Cyclotomic harmonic polylogarithms are related through the Mellin transform to the cyclotomic harmonic sums.

Iterated integrals also obey shuffle algebras [255] due to the identity

\begin{equation}
\Bigl(\int_0^x dy \, f(y) \Bigr) \Bigl(\int_0^x dy \, g(y) \Bigr) = \int_0^x dy \, f(y) \int_0^y dz \, g(z) + \int_0^x dy \, g(y) \int_0^y dz \, f(z) ,
\end{equation}

which allow to reduce the iterated integrals to a basis. Many algorithms pertaining to the symbolic manipulation of iterated integrals are available in HarmonicSums.

§.§.§ Mellin-Barnes integration

The Feynman parametrization is one of the standard ways to compute Feynman integrals. It consists in the repeated application to Feynman integrals of the identity

\begin{equation}
\frac{1}{A_1^{\nu_1}\cdots A_n^{\nu_n}} = \frac{\Gamma(\sum_{i=1}^n \nu_i)}{\Gamma(\nu_1)\cdots\Gamma(\nu_n)} \int_0^1 dx_1\cdots\int_0^1 dx_n \, \frac{x_1^{\nu_1-1}\cdots x_n^{\nu_n-1}}{(x_1 A_1+\cdots+x_n A_n)^{x_1+\cdots+x_n}} \,\delta\Big(1-\sum_{i=1}^n \nu_i\Big),
\end{equation}

which is valid for $A_i>0$, $\text{Re}(\nu_i)>0$. This identity is one of the methods which can be used to turn the integrals over the loop momenta into integrals over scalar Feynman parameters $x_i$. The evaluation of the integrals over $x_i$ is then the source of complexity in the Feynman representation, and the integrals will in general evaluate to special functions, in many cases yet unknown. One way to approach the integrals in the Feynman parametrization, which has been used in the calculations in this thesis, in part is through Mellin-Barnes integrals.

The Mellin-Barnes formula [272, 273, 274, 275]

\begin{equation}
\frac{1}{(A+B)^\lambda} = \frac{1}{2\pi i \Gamma(\lambda)} \int_{-i\infty}^{i\infty} d\sigma \,\Gamma(\sigma+\lambda)\Gamma(-\sigma) \frac{A^\sigma}{B^{\lambda+\sigma}}
\label{eq:mbFormula}
\end{equation}

is one important tool for the evaluation of Feynman integrals; early applications can be found in [276, 277], see also [278].
In (<ref>), the integration contour must separate the infinite sets of ascending and descending poles of the $\Gamma$-functions, and must otherwise stretch in the direction of the imaginary axis.
This formula is used to disentangle polynomials appearing from the Feynman parametrization. In general, the contour can be quite involved, and can necessitate the separate calculation of one or more residues, particularly if the integral is divergent in the dimensional parameter $\ep$. In any case, symbolic algebra packages exist to perform the analytic continuation and $\ep$-expansion of Mellin-Barnes integrals as well as for numerical evaluation. In this thesis we made use of MB and MBResolve [279, 280].

The evaluation of Mellin-Barnes integrals with these packages is not possible in general if the integrand depends on the further symbolic parameter $N$, as is the case in the calculation of OMEs. For this application, the packages were therefore used only to compute moments. In the general case, one can apply the residue theorem to evaluate the integral into a (nested) infinite sum, and turn to algorithms in summation theory from the package Sigma [202, 203, 204].

§ POLARIZED DEEP-INELASTIC SCATTERING

In the following, we present the calculation of the two-mass contributions to $A_{Qq}^{(3),\mathrm{PS}}$ and to $A_{gg,Q}^{(3)}$, which were performed in [196] and in [198] respectively.

§.§ The two-mass contribution to the polarized operator matrix element AQq(3),PS

The calculation presented here for the polarized $A_{Qq}^{(3),\mathrm{PS}}$ closely mirrors that of the corresponding unpolarized OME, which was performed in [171].

The Feynman rule used for the operator insertion, for the quark-quark-gluon vertex, is taken from [81]: one has

\begin{equation}
O_a^\mu(p,q)=-g T_a \Delta^\mu/\!\!\!p\gamma_5\sum_{i=0}^{N-2}(\Delta.p)^{N-i-2}(-\Delta.q)^i,
\end{equation}

where $\Delta^\mu$ is a light-like vector, $\Delta^2=0$ and $p$, $q$ are incoming quark momenta.

When using dimensional regularization in polarized physics, a choice must be made for the analytic continuation from four to $D$ dimensions of chiral quantities, such as $\gamma_5$, which are intrinsically four-dimensional. This calculation has been performed in the Larin scheme [197], which is the definition

\begin{equation}
\gamma_5 = \frac{i}{24}\gamma_\mu\gamma_\nu\gamma_\rho\gamma_\delta\varepsilon^{\mu\nu\rho\delta}.
\label{eq:gamma5Larin}
\end{equation}

The contraction of two Levi-Civita tensors is then

\begin{equation}
\varepsilon_{\alpha\beta\gamma\delta}\varepsilon^{\mu\nu\rho\sigma}=\left|\begin{array}{cccc}
g_{\alpha\mu} & g_{\alpha\nu} & g_{\alpha\rho} & g_{\alpha\sigma}\\
g_{\beta\mu} & g_{\beta\nu} & g_{\beta\rho} & g_{\beta\sigma}\\
g_{\gamma\mu} & g_{\gamma\nu} & g_{\gamma\rho} & g_{\gamma\sigma}\\
g_{\delta\mu} & g_{\delta\nu} & g_{\delta\rho} & g_{\delta\sigma}
\end{array}\right| ~.
\label{eq:LeviCivita}
\end{equation}

Other scheme choices have been made in the literature; in particular, the anomalous dimensions have been computed in the M-scheme, which was first defined in [86] and is obtained from the Larin scheme through a finite renormalization.

In order to compare with the literature, care must be taken to adopt a consistent scheme choice: see [85] for the relationship between the anomalous dimensions in the two schemes.

The pole structure for the OME can be derived from the renormalization structure; see [121] for this particular case. The predicted pole structure, which was used to check the calculation, is presented in Eq. (<ref>):

\begin{eqnarray}
\hat{\hat{\tilde{A}}}_{Qq}^{(3),\rm{PS, tm}} &=&
\frac{8}{3 \ep^3} \gamma_{gq}^{(0)} \hat{\gamma}_{qg}^{(0)} \beta_{0,Q}
2 \gamma_{gq}^{(0)} \hat{\gamma}_{qg}^{(0)} \beta_{0,Q} \left(L_1+L_2\right)
+\frac{1}{6} \hat{\gamma}_{qg}^{(0)} \hat{\gamma}_{gq}^{(1)}
-\frac{4}{3} \beta_{0,Q} \hat{\gamma}_{qq}^{\rm{PS},(1)}
\biggr]
\nonumber\\&&
\biggl[
\gamma_{gq}^{(0)} \hat{\gamma}_{qg}^{(0)} \beta_{0,Q} \left(L_1^2+L_1 L_2+L_2^2\right)
+\biggl\{\frac{1}{8} \hat{\gamma}_{qg}^{(0)} \hat{\gamma}_{gq}^{(1)}
- \beta_{0,Q} \hat{\gamma}_{qq}^{\rm{PS},(1)}
\biggr\} \left(L_2+L_1\right)
\nonumber\\&&
+\frac{1}{3} \hat{\tilde{\gamma}}_{qq}^{(2),\rm{PS}}
-8 a_{Qq}^{(2),\rm{PS}} \beta_{0,Q}
+ \hat{\gamma}_{qg}^{(0)} a_{gq}^{(2)}
\biggr]
\label{Ahhhqq3PSQ},
\end{eqnarray}

\begin{eqnarray}
\hat{\gamma}_{ij} &=& \gamma_{ij}(N_F+2) - \gamma_{ij}(N_F),
\\
\hat{\tilde{\gamma}}_{ij} &=& \frac{\gamma_{ij}(N_F+2)}{N_F+2}  - \frac{\gamma_{ij}(N_F)}{N_F} ,
\end{eqnarray}

and the notation $a_{ij}$, $\bar{a}_{ij}$ denote the respective $\mathcal{O}(\ep^0)$, $\mathcal{O}(\ep)$ terms of the OMEs, while $\zeta_k$ is the Riemann $\zeta$-function at integer values, Eq. (<ref>).
The quantity $\tilde{a}_{Qq}^{(3),\rm{PS}}\left(m_1^2,m_2^2,\mu^2\right)$ is the object of this calculation.

The diagrams for the two-mass contributions to $\tilde{A}_{Qq}^{(3), \rm PS}$. The dashed arrow line represents the external
massless quarks, while the thick solid arrow line represents a quark of mass $m_1$, and the thin arrow line a quark of mass $m_2$. We assume $m_1 > m_2$.

The contributing diagrams are shown in Figure <ref>. The unrenormalized OME is obtained from their sum by applying the projector [85, 178]

\begin{eqnarray}
\label{eq:PqNEW}
P_q \hat{G}_l^{ij} = - \delta_{ij} \frac{i (\Delta.p)^{-N-1}}{4 N_c (D-2)(D-3)} \ep_{\mu \nu p \Delta} {\rm tr} \left[p
\hspace*{-2mm}
\slash
\gamma^\mu \gamma^\nu \hat{G}_l^{ij}\right].
\end{eqnarray}

The numerator algebra was performed in Form [281] and in Mathematica.
Of the contributing diagrams, the nonzero ones are 9–12 and 13–16, which are related by a symmetry under the exchange of the two heavy quarks. One obtains for the unrenormalized OME

\begin{eqnarray}
A_{Qq}^{(3), {\rm PS, tm}}(N) &=&
2 \left[1 + (-1)^{N-1} \right] D_9(m_1,m_2,N) + 2 \left[1 + (-1)^{N-1} \right] D_9(m_2,m_1,N),
\label{eq:D9}
\end{eqnarray}

where we define the variable

\begin{equation}
\eta = \frac{m_2^2}{m_1^2}<1.
\end{equation}

The calculation made use of the following Mellin-Barnes representation for the massive bubbles of Figure <ref>:

Massive bubbles appearing in the Feynman diagrams shown in Figure <ref>.

\begin{eqnarray}
I^{\mu\nu,ab}_{a_1}(k) &=&
-\frac{8iT_F g_s^2}{(4\pi)^{D/2}}\delta_{ab}(k^2g^{\mu\nu}-k^\mu k^\nu)\int_0^1
I^{\mu\nu,ab}_{b_2}(k) &=&
\alpha_s T_F i e^{-\gamma_E \varepsilon/2}(k\cdot\Delta)^{N-1} (\mu^2)^{-\varepsilon/2}S_\varepsilon \epsilon^{\Delta k
\mu
\nu}\int_0^1 dx\, x^{N+D/2-1}(1-x)^{D/2-1}\nonumber\\
\times\Bigg\{\left(-k^2+\frac{m^2}{x(1-x)}\right)^{-2+D/2}2\Gamma(2-D/2)\left[(D-6)x^{-2}+(D+2N)x^{-1}\right]\\\nonumber&&
\end{eqnarray}

After the Feynman parametrization and the Mellin-Barnes decomposition are performed, one obtains for Diagram 9 the representation

\begin{eqnarray}
D_{9}(m_{1},m_{2},N) & = & C_F T_F^2 \alpha_s^3 S_\varepsilon^3 \frac{16}{2+\varepsilon}\Bigg\{ 4(2-\varepsilon)J_{1}-8\eta
\nonumber\\
& &
\times J_{5}
-8J_{6}-(\varepsilon-2)^{2}J_{7}+2(2-\varepsilon)\eta J_{8}+2(2-\varepsilon)(3+N)J_{9}-2(2-\varepsilon)J_{10}
\nonumber\\
&  & -2(2-\varepsilon)\bigg(2+\frac{\varepsilon}{2}+N\bigg)J_{11}+2(2-\varepsilon)J_{12}-8\eta J_{13}+2(2-\varepsilon)\eta
\end{eqnarray}

\begin{eqnarray}
J_{1} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-1+\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{2} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-1+\frac{\varepsilon }{2}+N}B_{3}\left(\frac{\eta}{x(1-x)}\right),\\
J_{3} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{4} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{\frac{\varepsilon }{2}+N}B_{2}\left(\frac{\eta}{x(1-x)}\right),\\
J_{5} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{1+\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{6} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{1+\frac{\varepsilon }{2}+N}B_{2}\left(\frac{\eta}{x(1-x)}\right),\\
J_{7} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-1+\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{8} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-1+\frac{\varepsilon }{2}+N}B_{3}\left(\frac{\eta}{x(1-x)}\right),\\
J_{9} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{10} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{\frac{\varepsilon }{2}+N}B_{2}\left(\frac{\eta}{x(1-x)}\right),\\
J_{11} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{1+\frac{\varepsilon }{2}+N}B_{1}\left(\frac{\eta}{x(1-x)}\right),\\
J_{12} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{1+\frac{\varepsilon }{2}+N}B_{2}\left(\frac{\eta}{x(1-x)}\right),\\
J_{13} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N)}{\Gamma\big(1+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-2+\frac{\varepsilon }{2}+N}B_{3}\left(\frac{\eta}{x(1-x)}\right),\\
J_{14} & = & \left(\frac{m_{1}^{2}}{\mu^{2}}\right)^{\frac{3}{2}\varepsilon}\frac{\Gamma(N+1)}{\Gamma\big(2+\frac{\varepsilon
}{2}+N\big)}\int_{0}^{1}dx\,(1-x)^{\frac{\varepsilon}{2}}x^{-2+\frac{\varepsilon }{2}+N}B_{3}\left(\frac{\eta}{x(1-x)}\right).
\end{eqnarray}

The functions $B_i$ are given by

\begin{eqnarray}
B_1(\xi) &=& \frac{1}{2 \pi i} \int_{-i \infty}^{i \infty} d\sigma \, \xi^{\sigma} \,
\Gamma(-\sigma) \Gamma(-\sigma+\varepsilon) \Gamma\left(\sigma-\frac{3 \varepsilon}{2}\right) \Gamma\left(\sigma
\frac{\Gamma^2(\sigma+2-\varepsilon)}{\Gamma(2 \sigma+4-2 \varepsilon)},
\label{B1}
\\
B_2(\xi) &=& \frac{1}{2 \pi i} \int_{-i \infty}^{i \infty} d\sigma \, \xi^{\sigma} \,
\Gamma(-\sigma) \Gamma(-\sigma+\varepsilon) \Gamma\left(\sigma-\frac{3 \varepsilon}{2}\right)
\Gamma\left(\sigma+1-\frac{\varepsilon}{2}\right)
\frac{\Gamma^2(\sigma+2-\varepsilon)}{\Gamma(2 \sigma+4-2 \varepsilon)},
\nonumber \\
\label{B2}
\\
B_3(\xi) &=& \frac{1}{2 \pi i} \int_{-i \infty}^{i \infty} d\sigma \, \xi^{\sigma} \,
\Gamma(-\sigma) \Gamma(-\sigma-1+\varepsilon) \Gamma\left(\sigma+1-\frac{3 \varepsilon}{2}\right)
\Gamma\left(\sigma+1-\frac{\varepsilon}{2}\right)
\nonumber \\ && \phantom{\frac{1}{2 \pi i} \int_{-i \infty}^{i \infty} d\sigma} \times
\frac{\Gamma^2(\sigma+3-\varepsilon)}{\Gamma(2 \sigma+6-2 \varepsilon)}.
\label{B3}
\end{eqnarray}

where we take

\begin{equation}
\xi = \frac{1}{\eta x(1-x)}
\end{equation}

for diagram 9.
It was possible to reduce this representation to the same class of functions as in the unpolarized calculation, as is expected from the fact that the two cases differ only in numerator structures. The calculation was then performed in $x$ space by taking residues of the functions $B_i$ in $\sigma$ and expanding in $\ep$. Because the functions $B_i$ are the same as those appearing in the unpolarized calculation, it was possible to refer to their $\mathcal{O}(\ep^0)$ behaviour as computed in [171], where this expansion was performed with the packages MB [279], MBresolve [280], Sigma [202, 203, 204], HarmonicSums [260, 249, 221], EvaluateMultiSums and SumProduction [282]. In this way, the integrals appearing in $B_i$ are evaluated in terms of sums involving harmonic sums.

The convergence of the integrals depends on the value of $\xi$. This implies that the interval $x\in [0,1]$ needs to be split into three intervals

\begin{eqnarray}
[0,\eta_-],~~~[\eta_-,\eta_+],~~~[\eta_+,1],~~~\text{with}~~\eta_\pm = \frac{1}{2}\left(1 \pm \sqrt{1-\eta}\right),
\end{eqnarray}

where the contours appearing in the functions $B_i$ are closed to the right (for the second region) or to the left (for the first and the third), giving rise to two different functional forms for the residue sums.
The constant parts in $\ep$, $B_i^{(0)}$, depend on harmonic sums [246, 247], defined in Eqs. (<ref>) and (<ref>).

Then, the prefactors appearing in $J_i$ can be expanded in $\ep$, giving rise, after partial fractioning, to denominators of the type

\begin{equation}
\frac{1}{N+l}, \quad {\rm with} \quad l \in \{0,1\},
\end{equation}

which can be absorbed inside an integral using the relation

\begin{eqnarray}
\frac{1}{N+l} \int_a^b dx \, x^{N-1} f(x) &=&
\frac{b^{N+l}}{N+l} \int_a^b dy \frac{f(y)}{y^{l+1}} - \int_a^b dx \, x^{N+l-1} \int_a^x dy \frac{f(y)}{y^{l+1}}
\label{absorbN1}
\\ &=&
\frac{a^{N+l}}{N+l} \int_a^b dy \frac{f(y)}{y^{l+1}} + \int_a^b dx \, x^{N+l-1} \int_x^b dy \frac{f(y)}{y^{l+1}}.
\label{absorbN2}
\end{eqnarray}

The purpose of this method is to ultimately obtain $\tilde{a}_{Qq}^{(3),\rm{PS}}$ in $x$-space, by leaving one of the Feynman parameters unintegrated.

The final result is expressed in terms of generalized iterated integrals, as defined in Eqs. (<ref>), (<ref>)
and harmonic polylogarithms [266] which can be considered a special case of iterated integrals and are defined as in (<ref>), (<ref>), (<ref>).

In this case, the letters appearing in the iterated integrals are

\begin{equation}
\frac{1}{\tau}, \quad \sqrt{4-\tau} \sqrt{\tau}, \quad \frac{\sqrt{1-4 \tau}}{\tau}.
\end{equation}

§.§.§ The x-space result

We obtain the following expression for the $O(\ep^0)$ term of the unrenormalized 3-loop two-mass pure singlet operator matrix element:

\begin{eqnarray}
\tilde{a}_{Qq}^{(3), \rm PS}(x) &=& C_F T_F^2 \Biggl\{
R_0(m_1,m_2,x) +\big(\theta(\eta_--x)+\theta(x-\eta_+)\big) x \, g_0(\eta,x)
\nonumber \\ &&
+\theta(\eta_+-x) \theta(x-\eta_-) \biggl[x \, f_0(\eta,x)
-\int_{\eta_-}^x dy \left(f_1(\eta,y)+\frac{x}{y} f_3(\eta,y)\right)\biggr]
\nonumber \\ &&
+\theta(\eta_--x) \int_x^{\eta_-} dy \left(g_1(\eta,y)+\frac{x}{y} g_3(\eta,y)\right)
\nonumber \\ &&
-\theta(x-\eta_+) \int_{\eta_+}^x dy \left(g_1(\eta,y)+\frac{x}{y} g_3(\eta,y)\right)
\nonumber \\ &&
+x \, h_0(\eta,x) +\int_x^1 dy \left(h_1(\eta,y)+\frac{x}{y} h_3(\eta,y)\right)
\nonumber \\ &&
+\theta(\eta_+-x) \int_{\eta_-}^{\eta_+} dy \left(f_1(\eta,y)+ \frac{x}{y} f_3(\eta,y)\right)
\nonumber \\ &&
+\int_{\eta_+}^1 dy \left(g_1(\eta,y)+\frac{x}{y} g_3(\eta,y)\right)
\Biggr\}.
\label{aQq}
\end{eqnarray}

Here we follow the notation used in Ref. [171]. Compared to that notation, in the present case no functions carrying the index 2 occur.
The functions $g_i(\eta,x)$ in Eq. (<ref>) shall not be confounded with polarized structure functions, also often
by $g_i$. Here $\theta(z)$ denotes the Heaviside function

\begin{eqnarray}
\theta(z) = \left\{\begin{array}{ll} 1 &~~z \geq 0\\
0 &~~z < 0
\end{array} \right.
\label{eq:Heaviside}
\end{eqnarray}

by which we divide the interval $x\in [0,1]$ as described earlier. We define for convenience

\begin{equation}
\end{equation}

\begin{equation}
L_1 = \ln\left(\frac{m_1^2}{\mu^2}\right), \quad L_2 = \ln\left(\frac{m_2^2}{\mu^2}\right),
\label{L1L2}
\end{equation}

with $\mu$ the renormalization scale. If in the following expressions the harmonic polylogarithms $H_{\vec{a}}$ are given without argument it is understood that their argument is $x$.
The functions appearing in Eq. (<ref>) are given by

\begin{eqnarray}
\bigg( L_1^3 +L_1  L_2  (L_1  +L_2  ) + L_2^3 \bigg)
\bigg[5 (-1+x) -2 (1+x) \HA_0 \bigg]\nonumber\\
& & +128 L_1  L_2  \bigg[
(x+1) \bigg(\frac{2}{3}\HA_{0,1}-\frac{10}{9}\HA_0-\frac{2}{3}\zeta_2\bigg)
\bigg]\nonumber\\
&&+32 \big(
+ L_2^2
\big)
\bigg[(x+1)\bigg(
\frac{2}{3}  \HA_{0,1}
-\frac{2}{3} \zeta_2
\bigg)
\frac{1}{9}
-\frac{5}{3} \HA_1
\bigg)\nonumber\\
&&        +\frac{1}{9} (17-37 x) \HA_0
\bigg]\nonumber\\
&&+64 (L_1 +L_2 ) \bigg[
(1+x) \bigg(
\Big(
2 \HA_{0,1}
-\frac{8 \zeta_2}{3}
\Big) \HA_0
-\frac{2}{9} \HA_0^3
-\frac{10}{3} \HA_{0,0,1}\nonumber\\
&&        -\frac{4}{3} \HA_{0,1,1}
+\frac{14}{3} \zeta_3
\bigg)
+(x-1) \bigg(
\frac{442}{27}
+\frac{5}{3} \HA_1^2
-\frac{5}{9} \HA_1 \Big(
1+9 \HA_0\Big)
\bigg)\nonumber\\
&&-\frac{2}{27} (56+137 x) \HA_0
+\frac{1}{9} (-5+4 x) \HA_0^2
+\frac{2}{9} (-17+28 x) \HA_{0,1}\nonumber\\
&&+\frac{2}{9} (-28+17 x) \zeta_2
\bigg]
+\frac{64}{1215} \Bigg[(1+x) \bigg(
\Big(
3240 \HA_{0,0,1}
+1620 \HA_{0,1,1}
\Big) \HA_ 0\nonumber\\
&&        +\Big(
-1620 \HA_{0,1}
+945 \zeta_2
\Big) \HA_ 0^2
+90 \HA_ 0^4
-1080 \HA_{0,0,0,1}\nonumber\\
&&        -2700 \HA_{0,0,1,1}
+540 \HA_{0,1,1,1}
+1296 \zeta_2^2
\bigg)+
(-1+x) \bigg(
20 (437+54 x)
&&        \Big(
1080 \HA_ 0
+4050 \HA_ 0^2
+2025 \zeta_2
\Big) \HA_ 1
-225 \HA_ 1^3\nonumber\\
&&        -45 \HA_ 1^2 \Big(
11+45 \HA_ 0\Big)
\bigg)
-10 \big(
-842+1111 x+81 x^2\big)\nonumber\\
&&        -540 (-7+11 x) \HA_{0,1}
-45 (-53+73 x) \zeta_2
-4860 (1+x) \zeta_3
\Big) \HA_ 0\nonumber\\
&&+165 (19+37 x) \HA_ 0^2
-30 (-19+8 x) \HA_ 0^3
+30 (-1+x) (157+27 x) \HA_ 1\nonumber\\
-30 (61+169 x)
-810 (1+x) \zeta_2
\Big) \HA_{0,1}
+180 (-11+25 x) \HA_{0,0,1}\nonumber\\
&&+180 (-14+13 x) \HA_{0,1,1}
+15 (131+329 x) \zeta_2
+90 (-55+29 x) \zeta_3
\Bigg],
\\
-\frac{32 (1-x)}{9} \Bigg[
-\frac{16 (-1+x) x}{\eta }
+18 \bigg(
-\frac{2 (\eta
-4 (-1+x) x
)^2}{9 \eta ^2}
+\frac{1}{3} \zeta_2
\bigg) \HA_0\big(
\nonumber\\
&&        +5 \HA_0^2\big(
+2 \bigg(
-4 (1-x) x
)^{3/2}}{\eta ^{3/2}}
\bigg) \zeta_2
\Bigg]\nonumber\\
&&-\frac{64 (1-x)}{9} \Bigg[
\bigg(
2 \frac{(\eta
-4 (1-x) x
)^{3/2}}{\eta ^{3/2}}
-3 \zeta_2
\bigg) G\left(
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)\Bigg]\nonumber\\
&&+\frac{64 (1-x)}{3} \left[
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)
\left\{\frac{\sqrt{1-4 \tau }}{\tau },\frac{\sqrt{1-4 \tau \
}}{\tau },\frac{1}{\tau }\right\},u\right)\right]\nonumber\\
-\frac{64 (1-x)}{9}
\frac{\sqrt{1-4 \tau }}{\tau },\frac{1}{\tau \
\frac{(\eta
-4 (1-x) x
)^{3/2}}{\eta ^{3/2}},\\
g_1(\eta,x)&=&\frac{64}{27 \eta ^2 x} \bigg[
-6 (1-x) \HA_ 0\big(
u\big) P_ 2
-8 \eta  (-1+x) (1+x) (7 \eta
+24 (1-x) x
&&        +3 \eta ^2 (-1+x) (-5+13 x) \HA_ 0^2\big(
-3 \eta ^2 (1-x) (-1+2 x) \HA_ 0^3\big(
&&        -(6 (1-x)) \bigg(
(1+x) \eta ^{3/2}
-4 \eta  (1+x) \sqrt{\eta
-4 (1-x) x
+2 (-1+x) x (1+10 x) \nonumber\\
&&           \times     \sqrt{\eta
-4 (1-x) x
\bigg) \sqrt{\eta } \zeta_2
\bigg]
+\frac{128 (-1+x)}{9 x} \bigg[
\Big(
-4 \frac{\sqrt{\eta
-4 (1-x) x
}}{\eta ^{3/2}} P_ 1\nonumber\\
&&                -3 (-1+2 x) \zeta_2
\Big) G\left(
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)\bigg]\nonumber\\
&&-\frac{128 (-1+x) (-1+2 x)}{3 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)\nonumber\\
&&-\frac{128 (-1+x) (-1+2 x)}{3 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau },\frac{\sqrt{1-4 \tau \
}}{\tau },\frac{1}{\tau }\right\},u\right)\nonumber\\
&&-\frac{256 (-1+x) P_ 1}{9 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau },\frac{1}{\tau \
\frac{\sqrt{\eta
-4 (1-x) x
}{\eta ^{3/2}},\\
g_3(\eta,x)&=&-\frac{32}{27 \eta ^2 x} \bigg[
8 \eta  (1-x) P_4
-6 (1-x) \HA_0\big(
u\big) P_5
+3 \eta ^2 (-1+x) (-5+8 x)\nonumber\\
&&         \times \HA_0^2\big(
+3 \eta ^2 (-1+x)^2 \HA_0^3\big(
-(6 (1-x)) \bigg(
(1+2 x) \eta ^{3/2}\nonumber\\
&&                -\eta  (4+5 x) \sqrt{\eta
-4 (1-x) x
+2 (-1+x) x (1+8 x) \sqrt{\eta
-4 (1-x) x
\bigg) \sqrt{\eta } \zeta_2
\bigg]\nonumber\\
&&-\frac{64 (1-x)}{9 x} \bigg[
\bigg(
2 \frac{\sqrt{\eta
-4 (1-x) x
}}{\eta ^{3/2}} P_3
+3 (-1+x) \zeta_2
\bigg) G\left(
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)\bigg]\nonumber\\
&&+\frac{64 (-1+x)^2}{3 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau }\right\},u\right)\nonumber\\
&&+\frac{64 (-1+x)^2}{3 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau },\frac{\sqrt{1-4 \tau \
}}{\tau },\frac{1}{\tau }\right\},u\right)\nonumber\\
&&-\frac{64 (1-x) P_3}{9 x}
\left\{\frac{\sqrt{1-4 \tau }}{\tau },\frac{1}{\tau \
}\right\},u\right) \frac{\sqrt{\eta
-4 (1-x) x
}}{\eta ^{3/2}},
\end{eqnarray}

with the polynomials

\begin{eqnarray}
P_1&=&2 \eta  (x+1)-10 x^3+9 x^2+x,\\
P_2&=&3 \eta ^2 (2 x \zeta_2+x-\zeta_2+1)+8 \eta  x \left(10 x^2-9 \
x-1\right)-16 (1-x)^2 x^2 (10 x+1),\\
P_3&=&\eta  (5 x+4)+2 x \left(-8 x^2+7 x+1\right),\\
P_4&=&7 \eta  (x+1)+6 x \left(-5 x^2+x+4\right),\\
P_5&=&\eta ^2 (x (3 \zeta_2+5)-3 \zeta_2+3)+8 \eta  x \left(8 x^2-7 \
x-1\right)-16 (x-1)^2 x^2 (8 x+1),
\end{eqnarray}

\begin{eqnarray}
\frac{16 (1-x)}{3}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau \
-\frac{4 P_6}{9 (-1+x) x^2} \Big[
+2 \HA_0\big(
v\big)\Big] \nonumber\\
&&                       \times \frac{(-\eta
+4 (1-x) x
)^{3/2}}{\eta ^{3/2}}
\bigg] G\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)\nonumber\\
&&+\frac{4 (1-x)}{3} \bigg[
\Big[
-1+2 \HA_0\big(
v\big)\Big] G^2\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)\bigg]\nonumber\\
&&+\frac{16 (1-x)}{3}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau },\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)
+\frac{1}{18} \bigg[
-1536 (1-x)\nonumber\\
&&        -\frac{9 \eta ^4}{(-1+x)^3 x^4}
-\frac{80 \eta ^3}{(-1+x)^2 x^3}
-\frac{104 \eta ^2}{(-1+x) x^2}
+\frac{576 \eta }{x}
+\frac{4  P_7}{(-1+x)^3 x^4}\HA_0\big(
&&        -320 (1-x) \HA_0^2\big(
+64 (1-x) \HA_0^3\big(
-\Big(128 (1-x)\Big) \Big(
-3 \HA_0\big(
v\big)\Big) \zeta_2\nonumber\\
&&        +768 (1-x) \zeta_3
\bigg]
-\frac{8 P_6}
{9 (1-x) x^2}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau
}\right\},v\right) \nonumber\\
+4 (1-x) x
)^{3/2}}{\eta ^{3/2}},\\
f_1(\eta,x)&=&\frac{1}{27 x^5} \Bigg[
-\frac{1}{(1-x)^3} \bigg[
6912 \eta  (-1+x)^3 x^4
-27 \eta ^4 (-1+2 x)
-240 \eta ^3 (-1+x) x \nonumber\\
&&                (-1+2 x)
+512 (-1+x)^4 x^4 (-16+11 x)
+12 \HA_ 0\big(
v\big) P_ 9\nonumber\\
&&                -24 \eta ^2 (-1+x)^2 x^2 (-25+2 x)
+192 (-1+x)^4 x^4 (-5+13 x) \HA_ 0^2\big(
&&                -192 (-1+x)^4 x^4 (-1+2 x) \HA_ 0^3\big(
\bigg]
+384 (1-x) x^4 \bigg(
-13 x\nonumber\\
&&                +3 (-1+2 x) \HA_ 0\big(
\bigg) \zeta_2
-2304 (-1+x) x^4 (-1+2 x) \zeta_3
\Bigg]
\frac{32 (-1+x) (-1+2 x)}{3 x}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau
&&        -\frac{8 P_ 8}{9 (1-x) x^3} \bigg(
-1+2 \HA_ 0\big(
v\big)\bigg) \frac{\sqrt{-
\eta
+4 (1-x) x
}}{\eta ^{3/2}}
\bigg] \nonumber\\
&& \times G\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)
-\frac{8 (-1+x) (-1+2 x)}{3 x} \bigg[
\bigg(
-1+2 \HA_ 0\big(
v\big)\bigg) \nonumber\\
&& \times G^2\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)\bigg]
+\frac{32 (1-x) (-1+2 x)}
{3 x} \nonumber\\
\times G\left(
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau },\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)
-\frac{16 P_ 8}{9 (-1+x) x^3} \nonumber\\
&&        \times G\left(
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau
}\right\},v\right) \frac{\sqrt{-\eta
+4 (1-x) x
}}{\eta ^{3/2}},\\
f_3(\eta,x)&=&\frac{1}{54 (-1+x)^2 x^5} \bigg[
27 \eta ^4
-240 \eta ^3 (1-x) x
-5184 \eta  (-1+x)^2 x^4
-1024 (-8+x) \nonumber\\
&&        \times (-1+x)^3 x^4
-12 \HA_ 0\big(
v\big) P_{10}
-24 \eta ^2 (-1+x) x^2 (25+11 x)\nonumber\\
&&        -192 (-1+x)^3 x^4 (-5+8 x) \HA_ 0^2\big(
+192 (-1+x)^4 x^4 \HA_ 0^3\big(
&&        +384 (-1+x)^3 x^4 \bigg(
-8 x
-3 (1-x) \HA_ 0\big(
\bigg) \zeta_2
+2304 (-1+x)^4 x^4 \zeta_3
\bigg]\nonumber\\
-\frac{16 (-1+x)^2}{3 x}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau
-\frac{4 P_{11}}{9 x^3} \bigg(
&&                +2 \HA_ 0\big(
v\big)\bigg) \frac{\sqrt{-
\eta
+4 (1-x) x
}}{\eta ^{3/2}}
\bigg] G\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)\nonumber\\
&&+\frac{4 (-1+x)^2}{3 x}
\bigg[
\bigg(
-1+2 \HA_ 0\big(
v\big)\bigg) G\left(
\left\{\sqrt{4-\tau } \sqrt{\tau }\right\},v\right)^2\bigg] \nonumber\\
&&+\frac{16 (-1+x)^2}{3 x}
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau },\sqrt{4-\tau } \sqrt{\tau }
\right\},v\right)
+\frac{8 P_{11}}{9 x^3} \nonumber\\
&& \times       G\left(
\left\{\frac{1}{\tau },\sqrt{4-\tau } \sqrt{\tau
}\right\},v\right) \frac{\sqrt{-\eta
+4 (1-x) x
}}{\eta ^{3/2}},
\end{eqnarray}

\begin{eqnarray}
P_6    &=&  3 \eta ^2+6 \eta  (1-x) x+4 (x-1)^2 x^2,\\
P_7    &=&  3 \eta ^4-24 \eta ^3 (1-x) x+20 \eta ^2 (x-1)^2 x^2-160 \eta  \
(x-1)^3   x^3-128 (x-1)^4 x^4,\\
P_8    &=&  \eta ^3 (6 x-3)+6 \eta ^2 x \left(2 x^2-3 x+1\right)-8 \eta  \
(x-1)^2 x^2 (8 x-1)\nonumber\\
& &  +8 (x-1)^3 x^3 (10 x+1),\\
P_9    &=&  \eta ^4 (6 x-3)+24 \eta ^3 x \left(2 x^2-3 x+1\right)-4 \eta ^2 (x-1)^2 x^2 (2 x+11)\nonumber\\
& &  -64 \eta  (x-1)^3 x^3 (8 x-1)+96 (x-1)^4 x^4 (x+4),\\
P_{10} &=&  3 \eta ^4-24 \eta ^3 (1-x) x-4 \eta ^2 x^2 \left(7 x^2+4 x-11\right)-32 \eta  (x-1)^2 x^3 (11 x-2)\nonumber\\
& &  +32 (x-1)^3 x^4 (7 x+12),\\
P_{11} &=&  3 \eta ^3-6 \eta ^2 (1-x) x-4 \eta  x^2 \left(11 x^2-13 x+2\right)+8 (1-x)^2 x^3 (8 x+1).
\end{eqnarray}

The functions $h_i$ are defined as follows

\begin{eqnarray}
h_i(\eta,x) = g_i\left(\frac{1}{\eta},x\right), \quad i=0,1,3.
\end{eqnarray}

[] The ratio of the 2-mass contributions to the massive OME $A_{Qq}^{{\sf PS},(3)}$ to all contributions to
$A_{Qq}^{{\sf PS},(3)}$
of $O(T_F^2)$ as a function of $x$ and $\mu^2$.
Dotted line (red): $\mu^2 = 30~\GeV^2$.
Dashed line (black): $\mu^2 = 50~\GeV^2$.
Dash-dotted line (blue): $\mu^2 = 100~\GeV^2$.
Full line (green): $\mu^2 = 1000~\GeV^2$. Here the on-shell heavy quark masses $m_c = 1.59~\GeV$ and $m_b = 4.78~\GeV$
[286, 287] have been used, from [196].

Figure <ref> shows the ratio of the two-mass contribution to the complete $O(T_F^2C_F)$ term. The two-mass correction grows relatively larger with $\mu^2$ and can reach the order of 40% of the total $O(T_F^2C_F)$ contribution.

The result is difficult to obtain in Mellin-$N$ space and we refrain from this. The main reason for this is the appearance of the Heaviside functions in (<ref>). In general one will obtain recursions not factorizing in first order.

§.§ The two-mass contribution to the polarized matrix element Agg,Q(3)

The calculation of the OME $A_{gg,Q}^{(3)}$ proceeds from the diagrams of Figure <ref>.

The 11 different topologies for $\tilde{A}_{gg,Q}^{(3)}$. Curly lines: gluons; thin arrow
lines: lighter massive quark; thick arrow lines: heavier massive quark; the symbol $\otimes$ represents the corresponding
local operator insertion, cf. [81] for the related Feynman rules.

The renormalization of the OME is performed in [121] and provides a check by predicting the pole terms of the full calculation.

We perform the calculation in $N$-space following the methods of the unpolarized calculation [172], which will be described in what follows. Then we obtain the $z$-space result by an inverse Mellin transform.
The Feynman rules applied to the operator insertion are found in [81], and to the sum of the Feynman diagrams we apply the projector [158, 175]

\begin{equation}\label{eq:projector}
\frac{1}{(D-2)(D-3)}
(\Delta\cdot p)^{-N-1}\ep^{\mu\nu\rho\sigma}
\Delta \hat{G}^{ab}_{Q,\mu\nu}
\Delta_{\rho}p_{\sigma}\,.
\end{equation}

We adopt the Larin scheme. The Feynman parametrization used for each diagram is chosen such that, after performing the Dirac algebra with FORM [281], numerator structures are not cancelled against denominators. At the price of having to deal with more complicated denominator structures, this allows us to reduce their number, which makes it easier to manipulate these terms.

We use the Feynman parametrization

\begin{eqnarray}
\Pi^{\mu\nu}_{ab} ( k ) = -i \frac{ 8 T_F g^2}{(4\pi)^{D/2}} \delta_{ab} ( k^2 g^{\mu\nu} - k^\mu k^\nu )
\int\limits_{0}^{1} \text{d}x \frac{ \Gamma(2-D/2) \left( x(1-x) \right)^{D/2-1} }
{ \left( -k^2 + \frac{m^2}{x(1-x)} \right)^{2-D/2} }
\end{eqnarray}

for massive quark bubbles. After the Feynman parametrization is obtained for the whole diagram, the numerator algebra is performed through the identities

\begin{eqnarray}
% \int \frac{\dd^d k}{(2\pi)^d} \frac{ 1 }{ k^2 + R^2 }  &=&  , \\
\int \frac{\dd^D k}{(2\pi)^D}  k_{\mu_1} k_{\mu_2} f(k^2)
&=&  \frac{g_{\mu_1\mu_2}}{D} \int \frac{\dd^D k}{(2\pi)^D}  k^2 f(k^2), \\
\int \frac{\dd^D k}{(2\pi)^D}  k_{\mu_1} k_{\mu_2} k_{\mu_3} k_{\mu_4} f(k^2)
&=& \frac{S_{\mu_1\mu_2\mu_3\mu_4}}{D(D+2)}  \int \frac{\dd^D k}{(2\pi)^D} (k^2)^2 f(k^2),
\nonumber\\
\\
\int \frac{\dd^D k}{(2\pi)^D}  k_{\mu_1} k_{\mu_2} k_{\mu_3} k_{\mu_4} k_{\mu_5} k_{\mu_6} f(k^2)
&=& \frac{S_{\mu_1\mu_2\mu_3\mu_4\mu_5\mu_6}}{D(D+2)(D+4)}
\int \frac{\dd^D k}{(2\pi)^D}  (k^2)^3 f(k^2),
\end{eqnarray}

with the symmetric tensors

\begin{eqnarray}
S_{\mu_1\mu_2\mu_3\mu_4} &=& g_{\mu_1\mu_2}g_{\mu_3\mu_4} + g_{\mu_1\mu_3}g_{\mu_2\mu_4}
+ g_{\mu_1\mu_4} g_{\mu_2\mu_3}
\\
S_{\mu_1\mu_2\mu_3\mu_4\mu_5\mu_6} &=&
\left[
g_{\mu_3\mu_4} g_{\mu_5\mu_6}
+ g_{\mu_3\mu_5} g_{\mu_4\mu_6}
+ g_{\mu_3\mu_6} g_{\mu_4\mu_5}\right]
\nonumber\\
\left[
g_{\mu_2\mu_4} g_{\mu_5\mu_6}
+ g_{\mu_2\mu_5} g_{\mu_4\mu_6}
+ g_{\mu_2\mu_6} g_{\mu_4\mu_5}\right]
\nonumber\\
\left[
g_{\mu_2\mu_3} g_{\mu_5\mu_6}
+ g_{\mu_2\mu_5} g_{\mu_3\mu_6}
+ g_{\mu_2\mu_6} g_{\mu_3\mu_5}\right]
\nonumber\\
\left[
g_{\mu_2\mu_3} g_{\mu_4\mu_6}
+ g_{\mu_2\mu_4} g_{\mu_3\mu_6}
+ g_{\mu_2\mu_6} g_{\mu_3\mu_4}\right]
\nonumber\\
\left[
g_{\mu_2\mu_3} g_{\mu_4\mu_5}
+ g_{\mu_2\mu_4} g_{\mu_3\mu_5}
+ g_{\mu_2\mu_5} g_{\mu_3\mu_4}\right].
\end{eqnarray}

The scalar integrals can then be performed using the relation

\begin{eqnarray}
\int \frac{\dd^D k}{(2\pi)^D} \frac{ (k^2)^m }{ ( k^2 + R^2 )^n }
&=& \frac{1}{(4\pi)^{D/2}} \frac{\Gamma(n-m-D/2)}{\Gamma(n)} \frac{\Gamma(m+D/2)}{\Gamma(D/2)}
\left( R^2 \right)^{m-n+D/2} .
\end{eqnarray}

After the integrals over the loop momenta are performed, only integrals over the Feynman parameters are left. They always appear in the form

\begin{equation}\label{eq:feynIntegrals}
\prod\limits_{i=1}^{j} \int\limits_0^1 \mathrm{d} x_i \ x_i^{a_i} (1-x_i)^{b_i} \ R_0^N \ \left[ R_1
\ m_1^2 + R_2 \ m_2^2 \right]^{-s},
\end{equation}

where $R_0$ is a polynomial in the Feynman parameters $x_i$, and $R_1$ and $R_2$ are rational functions
in $x_i$.

At this point the polynomial $R_0$ can be treated by applying the binomial theorem (multiple times if necessary)

\begin{equation}
(A + B)^N = \sum\limits_{i=0}^{N} \binom{N}{i} A^i B^{N-i}
\label{eq:I2}
\end{equation}

while the factor $\left[ R_1 \ m_1^2 + R_2 \ m_2^2 \right]^{-s}$ is treated via a Mellin-Barnes decomposition (<ref>).

The integrals over the Feynman parameters are then turned into infinite sums using the residue theorem. This sum representation is treated analytically by the algorithms of Sigma, HarmonicSums, EvaluateMultiSums and SumProduction, which reduce these nested sums into classes of hypergeometric sums which are shown, algorithmically, to be independent.

The target function space is, for the $N$-space result, that of harmonic sums, defined in Eqs. (<ref>) and (<ref>), generalized harmonic sums (<ref>),
cyclotomic sums [249] and binomial sums [252]. Harmonic polylogarithms [266] will also appear in the result.

The calculation was checked by computing fixed $N$-moments using MB and MBResolve and comparing the results with those obtained from the packages Q2E/EXP [288, 289]. In $N$-space, the result appears to exhibit spurious poles for $N=1/2,~3/2$. It has been verified analytically that the expression is actually regular at those points.

§.§.§ The N-space solution

We obtain, for the constant part in $\ep$ of the $N$-space $\tilde{A}_{gg,Q}^{(3)}$,

where the argument of the harmonic polylogarithms is implied to be $\eta$ where omitted, and that of the harmonic sums to be $N$ where omitted.

The polynomials $P_i$ are:

§.§.§ The result in momentum fraction $z$-space

From the $N$-space result, the Mellin inversion has been performed using the algorithms encoded in HarmonicSums. The result, which we reproduce here from [198], is reported as a one-dimensional integral, in a form amenable to numerical evaluation. It depends on iterated integrals of the type (<ref>) on an alphabet of root-valued letters. The iterated integrals are the same as those appearing in the unpolarized calculation, and can be found in Appendix D of [172].