††institutetext: Institut für Theoretische Teilchenphysik und Kosmologie, RWTH
Aachen University,
Sommerfeldstr. 16, 52056 Aachen, Germany

# Renormalization of the axial current operator in dimensional regularization
at four-loop in QCD

Long Chen and Michał Czakon<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We provide the renormalization constants of the axial current operators, both
singlet and non-singlet, in dimensional regularization up to four-loop order
in (QCD), determined using the off-shell Ward-Takahashi identity for an axial
current with a non-anticommuting $\gamma_{5}$. A possible application of the
result for the singlet axial current operator is the extraction of the non-
decoupling mass logarithms in the axial quark form factors.

††preprint: TTK-21-54, P3H-21-099

##  1 Introduction

Dimensional regularization (DR) tHooft:1972tcz ; Bollini:1972ui  is the
regularization framework that underlies most of the modern higher-order
perturbative calculations in Quantum Chromodynamics (QCD). It is well known
that special attention is needed in DR in the treatment of $\gamma_{5}$ – an
intrinsically 4-dimensional object. At the root of the issue is the
contradiction between a fully anticommuting (non-trivial) $\gamma_{5}$ in D
($\neq 4$) dimensions and the non-vanishing value of the trace of the product
of four $\gamma$ matrices and $\gamma_{5}$ in 4 dimensions. Nevertheless, the
anticommutativity of $\gamma_{5}$ is essential for the concept of chirality of
spinors in 4 dimensions. On the other hand, an anticommuting $\gamma_{5}$ in
DR, where the invariance of loop integrals under arbitrary loop-momentum
shifts is ensured, leads to the absence of the axial or Adler-Bell-Jackiw
(ABJ) anomaly Adler:1969gk ; Bell:1969ts . In order to overcome these issues,
various $\gamma_{5}$ prescriptions in DR have been developed in the literature
tHooft:1972tcz ; Akyeampong:1973xi ; Breitenlohner:1977hr ; Bardeen:1972vi ;
Chanowitz:1979zu ; Gottlieb:1979ix ; Ovrut:1981ne ; Espriu:1982bw ;
Buras:1989xd ; Kreimer:1989ke ; Korner:1991sx ; Larin:1991tj ; Larin:1993tq ;
Jegerlehner:2000dz ; Moch:2015usa ; Zerf:2019ynn  over a span of nearly 50
years, albeit each with its own pros and cons.

Given the fact that a local gauge theory with an internal gauge anomaly
suffers from severe theoretical problems, e.g. the loss of unitarity and all-
order (multiplicative) renormalizability, it is natural to wonder why one
would care about studying the anomalous singlet axial current at all. Putting
aside many related interesting theoretical developments, an external axial
anomaly (i.e. anomaly in the divergence of an external axial current which is
not coupled to the quantized gauge bosons in question) is not only allowed in
a gauge theory, but also crucial in understanding some important physical
observations, such as a low-energy theorem for $\pi^{0}\to\gamma\gamma$ decay
Sutherland:1967vf ; Veltman:1967 ; Adler:1969gk ; Bell:1969ts  and the
solution of the so-called $U(1)$ problem Weinberg:1975ui ; tHooft:1976snw ;
tHooft:1986ooh . Besides, there are practical applications of the
renormalization of the anomalous singlet axial current. One example is the
determination of the structure of the non-decoupling (heavy-quark) mass
logarithms in the quark form factors Collins:1978wz ; Chetyrkin:1993jm ;
Chetyrkin:1993ug ; Larin:1993ju ; Larin:1994va ; Chen:2021rft , which are
present in an anomaly-free theory such as the Standard Model. These non-
decoupling mass logarithms can be resummed by solving the corresponding
renormalization group equation where the anomalous dimension is precisely
given by that of the renormalized anomalous singlet axial current. Another
application concerns the calculation of polarized proton structure functions
(or parton distribution functions) and polarized splitting functions
Matiounine:1998re ; Kodaira:1998jn ; Vogt:2008yw ; Moch:2014sna ;
Behring:2019tus ; Tarasov:2020cwl ; Tarasov:2021yll ; Blumlein:2021ryt .

In this publication, we provide the renormalization constants of both the
singlet and non-singlet axial currents regularized with a non-anticommuting
$\gamma_{5}$ tHooft:1972tcz ; Akyeampong:1973xi ; Breitenlohner:1977hr  up to
the four-loop order in QCD. We employ, in particular, the non-anticommuting
$\gamma_{5}$ variant as prescribed in refs. Larin:1991tj ; Larin:1993tq . The
renormalization constant for a flavor non-singlet axial current was previously
determined in DR to $\mathcal{O}(\alpha_{s}^{3})$ in refs. Larin:1991tj ;
Larin:1993tq , whereas the finite part of the renormalization constant of the
flavor-singlet axial current was computed only to
$\mathcal{O}(\alpha_{s}^{2})$ due to the limitation of the approach in use
(see below). The missing $\mathcal{O}(\alpha_{s}^{3})$ finite part was
completed in ref. Ahmed:2021spj  where a derivative-free projector was
composed in order to efficiently project out the non-vanishing “anomalous”
form factor of the famous vector–vector–axial-vector (VVA)-amplitude
Adler:1969gk ; Adler:1969er , i.e. the matrix element of the axial current
operator between the vacuum and a pair of external gluons, at zero-momentum
(where the axial anomaly vanishes). The choice of the VVA-amplitude is
advantageous for determining the $\overline{\mathrm{MS}}$ renormalization
constant associated with the axial-anomaly operator (and subsequently
verifying its equality with that of the strong coupling constant
$\alpha_{s}$). This choice is, however, not very economical for studying axial
current renormalization.

The approach of refs. Gorishnii:1985xm ; Larin:1991tj ; Larin:1993tq  can be
used straightforwardly to obtain the renormalization constant for a non-
singlet axial current at $\mathcal{O}(\alpha_{s}^{4})$, provided the
propagator-type four-loop integrals can be efficiently evaluated which is
indeed the case thanks to Forcer Ruijl:2017cxj  and the analytic results for
the values of the relevant master integrals given in refs. Baikov:2010hf ;
Lee:2011jt . However, this approach does not allow to obtain the
renormalization constant for a singlet axial current at
$\mathcal{O}(\alpha_{s}^{4})$ via a four-loop calculation, because the axial
anomaly vanishes at zero momentum in the matrix element of the axial current
operator between the vacuum and a pair of external quarks. Alternatively, in
order to retain a non-vanishing anomaly in this matrix element, one has to
allow a non-zero momentum flowing through the singlet axial current. In
consequence, in order to deal with a propagator-type kinematic configuration
(and to avoid dealing with infrared/collinear divergences), one of the
external quarks must be off-shell, carrying away the same amount of input
momentum, while the other carries zero momentum. This naturally leads us to
propose to use the Ward-Takahashi identity for an axial current with off-shell
quark fields to determine the renormalization constant for a flavor singlet
(and non-singlet) axial current at $\mathcal{O}(\alpha_{s}^{4})$ via a four-
loop calculation. In this approach, the non-anticommuting $\gamma_{5}$ matrix
in the contact terms (due to off-shell external quark fields) does not require
any renormalization. We will address this in detail later in the article.

The article is organized as follows. In the next section, we recapitulate the
basic renormalization properties of the local composite operators involved in
our calculation, mainly to lay down the notations and conventions in use.
Section 3 is devoted to a detailed account of the core of our recipe, the
usage of the off-shell Ward-Takahashi identity for an axial current with a
non-anticommuting $\gamma_{5}$, for efficiently determining the
renormalization constants of axial currents, especially for the anomalous
singlet one. The perturbative results for the renormalization constants of
singlet and non-singlet axial currents determined in this way at four-loop
order in QCD are presented in section 4. We conclude in section 5.

##  2 Preliminaries

In this section, we specify the notations and conventions for the quantities
considered in the rest of this article. We use the non-anticommuting
definition of $\gamma_{5}$ in dimensional regularization, originally
introduced by ’t Hooft-Veltman tHooft:1972tcz  and Breitenlohner-Maison
Breitenlohner:1977hr

$\displaystyle\gamma_{5}=-\frac{i}{4!}\epsilon^{\mu\nu\rho\sigma}\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}\,,$
(1)

but with the Levi-Civita tensor111We use the convention
$\epsilon^{0123}=-\epsilon_{0123}=+1$. $\epsilon^{\mu\nu\rho\sigma}$ treated
according to refs. Larin:1991tj ; Zijlstra:1992kj ; Larin:1993tq , sometimes
called Larin’s prescription in the literature. As summarized in refs.
Larin:1991tj ; Larin:1993tq , the properly renormalized non-singlet and
singlet axial current operator with a non-anticommuting $\gamma_{5}$ in QCD
with $n_{f}$ massless quarks can be written as

$\displaystyle\left[J^{a,\mu}_{5,ns}\right]_{R}$ $\displaystyle=$
$\displaystyle Z_{ns}\,\bar{\psi}^{B}\,\gamma^{\mu}\gamma_{5}\,t^{a}\psi^{B}$
(2) $\displaystyle=$ $\displaystyle
Z^{f}_{ns}\,Z^{ms}_{ns}\,\bar{\psi}^{B}\,\frac{-i}{3!}\epsilon^{\mu\nu\rho\sigma}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}\,t^{a}\psi^{B}\,,$
$\displaystyle\left[J^{\mu}_{5,s}\right]_{R}$ $\displaystyle=$ $\displaystyle
Z_{s}\,\sum_{q}\,\bar{\psi}^{B}_{q}\,\gamma^{\mu}\gamma_{5}\,\psi^{B}_{q}$ (3)
$\displaystyle=$ $\displaystyle
Z^{f}_{s}\,Z^{ms}_{s}\,\sum_{q}\bar{\psi}^{B}_{q}\,\frac{-i}{3!}\epsilon^{\mu\nu\rho\sigma}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}\,\psi^{B}_{q}\,,$

where the subscript $R$ at a square bracket denotes operator renormalization.
In the non-singlet current eq.(2), $t^{a}$ denotes a (normalized) generator of
a flavor group in the representation spanned by the bare quark flavor-
multiplet $\psi^{B}$. Each $\psi^{B}_{q}$ in the singlet current eq.(3)
denotes a bare quark field of flavor $q$ with mass dimension222In ref.
Ahmed:2021spj , a global factor $\mu^{4-D}$ in the mass scale $\mu$ of
dimensional regularization was introduced in order for the mass dimension of
the r.h.s. operator be equal to the canonical dimension of the l.h.s. in four
dimensions. $(D-1)/2$, and the sum extends over all $n_{f}$ quark fields
active in the theory. In order to correctly define a hermitian axial current
with a non-anticommuting $\gamma_{5}$, the necessary “symmetrization” of the
matrix product,
$\gamma^{\mu}\gamma_{5}\rightarrow\frac{1}{2}\big{(}\gamma^{\mu}\gamma_{5}-\gamma_{5}\gamma^{\mu}\big{)}$,
before substituting eq. (1) is understood in the second equality in eqs. (2,
3).

The non-singlet axial current $J^{a,\mu}_{5,ns}$ is conserved, up to the
classically-expected mass term $2m\bar{\psi}i\gamma_{5}t^{a}\psi$ which
vanishes for massless quark fields, and hence is non-anomalous. This property
is sufficient to ensure that $J^{a,\mu}_{5,ns}$ should remain un-renormalized
if one had used an anticommuting $\gamma_{5}$, similarly to the non-
renormalization of the vector current
$\left[J^{\mu}\right]_{R}=\bar{\psi}^{B}\gamma^{\mu}\psi^{B}$. The appearance
of the ultraviolet (UV) renormalization constant $Z_{ns}\equiv
Z^{f}_{ns}\,Z^{ms}_{ns}$ is thus solely to amend the spurious terms
originating from the apparent loss of $\gamma_{5}$’s anticommutativity in the
chosen prescription (1), such that the following current conservation equation

$\displaystyle\partial_{\mu}\big{[}J^{a,\mu}_{5,s}\big{]}_{R}=0\,,$ (4)

holds for the renormalized non-singlet axial current. In this way, the
$\gamma_{5}$’s anticommutativity is effectively restored for the non-singlet
axial current Gorishnii:1985xm ; Larin:1991tj ; Larin:1993tq . This also
implies that the renormalization constant $Z_{ns}\equiv
Z^{f}_{ns}\,Z^{ms}_{ns}$, as a whole, has a vanishing anomalous dimension,
namely $\mu^{2}\frac{\mathrm{d}\ln Z_{ns}}{\mathrm{d}\mu^{2}}=0$.

On the other hand, the singlet axial current $J^{\mu}_{5,s}$ in eq. (3) is not
conserved, due to the quantum anomaly Adler:1969gk ; Bell:1969ts , and does
require genuine UV renormalization regardless of the $\gamma_{5}$ prescription
in use. Furthermore, it is known to renormalize multiplicatively Adler:1969gk
; Trueman:1979en . The factor $Z_{s}\equiv Z^{f}_{s}\,Z^{ms}_{s}$ denotes the
corresponding UV renormalization constant, conveniently parameterized as the
product of a pure $\overline{\mathrm{MS}}$-renormalization factor $Z^{ms}_{s}$
and an additional finite renormalization factor $Z^{f}_{s}$. The latter is
needed to restore the correct form of the axial Ward identity, or rather the
all-order axial-anomaly equation Adler:1969gk ; Adler:1969er , which in terms
of renormalized local composite operators reads

$\displaystyle\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{R}=a_{s}\,n_{f}\,{\displaystyle\mathrm{T}_{F}}\,\big{[}F\tilde{F}\big{]}_{R}\,,$
(5)

where ${\displaystyle\mathrm{T}_{F}}=1/2$,
$F\tilde{F}\equiv-\epsilon^{\mu\nu\rho\sigma}F^{a}_{\mu\nu}F^{a}_{\rho\sigma}=\epsilon_{\mu\nu\rho\sigma}F^{a}_{\mu\nu}F^{a}_{\rho\sigma}$
denotes the contraction of the field strength tensor
$F^{a}_{\mu\nu}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}+g_{s}\,f^{abc}A_{\mu}^{b}A_{\nu}^{c}$
of the gluon field $A_{\mu}^{a}$ with its dual form. We use the shorthand
notation $a_{s}\equiv\frac{\alpha_{s}}{4\pi}=\frac{g_{s}^{2}}{16\pi^{2}}$ for
the QCD coupling, and $f^{abc}$ denotes the structure constants of the non-
Abelian color group of QCD. The renormalization-group (RG) equation of $a_{s}$
in $D$ dimensions reads

$\displaystyle\mu^{2}\frac{\mathrm{d}\ln
a_{s}}{\mathrm{d}\mu^{2}}=-\epsilon-\mu^{2}\frac{\mathrm{d}\ln
Z_{a_{s}}}{\mathrm{d}\mu^{2}}\equiv-\epsilon+\beta\,,$ (6)

where $Z_{a_{s}}$ stands for the $\overline{\mathrm{MS}}$ coupling
renormalization constant and $\beta\equiv-\mu^{2}\frac{\mathrm{d}\ln
Z_{a_{s}}}{\mathrm{d}\mu^{2}}$ denotes the QCD beta function, i.e. the
anomalous dimension of the renormalized $a_{s}$ in 4 dimensions ($\mu$ denotes
the mass scale in dimensional regularization). In contrast to the l.h.s. of
eq. (5), the renormalization of the axial-anomaly operator $F\tilde{F}$ is not
strictly multiplicative (as known from ref. Adler:1969gk ), but involves
mixing with the divergence of the axial current operator Espriu:1982bw ;
Breitenlohner:1983pi ,

$\displaystyle\big{[}F\tilde{F}\big{]}_{R}=Z_{F\tilde{F}}\,\big{[}F\tilde{F}\big{]}_{B}\,+\,Z_{FJ}\,\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{B}\,,$
(7)

where the subscript $B$ implies that the fields in the local composite
operators are bare. In the computation of the matrix elements of the axial-
anomaly operator $F\tilde{F}$, we employ its equivalent form in terms of the
divergence of the Chern-Simons current $K^{\mu}$, namely

$\displaystyle F\tilde{F}$ $\displaystyle=$
$\displaystyle\partial_{\mu}K^{\mu}$ (8) $\displaystyle=$
$\displaystyle\partial_{\mu}\left(-4\,\epsilon^{\mu\nu\rho\sigma}\,\left(A_{\nu}^{a}\partial_{\rho}A_{\sigma}^{a}\,+\,g_{s}\,\frac{1}{3}f^{abc}A_{\nu}^{a}A_{\rho}^{b}A_{\sigma}^{c}\right)\right)\,,$

where, unlike $J^{\mu}_{5,s}$, the current
$K^{\mu}\equiv-4\,\epsilon^{\mu\nu\rho\sigma}\,\left(A_{\nu}^{a}\partial_{\rho}A_{\sigma}^{a}\,+\,g_{s}\,\frac{1}{3}f^{abc}A_{\nu}^{a}A_{\rho}^{b}A_{\sigma}^{c}\right)$
is not a gauge-invariant object Espriu:1982bw ; Breitenlohner:1983pi . The
Feynman rules for the r.h.s. of eq. (5) used in our calculation are directly
based on eq. (8). As verified explicitly to $\mathcal{O}(\alpha_{s}^{4})$ in
ref. Ahmed:2021spj , the equality $Z_{F\tilde{F}}=Z_{a_{s}}$ holds, which was
proved recently in ref. Luscher:2021bog . For the sake of self-explanatory
notations, we keep the symbol $Z_{F\tilde{F}}$ in the following discussion.

The renormalization of the operators $\partial_{\mu}J^{\mu}_{5}$ and
$F\tilde{F}$ specified above can be arranged into the following matrix form

$\displaystyle\begin{pmatrix}\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{R}\\\
\big{[}F\tilde{F}\big{]}_{R}\end{pmatrix}=\begin{pmatrix}Z_{s}&0\\\
Z_{FJ}&Z_{F\tilde{F}}\end{pmatrix}\cdot\begin{pmatrix}\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{B}\\\
\big{[}F\tilde{F}\big{]}_{B}\end{pmatrix}\,.$ (9)

The matrix of anomalous dimensions of these two renormalized operators is
defined by

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\,\ln\mu^{2}}\,\begin{pmatrix}\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{R}\\\
\big{[}F\tilde{F}\big{]}_{R}\end{pmatrix}=\begin{pmatrix}\gamma_{s}&0\\\
\gamma_{{\scriptscriptstyle FJ}}&\gamma_{{\scriptscriptstyle
F\tilde{F}}}\end{pmatrix}\cdot\begin{pmatrix}\big{[}\partial_{\mu}J^{\mu}_{5,s}\big{]}_{R}\\\
\big{[}F\tilde{F}\big{]}_{R}\end{pmatrix}\,.$ (10)

##  3 The Axial Ward-Takahashi Identity

Given the operator-level axial anomaly equation (5) in the canonical
quantization formalism, one can derive the following anomalous Ward-Takahashi
identity Adler:1969gk  in the momentum space:

$\displaystyle
q_{\mu}\,\mathrm{\Gamma}^{\mu}_{5,s}(p^{\prime},p)=-a_{s}\,n_{f}\,{\displaystyle\mathrm{T}_{F}}\,\mathrm{\Lambda}(p^{\prime},p)\,+\,\gamma_{5}\,\hat{S}^{-1}(p)\,+\,\hat{S}^{-1}(p^{\prime})\,\gamma_{5}\,,$
(11)

for the singlet axial current with massless quark fields333The trivial
identity color factor in the fundamental representation is suppressed; in case
of massive quark fields, there will be an additional contribution generated by
the classically-expected mass term $2m\bar{\psi}i\gamma_{5}\psi$..
Diagrammatically, this equation can be illustrated as in figure 1.

Figure 1: An illustration of the anomalous Ward-Takahashi identity (11). The
first diagram on the r.h.s. is only one of three types of contributions, while
the other two involve a local vertex with either three gluons or two quarks as
implied by eqs. (8) and (9).

$\mathrm{\Gamma}^{\mu}_{5,s}(p^{\prime},p)$ in eq. (11) is the amputated one-
particle irreducible (1PI) 3-point vertex function, computed using the
renormalized singlet axial current operator (i.e. eq. (3)) with a momentum
insertion $q=p^{\prime}-p$. $\hat{S}(p)$ in the contact terms on the r.h.s. is
defined by the full propagator of an off-shell massless quark,

$i\,\delta_{ij}\,\hat{S}(p)=i\,\delta_{ij}\,\frac{\,Z_{2}}{\not{p}-\mathrm{\Sigma}_{R}(p)}$
,

where $Z_{2}$ stands for the quark wavefunction renormalization constant,
$\mathrm{\Sigma}_{R}(p)$ is the renormalized quark self-energy function and
$\delta_{ij}$ is the identity color factor in the fundamental representation
of the color algebra. $\mathrm{\Lambda}(p^{\prime},p)$ denotes the amputated
1PI 3-point function

$\int d^{4}xd^{4}x^{\prime}\,e^{ip^{\prime}\cdot x^{\prime}-ip\cdot
x}\,\langle
0|\hat{\mathrm{T}}\left[\psi_{B}(x^{\prime})\,\bar{\psi}_{B}(x)\,\big{[}F\tilde{F}\big{]}_{R}(0)\right]|0\rangle|_{\mathrm{amp}}$
,

computed using the renormalized axial anomaly operator, i.e. the r.h.s. of eq.
(5).

The non-anomalous Ward-Takahashi identity for the renormalized non-singlet
axial current operator takes the usual form where the anomaly term in eq. (11)
is absent, namely,

$\displaystyle
q_{\mu}\,\mathrm{\Gamma}^{a,\mu}_{5,ns}(p^{\prime},p)=t^{a}\gamma_{5}\,\hat{S}^{-1}(p)\,+\,\hat{S}^{-1}(p^{\prime})\,\gamma_{5}t^{a}\,,$
(12)

which is very similar to the well-known Ward-Takahashi identity for a (color-
neutral) conserved vector current. Similarly,
$\mathrm{\Gamma}^{a,\mu}_{5,ns}(p^{\prime},p)$ in eq. (12) is the amputated
1PI vertex function computed using the renormalized non-singlet axial current
operator (i.e. eq. (2)) with a momentum insertion $q=p^{\prime}-p$.

One critical feature of eq.(11) and eq.(12), important for our approach,
deserves a comment: the $\gamma_{5}$ matrix in the r.h.s. contact terms
multiplies the quark self-energy functions from outside, hence never appears
within any loop correction, and clearly never on a closed fermion loop as
illustrated in figure 1. Indeed, the $\gamma_{5}$ matrix in the contact terms
of the axial Ward-Takahashi identities (2, 3) is always on an open quark line
and essentially “external”. Because of this, one can simply shift it, while
respecting anticommutativity, into the external projector used in the
computation, and subsequently combine it with the $\gamma_{5}$ present in the
said external projector. Therefore, one will not need any $\gamma_{5}$-related
renormalization constant, such as those in eqs. (2, 3), in the calculation of
the final renormalized contributions from the contact terms in the axial Ward-
Takahashi identities (11, 12). In particular, if one would combine these two
$\gamma_{5}$ matrices using the original definition by ’t Hooft-Veltman
tHooft:1972tcz  and Breitenlohner-Maison Breitenlohner:1977hr , one would have
$\gamma_{5}^{2}$ equal to an identity matrix in spinor space, and hence each
of the r.h.s. contact terms would contribute the same value as the self-energy
function of an off-shell quark without any $\gamma_{5}$ matrices involved.
Alternatively, if one proceeds literally according to the variant in refs.
Larin:1991tj ; Larin:1993tq , $\gamma_{5}^{2}$ would no longer be the identity
matrix, albeit, still proportional to it by a polynomial factor in $D$.
Consequently, the intermediate singular expressions on the r.h.s. would not be
identical to the self-energy function of an off-shell quark; but the same
finite remainder in the 4-dimensional limit will be produced (without the need
of any $\gamma_{5}$-related renormalization constant).

To determine the renormalization constants defined in eqs. (2, 3), we compute
first the bare expressions of both sides of the axial Ward-Takahashi
identities eq.(11) and eq.(12). The perturbative QCD corrections to these
matrix elements or correlation functions are computed in terms of Feynman
diagrams, which are manipulated in the usual way. The technical workflow
follows closely our previous work Ahmed:2021spj  where we refer for details,
in particular the tool chain in use. As mentioned in the introduction, to have
a non-vanishing divergence of the axial current, namely l.h.s of eqs. (2, 3),
one has to allow a non-zero momentum insertion through the axial current
operator. On the other hand, to keep a propagator-type kinematic
configuration, and also to avoid dealing with infrared/collinear divergences,
one of the external quarks must be off-shell, carrying away the same amount of
momentum inserted through the axial current operator. This then forces us to
take the kinematic configuration $q=-p$ with $p^{\prime}=0$. Consequently, the
quark self-energy term $\hat{S}^{-1}(p^{\prime})$ in the r.h.s. of eq. (11)
and eq. (12) vanishes completely in dimensional regularization. At this
kinematic limit, there is essentially just one Lorentz structure
$\not{q}\gamma_{5}$ (at least in 4 dimensions). We project the aforementioned
bare correlation functions of both sides onto this common Lorentz structure,
obtaining respectively a Lorentz scalar quantity. The pure
$\overline{\mathrm{MS}}$ parts of the renormalization constants in eqs. (2,
3), namely $Z_{s}^{ms}$ and $Z_{ns}^{ms}$, can be determined from the
remaining pole parts in these scalar projections after removing those that can
be accounted for by the wavefunction renormalization of the external quark
fields, by the $\alpha_{s}$ renormalization and by the renormalization of the
general-covariant-gauge fixing parameter $\xi$. The
non-$\overline{\mathrm{MS}}$ finite renormalization constants, $Z_{s}^{f}$ and
$Z_{ns}^{f}$, are determined by demanding the renormalized axial Ward-
Takahashi identities at $q=-p$ to hold order-by-order in the perturbative
coupling $a_{s}$ in the 4 dimensional limit.

##  4 Results

Below we present our final perturbative results for the renormalization
constants defined in eqs. (2, 3) for axial currents, as well as the $Z_{FJ}$
in eq.(9), at the four-loop order in QCD with $n_{f}$ massless quarks. The
perturbative QCD coupling $a_{s}\equiv\frac{\alpha_{s}}{4\pi}$ is renormalized
in the $\overline{\mathrm{MS}}$ scheme.

###  4.1 The non-singlet axial current

We start with the simple case, the renormalization constant $Z_{ns}\equiv
Z^{f}_{ns}\,Z^{ms}_{ns}$ for the non-anomalous non-singlet axial current. The
result for the $\overline{\mathrm{MS}}$ part reads:

$\displaystyle Z^{ms}_{ns}$ $\displaystyle=$ $\displaystyle
1\,+\,a_{s}^{2}\,\Big{\\{}\,C_{A}C_{F}\,\left(\frac{22}{3\epsilon}\right)\,+\,C_{F}n_{f}\,\left(-\frac{4}{3\epsilon}\right)\Big{\\}}$
(13) $\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(\frac{3578}{81\epsilon}-\frac{484}{27\epsilon^{2}}\right)\,+\,C_{A}C_{F}^{2}\,\left(-\frac{308}{9\epsilon}\right)\,+\,C_{A}C_{F}n_{f}\,\left(\frac{176}{27\epsilon^{2}}-\frac{832}{81\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}\,\left(\frac{32}{9\epsilon}\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{8}{81\epsilon}-\frac{16}{27\epsilon^{2}}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}\,C_{A}^{3}C_{F}\,\left(\frac{1331}{27\epsilon^{3}}-\frac{26411}{162\epsilon^{2}}+\frac{\frac{36607}{108}-154\zeta_{3}}{\epsilon}\right)\,+\,C_{A}^{2}C_{F}^{2}\,\left(\frac{121}{\epsilon^{2}}+\frac{440\zeta_{3}-\frac{29309}{54}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}n_{f}\,\left(-\frac{242}{9\epsilon^{3}}+\frac{631}{9\epsilon^{2}}+\frac{-\frac{4\zeta_{3}}{3}-\frac{8029}{81}}{\epsilon}\right)\,+\,C_{A}C_{F}^{3}\,\left(\frac{\frac{935}{6}-264\zeta_{3}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{2}n_{f}\,\left(\frac{\frac{2381}{27}-\frac{152\zeta_{3}}{3}}{\epsilon}-\frac{88}{3\epsilon^{2}}\right)\,+\,C_{A}C_{F}n_{f}^{2}\,\left(\frac{44}{9\epsilon^{3}}-\frac{206}{27\epsilon^{2}}+\frac{\frac{16\zeta_{3}}{3}+\frac{367}{162}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}n_{f}\,\left(\frac{48\zeta_{3}-\frac{40}{3}}{\epsilon}\right)\,+\,C_{F}^{2}n_{f}^{2}\,\left(\frac{4}{3\epsilon^{2}}+\frac{\frac{20}{27}-\frac{16\zeta_{3}}{3}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}n_{f}^{3}\,\left(-\frac{8}{27\epsilon^{3}}+\frac{4}{81\epsilon^{2}}+\frac{26}{81\epsilon}\right)\Big{\\}}\,.$

The definition of the quadratic Casimir color constants is as usual:
$C_{A}=N_{c}\,,\,C_{F}=(N_{c}^{2}-1)/(2N_{c})\,$ with $N_{c}=3$ in QCD and the
color-trace normalization factor ${\displaystyle\mathrm{T}_{F}}={1}/{2}$. The
result for the finite part reads:

$\displaystyle Z^{f}_{ns}$ $\displaystyle=$ $\displaystyle
1\,+\,a_{s}\,\Big{\\{}\,C_{F}\,\left(-4\right)\Big{\\}}\,+\,a_{s}^{2}\,\Big{\\{}C_{A}C_{F}\,\left(-\frac{107}{9}\right)\,+\,C_{F}^{2}\,\left(22\right)\,+\,C_{F}n_{f}\,\left(\frac{2}{9}\right)\Big{\\}}$
(14) $\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(56\zeta_{3}-\frac{2147}{27}\right)\,+\,C_{A}C_{F}^{2}\,\left(\frac{5834}{27}-160\zeta_{3}\right)\,+\,C_{A}C_{F}n_{f}\,\left(\frac{32\zeta_{3}}{3}+\frac{356}{81}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}\,\left(96\zeta_{3}-\frac{370}{3}\right)\,+\,C_{F}^{2}n_{f}\,\left(-\frac{32\zeta_{3}}{3}-\frac{62}{27}\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{52}{81}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}\,C_{A}^{3}C_{F}\,\left(\frac{10498\zeta_{3}}{27}-\frac{4120\zeta_{5}}{3}-\frac{77\pi^{4}}{30}-\frac{324575}{648}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}^{2}\,\left(-1570\zeta_{3}+\frac{17020\zeta_{5}}{3}+\frac{22\pi^{4}}{3}+\frac{10619}{6}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}n_{f}\,\left(\frac{1118\zeta_{3}}{9}-\frac{200\zeta_{5}}{3}-\frac{\pi^{4}}{45}+\frac{18841}{324}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{3}\,\left(\frac{3332\zeta_{3}}{3}-6760\zeta_{5}-\frac{22\pi^{4}}{5}-\frac{232949}{108}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{2}n_{f}\,\left(-\frac{1312\zeta_{3}}{9}-\frac{40\zeta_{5}}{3}-\frac{38\pi^{4}}{45}-\frac{1705}{81}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}^{2}\,\left(-\frac{88\zeta_{3}}{9}+\frac{4\pi^{4}}{45}+\frac{73}{54}\right)\,+\,C_{F}^{4}\,\left(564\zeta_{3}+1840\zeta_{5}+\frac{1553}{2}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}n_{f}\,\left(\frac{116\zeta_{3}}{3}+80\zeta_{5}+\frac{4\pi^{4}}{5}-\frac{850}{27}\right)\,+\,C_{F}^{2}n_{f}^{2}\,\left(\frac{88\zeta_{3}}{9}-\frac{4\pi^{4}}{45}-\frac{719}{81}\right)$
$\displaystyle+$ $\displaystyle
C_{F}n_{f}^{3}\,\left(\frac{5}{9}-\frac{16\zeta_{3}}{27}\right)\,+\,C_{1}C_{F}n_{f}\,\left(128\zeta_{3}-\frac{608}{3}\right)$
$\displaystyle+$ $\displaystyle
C_{2}C_{F}\,\left(1520\zeta_{3}+1920\zeta_{5}-\frac{32}{3}\right)\Big{\\}}\,,$

where the additional color constants are defined in terms of symmetric color
tensors444The symmetric tensor $d_{F}^{abcd}$ is defined by the color trace
$\frac{1}{6}\mbox{Tr}$
$\big{(}T^{a}T^{b}T^{c}T^{d}+T^{a}T^{b}T^{d}T^{c}+T^{a}T^{c}T^{b}T^{d}+T^{a}T^{c}T^{d}T^{b}+T^{a}T^{d}T^{b}T^{c}+T^{a}T^{d}T^{c}T^{b}\big{)}$
with $T^{a}$ the generators of the fundamental representation of the
SU($N_{c}$) group, and similarly $d_{A}^{abcd}$ for the adjoint
representation. as

$\displaystyle C_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{d_{F}^{abcd}d_{F}^{abcd}}{N_{c}^{2}-1}=\frac{N_{c}^{4}-6N_{c}^{2}+18}{96N_{c}^{2}}\,,\,$
$\displaystyle C_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{d_{F}^{abcd}d_{A}^{abcd}}{N_{c}^{2}-1}=\frac{N_{c}\left(N_{c}^{2}+6\right)}{48}\,.$
(15)

Up to three-loop order, we have reproduced the results in ref. Larin:1991tj .
We have checked that the anomalous dimension of the above $Z_{ns}$ vanishes up
to four-loop order (in the 4-dimensional limit), as expected from eq. (4). In
addition, we performed an alternative calculation of $Z_{ns}$ following the
recipe of refs. Gorishnii:1985xm ; Larin:1991tj ; Larin:1993tq  at four-loop
order, and we confirm that the same result follows.

###  4.2 The singlet axial current

Now we come to the more interesting and useful case of the anomalous singlet
axial current operator. The result for the $\overline{\mathrm{MS}}$ part
reads:

$\displaystyle Z^{ms}_{s}$ $\displaystyle=$ $\displaystyle
1\,+\,a_{s}^{2}\,\Big{\\{}\,C_{A}C_{F}\,\left(\frac{22}{3\epsilon}\right)\,+\,C_{F}n_{f}\,\left(\frac{5}{3\epsilon}\right)\Big{\\}}$
(16) $\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(\frac{3578}{81\epsilon}-\frac{484}{27\epsilon^{2}}\right)\,+\,C_{A}C_{F}^{2}\,\left(-\frac{308}{9\epsilon}\right)\,+\,C_{A}C_{F}n_{f}\,\left(\frac{149}{81\epsilon}-\frac{22}{27\epsilon^{2}}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}\,\left(-\frac{22}{9\epsilon}\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{20}{27\epsilon^{2}}+\frac{26}{81\epsilon}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}C_{A}^{3}C_{F}\,\left(\frac{1331}{27\epsilon^{3}}-\frac{26411}{162\epsilon^{2}}+\frac{\frac{36607}{108}-154\zeta_{3}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}^{2}\,\left(\frac{121}{\epsilon^{2}}+\frac{440\zeta_{3}-\frac{29309}{54}}{\epsilon}\right)\,+\,C_{A}^{2}C_{F}n_{f}\,\left(-\frac{121}{18\epsilon^{3}}+\frac{713}{36\epsilon^{2}}+\frac{-\frac{437\zeta_{3}}{6}-\frac{15593}{648}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{3}\,\left(\frac{\frac{935}{6}-264\zeta_{3}}{\epsilon}\right)\,+\,C_{A}C_{F}^{2}n_{f}\,\left(\frac{55}{6\epsilon^{2}}+\frac{\frac{46\zeta_{3}}{3}-\frac{1897}{108}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}^{2}\,\left(-\frac{22}{9\epsilon^{3}}+\frac{76}{27\epsilon^{2}}+\frac{-\frac{53\zeta_{3}}{3}-\frac{31}{81}}{\epsilon}\right)\,+\,C_{F}^{3}n_{f}\,\left(\frac{48\zeta_{3}+\frac{29}{12}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}^{2}\,\left(\frac{11}{6\epsilon^{2}}+\frac{\frac{56\zeta_{3}}{3}-\frac{851}{54}}{\epsilon}\right)\,+\,C_{F}n_{f}^{3}\,\left(\frac{10}{27\epsilon^{3}}+\frac{13}{81\epsilon^{2}}-\frac{35}{162\epsilon}\right)\Big{\\}}\,.$

The result for the finite part reads:

$\displaystyle Z^{f}_{s}$ $\displaystyle=$ $\displaystyle
1\,+\,a_{s}\,\Big{\\{}\,C_{F}\,\left(-4\right)\Big{\\}}\,+\,a_{s}^{2}\,\Big{\\{}\,C_{A}C_{F}\,\left(-\frac{107}{9}\right),C_{F}^{2}\,\left(22\right)\,+\,C_{F}n_{f}\,\left(\frac{31}{18}\right)\Big{\\}}$
(17) $\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(56\zeta_{3}-\frac{2147}{27}\right)\,+\,C_{A}C_{F}^{2}\,\left(\frac{5834}{27}-160\zeta_{3}\right)\,+\,C_{A}C_{F}n_{f}\,\left(\frac{110\zeta_{3}}{3}-\frac{133}{81}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}\,\left(96\zeta_{3}-\frac{370}{3}\right)\,+\,C_{F}^{2}n_{f}\,\left(\frac{497}{54}-\frac{104\zeta_{3}}{3}\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{316}{81}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}\,C_{A}^{3}C_{F}\,\left(\frac{10498\zeta_{3}}{27}-\frac{4120\zeta_{5}}{3}-\frac{77\pi^{4}}{30}-\frac{324575}{648}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}^{2}\,\left(-1570\zeta_{3}+\frac{17020\zeta_{5}}{3}+\frac{22\pi^{4}}{3}+\frac{10619}{6}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}n_{f}\,\left(\frac{3941\zeta_{3}}{9}-\frac{845\zeta_{5}}{4}-\frac{437\pi^{4}}{360}+\frac{98815}{1296}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{3}\,\left(\frac{3332\zeta_{3}}{3}-6760\zeta_{5}-\frac{22\pi^{4}}{5}-\frac{232949}{108}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{2}n_{f}\,\left(-\frac{8807\zeta_{3}}{18}-\frac{370\zeta_{5}}{3}+\frac{23\pi^{4}}{90}+\frac{10687}{162}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}^{2}\,\left(-\frac{1037\zeta_{3}}{18}+\frac{91\pi^{4}}{180}-\frac{235}{36}\right)\,+\,C_{F}^{4}\,\left(564\zeta_{3}+1840\zeta_{5}+\frac{1553}{2}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}n_{f}\,\left(\frac{443\zeta_{3}}{3}+260\zeta_{5}+\frac{4\pi^{4}}{5}-\frac{43025}{216}\right)\,+\,C_{F}^{2}n_{f}^{2}\,\left(\frac{484\zeta_{3}}{9}-\frac{22\pi^{4}}{45}-\frac{5971}{648}\right)$
$\displaystyle+$ $\displaystyle
C_{F}n_{f}^{3}\,\left(\frac{17}{4}-\frac{124\zeta_{3}}{27}\right)\,+\,C_{1}C_{F}n_{f}\,\left(264\zeta_{3}-520\zeta_{5}-\frac{320}{3}\right)$
$\displaystyle+$ $\displaystyle
C_{2}C_{F}\,\left(1520\zeta_{3}+1920\zeta_{5}-\frac{32}{3}\right)\Big{\\}}\,.$

Up to three-loop order, we have reproduced the results in ref. Larin:1993tq ;
Ahmed:2021spj  determined by computing the VVA-amplitudes. As can be easily
checked using these results, the differences $Z^{ms}_{s}-Z^{ms}_{ns}$ and
$Z^{f}_{s}-Z^{f}_{ns}$ all start from $\mathcal{O}(a_{s}^{2})$ and are
proportional to $n_{f}\,C_{F}$. As discussed in section 2, the renormalized
singlet axial current does have a non-zero anomalous dimension, which reads:

$\displaystyle\gamma_{s}$ $\displaystyle=$
$\displaystyle\mu^{2}\frac{\mathrm{d}\,\ln\big{(}Z^{f}_{s}\,Z^{ms}_{s}\big{)}}{\mathrm{d}\,\mu^{2}}$
$\displaystyle=$ $\displaystyle
a_{s}\,\Big{\\{}\,C_{F}\,\left(4\epsilon\right)\Big{\\}}\,+\,a_{s}^{2}\,\Big{\\{}C_{A}C_{F}\,\left(\frac{214\epsilon}{9}\right)\,+\,C_{F}^{2}\,\left(-28\epsilon\right)\,+\,C_{F}n_{f}\,\left(-\frac{31\epsilon}{9}-6\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(\left(\frac{2147}{9}-168\zeta_{3}\right)\epsilon\right)\,+\,C_{A}C_{F}^{2}\,\left(\left(480\zeta_{3}-\frac{4550}{9}\right)\epsilon\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}\,\left(\left(\frac{133}{27}-110\zeta_{3}\right)\epsilon-\frac{142}{3}\right)\,+\,C_{F}^{3}\,\left(\left(170-288\zeta_{3}\right)\epsilon\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}\,\left(\left(104\zeta_{3}-\frac{869}{18}\right)\epsilon+18\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{4}{3}-\frac{316\epsilon}{27}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}\,C_{A}^{3}C_{F}\,\left(\left(-\frac{41992\zeta_{3}}{27}+\frac{16480\zeta_{5}}{3}+\frac{154\pi^{4}}{15}+\frac{324575}{162}\right)\epsilon\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}^{2}\,\left(\left(5384\zeta_{3}-\frac{68080\zeta_{5}}{3}-\frac{88\pi^{4}}{3}-\frac{447472}{81}\right)\epsilon\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}n_{f}\,\left(\left(-\frac{15764\zeta_{3}}{9}+845\zeta_{5}+\frac{437\pi^{4}}{90}-\frac{98815}{324}\right)\epsilon-\frac{1607}{6}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{3}\,\left(\left(-\frac{5648\zeta_{3}}{3}+27040\zeta_{5}+\frac{88\pi^{4}}{5}+\frac{43967}{9}\right)\epsilon\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{2}n_{f}\,\left(\left(\frac{12334\zeta_{3}}{9}+\frac{1480\zeta_{5}}{3}-\frac{46\pi^{4}}{45}-\frac{25880}{81}\right)\epsilon+\frac{461}{2}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}^{2}\,\left(144\zeta_{3}+\left(\frac{2074\zeta_{3}}{9}-\frac{91\pi^{4}}{45}+\frac{235}{9}\right)\epsilon-\frac{82}{3}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{3}n_{f}\,\left(\left(-36\zeta_{3}-1040\zeta_{5}-\frac{16\pi^{4}}{5}+\frac{4145}{6}\right)\epsilon-63\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}^{2}\,\left(-144\zeta_{3}+\left(-\frac{1936\zeta_{3}}{9}+\frac{88\pi^{4}}{45}-\frac{530}{27}\right)\epsilon+107\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{4}\,\left(\left(-3792\zeta_{3}-7360\zeta_{5}-\frac{3950}{3}\right)\epsilon\right)\,+\,C_{F}n_{f}^{3}\,\left(\left(\frac{496\zeta_{3}}{27}-17\right)\epsilon+\frac{26}{3}\right)$
$\displaystyle+$ $\displaystyle
C_{1}C_{F}n_{f}\,\left(\left(-1056\zeta_{3}+2080\zeta_{5}+\frac{1280}{3}\right)\epsilon\right)\,+\,C_{2}C_{F}\,\left(\left(-6080\zeta_{3}-7680\zeta_{5}+\frac{128}{3}\right)\epsilon\right)\Big{\\}}\,.$

We note that due to the appearance of the finite renormalization $Z_{s}^{f}$,
the anomalous dimension of $Z_{s}$ given above contains terms linear in
$\epsilon$, due to the usage of the D-dimensional QCD beta function defined in
eq. (6).

As mentioned in the introduction, among the practical applications of the
knowledge of the renormalization of the anomalous singlet axial current is the
determination of the structure of the non-decoupling heavy-quark mass
logarithms in the axial quark form factors, which are present in an anomaly-
free theory such as the Standard Model. It is known Collins:1978wz ;
Chetyrkin:1993jm ; Chetyrkin:1993ug ; Larin:1993ju ; Larin:1994va  that the
effect of top-quark loops in axial quark form factors does not decouple in the
large top-quark mass or low-energy limit due to the presence of the axial-
anomaly type diagrams. A Wilson coefficient in front of the renormalized
singlet axial current operator can be introduced in the low-energy or heavy-
top effective Lagrangian, which encodes all non-decoupling top-quark mass
logarithms remaining in the total “physical” (non-anomalous) quark form
factors (See, for instance, refs. Chetyrkin:1993ug ; Chen:2021rft  for a
detailed discussion). The terms in this Wilson coefficient featuring the mass
logarithms at the four-loop order can be determined from the anomalous
dimension of the singlet axial current given in eq. (4.2) in 4 dimensions, by
solving the corresponding RG equation (e.g. in the form given by eq. (6.3) of
ref. Chen:2021rft ), provided the full knowledge of the lower-order results
available in refs. Chen:2021rft ; Ju:2021lah . Following the conventions of
ref. Chen:2021rft , the renormalized low-energy effective Lagrangian with the
axial current coupling to the Z-boson field $\text{Z}_{\mu}$ reads:

$\displaystyle\delta\mathcal{L}^{R}_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\mu^{\epsilon}\Big{(}Z_{ns}\,\sum_{i=1}^{n_{f}}a_{i}\,\bar{\psi}^{B}_{i}\,\gamma^{\mu}\gamma_{5}\,\psi^{B}_{i}\,+\,a_{b}\,Z_{ps}\,\big{[}J^{\mu}_{5,s}\big{]}_{B}$
(19)
$\displaystyle\,\,\,\,\quad+\,a_{t}\,C_{w}(a_{s},\mu/m_{t})\,\big{(}Z_{ns}+n_{f}\,Z_{ps}\big{)}\,\big{[}J^{\mu}_{5,s}\big{]}_{B}\Big{)}\,\text{Z}_{\mu}\,,$

where $C_{w}(a_{s},\mu/m_{t})$ with $m_{t}$ the on-shell top-quark mass is the
aforementioned Wilson coefficient. Furthermore, the axial electroweak coupling
of the quark $i$ is denoted by $a_{i}$. The renormalization constant
$Z_{ps}\equiv\frac{1}{n_{f}}\big{(}Z_{s}-Z_{ns}\big{)}$ is defined as the
difference between the non-singlet and singlet axial-current renormalization
constants defined respectively in eq. (2) and eq. (3), further normalized to
the case of a single quark flavor with axial coupling. Here, $n_{f}=5$ and
$a_{s}$ is to be taken after decoupling the top quark in the
$m_{t}\rightarrow\infty$ limit. Following from the RG invariance of the
effective Lagrangian eq. (19), the RG equation for $C_{w}(a_{s},\mu/m_{t})$
reads:

$\displaystyle\mu^{2}\frac{\mathrm{d}}{\mathrm{d}\mu^{2}}C_{w}(a_{s},\mu/m_{t})$
$\displaystyle=$
$\displaystyle\mu^{2}\frac{\partial}{\partial\mu^{2}}C_{w}(a_{s},\mu/m_{t})\,+\,\beta\,a_{s}\frac{\partial}{\partial
a_{s}}C_{w}(a_{s},\mu/m_{t})$ (20) $\displaystyle=$
$\displaystyle\frac{\gamma_{s}}{n_{f}}-\gamma_{s}\,C_{w}(a_{s},\mu/m_{t})\,,$

where only the 4-dimensional limit of $\gamma_{s}$ in eq. (4.2) is needed.
Solving eq. (20) with a perturbative ansatz for $C_{w}(a_{s},\mu/m_{t})$, the
logarithmic part of this Wilson coefficient can be explicitly reconstructed up
to $\mathcal{O}(\alpha_{s}^{4})$:

$\displaystyle C_{w}(a_{s},\mu/m_{t})|_{\text{log}}$ $\displaystyle=$
$\displaystyle a_{s}^{2}\,\Big{\\{}C_{F}\,\left(-6L_{\mu}\right)\Big{\\}}$
(21) $\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}C_{A}C_{F}\,\Big{(}-22L_{\mu}^{2}-\frac{76L_{\mu}}{3}\Big{)}\,+\,C_{F}^{2}\,\left(18L_{\mu}\right)\,+\,C_{F}n_{f}\,\Big{(}4L_{\mu}^{2}-\frac{8L_{\mu}}{3}\Big{)}\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}C_{A}^{2}C_{F}\,\Big{(}-\frac{242L_{\mu}^{3}}{3}-\frac{622L_{\mu}^{2}}{3}-\frac{10868L_{\mu}}{9}+924\zeta_{3}L_{\mu}\Big{)}\,+\,C_{F}^{3}\,\left(-63L_{\mu}\right)$
$\displaystyle+$ $\displaystyle
C_{F}n_{f}\,\Big{(}-\frac{328L_{\mu}}{9}\Big{)}\,+\,C_{A}C_{F}^{2}\,\Big{(}-792\zeta_{3}L_{\mu}+99L_{\mu}^{2}-50L_{\mu}\Big{)}$
$\displaystyle+$ $\displaystyle
C_{F}n_{f}^{2}\,\Big{(}-\frac{8L_{\mu}^{3}}{3}+\frac{8L_{\mu}^{2}}{3}-\frac{296L_{\mu}}{9}\Big{)}\,+\,C_{A}C_{F}\,\Big{(}\frac{1804L_{\mu}}{9}\Big{)}$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}\,\Big{(}\frac{88L_{\mu}^{3}}{3}+\frac{92L_{\mu}^{2}}{3}+\frac{3280L_{\mu}}{9}-24\zeta_{3}L_{\mu}\Big{)}$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}\,\Big{(}164L_{\mu}-24L_{\mu}^{2}\Big{)}\Big{\\}}\,+\,\mathcal{O}(\alpha_{s}^{5})\,,$

where $L_{\mu}\equiv\ln\frac{\mu^{2}}{m^{2}_{t}}$.555We note that the number
of the heavy decoupled quarks is set as 1 in eq. (21), and this is the reason
why there appear terms with factors like $C_{A}C_{F}$, $C_{F}n_{f}$ in the
$\mathcal{O}(\alpha_{s}^{4})$ result. The constant part of
$C_{w}(a_{s},\mu/m_{t})$, however, needs to be fixed by matching to an
explicit calculation in the heavy-top limit. We plan to come back to this in
the future.

###  4.3 The axial anomaly operator

As argued in the literature, e.g. refs. Breitenlohner:1983pi ; Bos:1992nd ;
Larin:1993tq ; Ahmed:2021spj ; Luscher:2021bog , it follows from the ABJ
equation with the proven relation $Z_{F\tilde{F}}=Z_{a_{s}}$, that the
equality

$\displaystyle\gamma_{s}=n_{f}\,{\displaystyle\mathrm{T}_{F}}\,a_{s}\,\gamma_{{\scriptscriptstyle
FJ}}$ (22)

with $\gamma_{s}$ and $\gamma_{{\scriptscriptstyle FJ}}$ defined in eq. (10)
must hold true, albeit only in the limit $\epsilon=0$. In principle, in the
matrix of renormalization constants in eq.(9) for the local composite
operators appearing in the ABJ equation, one only needs to compute $Z_{s}$.
This relation also provides a check of our 4-loop result eq. (4.2), as
$Z_{FJ}$ was already known to three-loop order Zoller:2013ixa ; Ahmed:2015qpa
.

The operator $\big{[}F\tilde{F}\big{]}_{R}$ is only required up to
$\mathcal{O}(\alpha_{s}^{3})$ in the calculation of $Z_{s}$ at four-loop order
by means of eq. (11). Nevertheless, it is possible to obtain $Z_{FJ}$ up to
$\mathcal{O}(\alpha_{s}^{4})$ using the same techniques as for the other
renormalisation constants. Indeed, we evaluate the matrix element of eq. (7)
corresponding to the first contribution on the r.h.s. in figure 1 at $q=-p$
with $p^{\prime}=0$ up to four-loop order as before. In this way we obtain:

$\displaystyle Z_{FJ}$ $\displaystyle=$ $\displaystyle
a_{s}\,\Big{\\{}C_{F}\,\left(\frac{12}{\epsilon}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{2}\,\Big{\\{}\,C_{A}C_{F}\,\left(\frac{142}{3\epsilon}-\frac{44}{\epsilon^{2}}\right)\,+\,C_{F}^{2}\,\left(-\frac{42}{\epsilon}\right)\,+\,C_{F}n_{f}\,\left(\frac{8}{\epsilon^{2}}-\frac{4}{3\epsilon}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{3}\,\Big{\\{}\,C_{A}^{2}C_{F}\,\left(\frac{484}{3\epsilon^{3}}-\frac{2378}{9\epsilon^{2}}+\frac{1607}{9\epsilon}\right)\,+\,C_{A}C_{F}^{2}\,\left(\frac{550}{3\epsilon^{2}}-\frac{2947}{9\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}\,\left(-\frac{176}{3\epsilon^{3}}+\frac{568}{9\epsilon^{2}}+\frac{\frac{164}{9}-96\zeta_{3}}{\epsilon}\right)\,+\,C_{F}^{3}\,\left(\frac{178}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}\,\left(\frac{96\zeta_{3}-\frac{548}{9}}{\epsilon}-\frac{16}{3\epsilon^{2}}\right)\,+\,C_{F}n_{f}^{2}\,\left(\frac{16}{3\epsilon^{3}}-\frac{8}{9\epsilon^{2}}-\frac{52}{9\epsilon}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
a_{s}^{4}\,\Big{\\{}\,C_{A}^{3}C_{F}\,\left(-\frac{5324}{9\epsilon^{4}}+\frac{36256}{27\epsilon^{3}}-\frac{37111}{27\epsilon^{2}}+\frac{\frac{8540\zeta_{3}}{9}-1760\zeta_{5}+\frac{110140}{81}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}^{2}\,\left(-\frac{726}{\epsilon^{3}}+\frac{50320}{27\epsilon^{2}}+\frac{-800\zeta_{3}+1760\zeta_{5}-\frac{200645}{81}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}^{2}C_{F}n_{f}\,\left(\frac{968}{3\epsilon^{4}}-\frac{5270}{9\epsilon^{3}}+\frac{352\zeta_{3}+\frac{3052}{9}}{\epsilon^{2}}+\frac{-\frac{5060\zeta_{3}}{3}+800\zeta_{5}+\frac{88\pi^{4}}{15}-\frac{35143}{162}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{2}n_{f}\,\left(\frac{220}{3\epsilon^{3}}+\frac{\frac{2143}{27}-352\zeta_{3}}{\epsilon^{2}}+\frac{\frac{5258\zeta_{3}}{3}+160\zeta_{5}-\frac{88\pi^{4}}{15}-\frac{77101}{162}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}n_{f}^{2}\,\left(-\frac{176}{3\epsilon^{4}}+\frac{196}{3\epsilon^{3}}+\frac{\frac{191}{9}-64\zeta_{3}}{\epsilon^{2}}+\frac{\frac{272\zeta_{3}}{3}-\frac{16\pi^{4}}{15}+\frac{1954}{81}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{4}\,\left(\frac{288\zeta_{3}-\frac{1397}{2}}{\epsilon}\right)\,+\,C_{F}^{3}n_{f}\,\left(\frac{22}{3\epsilon^{2}}+\frac{-\frac{280\zeta_{3}}{3}-960\zeta_{5}+\frac{5471}{9}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{F}^{2}n_{f}^{2}\,\left(\frac{32}{3\epsilon^{3}}+\frac{64\zeta_{3}-\frac{1352}{27}}{\epsilon^{2}}+\frac{-\frac{272\zeta_{3}}{3}+\frac{16\pi^{4}}{15}-\frac{1874}{81}}{\epsilon}\right)$
$\displaystyle+$ $\displaystyle
C_{A}C_{F}^{3}\,\left(\frac{\frac{18781}{9}-\frac{1088\zeta_{3}}{3}}{\epsilon}-\frac{2728}{3\epsilon^{2}}\right)\,+\,C_{F}n_{f}^{3}\,\left(\frac{32}{9\epsilon^{4}}-\frac{16}{27\epsilon^{3}}-\frac{104}{27\epsilon^{2}}+\frac{\frac{64\zeta_{3}}{9}-\frac{20}{3}}{\epsilon}\right)\Big{\\}}\,.$

##  5 Conclusion

We have computed the renormalization constants of the axial currents, both
singlet and non-singlet, in dimensional regularization up to four-loop order
in QCD, using the off-shell axial Ward-Takahashi identities with a non-
anticommuting $\gamma_{5}$. We determined, in addition, the mixing
renormalization constant $Z_{FJ}$ at the four-loop order by explicit
diagrammatic computation. With our recipe provided in this article, it seems
to even be possible to determine the renormalization constants for a flavor
singlet (and non-singlet) axial current at $\mathcal{O}(\alpha_{s}^{5})$ via a
five-loop calculation, provided the propagator-type five-loop integrals can be
evaluated.

## Acknowledgements

The work of L.C. and M.C. was supported by the Deutsche Forschungsgemeinschaft
under grant 396021762 – TRR 257.

## References

  * (1) G. ’t Hooft and M. J. G. Veltman, Regularization and Renormalization of Gauge Fields, Nucl. Phys. B44 (1972) 189–213.
  * (2) C. G. Bollini and J. J. Giambiagi, Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter, Nuovo Cim. B12 (1972) 20–26.
  * (3) S. L. Adler, Axial vector vertex in spinor electrodynamics, Phys. Rev. 177 (1969) 2426–2438.
  * (4) J. S. Bell and R. Jackiw, A PCAC puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$ model, Nuovo Cim. A60 (1969) 47–61.
  * (5) D. A. Akyeampong and R. Delbourgo, Dimensional regularization, abnormal amplitudes and anomalies, Nuovo Cim. A17 (1973) 578–586.
  * (6) P. Breitenlohner and D. Maison, Dimensional Renormalization and the Action Principle, Commun. Math. Phys. 52 (1977) 11–38.
  * (7) W. A. Bardeen, R. Gastmans, and B. E. Lautrup, Static quantities in Weinberg’s model of weak and electromagnetic interactions, Nucl. Phys. B46 (1972) 319–331.
  * (8) M. S. Chanowitz, M. Furman, and I. Hinchliffe, The Axial Current in Dimensional Regularization, Nucl. Phys. B159 (1979) 225–243.
  * (9) S. A. Gottlieb and J. T. Donohue, The Axial Vector Current and Dimensional Regularization, Phys. Rev. D20 (1979) 3378.
  * (10) B. A. Ovrut, Axial Vector Ward Identities and Dimensional Regularization, Nucl. Phys. B213 (1983) 241–265.
  * (11) D. Espriu and R. Tarrach, Renormalization of the Axial Anomaly Operators, Z. Phys. C16 (1982) 77.
  * (12) A. J. Buras and P. H. Weisz, QCD Nonleading Corrections to Weak Decays in Dimensional Regularization and ’t Hooft-Veltman Schemes, Nucl. Phys. B333 (1990) 66–99.
  * (13) D. Kreimer, The $\gamma_{5}$ Problem and Anomalies: A Clifford Algebra Approach, Phys. Lett. B237 (1990) 59–62.
  * (14) J. G. Körner, D. Kreimer, and K. Schilcher, A Practicable $\gamma_{5}$ scheme in dimensional regularization, Z. Phys. C54 (1992) 503–512.
  * (15) S. A. Larin and J. A. M. Vermaseren, The $\alpha_{s}^{3}$ corrections to the Bjorken sum rule for polarized electroproduction and to the Gross-Llewellyn Smith sum rule, Phys. Lett. B259 (1991) 345–352.
  * (16) S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B303 (1993) 113–118, arXiv:hep-ph/9302240 [hep-ph].
  * (17) F. Jegerlehner, Facts of life with $\gamma_{5}$, Eur. Phys. J. C18 (2001) 673–679, arXiv:hep-th/0005255 [hep-th].
  * (18) S. Moch, J. A. M. Vermaseren, and A. Vogt, On $\gamma_{5}$ in higher-order QCD calculations and the NNLO evolution of the polarized valence distribution, Phys. Lett. B748 (2015) 432–438, arXiv:1506.04517 [hep-ph].
  * (19) N. Zerf, Fermion Traces Without Evanescence, Phys. Rev. D101 no. 3, (2020) 036002, arXiv:1911.06345 [hep-ph].
  * (20) D. G. Sutherland, Current algebra and some nonstrong mesonic decays, Nucl. Phys. B 2 (1967) 433–440.
  * (21) M. Veltman, Theoretical aspects of high energy neutrino interactions, Proc. R. Soc. Lond. A 301 (1967) 107–112.
  * (22) S. Weinberg, The U(1) Problem, Phys. Rev. D 11 (1975) 3583–3593.
  * (23) G. ’t Hooft, Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle, Phys. Rev. D 14 (1976) 3432–3450. [Erratum: Phys.Rev.D 18, 2199 (1978)].
  * (24) G. ’t Hooft, How Instantons Solve the U(1) Problem, Phys. Rept. 142 (1986) 357–387.
  * (25) J. C. Collins, F. Wilczek, and A. Zee, Low-Energy Manifestations of Heavy Particles: Application to the Neutral Current, Phys. Rev. D 18 (1978) 242.
  * (26) K. G. Chetyrkin and J. H. Kuhn, Complete QCD corrections of order alpha-s**2 to the Z decay rate, Phys. Lett. B 308 (1993) 127–136.
  * (27) K. G. Chetyrkin and O. V. Tarasov, The alpha-s**3 corrections to the effective neutral current and to the Z decay rate in the heavy top quark limit, Phys. Lett. B 327 (1994) 114–122, arXiv:hep-ph/9312323.
  * (28) S. A. Larin, T. van Ritbergen, and J. A. M. Vermaseren, The alpha(s)**3 correction to Gamma (Z0 —$>$ hadrons), Phys. Lett. B 320 (1994) 159–164, arXiv:hep-ph/9310378.
  * (29) S. A. Larin, T. van Ritbergen, and J. A. M. Vermaseren, The Large quark mass expansion of Gamma (Z0 —$>$ hadrons) and Gamma (tau- —$>$ tau-neutrino + hadrons) in the order alpha-s**3, Nucl. Phys. B 438 (1995) 278–306, arXiv:hep-ph/9411260.
  * (30) L. Chen, M. Czakon, and M. Niggetiedt, The complete singlet contribution to the massless quark form factor at three loops in QCD, arXiv:2109.01917 [hep-ph].
  * (31) Y. Matiounine, J. Smith, and W. L. van Neerven, Two loop operator matrix elements calculated up to finite terms for polarized deep inelastic lepton - hadron scattering, Phys. Rev. D 58 (1998) 076002, arXiv:hep-ph/9803439.
  * (32) J. Kodaira and K. Tanaka, Polarized structure functions in QCD, Prog. Theor. Phys. 101 (1999) 191–242, arXiv:hep-ph/9812449.
  * (33) A. Vogt, S. Moch, M. Rogal, and J. A. M. Vermaseren, Towards the NNLO evolution of polarised parton distributions, Nucl. Phys. B Proc. Suppl. 183 (2008) 155–161, arXiv:0807.1238 [hep-ph].
  * (34) S. Moch, J. A. M. Vermaseren, and A. Vogt, The Three-Loop Splitting Functions in QCD: The Helicity-Dependent Case, Nucl. Phys. B 889 (2014) 351–400, arXiv:1409.5131 [hep-ph].
  * (35) A. Behring, J. Blümlein, A. De Freitas, A. Goedicke, S. Klein, A. von Manteuffel, C. Schneider, and K. Schönwald, The Polarized Three-Loop Anomalous Dimensions from On-Shell Massive Operator Matrix Elements, Nucl. Phys. B 948 (2019) 114753, arXiv:1908.03779 [hep-ph].
  * (36) A. Tarasov and R. Venugopalan, Role of the chiral anomaly in polarized deeply inelastic scattering: Finding the triangle graph inside the box diagram in Bjorken and Regge asymptotics, Phys. Rev. D 102 no. 11, (2020) 114022, arXiv:2008.08104 [hep-ph].
  * (37) A. Tarasov and R. Venugopalan, The role of the chiral anomaly in polarized deeply inelastic scattering II: Topological screening and transitions from emergent axion-like dynamics, arXiv:2109.10370 [hep-ph].
  * (38) J. Blümlein, P. Marquard, C. Schneider, and K. Schönwald, The three-loop polarized singlet anomalous dimensions from off-shell operator matrix elements, arXiv:2111.12401 [hep-ph].
  * (39) T. Ahmed, L. Chen, and M. Czakon, Renormalization of the flavor-singlet axial-vector current and its anomaly in dimensional regularization, JHEP 05 (2021) 087, arXiv:2101.09479 [hep-ph].
  * (40) S. L. Adler and W. A. Bardeen, Absence of higher order corrections in the anomalous axial vector divergence equation, Phys. Rev. 182 (1969) 1517–1536.
  * (41) S. G. Gorishnii and S. A. Larin, QCD Corrections to the Parton Model Rules for Structure Functions of Deep Inelastic Scattering, Phys. Lett. B 172 (1986) 109–112.
  * (42) B. Ruijl, T. Ueda, and J. Vermaseren, Forcer, a FORM program for the parametric reduction of four-loop massless propagator diagrams, Comput. Phys. Commun. 253 (2020) 107198, arXiv:1704.06650 [hep-ph].
  * (43) P. Baikov and K. Chetyrkin, Four Loop Massless Propagators: An Algebraic Evaluation of All Master Integrals, Nucl. Phys. B 837 (2010) 186–220, arXiv:1004.1153 [hep-ph].
  * (44) R. Lee, A. Smirnov, and V. Smirnov, Master Integrals for Four-Loop Massless Propagators up to Transcendentality Weight Twelve, Nucl. Phys. B 856 (2012) 95–110, arXiv:1108.0732 [hep-th].
  * (45) E. B. Zijlstra and W. L. van Neerven, ${\cal O}(\alpha_{s}^{2})$ correction to the structure function $F_{3}$ ($x,Q^{2}$) in deep inelastic neutrino - hadron scattering, Phys. Lett. B297 (1992) 377–384.
  * (46) T. L. Trueman, Chiral Symmetry in Perturbative QCD, Phys. Lett. 88B (1979) 331–334.
  * (47) P. Breitenlohner, D. Maison, and K. S. Stelle, Anomalous Dimensions and the Adler-bardeen Theorem in Supersymmetric Yang-Mills Theories, Phys. Lett. 134B (1984) 63–66.
  * (48) M. Lüscher and P. Weisz, Renormalization of the topological charge density in QCD with dimensional regularization, arXiv:2103.15440 [hep-ph].
  * (49) W.-L. Ju and M. Schönherr, The qT and $\Delta$$\phi$ spectra in W and Z production at the LHC at N3LL’+N2LO, JHEP 10 (2021) 088, arXiv:2106.11260 [hep-ph].
  * (50) M. Bos, Explicit calculation of the renormalized singlet axial anomaly, Nucl. Phys. B 404 (1993) 215–244, arXiv:hep-ph/9211319.
  * (51) M. Zoller, OPE of the pseudoscalar gluonium correlator in massless QCD to three-loop order, JHEP 07 (2013) 040, arXiv:1304.2232 [hep-ph].
  * (52) T. Ahmed, T. Gehrmann, P. Mathews, N. Rana, and V. Ravindran, Pseudo-scalar Form Factors at Three Loops in QCD, JHEP 11 (2015) 169, arXiv:1510.01715 [hep-ph].