

\section{Introduction}
\label{sec:1}

\vspace{1mm}
\noindent
The heavy flavor corrections to deeply inelastic scattering exhibit different scaling violations
if compared to the light flavor contributions. At the present experimental accuracy of the
deep--inelastic data, precision determinations of the strong coupling constant $\alpha_s(M_Z^2)$ 
\cite{ALPHA} and the determination of the parton distribution functions \cite{PDF} require their 
calculation to next-to-next-to leading order (NNLO). If the virtuality, $Q^2$, is significantly larger 
than the heavy quark mass squared, $m_Q^2$, i.e. $Q^2/m_Q^2\raisebox{-0.07cm    10$, one may calculate the heavy 
flavor corrections analytically \cite{Buza:1995ie}. 
Here $m_Q$ denotes the heavy quark mass. This corresponds 
to $Q^2$ values above about $25~\rm GeV^2$ in the case of charm. In Refs.~\cite{Bierenbaum:2009mv,
Ablinger:2017err,Behring:2014eya} it has been shown to NNLO in the single-- and double--mass cases, 
how to express the massive coefficient functions for the structure functions in terms of 
massive operator matrix elements (OMEs) $A_{ij}$ and the massless coefficient functions 
\cite{WILS3} 
at large scales, as well as the corresponding relations for the parton distribution functions 
in the variable flavor number scheme (VFNS)\footnote{The two--loop corrections can be found in 
Refs.~\cite{Buza:1995ie,TWOL,DELG2,Buza:1996wv}.}. From NNLO onward, five of the OMEs 
\cite{THREE,Ablinger:2014nga} are needed to express the deep--inelastic structure functions \cite{SF}. 
The additional OMEs, $A_{gg,Q}$ and $A_{gq,Q}$, at NNLO contribute only to the relations in the 
VFNS. 
The three--loop corrections to $A_{gq,Q}$ have been computed in Ref.~\cite{Ablinger:2014lka}.

In this  paper we calculate the massive OMEs $A_{gg,Q}$ and $\Delta A_{gg,Q}$ to three--loop order. 
In Refs.~\cite{Blumlein:2012vq,Ablinger:2014uka} the $O(T_F^2 N_F)$ and $O(T_F^2)$ contributions to 
the operator matrix element  have been calculated in the unpolarized case already. All logarithmic 
contributions have 
been derived in Refs.~\cite{Behring:2014eya,Blumlein:2021xlc}. A series of moments and scalar 
integrals in the double--mass case were obtained in Refs.~\cite{Ablinger:2017err,TWOMASS}. 
The $O(T_F)$ contribution to the three--loop anomalous dimensions 
$(\Delta)\gamma_{gg,Q}^{(2)}$, resulting 
from the single pole term of the unrenormalized OME, has been 
calculated in 
Refs.~\cite{ Ablinger:2017tan,Ablinger:2019etw,Blumlein:2021ryt}. The double--mass 
three--loop contributions were computed in Refs.~\cite{Ablinger:2018brx,Ablinger:2020snj}. Here both charm and 
bottom quark lines are contained in the corresponding Feynman diagrams. Contributions of this kind emerge
first at two--loop order due to reducible diagrams \cite{Blumlein:2018jfm,DELG2}.

We calculate the massive OME using the techniques described in Ref.~\cite{Ablinger:2015tua}.
These include the method of (generalized) hypergeometric functions \cite{HYP}, the Mellin--Barnes method 
\cite{MB}, the method of ordinary differential equations \cite{DIFF,Ablinger:2018zwz}, and the 
(multivartiate) Almkvist-Zeilberger algorithm \cite{AZ,Ablinger:21}. All these methods finally map the 
problem to multiply nested sums, which are solved using the packages {\tt Sigma} \cite{SIG1,SIG2}, 
{\tt EvaluateMultisums} and {\tt SumProduction} \cite{EMSSU}, based on difference--ring theory 
\cite{Karr:81,Schneider:01,Schneider:05a,Schneider:07d,Schneider:10b,Schneider:10c,Schneider:15a,
Schneider:08c}, as well as the packages {\tt OreSys} \cite{ORESYS}, and {\tt MultiIntegrate} \cite{MI,
Ablinger:21}.\footnote{For surveys on the computation methods and the related 
function and number spaces, see Refs.~\cite{Blumlein:2018cms,Blumlein:2022zkr}.} For an efficient 
treatment of the 
different sum- and function algebras, such as the harmonic sums \cite{Vermaseren:1998uu,Blumlein:1998if} 
and harmonic polylogarithms \cite{Remiddi:1999ew}, generalized harmonic sums and polylogarithms 
\cite{Moch:2001zr,Ablinger:2013cf}, cyclotomic harmonic sums and polylogarithms 
\cite{Ablinger:2011te}, iterated integrals induced by quadratic forms \cite{Ablinger:2021fnc}, 
and finite 
binomial sums and square root--valued iterated integrals \cite{Ablinger:2014bra}, 
we use the package {\tt HarmonicSums} \cite{HARMSU,Blumlein:2009fz,
Blumlein:2003gb,Ablinger:2011te,Ablinger:2013cf,Ablinger:2014bra,
Vermaseren:1998uu,Blumlein:1998if,Remiddi:1999ew,Blumlein:2009ta,Ablinger:2021fnc}.
We mention, in particular, that our solution methods do not require to refer to special bases of 
master integrals.

One obtains a general expression for the constant part of the unrenormalized 
massive OMEs $(\Delta) A_{gg,Q}^{(3)}$, $(\Delta) a_{gg,Q}^{(3)}$, at even (odd) integer values of the 
Mellin variable $N$, according to the light--cone expansion in the unpolarized and polarized cases, 
cf.~\cite{Blumlein:2012bf}. The representation is given in terms of synchronized sum--product 
representations. In the present case it turned out that all but two diagrams could be calculated in 
this way and for the latter cases the problem had to be solved in $x$--space, deriving the $N$--space 
representation afterwards.

We consider the case of a single heavy quark. With the present results, both the unpolarized 
and polarized gluon distributions can be mapped in the single--mass VFNS \cite{Buza:1996wv} 
and the double--mass VFNS \cite{Ablinger:2017err} to three--loop order.

The paper is organized as follows. In Section~\ref{sec:2} we present the basic formalism and 
describe the new functional aspects in the present case, both in Mellin $N$ and $x$--space.
In Section~\ref{sec:3} the results in Mellin $N$-space are presented for 
$a_{gg,Q}^{(3)}$ and $\Delta a_{gg,Q}^{(3)}$. The characteristic aspects of the $x$--space 
results are discussed in Section~\ref{sec:4}. Because the full $x$--space results 
are rather lengthy, we
describe its principal structure and  
present it in full form in a  computer--readable ancillary file to this paper.
Numerical results  are discussed in Section~\ref{sec:5} deriving a fast and precise numerical 
representation,
and Section~\ref{sec:6} contains 
the conclusions. A series of technical aspects are presented in the Appendices \ref{sec:A}--\ref{sec:D}.
\section{Basic Formalism and the computation method}
\label{sec:2}

\vspace{1mm}
\noindent
The unrenormalized massive OMEs $A_{gg,Q}^{(k)}$ from one-- to three--loop order for a 
single heavy quark 
and 
$N_F$ massless quarks, cf.~Ref.~\cite{Bierenbaum:2009mv}, is given by
\begin{eqnarray}
    \hat{\hspace*{-1mm}\hat{A}}_{gg,Q}^{(1)}&=& 
             \Bigl(\frac{\hat{m}^2}{\mu^2}\Bigr)^{\varepsilon/2}\Biggl[
                          \frac{\hat{\gamma}_{gg}^{(0)}}{\varepsilon}
                         +a_{gg,Q}^{(1)}
                         +\varepsilon\overline{a}_{gg,Q}^{(1)}
                         +\varepsilon^2\overline{\overline{a}}_{gg,Q}^{(1)}
                        \Biggr] + O(\varepsilon^3),
\\
    \hat{\hspace*{-1mm}\hat{A}}_{gg,Q}^{(2)}&=&
             \Bigl(\frac{\hat{m}^2}{\mu^2}\Bigr)^{\varepsilon}
                  \Biggl[
                          \frac{1}{\varepsilon^2} c_{gg,Q,(2)}^{(-2)}
                         + \frac{1}{\varepsilon} c_{gg,Q,(2)}^{(-1)}
                         +  c_{gg,Q,(2)}^{(0)}
                         + \varepsilon c_{gg,Q,(2)}^{(1)} \Biggr]+ O(\varepsilon^2),
\\
    \hat{\hspace*{-1mm}\hat{A}}_{gg,Q}^{(3)}&=&
             \Bigl(\frac{\hat{m}^2}{\mu^2}\Bigr)^{3\varepsilon/2}
                  \Biggl[
                          \frac{1}{\varepsilon^3} c_{gg,Q,(3)}^{(-3)}
                         + \frac{1}{\varepsilon^2} c_{gg,Q,(3)}^{(-2)}
                         + \frac{1}{\varepsilon} c_{gg,Q,(3)}^{(-1)}
                         + a_{gg,Q}^{(3)}\Biggr] + O(\varepsilon)~.
 \end{eqnarray}
Here and in the following $\varepsilon = D - 4$ denotes the dimensional parameter, 
$c_{gg,Q,(k)}^{(-l)}$ are the expansion coefficients of the unrenormalized OME $\hat{\hspace*{-1mm}\hat{A}}_{gg,Q}$,
$\hat{m}$ denotes the unrenormalized heavy quark 
mass, $\mu$ is both the factorization and
renormalization scale, $\zeta_i$ is Riemann's $\zeta$-function at integer values $i \geq 2$, 
$\beta_k$ and $\beta_{k,Q}$ are expansion coefficients of the  $\beta$-function in Quantum Chromodynamics 
(QCD) in different 
schemes, $\gamma_{ij}^{(k)}$ are anomalous 
dimensions ~\cite{ANOM,Ablinger:2014nga,Ablinger:2017tan,Blumlein:2021ryt,WILS3}
and $m_k^{(l)}$ are the expansion coefficients
of the heavy mass, while $a_{ij}^{(k)}$ and $\bar{a}_{ij}^{(k)}$  
are the expansion coefficients in $\varepsilon$ to $O(\varepsilon^0)$ and $O(\varepsilon)$ of the different 
OMEs. All explanations have been given in Ref.~\cite{Bierenbaum:2009mv}, to which we refer.
Analogous expressions hold for $\Delta A_{gg,Q}^{(k)}$. Here the anomalous dimensions have to be replaced
by the polarized ones, and the coefficients $a_{ij}^{(k)} (\bar{a}_{ij}^{(k)})$ by $\Delta a_{ij}^{(k)}
(\Delta \bar{a}_{ij}^{(k)})$.

The renormalization of $\hat{\hspace*{-1mm}\hat{A}}_{gg,Q}^{(3)}$ includes mass and coupling renormalization, as well as 
the renormalization of the local operator, and the subtraction of the collinear
singularities, see~\cite{Bierenbaum:2009mv}. It is essential to work in a MOM-scheme
for charge renormalization first and then transform to the $\overline{\rm MS}$ scheme.
The renormalized OME from first to third order are then obtained by
\begin{eqnarray}
    A_{gg,Q}^{(1), \overline{{\sf MS}}}&=& - \beta_{0,Q} 
\ln\left(\frac{m^2}{\mu^2}\right)~,\\
    A_{gg,Q}^{(2), \overline{{\sf MS}}}&=&
      \frac{1}{8}\Biggl\{
                          2\beta_{0,Q}
                             \Bigl(
                                    \gamma_{gg}^{(0)}
                                   +2\beta_0 
                             \Bigr)
                         +\gamma_{gq}^{(0)}\hat{\gamma}_{qg}^{(0)}
                         +8\beta_{0,Q}^2
                 \Biggr\}
                    \ln^2 \Bigl(\frac{m^2}{\mu^2}\Bigr)
                +\frac{\hat{\gamma}_{gg}^{(1)}}{2}
                     \ln \Bigl(\frac{m^2}{\mu^2}\Bigr)
\nonumber\\ &&
                -\frac{\zeta_2}{8}\Bigl[
                                         2\beta_{0,Q}
                                            \Bigl(
                                                   \gamma_{gg}^{(0)}
                                                  +2\beta_0
                                            \Bigr)
                                        +\gamma_{gq}^{(0)}
                                            \hat{\gamma}_{qg}^{(0)}
                                  \Bigr]
                +a_{gg,Q}^{(2)}~,   \label{AggQ2MSren} 
\\
  A_{gg,Q}^{(3), \overline{{\sf MS}}}&=&
                    \frac{1}{48}\Biggl\{
                            \gamma_{gq}^{(0)}\hat{\gamma}_{qg}^{(0)}
                                \Bigl(
                                        \gamma_{qq}^{(0)}
                                       -\gamma_{gg}^{(0)}
                                       -6\beta_0
                                       -4n_f\beta_{0,Q}
                                       -10\beta_{0,Q}
                                \Bigr)
                           -4
                                \Bigl(
                                        \gamma_{gg}^{(0)}\Bigl[
                                            2\beta_0
                                           +7\beta_{0,Q}
                                                         \Bigr]
\nonumber\\ &&
                                       +4\beta_0^2
                                       +14\beta_{0,Q}\beta_0
                                       +12\beta_{0,Q}^2
                                \Bigr)\beta_{0,Q}
                     \Biggr\}
                     \ln^3 \Bigl(\frac{m^2}{\mu^2}\Bigr)
                    +\frac{1}{8}\Biggl\{
                            \hat{\gamma}_{qg}^{(0)}
                                \Bigl(
                                        \gamma_{gq}^{(1)}
                                       +(1-n_f)\hat{\gamma}_{gq}^{(1)}
                                \Bigr)
\nonumber\\ &&
                           +\gamma_{gq}^{(0)}\hat{\gamma}_{qg}^{(1)}
                           +4\gamma_{gg}^{(1)}\beta_{0,Q}
                           -4\hat{\gamma}_{gg}^{(1)}[\beta_0+2\beta_{0,Q}]
                           +2\gamma_{gg}^{(0)}\beta_{1,Q}
                           +4[\beta_1+\beta_{1,Q}]\beta_{0,Q}
                     \Biggr\}
                     \ln^2 \Bigl(\frac{m^2}{\mu^2}\Bigr)
\nonumber\\ &&
                    +\frac{1}{16}\Biggl\{
                            8\hat{\gamma}_{gg}^{(2)}
                           -8n_fa_{gq,Q}^{(2)}\hat{\gamma}_{qg}^{(0)}
                           -16a_{gg,Q}^{(2)}(2\beta_0+3\beta_{0,Q})
                           +8\gamma_{gq}^{(0)}a_{Qg}^{(2)}
                           +8\gamma_{gg}^{(0)}\beta_{1,Q}^{(1)}
\nonumber\\ &&
                   +\gamma_{gq}^{(0)}\hat{\gamma}_{qg}^{(0)}\zeta_2
                                \Bigl(
                                        \gamma_{gg}^{(0)}
                                       -\gamma_{qq}^{(0)}
                                       +6\beta_0
                                       +4n_f\beta_{0,Q}
                                       +6\beta_{0,Q}
                                \Bigr)
\nonumber\\ &&
                   +4\beta_{0,Q}\zeta_2
                                \Bigl( 
                                       \gamma_{gg}^{(0)}
                                      +2\beta_0
                                \Bigr)
                                \Bigl(
                                       2\beta_0
                                      +3\beta_{0,Q}
                                \Bigr)
                     \Biggr\}
                     \ln \Bigl(\frac{m^2}{\mu^2}\Bigr)
                   +2(2\beta_0+3\beta_{0,Q})\overline{a}_{gg,Q}^{(2)}
\nonumber\\ &&
                   +n_f\hat{\gamma}_{qg}^{(0)}\overline{a}_{gq,Q}^{(2)}
                   -\gamma_{gq}^{(0)}\overline{a}_{Qg}^{(2)}
                   -\beta_{1,Q}^{(2)} \gamma_{gg}^{(0)}
                   +\frac{\gamma_{gq}^{(0)}\hat{\gamma}_{qg}^{(0)}\zeta_3}{48}
                                \Bigl(
                                        \gamma_{qq}^{(0)}
                                       -\gamma_{gg}^{(0)}
                                       -2[2n_f+1]\beta_{0,Q}
\nonumber\\ &&
                                       -6\beta_0
                                \Bigr)
                   +\frac{\beta_{0,Q}\zeta_3}{12}
                                \Bigl(
                                        [\beta_{0,Q}-2\beta_0]\gamma_{gg}^{(0)}
                                       +2[\beta_0+6\beta_{0,Q}]\beta_{0,Q}
                                       -4\beta_0^2
                                \Bigr)
\nonumber\\ &&
                   -\frac{\hat{\gamma}_{qg}^{(0)}\zeta_2}{16}
                                \Bigl(
                                        \gamma_{gq}^{(1)}
                                       +\hat{\gamma}_{gq}^{(1)}
                                \Bigr)
                   +\frac{\beta_{0,Q}\zeta_2}{8}
                                \Bigl(
                                        \hat{\gamma}_{gg}^{(1)}
                                      -2\gamma_{gg}^{(1)}
                                      -2\beta_1
                                      -2\beta_{1,Q}
                                \Bigr)
                           +\frac{\delta m_1^{(-1)}}{4}
                                \Bigl(
                                     8 a_{gg,Q}^{(2)}
\nonumber\\ &&
                                    +24 \delta m_1^{(0)} \beta_{0,Q}
                                    +8 \delta m_1^{(1)} \beta_{0,Q} 
                                    +\zeta_2 \beta_{0,Q} \beta_0
                                    +9 \zeta_2 \beta_{0,Q}^2
                                \Bigr)
                           +\delta m_1^{(0)}
                                \Bigl(
                                     \beta_{0,Q} \delta m_1^{(0)}
                                    +\hat{\gamma}_{gg}^{(1)}
                                \Bigr)
\nonumber\\ &&
                           +\delta m_1^{(1)}
                                \Bigl(
                                     \hat{\gamma}_{qg}^{(0)} \gamma_{gq}^{(0)}
                                    +2 \beta_{0,Q} \gamma_{gg}^{(0)}
                                    +4 \beta_{0,Q} \beta_0
                                    +8 \beta_{0,Q}^2
                                \Bigr)
                           -2 \delta m_2^{(0)} \beta_{0,Q}
                 +a_{gg,Q}^{(3)}~. \label{Agg3QMSren}
   \end{eqnarray}
Here the heavy quark mass is renormalized in the on--shell scheme (OMS) and the strong coupling constant in 
the $\overline{\sf MS}$ scheme. The transformation from the OMS to the $\overline{\sf MS}$ scheme for the 
quark mass is described e.g. in Eq.~(5.13) of Ref.~\cite{Klein:2009ig}, see also 
Refs.~\cite{MASS}.

The Feynman diagrams are generated by using {\tt QGRAF} \cite{Nogueira:1991ex,Bierenbaum:2009mv}, the 
spinor and 
Lorentz--algebra 
is performed using {\tt FORM} \cite{FORM}, and the color algebra by using {\tt Color} 
\cite{vanRitbergen:1998pn}. The 
operator insertions are resummed into linear propagators as has been 
discussed in Ref.~\cite{Ablinger:2014yaa}, e.g. by 
\begin{eqnarray}
\sum_{k = 0}^\infty (\Delta. p)^k t^k = \frac{1}{1 - t \Delta.p},
\label{eq:rest}
\end{eqnarray}
where $t$ denotes an auxiliary parameter. Similar relations hold for the more complicated operator insertions. 
Taking the $k$th moment of (\ref{eq:rest}),
i.e. the coefficient of $t^k$, one obtains the contribution due to $(\Delta. p)^k$. Here $\Delta$ denotes 
a light--like vector.

In total, 642 irreducible Feynman diagrams contribute to the OMEs. The reduction to 
master 
integrals using the integration-by-parts relations
\cite{IBP} has been performed using {\tt Reduze~2} \cite{Studerus:2009ye,vonManteuffel:2012np}.

Unlike the case in later computations starting in 2017, we did not use the method of arbitrary high moments
\cite{Blumlein:2017dxp}, establishing the corresponding difference equations by the method of guessing \cite{GUESS} 
and solving them using the package {\tt Sigma} \cite{SIG1,SIG2}. Instead we calculated the master integrals directly in 
Mellin $N$--space using the different 
methods mentioned in Section~\ref{sec:1}, from which the individual Feynman diagrams were calculated. 
Their first few 
moments were compared to Mellin moments computed using {\tt Matad} 
\cite{Steinhauser:2000ry}, 
cf.~Ref.~\cite{Bierenbaum:2009mv} in the unpolarized case. In the case of two Feynman 
diagrams,
which are related to each other by reversal of the internal fermion line and which give
the same result, we had
to use a different computation method. Here we applied the method of differential equations 
\cite{Ablinger:2018zwz} to 
the master integrals working in the auxiliary variable $t$ used for the operator 
resummation. 
The transition to 
momentum fraction $x$--space by an analytic continuation is described in detail in Ref.~\cite{BS22}. This method avoids
going to Mellin $N$--space. The result can be Mellin transformed analytically at the end of the calculation, 
cf.~Section~\ref{sec:4}. Working in Mellin $N$--space implies 
that one decides for expressions at either even or odd integer values of $N$ starting with a value $N_0$ implied by 
the crossing relations \cite{Politzer:1974fr,Blumlein:1996vs}. If one would like to extract the small $x$ behavior 
of the OME, one has to do this in $x$--space, since the corresponding poles are situated at $N=1$ in the unpolarized 
case and at $N=0$ in the polarized case, for which no $N$--space representation is available.

In Mellin $N$--space, the present OMEs exhibit nested finite binomial sums \cite{Ablinger:2014bra}. 
Transforming to $x$--space, 
these structures lead to ${\rm G}$--functions \cite{Ablinger:2014bra} containing also 
square root--valued 
letters. Along with these very many new constants ${\rm G}(\{a_i\},1)$ occur in intermediary steps. They can be 
rationalized by procedures contained in the package {\tt 
HarmonicSums}, leading 
to tabulated cyclotomic constants \cite{Ablinger:2011te}. The latter turn out to 
further reduce to multiple zeta values \cite{Blumlein:2009cf} in all cases. We show a series of examples in
Appendix~\ref{sec:B}. Also the iterated integrals ${\rm G}(\{a_i\},x)$ containing root--valued letters can be 
rationalized, leading to cyclotomic harmonic polylogarithms \cite{Ablinger:2011te}. The latter 
representation can in principle be used in the numerical representation, cf.~Section~\ref{sec:5}.

Let us now discuss the new structures appearing in the present OMEs, which are nested finite binomial sums in Mellin 
$N$--space.
Their Mellin inversion can be expressed by iterative ${\rm G}$--functions and 
usual harmonic polylogarithms over the alphabet $\mathfrak{A}$
\begin{eqnarray}
\mathfrak{A} = \left. \left\{f_{k}(x)\right\}\right|_{k=1..6} =
\Biggl\{\frac{1}{x}, \frac{1}{1-x}, \frac{1}{1+x}, 
\frac{\sqrt{1-x}}{x}, \sqrt{x(1-x)}, \frac{1}{\sqrt{1-x}}\Biggr\}.
\label{eq:alph}
\end{eqnarray}
The ${\rm G}$--functions are defined by \cite{Ablinger:2014bra}
\begin{eqnarray}
{\rm G}(\{b, \vec{a}\},x) = \int_0^x dy f_b(y) {\rm G}(\{\vec{a}\},y),~~~{\rm G}(\{\emptyset\},y) = 1,~~~ f_b(x), 
f_{a_i}(x) \in 
\mathfrak{A}.
\end{eqnarray}
In $N$--space we consider the objects
\begin{eqnarray} 
{\sf BS}_0(N) &=& \frac{1}{2N - (2l+1)},~~~l \in \mathbb{N},
\\
{\sf BS}_1(N) &=& 4^N \frac{(N!)^2}{(2N)!},
\\
{\sf BS}_2(N) &=& \frac{1}{4^N} \frac{(2N)!}{(N!)^2},
\\
{\sf BS}_3(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{-\tau_1} \big(
        2 \tau_1\big)!}{\big(
        \tau_1!\big)^2 \tau_1},
\\
{\sf BS}_4(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{\tau_1} \big(
        \tau_1!\big)^2}{\big(
        2 \tau_1\big)! \tau_1^2},
\\
{\sf BS}_5(N) &=& 
\sum_{\tau_1=1}^N \frac{4^{\tau_1} \big(
        \tau_1!\big)^2}{\big(
        2 \tau_1\big)! \tau_1^3},
\\
{\sf BS}_6(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{-\tau_1} \big(
        2 \tau_1\big)! 
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{\tau_2} \big(
        \tau_2!\big)^2}{\big(
        2 \tau_2\big)! \tau_2^2}}{\big(
        \tau_1!\big)^2 \tau_1},
\\
{\sf BS}_7(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{-\tau_1} \big(
        2 \tau_1\big)! 
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{\tau_2} \big(
        \tau_2!\big)^2}{\big(
        2 \tau_2\big)! \tau_2^3}}{\big(
        \tau_1!\big)^2 \tau_1},
\\
{\sf BS}_8(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{\tau_2} \big(
        \tau_2!\big)^2}{\big(
        2 \tau_2\big)! \tau_2^2}}{\tau_1},
\\
{\sf BS}_9(N) &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{-\tau_1} \big(
        2 \tau_1\big)! 
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{\tau_2} \big(
        \tau_2!\big)^2 
\sum_{\tau_3=1}^{\tau_2} \displaystyle \frac{1}{\tau_3}}{\big(
        2 \tau_2\big)! \tau_2^2}}{\big(
        \tau_1!\big)^2 \tau_1},
\\
{\sf BS}_{10}(N) &=& \sum_{\tau_1=1}^N \frac{4^{\tau_1}}{\displaystyle \binom{2 \tau_1}{\tau_1}} 
\frac{1}{\tau_1^2} S_1(\tau_1),
\end{eqnarray}
where $S_{\vec{a}} \equiv S_{\vec{a}}(N)$ denote the nested harmonic 
sums \cite{Vermaseren:1998uu,Blumlein:1998if}.

The above finite binomial sums obey the following recursion relations
\begin{eqnarray}
{\sf BS}_3(N) - {\sf BS}_3(N-1)  &=& \frac{1}{N} {\sf BS}_2(N),
\\
{\sf BS}_4(N) - {\sf BS}_4(N-1)  &=& \frac{1}{N^2} {\sf BS}_1(N),
\\
{\sf BS}_5(N) - {\sf BS}_5(N-1)  &=& \frac{1}{N^3} {\sf BS}_1(N),
\\
{\sf BS}_6(N) - {\sf BS}_6(N-1)  &=& \frac{1}{N} {\sf BS}_2(N) {\sf BS}_4(N),
\\
{\sf BS}_7(N) - {\sf BS}_7(N-1)  &=& \frac{1}{N} {\sf BS}_2(N) {\sf BS}_5(N),
\\
{\sf BS}_8(N) - {\sf BS}_8(N-1)  &=& \frac{1}{N} {\sf BS}_4(N),
\\
{\sf BS}_9(N) - {\sf BS}_9(N-1)  &=& \frac{1}{N} {\sf BS}_2(N) {\sf BS}_{10}(N),
\\
{\sf BS}_{10}(N) - {\sf BS}_{10}(N-1)  &=& \frac{1}{N^2} {\sf BS}_1(N) S_1.
\end{eqnarray}
In Appendix~\ref{sec:D} we will also calculate representations of $a_{gg,Q}^{(3)}(N)$ and $\Delta 
a_{gg,Q}^{(3)}(N)$ in the asymptotic
region $|N| \gg 1$  up to $O(1/N^{10})$. In the analyticity region this and 
the recurrences allow one to compute $(\Delta)a_{gg,Q}^{(3)}(N)$
for $N \in \mathbb{C}$, see Ref.~\cite{Blumlein:2009ta}.
The asymptotic expansion of $(\Delta) a_{gg,Q}^{(3)}$ can be obtained from the asymptotic expansions of its building 
blocks. The asymptotic representations of the finite binomial sums ${\sf BS}_k(N),~k = 
0 \dots 10$ are also given in Appendix~\ref{sec:D}.
The constants contributing to the above sums can be calculated using infinite binomial and inverse binomial sums
\cite{Davydychev:2003mv,Weinzierl:2004bn,Ablinger:2014bra}.
In the calculation, first partly different binomial sums appear after the sum reduction performed by {\tt 
Sigma}, which have the following relations
\begin{eqnarray}
{\sf \overline{BS}}_4 &=& \sum_{\tau_1=1}^N \frac{\displaystyle 4^{-\tau_1} \big(
        2 \tau_1\big)!}{\big(
        \tau_1!\big)^2 \big(
        1+\tau_1\big)^2}
=
3
-\frac{(1+2 N) (3+2 N)}{(1+N)^2}
\frac{1}{4^{N}} \binom{2 N}{N}
-2 {\sf BS}_3,
\\
{\sf \overline{BS}}_6 &=& 
\sum_{\tau_1=1}^N \frac{\displaystyle 4^{\tau_1} \big(
        \tau_1!\big)^2 
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{-\tau_2} \big(
        2 \tau_2\big)!}{\big(
        \tau_2!\big)^2 \big(
        1+\tau_2\big)^2}}{\big(
        2 \tau_1\big)! \tau_1^2}
= -2 S_2 - 2 S_3 
-\frac{N (4+3 N)}{(1+N)^2}
\nonumber\\ &&
+3 {\sf BS}_4 - 2 {\sf BS}_3 {\sf BS}_4 + 2 {\sf BS}_6,
\\
{\sf \overline{BS}}_9 &=& \sum_{\tau_1=1}^N \frac{\displaystyle
\sum_{\tau_2=1}^{\tau_1} \frac{\displaystyle 4^{\tau_2} \big(
        \tau_2!\big)^2 
\sum_{\tau_3=1}^{\tau_2} \frac{\displaystyle 4^{-\tau_3} \big(
        2 \tau_3\big)!}{\big(
        \tau_3!\big)^2 \big(
        1+\tau_3\big)^2}}{\big(
        2 \tau_2\big)! \tau_2^2}}{\tau_1} =
 2 {\sf BS}_7
+3 {\sf BS}_8
-2 {\sf BS}_3 {\sf BS}_8
-2 {\sf BS}_9
\nonumber\\ && 
+\frac{N (5+4 N)}{(1+N)^2}
+2 {\sf BS}_6 S_1
-\big(
        3
        +2 S_2
        +2 S_3
\big) S_1
-S_2
-2 S_3
-2 S_4
+2 S_{2,1}
\nonumber\\ &&
+2 S_{3,1},
\\
{\sf \overline{BS}}_{11} &=& 
\sum_{i_1=3}^N \frac{2^{2 i_1} \big(
        \big(
                -2+i_1\big)!\big)^2}{\big(
        2 i_1\big)!} = -\frac{44}{3}
+3 {\sf BS}_4
+\frac{(-1+3 N)}{N^2} 2^{1+2 N} \binom{2 N}{N}^{-1},
\\
{\sf \overline{BS}}_{12} &=& 
\sum_{i_1=3}^N \frac{2^{2 i_1} \big(
        \big(
                -2+i_1\big)!\big)^2}{\big(
        2 i_1\big)! i_1} = -\frac{67}{3}
+ 4 {\sf BS}_4
+  {\sf BS}_5
+\frac{(-1+4 N)}{N^2} 2^{1+2 N} \binom{2 N}{N}^{-1},
\\
{\sf \overline{BS}}_{13} &=& 
\sum_{i_1=3}^N \frac{\displaystyle 2^{2 i_1} \big(
        \big(
                -2+i_1\big)!\big)^2 S_1\big({i_1}\big)}{\big(
        2 i_1\big)!}
=
-23
+3 {\sf BS}_5
-3 {\sf BS}_8
+ {\sf BS}_4 \big(
        -2
        +3 S_1
\big)
\nonumber\\ &&
+\frac{2(-1+4 N)}{N^2} 4^N \binom{2 N}{N}^{-1}
+\frac{2 (-1+3 N)}{N^2} 4^N  \binom{2 N}{N}^{-1} S_1.
\end{eqnarray}

The above functions emerge in products in part and one has to calculate the corresponding 
Mellin convolutions, also with harmonic sums. Besides the above binomial sums, the following  34 
harmonic sums up to weight {\sf w = 5} contribute
\begin{eqnarray}
&&
\{
S_1,
S_{-1},
S_2,
S_{-2},
S_3,
S_{-3},
S_4,
S_{-4},
S_{-5},
S_5,
S_{2,1},
S_{-2,1},
S_{2,-1},
S_{-2,-1},
S_{-2,2},
S_{3,1},
S_{-3,1},
S_{4,1},
\nonumber\\ &&
S_{-4,1},
S_{2,3},
S_{2,-3},
S_{-2,3},
S_{-2,-3},
S_{3,1,1},
S_{-3,1,1},
S_{2,1,1},
S_{-2,1,1},
S_{2,2,1},
S_{2,1,-2},
S_{-2,1,-2},
S_{-2,2,1},
\nonumber\\ &&
S_{-2,1,1,1},
S_{2,1,1,1}\big\}
\end{eqnarray}
after algebraic reduction \cite{Blumlein:2003gb}. Here we suppressed
the argument $N$ of the harmonic sums for brevity.

All logarithmic terms to the OMEs and the contributions to the constant term implied by renormalization 
were given in Refs.~\cite{Behring:2014eya,Blumlein:2021xlc} before. 
\section{\boldmath Results for $a_{gg,Q}^{(3)}(N)$ and $\Delta 
a_{gg,Q}^{(3)}(N)$}
\label{sec:3}

\vspace{1mm}
\noindent
We now present the constant contributions of the unrenormalized 
massive OMEs $(\Delta) 
\hat{\hspace*{-1mm}\hat{A}}_{gg,Q}^{(3)}$, denoted by 
$a_{gg,Q}^{(3)}$ and $\Delta a_{gg,Q}^{(3)}$. The expression for 
$a_{gg,Q}^{(3)}(N)$, valid for even moments $N \in 
\mathbb{N}, N \geq 2$, reads 
\input{aggQ3.tex}

\noindent
with
\begin{eqnarray}
{\sf B_4} = - 4 \zeta_2 \ln^2(2) + \frac{2}{3} \ln^4(2) - \frac{13}{2} \zeta_4 + 16 
{\rm Li}_4\left(\frac{1}{2}\right).
\end{eqnarray}
The polynomials $P_i$ are listed in Appendix~\ref{sec:A}.
Eq.~(\ref{aggQUNP}) possesses a removable pole at $N=2$. By a series 
expansion one obtains 
Eq.~(8.67) of Ref.~\cite{Bierenbaum:2009mv}. Furthermore, also the even moments for $N=4$ to $N=10$ agree with 
the result in Ref.~\cite{Bierenbaum:2009mv}. There are removable poles also at $N = 3, N = 1/2$ 
and $N = 
3/2$, 
see also Ref.~\cite{Ablinger:2014uka}. To see their cancellation, one has to expand the Mellin inversion to 
$x$--space around $x=0$. In the present case one finds the most singular terms $\propto 
\ln(x)/x$ and 
$\propto 1/x$, i.e. it is 
proven that rightmost pole is situated at $N=1$, as expected in the unpolarized gluonic case. Note that in 
general one may not expect to cancel the above poles by just using the $N$--space 
representation.

It is also instructive to see how accurate the asymptotic expansion of $a_{gg,Q}^{(3)}(N)$ 
represents higher moments. We compare the expressions for the $N_F$--independent part only, 
since the $N_F$--dependent 
part solely consists of harmonic sums. One may represent $a_{gg,Q}^{(3)}(N)$ by
\begin{eqnarray}
a_{gg,Q}^{(3)}(N) &=& a_{gg,Q,\delta}^{(3)} +  a_{gg,Q,\rm pl}^{(3)}(N) +
\tilde{a}_{gg,Q}^{(3)}(N).
\label{eq:aggQtil}
\end{eqnarray}
Here $a_{gg,Q,\delta}^{(3)}$ denotes the $N$--independent part,~cf.~(\ref{eq:4.6}),
$a_{gg,Q,\rm pl}^{(3)}(N)$ the part $\propto L$ with
\begin{eqnarray}
L = \ln(N) + \gamma_E, 
\label{eq:LG}
\end{eqnarray}
where $\gamma_E$ denotes the Euler--Mascheroni constant. 
The explicit 
expressions for the asymptotic expansion of 
$(\Delta) \tilde{a}_{gg,Q}^{(3)}(N)$ for the $N_F$--independent part are given in Appendix~\ref{sec:D}.
The asymptotic representation for $\tilde{a}_{gg,Q}^{(3)}(N)$ for positive even integer values converges 
very quickly, as shown in Table~1. Already for $N = 4$ a reasonable approximation is obtained by 
expanding to $O(1/N^{10})$.
\begin{table}[H]
\centering
\begin{tabular}{|r||r|r|}
\hline \hline
$N$  & complete expression & rel. asymp. accuracy \\
\hline \hline
4     & --430.337594532836914 & $ -2.10169$D--3  \\
6     & --261.554324759832203 & $ -3.13097$D--5   \\
8     & --193.698203673549029 & $ -1.29985$D--6  \\
10    & --156.572494406521072 & $ -9.09255$D--8  \\
12    & --132.845604857074259 & $ -7.66447$D--9  \\
22    & --79.9576964391278831 & $  3.81320$D--11   \\
42    & --48.0359673792636099 & $  1.41090$D--13  \\
102   & --24.2090141858051135 & $  3.14898$D--17  \\
\hline
\hline
\end{tabular}
\label{TAB1}
\caption[]{\sf 
Numerical comparison of $\tilde{a}_{gg,Q}^{(3)}(N)$ in QCD              with
its asymptotic
representation for $N_F = 0$ retaining 10 terms of the asymptotic expansion.}
\end{table}

We turn now to $\Delta a_{gg,Q}(N)$ which is given in the Larin scheme  \cite{Larin:1993tq} by
\input{POLaggQ3.tex}

\noindent
One may represent $\Delta a_{gg,Q}^{(3)}(N)$ in an analogous way to Eq.~(\ref{eq:aggQtil}). As in the 
unpolarized case, the expansion 
converges very quickly for the positive odd integer values, cf.~Table~2, with a reasonable description 
down to $N = 3$.
\begin{table}[H]
\centering
\begin{tabular}{|r||r|r|}
\hline \hline
$N$  & complete expression & rel. asymp. accuracy \\
\hline \hline
3     & --429.345090408771279 & $ -9.39608$D--4  \\
5     & --264.713704879676430 & $ -7.20184$D--6  \\
7     & --195.879611926179533 & $ -2.77043$D--7  \\
9     & --158.086523949063478 & $ -2.36920$D--8  \\
11    & --133.968747075663141 & $ -3.28594$D--9  \\
21    & --80.4139205966345783 & $ -5.45018$D--12  \\
41    & --48.2243776140486971 & $ -7.23295$D--15  \\
101   & --24.2616935237270326 & $ -1.03115$D--18  \\
\hline
\hline
\end{tabular}
\label{TAB2}
\caption[]{\sf Numerical comparison of $\Delta \tilde{a}_{gg,Q}^{(3)}(N)$ in QCD 
with 
its asymptotic 
representation for $N_F = 0$ retaining 10 terms of the asymptotic expansion.}
\end{table}
\section{The \boldmath $x$-space representation}
\label{sec:4}

\vspace{1mm}
\noindent
We perform an analytic inverse Mellin transform to $x$--space by using algorithms 
implemented in {\tt 
HarmonicSums}.
In momentum fraction space, the quantities $(\Delta) a_{gg,Q}(x)$ depend besides  
harmonic polylogarithms, ${\rm H}_{\vec{a}}(x)$, on ${\rm G}$--functions up to weight {\sf w = 5} 
at arguments $x$ and $1$ over the alphabet (\ref{eq:alph}), e.g. 
\begin{eqnarray}
{\rm G}\left(\left\{
\sqrt{(1-y) y}, 
\frac{1}{y},
\sqrt{(1-y) y}, 
\frac{1}{y},\frac{1}{1-y}\right\},x\right) = {\rm G}(\{5,1,5,1,2\},x),
\end{eqnarray}
and 17 similar functions both in the unpolarized and polarized case. The appearing 
constants can all be calculated analytically by using {\tt HarmonicSums} and only 
multiple zeta values remain at the end. In addition to this, the following 48 harmonic 
polylogarithms contribute
\begin{eqnarray}
&&\bigl\{
{\rm H}_0, {\rm H}_{-1}, 
{\rm H}_1, 
{\rm H}_{-1,1}, 
{\rm H}_{0,-1},{\rm H}_{0,1},{\rm H}_{0,-1,-1},
{\rm H}_{0,-1,1},{\rm H}_{0,0,-1},
{\rm H}_{0,0,1},{\rm H}_{0,1,-1},{\rm H}_{0,1,1},{\rm H}_{0,-1,-1,-1},
\nonumber\\ && 
{\rm H}_{0,-1,-1,1},{\rm H}_{0,-1,0,1},
{\rm H}_{0,-1,1,-1},{\rm H}_{0,-1,1,1},{\rm H}_{0,0,-1,-1},{\rm H}_{0,0,-1,1},{\rm H}_{0,0,0,-1},
{\rm H}_{0,0,0,1},{\rm H}_{0,0,1,-1},{\rm H}_{0,0,1,1},
\nonumber\\ && 
{\rm H}_{0,1,-1,-1},{\rm H}_{0,1,-1,1},{\rm H}_{0,1,1,-1},
{\rm H}_{0,1,1,1},{\rm H}_{0,-1,-1,0,1},{\rm H}_{0,-1,0,-1,-1},{\rm H}_{0,-1,0,1,1},
{\rm H}_{0,-1,1,0,1},
\nonumber\\ && 
{\rm H}_{0,0,-1,-1,-1},
{\rm H}_{0,0,-1,0,-1},{\rm H}_{0,0,-1,0,1},{\rm H}_{0,0,-1,1,1},{\rm H}_{0,0,0,-1,-1},{\rm H}_{0,0,0,-1,1},
{\rm H}_{0,0,0,0,-1},{\rm H}_{0,0,0,0,1},
\nonumber\\ && 
{\rm H}_{0,0,0,1,-1},{\rm H}_{0,0,0,1,1},
{\rm H}_{0,0,1,0,-1},
{\rm H}_{0,0,1,0,1},{\rm H}_{0,0,1,1,-1},
{\rm H}_{0,0,1,1,1},{\rm H}_{0,1,0,1,-1},
{\rm H}_{0,1,0,1,1},{\rm H}_{0,1,1,1,1}\Bigr\},
\end{eqnarray}
both after algebraic reduction for the ${\rm G}$- and ${\rm H}$-functions \cite{Blumlein:2003gb}, where we 
suppressed
the argument $x$ of the harmonic polylogarithms.

Furthermore, denominator structures of the kind
\begin{eqnarray}
\frac{1}{(1 \pm x)^k},~~~k = 2,3 
\end{eqnarray}
appear, which are also known from other massive calculations 
\cite{Berends:1987ab,Blumlein:2020jrf}. Since the 
corresponding expressions are very lengthy, we present them only in an ancillary file in 
computer--readable form. Here we discuss their principal structure. The expressions for $(\Delta) 
a_{gg,Q}^{(3)}(x)$ have the following form
\begin{eqnarray}
(\Delta) a_{gg,Q}^{(3)}(x) = (\Delta) a_{gg,Q,\delta}^{(3)} \delta(1-x) + [(\Delta) a_{gg,Q,\rm 
plus}^{(3)}(x)]_+
+ (\Delta) a_{gg,Q,\rm reg}^{(3)}(x)],
\end{eqnarray}
where 
\begin{eqnarray}
\int_0^1 dx g(x) [f(x)]_+ := \int_0^1 [g(x) - g(1)] f(x).
\end{eqnarray}
For the $\delta(1-x)$ and $+$-function contributions we obtain 
\begin{eqnarray}
\label{eq:4.6}
(\Delta) a_{gg,Q,\delta}^{(3)} &=& 
\textcolor{blue}{T_F} \Biggl\{
        \textcolor{blue}{C_F} \Biggl[
                \textcolor{blue}{C_A} \Biggl(
                        \frac{16541}{162}
                        -\frac{64 {\sf B}_4}{3}
                        +\frac{128 \zeta_4}{3}
                        +52 \zeta_2
                        -\frac{2617 \zeta_3}{12}
                \Biggr)
                +\textcolor{blue}{T_F} \Biggl(
                        -\frac{1478}{81}
\nonumber\\ &&                   
      +\textcolor{blue}{N_F} \Biggl(
                                -\frac{1942}{81}
                                -\frac{20 \zeta_2}{3}
                        \Biggr)
                        -\frac{88 \zeta_2}{3}
                        -7 \zeta_3
                \Biggr)
        \Biggr]
        +\textcolor{blue}{C_A^2} \Biggl[
                \frac{34315}{324}
                +\frac{32 {\sf B}_4}{3}
                -\frac{3778 \zeta_4}{27}
\nonumber\\ &&              
   +\frac{992}{27} \zeta_2
                +\Biggl(
                        \frac{20435}{216}
                        +24 \zeta_2
                \Biggr) \zeta_3
                -\frac{304}{9} \zeta_5
        \Biggr]
        +\textcolor{blue}{C_A T_F} \Biggl[
                \frac{2587}{135}
                +\textcolor{blue}{N_F} \Biggl(
                        -\frac{178}{9}
                        +\frac{196 \zeta_2}{27}
                \Biggr)
\nonumber\\ &&            
     +\frac{572 \zeta_2}{27}
                -\frac{291 \zeta_3}{10}
        \Biggr]
        +\textcolor{blue}{C_F^2} \Biggl[
                \frac{274}{9}
                +\frac{95 \zeta_3}{3}
        \Biggr]
        +\frac{64}{27} \textcolor{blue}{T_F^2} \zeta_3
\Biggr\},
\\
(\Delta) a_{gg,Q,\rm plus}^{(3)} &=& 
\frac{\textcolor{blue}{T_F}}{1-x} \Biggl\{
        \textcolor{blue}{C_A T_F}  \Biggl[
                \frac{35168}{729}
                +\textcolor{blue}{N_F} \Biggl(
                        \frac{55552}{729}
                        +\frac{160 \zeta_2}{27}
                        -\frac{448 \zeta_3}{27}
                \Biggr)
                +\frac{560}{27} \zeta_2
                +\frac{1120}{27} \zeta_3
        \Biggr]
\nonumber\\ && 
        +\textcolor{blue}{C_A^2} \Biggl[
                -\frac{32564}{729}
                -\frac{32 {\sf B}_4}{3}
                +104 \zeta_4
                -\frac{3248 \zeta_2}{81}
                -\frac{1796 \zeta_3}{27}
        \Biggr]
        +\textcolor{blue}{C_A C_F}  \Biggl[
                -\frac{6152}{27}
                +\frac{64 {\sf B}_4}{3}
\nonumber\\ &&               
  -96 \zeta_4
                -40 \zeta_2
                +\frac{1208 \zeta_3}{9}
        \Biggr]
\Biggr\}.
\end{eqnarray}
They are the same in both cases.
The regular parts $(\Delta) a_{gg,Q,\rm reg}^{(3)}(x)$ are given in an ancillary file 
to the present paper.

We further expand $(\Delta) a_{gg,Q}(x)$  in the small $x$ and large $x$ regions, where
much simpler structures ruled by logarithms are obtained. 
The principal $x$--space structure to higher powers in $x$ is
\begin{eqnarray}
(\Delta) a_{gg,Q}^{x \rightarrow 0}(x) &=& 
c_1 \frac{\ln(x)}{x} + c_2 \frac{1}{x} +
\sum_{k=0}^\infty \Biggl[c_{3,k} + c_{4,k} \ln(x) 
+ c_{5,k} \ln^6(x) + c_{7,k} \ln^3(x)
+ c_{8,k} \ln^4(x) 
\nonumber\\ &&
+ c_{9,k} \ln^5(x)\Biggr] x^k
\\
(\Delta) a_{gg,Q}^{x \rightarrow 1}(x) &=& d_1 \delta(1-x) + \frac{d_2}{1-x}  + 
\sum_{k=0}^\infty
\Biggl[d_{3,k} + d_{4,k} \ln(1-x) + d_{5,k} \ln^2(1-x) 
\nonumber\\ &&
+ d_{6,k} \ln^3(1-x)+ d_{7,k} \ln^4(1-x)\Biggr] x^k,
\end{eqnarray}
where some of the coefficients can be zero.
The leading behavior for the two expansions around $x=0$ and $x=1$ in the unpolarized case
are given by
\begin{eqnarray}
\lefteqn{ a_{gg,Q}^{x \rightarrow 0}(x) \propto} \nonumber\\ && 
\frac{1}{x}
\Biggl\{
\ln(x) 
\Biggl[
        \textcolor{blue}{C_A^2 T_F}
        \biggl(
                -\frac{11488}{81}
                +\frac{224 \zeta_2}{27}
                +\frac{256 \zeta_3}{3}
        \biggr) 
        + \textcolor{blue}{C_A C_F T_F}
        \biggl(
                -\frac{15040}{243}
                -\frac{1408 \zeta_2}{27}
\nonumber\\ &&      
           -\frac{1088 \zeta_3}{9}
        \biggr)
\Biggr]
+\textcolor{blue}{C_A T_F^2} 
\Biggl[
        \frac{112016}{729}
        +\frac{1288}{27} \zeta_2
        +\frac{1120}{27} \zeta_3
        +\biggl(
                \frac{108256}{729}
                +\frac{368 \zeta_2}{27}
                -\frac{448 \zeta_3}{27}
        \biggr) 
\nonumber\\ && \times
\textcolor{blue}{N_F}
\Biggr]
+ \textcolor{blue}{C_F} 
\Biggl[
        \textcolor{blue}{T_F^2} 
        \Biggl(
                -\frac{107488}{729}
                -\frac{656}{27} \zeta_2
                +\frac{3904}{27} \zeta_3
                +\biggl(
                        \frac{116800}{729}
                        +\frac{224 \zeta_2}{27}
                        -\frac{1792 \zeta_3}{27}
                \biggr) \textcolor{blue}{N_F}
        \Biggr)
\nonumber\\ && 
        + \textcolor{blue}{C_A T_F}
        \Biggl(
                -\frac{5538448}{3645}
                +\frac{1664 {\sf B}_4}{3}
                -\frac{43024 \zeta_4}{9}
                +\frac{12208}{27} \zeta_2
                +\frac{211504}{45} \zeta_3
        \Biggr) 
\Biggr]
\nonumber\\ && 
+ \textcolor{blue}{C_A^2 T_F}
\biggl(
        -\frac{4849484}{3645}
        -\frac{352 {\sf B}_4}{3}
        +\frac{11056 \zeta_4}{9}
        -\frac{1088}{81} \zeta_2
        -\frac{84764}{135} \zeta_3
\biggr) 
\nonumber\\ && 
+ \textcolor{blue}{C_F^2 T_F}
\biggl(
        \frac{10048}{5}
        -640 {\sf B}_4
        +\frac{51104 \zeta_4}{9}
        -\frac{10096}{9} \zeta_2
        -\frac{280016}{45} \zeta_3
\biggr) 
\Biggr\}
\nonumber\\ &&
+
\Biggl[
        -\frac{4}{3}  \textcolor{blue}{C_F C_A T_F} 
        +\frac{2}{15} \textcolor{blue}{C_F^2 T_F}
\Biggr] \ln^5(x)
+ \Biggl[
-\frac{40}{27} \textcolor{blue}{C_A^2 T_F} 
+\frac{4}{9} \textcolor{blue}{C_F^2 T_F}
+\textcolor{blue}{C_F} \Biggl(
        -\frac{296}{27} \textcolor{blue}{C_A T_F} 
\nonumber\\ && 
        +\Biggl(
                \frac{28}{27}
+\frac{56}{27} \textcolor{blue}{N_F}\Biggr) \textcolor{blue}{T_F^2}
\Biggr)
\Biggr] \ln^4(x)
+ \Biggl[
\frac{112}{81} \textcolor{blue}{C_A} (1+2 \textcolor{blue}{N_F}) \textcolor{blue}{T_F^2}
+\textcolor{blue}{C_F} \Biggl(
        \Biggl(
                \frac{1016}{81}+\frac{496}{81} \textcolor{blue}{N_F}\Biggr) \textcolor{blue}{T_F^2}
\nonumber\\ &&            
        +\textcolor{blue}{C_A T_F} \Biggl(
     -\frac{10372}{81}
                -\frac{328 \zeta_2}{9}
        \Biggr)
\Biggr)
+\textcolor{blue}{C_F^2 T_F} \Biggl[
        -\frac{2}{3}
        +\frac{4 \zeta_2}{9}
\Biggr]
+\textcolor{blue}{C_A^2 T_F} \Biggl[
        -\frac{1672}{81}
        +8 \zeta_2
\Biggr]
\Biggr] \ln^3(x)
\nonumber\\ && 
+ \Biggl[
\frac{8}{81} \textcolor{blue}{C_A} (155+118 \textcolor{blue}{N_F}) \textcolor{blue}{T_F^2}
+\textcolor{blue}{C_F} \Biggl[
        \textcolor{blue}{T_F^2} \Biggl(
                -\frac{32}{81}
                +\textcolor{blue}{N_F} \Biggl(
                        \frac{3872}{81}
                        -\frac{16 \zeta_2}{9}
                \Biggr)
                +\frac{232 \zeta_2}{9}
        \Biggr)
\nonumber\\ && 
        +\textcolor{blue}{C_A T_F}  \Biggl(
                -\frac{70304}{81}
                -\frac{680 \zeta_2}{9}
                +\frac{80 \zeta_3}{3}
        \Biggr)
\Biggr]
+\textcolor{blue}{C_A^2 T_F} \Biggl[
        \frac{4684}{81}
        +\frac{20 \zeta_2}{3}
\Biggr]
+\textcolor{blue}{C_F^2 T_F} \Biggl[
        56
\nonumber\\ &&         
+\frac{8 \zeta_2}{3}
        -40 \zeta_3
\Biggr]
\Biggr] \ln^2(x)
+ \Biggl[
\textcolor{blue}{C_F} \Biggl[
        \textcolor{blue}{T_F^2} \Biggl(
                \frac{140992}{243}
                +\textcolor{blue}{N_F} \Biggl(
                        \frac{182528}{243}
                        -\frac{400 \zeta_2}{27}
                        -\frac{640 \zeta_3}{9}
                \Biggr)
\nonumber\\ &&                
 -\frac{728}{27} \zeta_2
                -\frac{224}{9} \zeta_3
        \Biggr)
        +\textcolor{blue}{C_A T_F} \Biggl(
                -\frac{514952}{243}
                +\frac{152 \zeta_4}{3}
                -\frac{21140 \zeta_2}{27}
                -\frac{2576 \zeta_3}{9}
        \Biggr)
\Biggr]
\nonumber\\ && 
+\textcolor{blue}{C_A T_F^2}  \Biggl[
        \frac{184}{27}
        +\textcolor{blue}{N_F} \Biggl(
                \frac{656}{27}
                -\frac{32 \zeta_2}{27}
        \Biggr)
        +\frac{464 \zeta_2}{27}
\Biggr]
+\textcolor{blue}{C_A^2 T_F} \Biggl[
        -\frac{42476}{81}
        -92 \zeta_4
        +\frac{4504 \zeta_2}{27}
\nonumber\\ &&  
       +\frac{64 \zeta_3}{3}
\Biggr]
+\textcolor{blue}{C_F^2 T_F} \Biggl[
        -\frac{1036}{3}
        -\frac{976 \zeta_4}{3}
        -\frac{58 \zeta_2}{3}
        +\frac{416 \zeta_3}{3}
\Biggr]
\Biggr] \ln(x),
\end{eqnarray}
and 
\begin{eqnarray}
a_{gg,Q}^{(3),x \rightarrow 1}(x) &\propto& 
a_{gg,Q,\delta}^{(3)} \delta(1-x)
+ a_{gg,Q,\rm plus}^{(3)}(x)
+ \Biggl[
        -\frac{32}{27} \textcolor{blue}{C_A T_F^2} (17+12 \textcolor{blue}{N_F}) 
        +\textcolor{blue}{C_A C_F T_F} \Biggl(
                56
                -\frac{32 \zeta_2}{3}
        \Biggr)
\nonumber\\ && 
        +\textcolor{blue}{C_A^2 T_F} \Biggl(
                \frac{9238}{81}
                -\frac{104 \zeta_2}{9}
                +16 \zeta_3
        \Biggr)
\Biggr] \ln(1-x)
+\Biggl[
        -\frac{8}{27} \textcolor{blue}{C_A T_F^2} (7+8 \textcolor{blue}{N_F}) 
\nonumber\\ &&       
 +\textcolor{blue}{C_A^2 T_F} \Biggl(
                \frac{314}{27}
                -\frac{4 \zeta_2}{3}
        \Biggr)
\Biggr] \ln^2(1-x)
+\frac{32}{27} \textcolor{blue}{C_A^2 T_F} \ln^3(1-x).
\end{eqnarray}
In the polarized case one has
\begin{eqnarray}
\Delta a_{gg,Q}^{x \rightarrow 0}(x) &\propto& 
\Biggl[
        -\frac{4}{15} \textcolor{blue}{C_A^2} \textcolor{blue}{T_F}
        +\frac{12}{5} \textcolor{blue}{C_A} \textcolor{blue}{C_F} \textcolor{blue}{T_F}
        +\frac{2}{15} \textcolor{blue}{C_F^2} \textcolor{blue}{T_F}
\Biggr] \ln^5(x)
+\Biggl[
        \frac{34}{27} \textcolor{blue}{C_A^2} \textcolor{blue}{T_F}
        +\frac{14}{9} \textcolor{blue}{C_F^2} \textcolor{blue}{T_F}
        \nonumber \\ &&
        +\textcolor{blue}{C_F} \Biggl(
                \frac{760}{27} \textcolor{blue}{C_A} \textcolor{blue}{T_F}
                + \textcolor{blue}{T_F^2}
                \biggl(
                        \frac{28}{27}
                        +\frac{56 \textcolor{blue}{N_F}}{27}
                \biggr)
        \Biggr)
\Biggr] \ln^4(x)
+\Biggl[
        \frac{112}{81} \textcolor{blue}{C_A} (1+2 \textcolor{blue}{N_F}) \textcolor{blue}{T_F^2}
        \nonumber \\ &&
        +\textcolor{blue}{C_F} 
        \Biggl(
                \textcolor{blue}{T_F^2} 
                \biggl(
                        \frac{968}{81}
                        +\frac{1552}{81} \textcolor{blue}{N_F}
                \biggr) 
                +\textcolor{blue}{C_A} \textcolor{blue}{T_F} 
                \biggl(
                        \frac{13484}{81}
                        -\frac{184 \zeta_2}{9}
                \biggr)
        \Biggr)
        +\textcolor{blue}{C_F^2} \textcolor{blue}{T_F} 
        \Biggl(
                -\frac{70}{9}
                +\frac{4 \zeta_2}{9}
        \Biggr)
        \nonumber \\ &&
        +\textcolor{blue}{C_A^2} \textcolor{blue}{T_F} 
        \Biggl(
                \frac{2848}{81}
                +8 \zeta_2
        \Biggr)
\Biggr] \ln^3(x)
+\Biggl[
        \frac{16}{81} \textcolor{blue}{C_A} (85+146 \textcolor{blue}{N_F}) \textcolor{blue}{T_F^2}
        +\textcolor{blue}{C_F} 
        \Biggl(
                \textcolor{blue}{T_F^2} 
                \biggl(
                        \frac{2680}{81}
                        \nonumber \\ &&
                        +\textcolor{blue}{N_F} 
                        \biggl(
                                \frac{6704}{81}
                                -\frac{16 \zeta_2}{9}
                        \biggr)
                        +\frac{232 \zeta_2}{9}
                \biggr)
                +\textcolor{blue}{C_A} \textcolor{blue}{T_F} 
                \biggl(
                        \frac{44476}{81}
                        -\frac{1184 \zeta_2}{9}
                        -\frac{544 \zeta_3}{3}
                \biggr)
        \Biggr)
        \nonumber \\ &&
        +\textcolor{blue}{C_F^2} \textcolor{blue}{T_F} 
        \Biggl(
                -\frac{358}{3}
                -\frac{412 \zeta_2}{3}
                +\frac{8 \zeta_3}{3}
        \Biggr)
        +\textcolor{blue}{C_A^2} \textcolor{blue}{T_F} 
        \Biggl(
                \frac{2588}{27}
                +\frac{244 \zeta_2}{9}
                +112 \zeta_3
        \Biggr)
\Biggr] \ln^2(x)
\nonumber \\ &&
+ \Biggr[
        \textcolor{blue}{C_F} 
        \Biggl(
                \textcolor{blue}{T_F^2} 
                \biggl(
                        \frac{53920}{243}
                        +\textcolor{blue}{N_F} \biggl(
                                \frac{121280}{243}
                                -\frac{304 \zeta_2}{27}
                                -\frac{640 \zeta_3}{9}
                        \biggr)
                        +\frac{2776}{27} \zeta_2
                        -\frac{224}{9} \zeta_3
                \biggr)
                \nonumber \\ &&
                +\textcolor{blue}{C_A} \textcolor{blue}{T_F} 
                \biggl(
                        -\frac{638672}{243}
                        +\frac{4400 \zeta_4}{3}
                        -\frac{16484 \zeta_2}{27}
                        +\frac{5968 \zeta_3}{9}
                \biggr)
        \biggr)
        +\textcolor{blue}{C_A} \textcolor{blue}{T_F^2} 
        \Biggl(
                \frac{4880}{81}
                \nonumber \\ &&
                +\textcolor{blue}{N_F} 
                \biggl(
                        \frac{9760}{81}
                        -\frac{32 \zeta_2}{27}
                \biggr)
                +\frac{464 \zeta_2}{27}
        \Biggr)
        +\textcolor{blue}{C_F^2} \textcolor{blue}{T_F} 
        \Biggr(
                -\frac{2152}{3}
                -\frac{1744 \zeta_4}{3}
                -\frac{530 \zeta_2}{3}
                -\frac{3136 \zeta_3}{3}
        \Biggl)
        \nonumber \\ &&
        +\textcolor{blue}{C_A^2} \textcolor{blue}{T_F} 
        \Biggl(
                \frac{964}{9}
                -\frac{452 \zeta_4}{3}
                +\frac{2240 \zeta_2}{27}
                +\frac{1904 \zeta_3}{9}
        \Biggr)
\Biggr] \ln(x),
\end{eqnarray}
and 
\begin{eqnarray}
\Delta a_{gg,Q}^{x \rightarrow 1}(x) =
a_{gg,Q}^{x \rightarrow 1}(x)
\end{eqnarray}
for all terms up to $\propto \ln(1-x)$.
There are no predictions for these limits in the literature. As will be shown in Section~\ref{sec:5}, 
both the 
formally most leading small $x$
and large $x$ contributions are insufficient approximations to $(\Delta) a_{gg,Q}(x)$, as expected from other 
results
\cite{Ablinger:2014nga,SMX}. This also applies to the leading large $x$ result due to its limited reach.
\section{Numerical Results}
\label{sec:5}

\vspace{1mm}
\noindent
For the numerical representations one needs to calculate the ${\rm G}$--functions in a new way, while for 
the harmonic polylogarithms different numerical codes exist, 
cf.~e.g.~Refs.~\cite{Gehrmann:2001pz,Vollinga:2004sn,Ablinger:2018sat}. One way to 
proceed in the present case would be to 
rationalize the letters of the ${\rm G}$--functions by a formalism available in the package {\tt HarmonicSums} leading 
to the cyclotomic letters \cite{Ablinger:2011te}
\begin{eqnarray} 
\Biggl\{
\frac{1}{1+t^2}, \frac{t}{1+t^2} \Biggr\}. 
\end{eqnarray}
The corresponding functions can also be decomposed into letters belonging to generalized harmonic polylogarithms 
\cite{Moch:2001zr,Ablinger:2013cf},
\begin{eqnarray} 
\Biggl\{\frac{1}{t}, \frac{1}{1-t}, \frac{1}{1+t},
\frac{1}{1 + i t}, \frac{1}{1 - i t} \Biggr\}. 
\end{eqnarray}
Iterated integrals over the latter letters can be calculated numerically using the H\"older convolution 
\cite{Borwein:1999js}\footnote{This is a particular convolution in which the H\"older condition \cite{NASM} and  
H\"older mean \cite{HOELDER} are used.}
implemented in Ref.~\cite{Vollinga:2004sn}.

Since the above representation is still somewhat slow, we have alternatively expanded the final expressions for
$a_{gg,Q}^{(3)}(x)$ and
$\Delta a_{gg,Q}^{(3)}(x)$ in series around $x=0$ and $x=1$ analytically to 50 terms. Both 
representations are matched
in the middle of the $x$ range $]0,1]$ at an accuracy of $5 \cdot 10^{-15}$. The corresponding routines are found in an
attachment to the present paper. Since our expansions are fully analytic, they can be extended to an even higher 
accuracy if needed.

\begin{figure}[H]   
\centering
\includegraphics[width=0.7\textwidth]{aaggQ3.pdf}
\caption{\sf 
The non--$N_F$ terms of $a_{gg,Q}^{(3)}(N)$ (rescaled) as a function of $x$. 
Full line (black): complete result;
upper dotted line (red): term $\propto \ln(x)/x$;
lower dashed line (cyan): small $x$ terms $\propto 1/x$;
lower dotted line (blue): small $x$ terms including all $\ln(x)$ terms
up to the constant term; 
upper dashed line (green): large $x$ contribution up to the constant 
term;
dash-dotted line (brown): complete large $x$ contribution.}
\label{fig:1} 
\end{figure}
\begin{figure}[H]   
\centering
\includegraphics[width=0.7\textwidth]{aggQ3NF.pdf}
\caption{\sf 
The $N_F$ terms of $a_{gg,Q}^{(3)}(N)$ (rescaled) as a function of $x$. 
Full line (black): complete result; 
upper dotted line (red): term $\propto \ln(x)/x$;
upper dashed line (cyan): small $x$ terms $\propto 1/x$;
lower dotted line (blue): small $x$ terms including all $\ln(x)$ terms
up the constant term;
lower dashed line (green): large $x$ contribution
up to the constant term;
dash-dotted line (brown): full large $x$ contribution.}
\label{fig:2} 
\end{figure}

We split the discussion of $(\Delta) a_{gg,Q}^{(3)}(x)$ into the terms free of $N_F$ and the linear $N_F$ term, since
the parameter $N_F$ is arbitrary. In the unpolarized case
the corresponding results are given in Figures~\ref{fig:1} and \ref{fig:2}, where 
we have used the rescaling factor $x (1-x)^2$ for better visibility.

The formally leading small $x$ term $\propto \ln(x)/x$ of the $N_F$--independent part deviates from the 
complete result everywhere, cf.~Figure~\ref{fig:1}.\footnote{This has also been observed before for other 
OMEs in Refs.~\cite{Blumlein:1998mg,Ablinger:2014nga}.} Adding the $1/x$ allows to describe the region 
$x < 0.001$. Adding also the remaining terms of the expansion around $x=0$ up to the constant term one 
covers the region $x < 0.1$. The large $x$ singular and logarithmic contributions $\propto \ln^l(1-x)$ 
up to the constant term stop to agree with the complete result at $x \sim 0.7$. 
The matching between the small $x$ and large $x$ expansions, taking into account 50 expansion
terms in both cases, can be performed with a
relative accuracy of $O(10^{-14})$ and finally allows to describe the whole region $x \in ]0,1[$.

Both the $N_F$--dependent terms in the unpolarized and polarized cases are smaller than the corresponding 
$N_F$--independent parts. The unpolarized $N_F$--dependent part of $a_{gg,Q}^{(3)}(x)$ behaves like 
$1/x$ in the small $x$ limit, as shown in Figure~\ref{fig:2}. The expansions around $x=0$ and $x=1$ to 
$O(x^{50})$ can be matched at $x \sim 0.5$ at a relative accuracy of $O(10^{-14})$. Just keeping 
the divergent and $\ln^l(1-x)$ up to the constant term in the large $x$ region yields
agreement with the complete result for $x \raisebox{-0.07cm    0.7$. The leading small $x$ term, which in view of the 
complete expression for $a_{gg,Q}^{(3)}(x)$ is subleading, starts to deviate for values $x > 0.002$ 
from the complete expression. Accounting in addition for the next term $\propto \ln^4(x)$ the range is 
even only $x < 0.0002$. Including all small $x$ logarithmic terms up to the constant term cover the 
region $x < 0.02$.

In the polarized case, we rescale $\Delta a_{gg,Q}^{(3)}(N)$ by the factor $\sqrt{x} (1-x)^2$ for 
better visibility, which is different from the one in the unpolarized case, cf.~Figures~3,4.
\begin{figure}[H] 
\centering 
\includegraphics[width=0.7\textwidth]{DaggQ3.pdf} 
\caption{\sf The 
non--$N_F$ terms of $\Delta a_{gg,Q}^{(3)}(N)$ (rescaled) as a function of $x$. 
Full line (black): complete result; 
lower dotted line (red): term $\ln^5(x)$; 
upper dotted line (blue): small $x$ terms $\propto \ln^5(x)$ and $\ln^4(x)$;
upper dashed line (cyan): small $x$ terms including all $\ln(x)$ terms up to the constant term; 
lower dash-dotted line (green): large $x$ contribution up to the constant term; 
dash-dotted line (brown): full large $x$ contribution.} 
\label{fig:3} 
\end{figure}
In the $N_F$--independent case, the leading small $x$ 
result $\propto \ln^5(x)$ deviates from the complete result significantly, which is also the case adding one subleading 
term. Taking all small $x$ logarithms into account 
one covers the region $x < 0.002$. $x \sim 0.5$ is 
the ideal matching point of the 50--term  small $x$ and large $x$ expansions at a relative accuracy of $ \raisebox{-0.07cm    
10^{-14}$. The 
large $x$ expansion of the leading terms $\propto 1/(1-x)$ and $\ln^k(1-x)$ down to the constant term range to $0.7 < x$, 
while the 50--term 
expansion covers $0.02 < x$ and does then deviate from the complete result.
\begin{figure}[H]   
\centering
\includegraphics[width=0.7\textwidth]{DaggQ3NF.pdf}
\caption{\sf The $\Delta a_{gg,Q}^{(3)}(N)$  term  $\propto N_F$ (rescaled) as a 
function of $x$.
Full line (black): complete result;
upper dotted line (red): term $\propto \ln^4(x)$;
lower dotted line (blue): small $x$ terms $\propto \ln^4(x)$ and $\ln^3(x)$;
upper dashed line (cyan): small $x$ terms including all $\ln(x)$ terms up to the constant term;
lower dashed line (green): large $x$ contribution up to the constant term;
dash-dotted line (brown): full large $x$ contribution.}
\label{fig:4} 
\end{figure}


A similar behavior of the different curves in the $N_F$--dependent part in the polarized case to 
the $N_F$--independent part is observed in Figure~\ref{fig:4}. Again, the leading small $x$ term 
$\propto \ln^4(x)$ does not describe the complete expression anywhere, 
as well as taking into account the next logarithmic order,
while the complete set of the small 
$x$ logarithms $\ln^k(x)$ covers the range $x  < 0.008$. The 50--term large $x$ expansion describes the 
complete result for $x \raisebox{-0.07cm    0.3$, while the leading large $x$ singular contributions agree with it only 
for $x > 0.8$. 

For the numerical representation used we present the regular parts expanded in terms of power series 
to $O(x^{50})$ around $x=0$ and $x=1$ in a {\tt Mathematica} .m file attached to this paper. 
Since the corresponding 
expansions are performed analytically, it can be extended to work at an even higher precision, if 
needed.
\section{Conclusions}
\label{sec:6}

\vspace*{1mm}
\noindent
We have calculated the single heavy quark mass OMEs $A_{gg,Q}^{(3)}$ and $\Delta A_{gg,Q}^{(3)}$ in the 
unpolarized and polarized cases, and the previously known logarithmic contributions are now supplemented 
by the functions $a_{gg,Q}^{(3)}$ and $\Delta a_{gg,Q}^{(3)}$. In Mellin $N$--space, they contain besides 
nested harmonic sums also nested finite binomial sums. In momentum fraction $x$--space, theses 
quantities are represented by iterated integrals, the ${\rm G}$--functions, over an alphabet containing also 
square root--valued letters. We provided recursions and asymptotic expansions in 
Mellin $N$--space for 
these quantities, allowing also the analytic continuation of $a_{gg,Q}^{(3)}$ and 
$\Delta 
a_{gg,Q}^{(3)}$ from the even (odd) moments to $N \in \mathbb{C}$. In $x$--space it is convenient to
use highly accurate Taylor series expansions. Large expressions both in $N$-- and $x$--space are 
provided in computer--readable ancillary files to the present paper. We find that leading small and 
large $x$ expansions of $a_{gg,Q}^{(3)}$ and $\Delta  a_{gg,Q}^{(3)}$ are of limited use, since they 
cover rather small regions in $x$ only.

The present results complete the transition matrix elements in the single-- and double--mass
variable flavor number scheme for the gluon distributions in the unpolarized and polarized cases at 
three--loop order.



