Renormalization of the flavor-singlet axial-vector current and its anomaly in dimensional regularization

The renormalization constant ZJ of the flavor-singlet axial-vector current with a non-anticommuting γ5 in dimensional regularization is determined to order αs3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_s^3 $$\end{document} in QCD with massless quarks. The result is obtained by computing the matrix elements of the operators appearing in the axial-anomaly equation ∂μJ5μR=αs4πnfTFFF˜R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\left[{\partial}_{\mu }{J}_5^{\mu}\right]}_R=\frac{\alpha_s}{4\pi }{n}_f{\mathrm{T}}_F{\left[F\tilde{F}\right]}_R $$\end{document} between the vacuum and a state of two (off-shell) gluons to 4-loop order. Furthermore, through this computation, the equality between the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} renormalization constant ZFF˜\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {Z}_{F\tilde{F}} $$\end{document} associated with the operator FF˜R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\left[F\tilde{F}\right]}_R $$\end{document} and that of αs is verified explicitly to hold true at 4-loop order. This equality automatically ensures a relation between the respective anomalous dimensions, γJ=αs4πnfTFγFJ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\gamma}_J=\frac{\alpha_s}{4\pi }{n}_f{\mathrm{T}}_F{\gamma}_{FJ} $$\end{document}, at order αs4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_s^4 $$\end{document} given the validity of the axial-anomaly equation which was used to determine the non-MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} piece of ZJ for the particular γ5 prescription in use.


Introduction
Dimensional regularization (DR) [1,2] has become the choice of regularization framework that underlies most of the modern higher order perturbative calculations in Quantum Chromodynamics (QCD). This is due to its many celebrated features. For example, it manifestly preserves the gauge symmetry and Lorentz invariance of QCD, renders all loop integrals invariant under arbitrary loop momentum shifts, and allows one to handle both ultraviolet (UV) and infrared/collinear (IR) divergences in the same manner without introducing separate mass scales.
On the other hand, ever since the introduction of DR, it has been recognized that special attention is required in the treatment of γ 5 , an intrinsically 4-dimensional object. At the root of the issue is that a fully anticommuting γ 5 is algebraically incompatible with the Dirac algebra in D ( = 4) dimensions. However, the anticommutativity of γ 5 is essential for the concept of chirality of spinors in 4 dimensions and (non-anomalous) chiral symmetries in quantum field theory (QFT). Furthermore, there is one well-known subtlety related to γ 5 and chiral symmetries in QFT, namely, that the flavor-singlet axial-vector current, 1 defined with γ 5 , exhibits a quantum anomaly, the axial or Adler-Bell-Jackiw (ABJ) anomaly [3,4], in its divergence. An anticommuting γ 5 together with the invariance of loop integrals under arbitrary loop-momentum shifts, however, leads to the absence of the axial anomaly, which can be easily checked at one-loop order which subsequently holds to all orders, due to the Adler-Bardeen theorem [5].
A proper and practical treatment of γ 5 in DR must reproduce the known physical properties of γ 5 -related objects, established in a more traditional or physical regularization JHEP05(2021)087 scheme (e.g. the Pauli-Villars regularization [6]). To be more specific, the following criteria characterise a consistent γ 5 prescription in dimensional regularization: (i) in the case of axial-anomaly-free Feynman diagrams (e.g. with all fermion lines having an even number of axial currents), the anticommutativity of γ 5 must be effectively ensured if not directly asserted from the outset; (ii) in the case of Feynman diagrams with a genuine axial anomaly, the all-order axialanomaly equation [3,5] must be preserved.
Despite the aforementioned difficulties, various γ 5 prescriptions in DR have been developed in the literature [1,[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] over a span of nearly 50 years, albeit each with its own pros and cons. In spite of the amount of work, the subject is still capable of holding some amusing surprises, as recently demonstrated in ref. [22], where a set of individually finite Feynman diagrams with an axial current regularized with a non-anticommuting γ 5 in DR, derived from the heavy-top effective Lagrangian for qq → ZH, calls for additional manual corrections. 2 In this article, we use a particular variant of the original γ 5 prescription (and γ µ γ 5 ) [1,7,8] in DR, as constructively defined in refs. [17,18], where γ 5 does not anticommute with all γ µ in D dimensions. We compute the renormalization constants of the local composite operators in the axial-anomaly equation with this γ 5 prescription within QCD. 3 Closely related to the axial anomaly, it was known already since the early work [3] that these operators, namely, the divergence of the axial current operator and the axial-anomaly operator, require additional non-trivial renormalizations in an Abelian gauge theory even with Pauli-Villars regularization (where γ 5 is just the standard 4-dimensional anticommuting Dirac matrix). The situation becomes only more involved with a non-anticommuting γ 5 prescription in DR due to the appearance of spurious terms originating from the apparent loss of γ 5 's anticommutativity. With a constructive non-anticommuting γ 5 in DR, every Feynman diagram derived from the Feynman rules of the theory under consideration is mathematically unambiguously determined. However, the price one has to pay in exchange for this mathematical clarity, besides the increment in the cost of computing traces with a large number of γ matrices, is that one needs a non-trivial renormalization to eliminate some spurious terms whose presence spoils the correct form of chiral Ward identities, both for the anomalous and non-anomalous axial currents. For this reason, a pure MS renormalization [26] is typically insufficient for γ 5 -related objects in QFT [10,11,13,17,18,27,28].
As far as the work of this article is concerned, the renormalization constants for a flavor non-singlet axial current and a pseudoscalar current were determined in DR to O(α 3 s ) in refs. [17,18], whereas that of the flavor-singlet axial current was computed only to O(α 2 s ). In addition, these renormalization constants including the -suppressed terms at O(α 2 s ) were determined and discussed in the calculation of polarized splitting functions [29][30][31][32].

JHEP05(2021)087
The renormalization constants needed for the axial-anomaly operator were known to O(α 2 s ) from refs. [18,28] and were computed to O(α 3 s ) in refs. [33,34]. In this article, we extend the knowledge of the renormalization constant of the singlet axial current to O(α 3 s ) by computing the matrix elements of the operators in the axial-anomaly equation between the vacuum and a pair of (off-shell) gluons to 4-loop order. With the same computational set-up, the MS renormalization constant associated with the axial-anomaly operator is determined to O(α 4 s ) and subsequently shown to coincide with that of α s , a result claimed to hold in ref. [35] using the background field method [36]. The consequence of this equality in QCD in the low-energy region was discussed to O(a 2 s ) in ref. [37] in the context of a heavy-top effective Lagrangian describing the interactions of a CP-odd scalar with gluons.
The article is organized as follows. In section 2, we specify the prescription for γ 5 in general and the (singlet) axial current in particular. Furthermore, we describe the renormalization of operators in the axial-anomaly equation. In section 3, we explain in detail the simple projector employed to project out the matrix elements of operators involved with the chosen special kinematics needed to extract the relevant renormalization constants. The computational details are exposed in section 4, including the renormalization procedure, followed by our results on the renormalization constants of the singlet axial-current, the axial-anomaly operator, with a discussion of their implications in section 5. We conclude in section 6.

Renormalization of the singlet axial current
Multi-loop calculations in QFT involving axial currents in dimensional regularization [1,2] face the issue of defining inherently 4-dimensional objects, e.g. Dirac's γ 5 (and the closely related Levi-Civita tensor µνρσ ), properly in D dimensions. In this article, we use the non-anticommuting definition of γ 5 in dimensional regularization, originally introduced by 't Hooft-Veltman [1] and Breitenlohner-Maison [8] but with the Levi-Civita tensor 4 µνρσ treated according to refs. [17,18,38]. In particular, the contraction of the µνρσ tensor in eq. (2.1) with the one from an external projector will be done according to the usual mathematical identity in 4 dimensions but with the Lorentz indices of the resulting spacetime metric tensors all taken as D-dimensional. Due to this different treatment of µνρσ in the above constructive definition of γ 5 in DR (compared to the original), this prescription is sometimes called Larin's prescription in the literature. The spacetime dimension D is parameterized conveniently as D ≡ 4 − 2 .

JHEP05(2021)087
apparent violation of the anticommutativity of γ 5 with γ µ in the "evanescent" 4 − D = 2 dimensional space gives rise to some additional spurious UV-divergence-originated terms in dimensionally regularized amplitudes, which typically leads to the loss of the correct form of chiral Ward identities for both non-singlet and singlet axial currents. This calls for non-trivial UV renormalization of γ 5 -related objects [10,13,17,18,27,28,39,40]. Due to the aforementioned technical origin of these spurious terms, the amendment is typically realized in the form of some additional UV renormalization constants on top of the pure MS renormalization constants, introduced specifically for the axial or γ 5 -dependent part of the amplitudes.
In this article, we work in QCD with n f massless quarks. As nicely summarized in refs. [17,18], the properly renormalized singlet axial current, with γ µ γ 5 defined in eq. (2.2), can be written as where ψ B denotes a bare quark field with mass dimension (D − 1)/2 and the subscript R at a square bracket denotes operator renormalization. A factor µ 4−D in the mass scale µ of dimensional regularization has been inserted in eq. (2.3) 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. The sum extends over all n f quark fields. Here and below J µ 5 denotes the bare flavor-singlet axial current. It is known to renormalize multiplicatively [3,27], as it is the only local composite current operator in the context of QCD that has the correct mass dimension and conserved quantum numbers (which are preserved under renormalization). The factor Z J ≡ Z f 5 Z ms 5 denotes the UV renormalization constant of the current, conveniently parameterized as the product of a pure MS-renormalization part Z ms 5 and an additional finite renormalization factor Z f 5 . The latter is needed to restore the correct form of the axial Ward identity, namely, the all-order axial-anomaly equation [3,5], which reads in terms of renormalized local composite operators 5 where the gluon field A a µ with its dual form. We use the shorthand notation a s ≡ αs 4π = g 2 s 16π 2 for the QCD coupling, and f abc denotes the structure constants of the non-Abelian color group of QCD. In contrast to the l.h.s. of (2.4), the renormalization of the axial-anomaly operator FF is not strictly multiplicative (as known from ref. [3]), but involves mixing with the divergence of the axial current operator [13,35], where the subscript B implies that the fields in the local composite operators are bare, and a mass-dimension compensation factor µ 4−D has been introduced as in eq. (2.3). In the computation of the matrix elements of the axial-anomaly operator FF , we employ its equivalent form in terms of the divergence of the "axial gluon current" K µ (a Chern-Simons 3-form), namely [13,35]. The practical advantage of using this alternative form of the axial-anomaly operator in perturbative computations is already clear from the fact that the number of gluon fields therein is one less than in its original defining form. We derive the Feynman rules for the r.h.s. of eq. (2.4) directly based on eq. (2.6). These rules are listed in figure 1, and are used in our computations to be presented in the following sections.
The renormalization of the operators ∂ µ J µ 5 and FF specified, respectively, in eq. (2.3) and eq. (2.5) can be arranged into the following matrix form The matrix of anomalous dimensions of these two renormalized operators is defined by Inserting eq. (2.7) into eq. (2.8) leads to the following expressions for anomalous dimensions JHEP05(2021)087 in terms of the renormalization constants: (2.9) Before we plunge into the technical details of the computation, let us quickly comment on our choice of computing the matrix element of the flavor-singlet axial current J µ 5 , or rather its divergence ∂ µ J µ 5 , between the vacuum and a pair of off-shell gluons. The factor in QCD is an overall UV-renormalization constant for the axial-current operator, which by itself carries no Lorentz index and does not distinguish the components of the axial current. The very same universal UV-renormalization constant for the axialcurrent operator, at least the pure MS part Z ms 5 , can be determined by computing the matrix elements of J µ 5 , or alternatively those of its "scalar component" ∂ µ J µ 5 , between the vacuum and a pair of external quarks, or any other external states, provided these matrix elements are not vanishing with the chosen kinematics. To further fix the additional non-MS finite Z f 5 for an anomalous singlet-axial current, one could compute the matrix elements of both sides of the axial-anomaly equation (2.4) between the same external states, and subsequently demand the matching between the l.h.s. and r.h.s. matrix elements order by order in perturbation theory. Having in mind that, besides determining Z J ≡ Z f 5 Z ms 5 , we would also like to verify the equality between the MS renormalization constant Z FF and that of a s at O(α 4 s ), we choose to only compute the matrix elements with two external gluons. However, to determine Z f 5 to O(a N s ) from this choice of external states, an (N+1)loop computation is needed, because in this case the matrix element of ∂ µ J µ 5 starts at oneloop order (i.e. the anomalous fermion-triangle diagram). In particular, Z f 5 was determined to O(a 2 s ) from a 3-loop computation of these matrix elements in ref. [18], while in this article we will extend this result to one order higher in a s .

The axial-anomaly projector
We choose to determine the renormalization constant Z J ≡ Z f 5 Z ms 5 of the singlet axial current to O(a 3 s ), and Z FF to O(a 4 s ), in DR by computing matrix elements of operators appearing in the axial-anomaly equation between the vacuum and a pair of off-shell gluons evaluated at a specifically chosen single-scale kinematic configuration [18,28], as illustrated in figure 2.
Let us denote by 0|T | amp the amputated one-particle irreducible (1PI) vacuum expectation value of the time-ordered covariant product of the (singlet) axial current and two gluon fields in coordinate space with open Lorentz indices. Subsequently, we introduce the following rank-3 matrix element Γ µµ 1 µ 2 lhs (p 1 , p 2 ) in momentum space, defined by where translation invariance has been used to shift the coordinate of one gluon field to the origin, and the resulting total momentum conservation factor, ensuring p 2 = −q − p 1 , is implicit. At leading order, this is just the matrix element corresponding to the famous one-loop fermion-triangle diagram [3,4] with the polarization vectors of the external gluons stripped off and without imposing on-shell constraints on the incoming gluon momenta. The higher order corrections consist of all 1PI diagrams with amputated external gluon legs to the specific loop order in question. However, rather than performing the projection for the axial anomaly literally as devised in eq. (19) of ref. [18] where a derivative w.r.t. the momentum q is taken before going to the limit q µ → 0 (see also ref. [28]), we use instead the following projector and directly project Γ µµ 1 µ 2 (p 1 , −p 1 ) onto this structure right at q = 0 (i.e. p 2 = −p 1 ). This leads to the scalar mass-dimensionless matrix element In fact, eq. (3.2) encodes the one and only "physical" Lorentz structure of Γ µµ 1 µ 2 lhs (p 1 , −p 1 ) surviving at the chosen kinematics q = 0, p 2 = −p 1 under the condition of having one Levi-Civita tensor and being Bose-symmetric. Consequently, the projector P µµ 1 µ 2 projects out the form factor in front of this unique structure, which is not vanishing in the limit q = 0. Strictly speaking, with all Lorentz indices taken to be D-dimensional as stated below eq. (2.1), the squared norm of the Lorentz structure µµ 1 µ 2 ν p ν 1 is equal to (6-11D+6D 2 -D 3 )p 1 · p 1 , but it is known [41,42] that the parameter D here can be safely set to 4 consistently throughout the computation in DR without problem. 6 In fact, it is also convenient to do so in our computational set-up, to be discussed in the next section, as the JHEP05 (2021)087 projection of M lhs is separated from the evaluation and substitution of each individual loop integral therein.
The scalar quantity M lhs is projected out and evaluated right at the limit of zero momentum q = −(p 1 +p 2 ) flowing through the inserted operator J µ 5 (y), but with off-shell gluon momentum p 2 1 = 0. The reasons for choosing such a special single-scale kinematic configuration, common in computations of QCD UV renormalization constants in MS (e.g. in refs. [43][44][45][46][47][48][49][50]), are the following: • the possible IR divergences in Γ µµ 1 µ 2 lhs (p 1 , p 2 ) are nullified and all the poles in evaluated with off-shell momentum p 1 are of UV origin owing to the IR-rearrangement [51]; • all loop integrals involved are reduced to massless propagator-type integrals which are well studied and known numerically [52] as well as analytically [53,54] to 4-loop order.
Because p 2 1 = 0 is the only physical scale involved in the matrix elements under consideration, we set p 2 1 = 1 from the outset in our computations without loss of generality. Although, at first sight, the projector P µµ 1 µ 2 in eq. (3.2) does not seem to have anything to do with the axial anomaly, it can be shown that with the appropriate regularity condition it is indeed equivalent to the operation devised for projecting out the anomaly in eq. (19) of ref. [18]; see appendix A. The very same projector P µµ 1 µ 2 is also used in extracting a scalar quantity from the matrix element of FF R between the vacuum and the same external (off-shell) gluon state. In order to be able to apply eq. (3.2), it is crucial to use the equivalent form of the axial-anomaly operator FF = ∂ µ K µ in terms of the divergence of the current K µ . To be specific, we define in analogy to eq. (3.1): Comments below eq. (3.1) apply here as well. Subsequently, one contracts Γ µµ 1 µ 2 rhs (p 1 , −p 1 ) with P µµ 1 µ 2 , yielding a r.h.s. counterpart to eq. (3.3): The perturbative corrections to M rhs in terms of Feynman diagrams must be computed using the Feynman rules involving K µ listed in figure 1.

Details of the computation
We are now in the position to discuss the technical aspects of the computation of the scalar projections M lhs and M rhs to 4-loop order in DR out of which Z J ≡ Z f 5 Z ms 5 is determined to order α 3 s , and Z FF to order α 4 s . FORM [56]. The computations are done in a general Lorentz-covariant gauge for the gluon field, with the general-covariant-gauge fixing parameter ξ defined through the gluon propagator i , k being the momentum of the gluon. With this parametrization, the Feynman gauge corresponds to ξ = 0. As will become clear later, keeping the ξ dependence in the results for M lhs and M rhs to at least 3-loop order, does not only serve as a check of ξ-independence of the renormalization constants extracted to this order, but is also necessary to achieve the full UV renormalization of M lhs and M rhs at the 4-loop order.
Up to 3-loop order, the reduction of loop integrals in the bare M lhs and M rhs to master integrals, by means of integration-by-parts identities (IBP) [57,58], is done using IdSolver [55], while analytical expressions for master integrals computed in refs. [53,54] are used. 8 For the reduction of 4-loop (massless propagator-type) integrals in M lhs and M rhs , we employ the program Forcer [59] which is specifically designed for the parametric reduction of this type of integrals with high performance. To this end, we first extract the list of unreduced loop integrals appearing in our unreduced M lhs and M rhs , the numbers of which are about 6 × 10 4 and 9 × 10 4 respectively. We then prepare for each loop integral the Forcer input form (which is essentially a compact notation encoding the incidence matrix of the integral's graph representation 9 ), and subsequently make a parallel 7 The C++ library DiaGen provides, besides diagram generation for arbitrary Feynman rules, topological analysis tools and an interface to the C++ library IdSolver that allows to directly apply integration-byparts identities to the integrals occurring in the generated diagrams. IdSolver has been originally written for the calculation of ref. [46], while DiaGen predates this software. 8 In particular, we take the analytic expressions as provided in the ancillary files of ref. [53]. 9 Forcer does not accept graphs with vertices with degree higher than 4, which are typically associated JHEP05(2021)087  In a parallel set-up serving as a cross-check, the contributing QCD diagrams are generated symbolically using the diagram generator QGRAF [61], with the number of diagrams listed in the fourth and fifth column in table 1. The difference in the number of diagrams explicitly taken into account between the two set-ups is mainly due to the diagrams with sub loop-graphs rendered scaleless because of q = 0, which are not explicitly excluded in the set-up with QGRAF. The diagrams are subsequently processed through a series of in-house routines based on FORM in order to apply the Feynman rules, perform SU(N c ) color, D-dimensional spinor and Lorentz algebras. With the projector of eq. (3.2), M lhs (and M rhs ) is expressed as a linear combination of a large number of scalar Feynman integrals belonging to the family of massless propagator-type loop integrals at each loop order, and REDUZE2 [62,63] is used to find the appropriate loop momentum shifts. The scalar with unreduced 4-loop integrals of sub-sectors resulting from canceling certain loop propagators between denominator and numerator. We overcome this by introducing auxiliary graphs with only degree-3 and -4 vertices where the contracted graphs representing the sub-sector integrals can be embedded. integrals present in the projected matrix elements are reduced to linear combinations of 28 master integrals using FIRE6 [64,65] combined with LiteRed [66,67], and the analytical results of these master integrals are taken from ref. [53]. As expected, the time for performing the IBP reductions using the Laporta-algorithm [68] based tools is significantly longer than that needed by Forcer. Insertion of the IBP table (in terms of symbolic masters) into M lhs and M rhs , and subsequent simplification of the final analytical results, are performed using Mathematica and fermat.

JHEP05(2021)087
Despite the superficial difference in the numbers of Feynman diagrams generated between these two set-ups, as well as the difference in the implemented technical treatments, the analytical results for the bare M lhs and M rhs were found to be identical, to the order needed for determining the renormalization constants presented below.
UV renormalization. In order to renormalize the bare expressions of M lhs and M rhs within our computational set-ups as described above, we make use of the multiplicative renormalizability of the QCD Lagrangian and of the additional local composite operators (apart from the operator mixing discussed in section 2). To be definite about the notations and conventions in use, let us start with the renormalization of the bare QCD couplingâ s , where the dimensionless renormalized coupling a s ≡ αs 4π = g 2 s 16π 2 was already introduced below eq. (2.4). The bare couplingâ s has mass-dimension 2 as is exhibited on the r.h.s. of (4.1). In the MS scheme for dimensionally-regularized loop integrals, S = (4π) e − γ E (with γ E the Euler constant). All UV poles on the r.h.s. of eq. (4.1) are explicitly encoded in JHEP05(2021)087 Z as whose dependence on the scale µ is implicit and enters solely through the renormalized coupling a s . (The dependence on µ of these quantities is suppressed from here onward whenever there is no confusion.) The independence ofâ s on µ implies the renormalizationgroup (RG) equation of a s in D dimensions: where β ≡ −µ 2 d ln Za s dµ 2 denotes the QCD beta function, i.e., the anomalous dimension of the renormalized a s in 4 dimensions. Below we set µ = 1 and eq. (4.1) reduces toâ s = Z as a s . All loop integrals are evaluated in the MS scheme.
The multiplicative renormalizability of the QCD Lagrangian and of the axial current operator J µ 5 in eq. (2.3) implies that M lhs can be renormalized as where we have suppressed the kinematic dependence of M lhs on p 1 · p 1 (set to 1), and as in eq. (4.1) we used symbols with a hat to denote the unrenormalized quantities while the hat is omitted in the renormalized counterparts. Furthermore, in the last line of the above equation, we have introduced a shorthand notationM lhs ≡ Z ms 5 Z 3Mlhs (Z as a s , 1 − Z 3 + Z 3 ξ) denoting the purely MS renormalized form of the matrix element, which, after an additional finite renormalization by Z f 5 , gives the final properly renormalized M lhs . There are two more points in the formula eq. (4.3) that require comments. First, since the external gluons are set off-shell, we need to include the corresponding MS wavefunction renormalization constant Z 3 . Second, the matrix element M lhs projected out using the projector eq. (3.2) at the chosen kinematics is actually dependent on the covariant-gauge fixing parameter ξ, which by itself requires renormalization in QCD. With ξ defined through the gluon propagator i k 2 −g µν + ξ k µ k µ k 2 (in order to have less terms when inserted into diagrams with many gluons), its renormalization readsξ = 1 − Z 3 + Z 3 ξ where one has used the fact that 1 −ξ renormalizes multiplicatively with the renormalization constant Z ξ = Z 3 (see, e.g. [46,47]). This is why keeping the explicit ξ-dependence inM lhs , at least to 3-loop order, is necessary to achieve its complete UV renormalization at 4-loop order. After inserting the explicit expressions of the needed renormalization constants into eq. (4.3), its r.h.s. will be expanded and subsequently truncated to O(a 4 s ). Taking the operator mixing described in eq. (2.5) into account, the UV renormalization of M rhs proceeds as follows: Comments given above apply here as well. We note that the renormalization of M rhs chosen in eq. With the analytical expressions of these bare matrix elements to 4-loop order (provided as supplementary material associated with this article) at hand, we present below our results for Z J ≡ Z ms 5 Z f 5 and Z FF .

Results and discussions
Whereas the off-shell matrix elements M lhs and M rhs depend on the gluon-field gaugefixing parameter ξ, the renormalization constant Z J ≡ Z ms 5 Z f 5 of the gauge-invariant axial-current operator J µ 5 is independent of ξ. As shown in eq. (4.3) and eq. (4.4), the ξ dependence of these matrix elements should be kept to 3-loop order 10 in order to renormalize them at 4-loop order. We only computed the 4-loop matrix elements in the Feynman gauge ξ = 0.
By using eq. (4.3), one can extract the MS renormalization constant Z ms 5 to O(a 3 s ) from the poles of the 4-loop expression ofM lhs . We obtain the same result as the one given in refs. [18,34]. For reader's convenience, we document here the analytic result: The definition of the quadratic Casimir color constants is as usual: and in the same gauge, the MS renormalized r.h.s. is given to O(a 4 s ) by lhs andM rhs with full ξ-dependence to 3-loop order can be found in the supplementary material associated with this article.

JHEP05(2021)087
It is obtained by perturbatively expanding the ratio of the finiteM lhs and a s n f T F M rhs to O(a 3 s ). The first two orders of eq. (5.4) agree with ref. [18], while the third order terms are our new result. Furthermore,M lhs to 3-loop order and M rhs to 2-loop order (with full ξ dependence) are checked to be the same as those given in ref. [18], which serves as an indirect confirmation of the statement made in section 3: our projector eq. (3.2) is equivalent to the projection operation devised in ref. [18]. In addition, we have verified that the ratio determined above is independent of µ by inserting expressions of M lhs and M rhs with explicitly restored µ dependence. Because Z 3 appears as a global factor common on both sides, its incorporation in the calculation is not necessary if one is just interested in determining Z f 5 . However, without including Z 3 , the matrix elements of both sides are not fully UV renormalized, and one needs to be careful with truncating these UV-singular expressions to appropriate orders in . The result for Z f 5 in Quantum Electrodynamics (QED) can be obtained from eq. (5.4) by the known substitution rules: where α denotes the QED fine-structure constant and Finally, we note that the difference between Z f 5 in eq. (5.4) and the additional finite renormalization constant for the non-singlet axial current computed to O(a 3 s ) in ref. [17] starts at O(a 2 s ) and is proportional to n f C F . With both Z f 5 and Z ms 5 known to O(a 3 s ), we can now compute the full anomalous dimension γ J according to eq. (2.9), The 4-dimensional limit of γ J was determined in ref. [18], with which we agree. Thedependent parts in eq. (5.5) to O(a 2 s ) agree with those given in ref. [34], apart from a pure term at O(a 0 s ) due to the µ 2 factor introduced in eq. (2.7). We stress that the full form of γ J given by eq. (5.5), including its -dependent terms, should be used when trying to solve eq. (2.9) for renormalization constants. As far as γ J at O(a 3 s ) is concerned, the newly found O(a 3 s ) terms of Z f 5 contribute only to its -dependent part. However, when γ J will be computed at O(a 4 s ) these terms will affect its 4-dimensional part as well. The matrix element M rhs is only needed to 3-loop order in eq. . This equality was shown in ref. [35] to ensure the extension of the Adler-Bardeen theorem to a non-abelian gauge theory. The equality was further claimed to hold true in the same reference by showing that, apart from the mixing with the divergence of the axial-current operator, the operatorâ s FF does not need any multiplicative renormalization in a gauge-invariant regularization scheme using the background field method [36]. We refrain from reproducing here the 4-loop expression of Z FF , which contains three non-quadratic color constants, as it is checked to be in perfect agreement with Z as (whose expression can be found in refs. [44][45][46][47]). Curious readers are invited to repeat this check in accordance with eq. (4.4) using the bare matrix elements provided in the accompanying supplemental material. We note that among the three non-quadratic color constants involved in this calculation, the color factor 12 d abcd is absent in the 1 pole of the bare 4-loop M rhs . It appears, however, both in Z FF and Z 3 , but cancels in the product Z FF Z 3 . The Abelian version of eq. (5.6) was known already since the early reference [3]. This relation was first verified explicitly in QCD to O(a s ) in ref. [13], to O(a 2 s ) in refs. [18,28] and subsequently to O(a 3 s ) in refs. [33,34], and now finally in this article to O(a 4 s ) where non-quadratic color Casimirs start to appear. Before we move on, it is interesting to note that with eq. (5.6) the renormalization constant in front ofM rhs in eq. (4.4) reduces to unity in QED owing to the special allorder equality between the fermion field renormalization constant Z 2 and the QED-vertex renormalization constant Z 1 . To be more specific, in QED it reads In consequence, the only UV counterterms needed in this case to renormalizê M rhs are due to the operator mixing of FF with ∂ µ J µ 5 . Subjecting both sides of the axial-anomaly equation (2.4) to the logarithmic derivative µ 2 d d µ 2 , one obtains [18] ( with anomalous dimensions defined in eq. (2.9), which actually holds only in the 4dimensional limit = 0. The equality between the MS renormalization constants, Z FF and Z as , explicitly verified to 4-loop order, implies the equality γ FF = −β + = −µ 2 d ln as dµ 2 , to the same perturbative order. Clearly, eq. (5.7), which is a natural consequence of the allorder axial-anomaly equation (2.4), combined with the equality γ FF = −β + necessarily implies that must hold true, albeit only in the limit = 0. This relation has been checked by explicit calculations to O(a 2 s ) in ref. [18] and O(a 3 s ) in refs. [33,34], and discussed in refs. [29][30][31][32]. This relation was also shown in refs. [28,35] to be natural in a non-Abelian extension of the Adler-Bardeen theorem. 12 The symmetric tensor d abcd F is defined by the color trace 1 6 Tr with T a the generators of the fundamental representation of the SU(Nc) group.

JHEP05(2021)087
Reformulating the axial-anomaly equation (2.4) in terms of the bare fields, one arrives at where we note the appearance of the factor µ 2 , as well as the RG-invariant bare couplinĝ a s , and all terms have their canonical mass-dimensions as in four dimensions. Applying the logarithmic derivative in µ 2 to the axial current on the l.h.s. of eq. (5.9), we obtain which is explicitly satisfied with the anomalous dimensions determined in the present publication but only in the 4-dimensional limit owing to eq. (5.8). Clearly, the RG-invariance of the operator on the r.h.s. of eq. (5.9) is exact in D dimensions as long as the Levi-Civita tensor in the definition of FF is algebraically consistently defined. The RG-invariance of the l.h.s. operator, on the other hand, is only conditionally true in the 4-dimensional limit. This reflects the fact that eq. (2.4), as well as eq. (5.9), should be regarded as a relation valid in the 4-dimensional limit, rather than a D-dimensional identity. It was shown by explicit calculations in ref. [13] that the coefficient Z J − n f T F a s Z F J in front of the bare axial current J µ 5 B on the l.h.s. of eq. (5.9) reduces to unity up to O(a 2 s ) in QCD with the γ 5 prescription of ref [10]. Apparently, with Larin's prescription of γ 5 [17,18], this coefficient is no longer unity, but the axial current in eq. (5.9) remains RG-invariant, albeit only in the 4-dimensional limit.

Conclusion
We have extended the knowledge of the renormalization constant Z J , and in particular, of the finite non-MS factor Z f 5 , of the flavor-singlet axial-vector current regularized with Larin's prescription of γ 5 to O(α 3 s ) in QCD through computations of matrix elements of operators appearing in the axial-anomaly equation ∂ µ J µ 5 R = a s n f T F FF R between the vacuum and a pair of (off-shell) gluons to 4-loop order. Furthermore, we have verified explicitly up to 4-loop order the equality between the MS renormalization constant Z FF associated with the operator FF R and that of α s . This equality automatically ensures that the relation γ J = a s n f T F γ F J remains valid to O(α 4 s ) given the correctness of the all-order axial-anomaly equation which was used to determine the non-MS piece of Z J in Larin's prescription of γ 5 .
Note added. While this work was under review, a proof of the absence of a (divergent) multiplicative renormalization of the axial-anomaly operator in dimensionally regularized QCD was recently completed in ref. [69] in Becchi-Rouet-Stora quantization [70,71]. In a private communication, the authors of ref. [69] pointed to us that the proof in ref. [35] is incomplete.

JHEP05(2021)087
Forschungsgemeinschaft under grant 396021762 -TRR 257. The work of T.A. received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 772099. We gratefully acknowledge the computing resources provided by Max-Planck-Institute for Physics and by RWTH Aachen University where some of the computations have been completed.

JHEP05(2021)087
With all UV divergences in Γ µµ 1 µ 2 lhs (p 1 , p 2 ) regularized dimensionally, which can all be factorized and subtracted by momentum/mass-independent MS renormalization constants, the coefficients of its Laurent expansion series in are expected to behave regularly, i.e. not be divergent, when q → 0 with p 2 1 fixed at an off-shell point (as no IR divergence is present here). Under this regularity condition on Γ µµ 1 µ 2 lhs (p 1 , p 2 ) regarding its asymptotic behavior at q → 0, we see that the projector P µµ 1 µ 2 in eq. (3.2) should project out the same object as the operation in eq. (A.1). The fact that we obtained the sameM lhs to 3-loop order and M rhs to 2-loop order (with full ξ dependence) as those given in ref. [18], as discussed in section 5, indirectly confirms this point.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.