\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.