The complete singlet contribution to the massless quark form factor at three loops in QCD

It is well known that the effect of top quark loop corrections in the axial part of quark form factors (FF) does not decouple in the large top mass or low energy limit due to the presence of the axial-anomaly type diagrams. The top-loop induced singlet-type contribution should be included in addition to the purely massless result for quark FFs when applied to physics in the low energy region, both for the non-decoupling mass logarithms and for an appropriate renormalization scale dependence. In this work, we have numerically computed the so-called singlet contribution to quark FFs with the exact top quark mass dependence over the full kinematic range. We discuss in detail the renormalization formulae of the individual subsets of the singlet contribution to an axial quark FF with a particular flavor, as well as the renormalization group equations that govern their individual scale dependence. Finally we have extracted the 3-loop Wilson coefficient in the low energy effective Lagrangian, renormalized in a non-$\overline{\mathrm{MS}}$ scheme and constructed to encode the leading large mass approximation of our exact results for singlet quark FFs. We have also examined the accuracy of the approximation in the low energy region.


Introduction
The form factors (FF) of the vertices that couple an external color-neutral boson, such as a Higgs or an electroweak gauge boson, to a pair of quarks or gluons are important ingredients for calculating a number of phenomenologically interesting processes. The knowledge of high order perturbative corrections to these vertex FFs in Quantum Chromodynamics (QCD) is essential to make precision predictions for collider processes such as quark pair production in electron-position collisions, the Drell-Yan processes, hadronic production and decay of the Higgs boson and massive electroweak bosons. Furthermore, these vertex FFs constitute simple yet important objects of which the high-order QCD corrections can be used to extract certain universal QCD quantities of particular theoretical interest, such as the cusp anomalous dimensions [1,2] and the collinear quark and gluon anomalous dimensions (see, e.g., ref. [3][4][5][6][7][8][9]).
Due to the aforementioned importance, there has been a great amount of work on these objects in the literature. In this article, we are concerned with the so-called singlet type contribution to quark FFs describing the coupling of a pair of massless quarks to an external (axial) vector current to three loops in QCD including effects of a massive top quark. The QCD virtual corrections to quark FFs can be conveniently divided into two classes depending on whether the external color-neutral boson couples directly to the external quarks. Such a separation of QCD corrections is convenient for a number of practical reasons, such as allowing one to apply different γ 5 prescriptions [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] in the case of an axial current, as well as simplification of the calculation of loop integrals if the internal fermion loops are massive. Limited to purely massless QCD corrections, the 3loop results for vector [3,[26][27][28], scalar [29] and pseudo-scalar [30] part of quark FFs were derived in the literature, and there was recently important progress towards the 4-loop corrections to the vector quark FFs [8,9,31,32]. 1 The three-loop singlet contribution to the axial part of quark FFs in purely massless QCD was determined only very recently in ref. [42].
However, for physical application of the result for the axial quark FF, such as for theoretical predictions of the Z-mediated Drell-Yan processes to the third order in QCD coupling α s , it is necessary to incorporate the singlet QCD contribution with top quark loops, for at least two reasons both related to the presence of the axial-anomaly type diagrams [43,44]. First, in the absence of the top-loop contribution, the purely massless contribution to the axial FFs contains an explicit logarithmic renormalization scale dependence beyond that expected from the running of the MS renormalized α s , which is related to the non-vanishing anomalous dimension of the singlet axial current (e.g., determined in refs. [21,22,45]). Second, as well known in refs. [46][47][48][49][50], the top quark contribution to the axial FFs does not actually decouple in the large top mass or low energy limit, in contrast to the case of vector FFs. In particular, the singlet-type QCD contribution to the inclusive Z boson decay rate has been investigated in detail in the large top mass limit to O(α 3 s ) [48][49][50] and to O(α 4 s ) [51], and was found to be considerable. Here we compute the exact top quark contribution to massless quark FFs to 3-loop order in QCD, especially for the axial part. Most of the top mass dependent master integrals involved can be mapped to those in the 3-loop Higgs-gluon FF determined in ref. [52], and the additional ones are computed analytically and verified numerically by the same technique through this work. With the complete singlet axial current renormalization constant determined to O(a 3 s ) in ref. [45], including the non-MS finite piece, we are able to properly treat the individual subsets of singlet diagrams separated according to the flavor of the internal quark coupled to the Z boson, namely each separated flavor subset is mathematically consistent on its own. Not only an interesting theoretical question on its own, the UV renormalization of the anomalous top-quark loop contribution to the axial quark form factor determines the structure of the non-decoupling mass logarithms as well. We will therefore discuss the relevant renormalization formula in detail later in the article. We note that as long as one is only concerned with the anomaly-free sum of all singlet-type QCD diagrams from each electroweak doublet, e.g., in refs. [47][48][49][50][51], it is not necessary to include the non-MS part in the renormalization of the (singlet) axial current. The result presented here provides one of the missing ingredients needed to push the theoretical predictions of Z-mediated Drell-Yan processes to the third order in α s , such as done recently for those mediated by a virtual photon [53,54] or a W boson [55].
The article is organized as follows. In the next section, we introduce our conventions and notations for the quark FFs, and subsequently discuss the technicalities of their computation in section 3. In section 4 we discuss in detail the ultraviolet (UV) renormalization formulae for the individual subsets of singlet contributions, as well as the renormalizationgroup (RG) equations that govern the scale dependence of their finite remainders defined after subtraction of infrared (IR) divergences. We then present our exact numerical results for these finite remainders in section 5, and examine the quality of the large mass expansion results. In section 6 we extract the Wilson coefficient in front of the axial current of massless quarks in the low energy effective QCD, which can be conveniently used to approximate the leading behavior of the singlet contribution in the large top mass limit. We conclude in section 7.

Preliminaries
We consider the one-particle-irreducible corrections to the 3-point vertex function of an external (off-shell) Z boson and a pair of massless quarks of flavor q with on-shell outgoing momenta p 1 and p 2 , in QCD with n f = n l + 1 = 6 flavors and only the top quark kept massive. This vertex function admits the following Lorentz tensor decomposition where δ ij denotes the color factor, and v q and a q are respectively the vector and axial vector couplings of the external quark q to the Z boson. In eq.(2.1), we have used the fact that limited to gauge interactions, there are only two Lorentz structures, one parity-even and the other parity-odd, sandwiched between the two on-shell massless spinors which are linearly independent in 4 dimensions. The Lorentz-invariant coefficients F V and F V are, respectively, the vector and axial FF of the massless quark q, which are functions of s = (p 1 + p 2 ) 2 = 2 p 1 · p 2 as well as the top mass m t (when the top quark loop contributes).
The normalization is such that the tree-level values of these FFs read (in 4 dimensions): The two Lorentz structures in eq.(2.1) are orthogonal to each other. Pulling out the color factor and putting v q = a q = 1, the dimensionless F V and F A can be projected out in the following way: with the Dirac algebra and trace Tr done in D = 4 − 2 dimensions, and the γ 5 treated as anticommuting in the second line. In our calculations of the singlet contribution to F A , we used a non-anticommuting γ 5 definition [10,12] in the variant as prescribed in refs. [21,22]. Notice that the same projection applies, even though the form of the projection has been determined assuming an anti-commuting γ 5 . In fact, as long as one is only concerned with the finite remainders of these FFs in 4 dimensions, to be discussed in the following, one can set the -parameter in these projectors to be 0 from the outset [56,57], which is what we actually did regarding the axial FF projector.
The QCD virtual corrections to F V and F A can be conveniently classified into two parts, the non-singlet and singlet part, depending on whether the external Z boson couples directly to the external quarks or not. For the sake of later convenience, we have pulled out the Z boson couplings from the respective singlet contribution. In the remainder of this article, we adopt the convention regarding the terminology for the singlet and non-singlet type QCD corrections to F V (A) where the classification is solely based on the topology of the contributing Feynman diagrams. The non-singlet QCD corrections have the Z boson coupled directly to the open fermion line of the external quark q, which starts from the tree level. It thus depends only on the electroweak coupling of the external quark q. With q massless and an anticommuting γ 5 (which is straightforward to apply here), one has F A ns = F V ns to all orders in QCD owing to chirality conservation.
On the other hand, the singlet contribution F internal quarks running in the loops, which are normalized w.r.t that of external quark q as defined in eq.(2.1). In the Standard Model, quarks in a weak doublet couple with opposite sign to the Z boson in the axial part of the neutral current, and hence axial contributions from doublets add up to zero in the massless limit for singlet diagrams. Therefore in the usual approximation taken here, the only non-zero axial contribution comes from the top-bottom doublet due to the large mass difference. Specifically, we denote with λ q ≡ a b aq equal to ±1 depending on whether the external q is an upper or lower quark (i.e., having the same weak isospin as the bottom quark). The full QCD corrections to F A s were determined to 2-loop order in refs. [35] for both massless and massive external quarks. Very recently the 3-loop bottom contribution F A s,b was derived in effective QCD with n l massless quarks in ref. [42]. In this work, we provide the result for the top-loop contribution F A s,t to 3-loop order with exact top mass dependence, and also the part F A s,b that contains top quark loops inserted through gluon self-energy corrections.
Concerning the vector part of the singlet contribution, F V s vanishes at the 2-loop order due to the same reason that underlies the Furry theorem, and starts to contribute only from the 3-loop order with sample diagrams shown in fig. 2. The leading 3-loop result is completely UV and IR finite, as computed in refs. [3,26,28] but with only massless quarks included. The previously-missing top-loop induced contribution, i.e., the diagram with thick lines in fig. 2, is computed in this work for completeness. It is known to be power suppressed in the low energy limit.

Calculation of the bare singlet form factors
The bare quark FFs can be expanded formally in the bare QCD coupling constantâ s ≡α s 4π : where the perturbative expansion starts from the 2-loop order. For the calculation of these QCD loop corrections, we use dimensional regularization [10,58] with a non-anticommuting γ 5 [10,12] in a variant as prescribed in refs. [21,22]. The techniques employed closely follow those used in the computation of Higgs-gluon FF in ref. [52]. Here we merely sketch the main points specific to the amplitude in question. Symbolic expressions of the contributing Feynman diagrams to the 3-loop order are generated by the C++ diagram generator DiaGen [59]. There are 2 at the 2-loop order, and 57 at the 3-loop order, for both F A s,b and F A s,t . There are 4 diagrams contributing to F A,3 s,b with top-quark loops (see fig. 7 for an example), while the remaining ones are completely massless and have been determined in ref. [42], which we include here only for completeness.
The loop integrals in all contributing diagrams are then reduced to a much smaller set of master integrals via Integration-By-Parts (IBP) identities [60,61] with the help of a C++ implementation of the Laporta algorithm [62]. All massless master integrals involved in F A s,b are known analytically [3,[26][27][28]. The master integrals in F A s,t can all be mapped to those solved numerically in ref. [52]. There are, however, 15 top mass dependent masters M i ( , s m 2 t ) (with i = 1, 2, · · · , 15) appearing in F A,3 s,b , which we have to solve in addition in the present work. They are, fortunately, simple enough to be solved analytically using the differential equation approach [63,64]. We first derive the system of first-order homogeneous linear differential equations in the variable x ≡ s m 2 t for these 15 masters: by IBP reducing their derivatives back to themselves. The coefficients A ij ( , x) are rational functions in x and by the virtue of IBP identities. After performing a change of variable 2 x = 2−y− 1 y , the differential equation system (3.2) is then fed to the package CANONICA [65,66], which finds an -form [67] together with the rational transformation of the basis of master integrals 3 . The letters involved in this -form differential equation are {y, y + 1, y − 1}, and it is then straightforward to read off the solutions in terms of Harmonic polylogarithms (HPL) [68] with certain integration constants to be fixed by boundary conditions. We determine the boundary conditions by computing these master integrals in the large mass limit, which is located at x → 0 (or y → 1). In this way, all 15 top mass dependent masters in F A,3 s,b are solved analytically in terms of HPLs, expanded in to the orders needed to obtain the F A, 3 s,b at O( 0 ). In the supplemental material associated with this article, we attach the analytical result for the UV renormalized top-quark loop contribution to F A,3 s,b .

The renormalization formulae and RG equations
With the bare results for F A s,b and F A s,t at hand, we are now ready to perform the UV renormalization and define the finite remainders after IR subtraction. Although the nonsinglet and singlet part of the axial quark FF contribute to physical observables in a coherent physically-indistinguishable way, they do depend on (potentially) different electroweak couplings. Therefore, as far as the QCD corrections are concerned, they can be treated independently. In particular, one can derive RG equations for them separately.
To set up the notations and conventions in use, let us start with the renormalization of the bare QCD couplingâ s ,â with the dimensionless renormalized coupling a s ≡ αs 4π = g 2 s 16π 2 . The bare couplingâ s has mass-dimension 2 as is exhibited on the r.h.s. of (4.1). We work 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 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 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. When discussing the extraction of the Wilson coefficient in the low energy effective theory in section 6, we will further perform an additional finite renormalization of a s to decouple the top quark effect in the gluon self-energy correction.
By the virtue of the multiplicative renormalizablity of the QCD Lagrangian and the definition of the renormalized axial currents with a non-anticommuting γ 5 as summarized in ref. [45], we derive the following renormalization formulae for the a b -and a t -dependent singlet contributions individually: where theâ s on the r.h.s will be substituted according to eq.(4.1) and the bare massm t is to be renormalized on-shell bym t = Z m m t . The dependence on s is suppressed in order not to overload the notations. It is understood that n f i=1 is a shorthand notation for i∈{u,··· ,t} , and similarly n l i=1 (to appear below) refers to i∈{u,··· ,b} . The Z 2 is the on-shell wavefunction renormalization constant of the external light quark, which differs from one due to the presence of massive top loops starting from 2-loop order. The Z s ≡ 1 n f Z S − Z ns is the difference between the usual singlet and non-singlet axial current renormalization constants determined in refs. [21,45], further normalized to the case of having just one single flavor, e.g., either the bottom or top quark, coupled to the Z boson. To be more specific, the constant Z S is defined in the following renormalized singlet axial current with a non-anticommuting γ 5 , and is given by the product of eq.(5.1) and eq.(5.4) of ref. [45]. And to distinguish the different definitions of the terminology "singlet" we have given it a capitalized subscript. The constant Z ns is the one needed to renormalize the non-singlet axial current J µ ns,5 = i with a non-anticommuting γ 5 and can be obtained by taking the product of eq.(8) and eq.(11) of ref. [21]. For reader's convenience, we reproduce here the result for Z s , which reads The definition of the quadratic Casimir color constants is as usual: Notice that Z ns is scale independent, since the renormalized non-singlet axial and vector current are related by multiplication of γ 5 and the vector current is not renormalized. It is worthy to emphasize that what appears in the r.h.s of eq. (4.6) is Z ns + n f Z s rather than just Z s . Furthermore, the F A ns (â s ,m t ) in eq.(4.3) must be computed using the same non-anticommuting γ 5 prescription as used in calculating the singlet FFs. Therefore, this renormalization formula shows that as soon as one applies a non-anticommuting γ 5 prescription in the calculation of an isolated individual subset of singlet diagrams, one is forced to have the knowledge of the purely non-singlet diagrams in the same convention for the sake of renormalization, albeit only up to an order less by two loops.
Note that our formula eq.(4.3) is more general than what is needed here: once expanded up to O(a 3 s ) in question, only actually contribute, because the singlet quantities, Z s and F A s , all start from O(a 2 s ). For the same reason, the on-shell Z 2 does not contribute neither in our 3-loop results. The remaining terms, which are quite non-trivial due to mixing with singlet diagrams featuring quarks of other flavors (with potentially different masses), should get involved but only starting from 4-loop order. We note also that the appearance of them t dependence in various pieces in eq.(4.3) starts typically from 2-loop order, as the external quark q is massless. Furthermore, the difference between the two subsets of singlet contributions in eq.(4.3), which is essentially eq.(2.4), requires only the non-singlet axial current renormalization: as expected.
In the case of the renormalization of the singlet contribution to the vector counterpart, F V s,f , it is well known that one needs only the overall (on-shell) wavefunction renormalization constant Z 2 in addition to the renormalization ofâ s and top mass. Furthermore, up to the O(a 3 s ) considered in the present work, it is completely UV (and IR) finite.
Starting with eq.(4.3) and noting the non-zero anomalous dimension of the singlet axial-current operator eq.(4.6), one can then derive the RG equations for the renormalized a b -and a t -dependent singlet contribution, which governs their respective µ-dependence. With the external quark field and top mass renormalized on-shell while the a s in the MS scheme, the RG equations read where all FFs on both sides are the UV renormalized ones. Again once expanded and truncated up to O(a 3 s ), only the term with F A ns (a s , m t , µ) in the r.h.s. contributes to the 3loop calculations considered in the present work. Just like the pure non-singlet contribution F A ns (a s , m t , µ), one sees that the "physical" combination F A s,b (a s , m t , µ) − F A s,t (a s , m t , µ) has a zero anomalous dimension as a direct consequence of eq.(4.7): which is also clear from eq.(4.8). Therefore, the net µ dependence in the a t -dependent singlet contribution F A s,t (a s , m t , µ) is necessary to cancel that of F A s,b (a s , m t , µ), such that the remaining explicit µ dependence is related to the MS renormalization of a s in the usual way. Based on this, one anticipates already that the top-quark contribution F A s,t (a s , m t , µ) cannot completely decouple in the naive sense in the large top mass or low energy limit, because the "massless" contribution F A s,b (a s , m t , µ) still has a non-zero anomalous dimension to be compensated.
The UV-renormalized singlet FFs still contain IR divergences, starting from 3-loop order, from exchange of virtual soft and/or collinear particles, regularized as poles in the dimensional regulator . By factorizing out the IR singularities, we define the following finite remainders of F A s,b(t) (a s , m t , µ): where the dependence on s is suppressed as before, and the I qq denotes the IR-singular factor determined in ref. [69]. For the present application, the I qq needed reads With this I qq operator, which has a vanishing anomalous dimension, the same form of the RG equations derived in eq.(4.8) simply carries over to the cases of finite remainders defined above, namely, where the F A ns and F A s,i denote the finite remainders of the corresponding UV-renormalized FFs defined in the same way as in eq.(4.10), in particular by the same I qq operator. Once again, the combination F A s (a s , m t , µ) ≡ F A s,b (a s , m t , µ) − F A s,t (a s , m t , µ) has no anomalous dimension, just like the pure non-singlet contribution F A ns (a s , m t , µ). The top mass dependence in various parts above typically enter starting from O(a 2 s ), except in F A s,b where the dependence starts from the 3-loop order (see, e.g. fig. 7). Since all IR-subtracted FFs involved in eq.(4.12) are finite in 4 dimensions, one can simply take the 4-dimensional expressions of the anomalous dimensions involved, in particular the γ s as defined in eq. (4.6).
Alternatively, one could also consider finite remainders defined in the MS factorization scheme [70] where the IR factors contain only poles. Denoting the MS IR factor in need as Z qq , and by F A s,b (a s , m t , µ) the MS counterpart of F A s,b (a s , m t , µ), one has where the transformation factor I qq Z qq is free of poles and one can take its 4-dimensional expression. To be more specific, to the perturbative order needed here, it reads in the 4-dimensional limit where the a s -coefficient is the O( 0 ) term of the a s -coefficient of I qq in eq.(4.11). The same transformation factor holds for F A s,t (a s , m t , µ) as well. While the RG equations of these MS finite remainders would develop additional terms in the r.h.s. compared to those of eq. (4.12), due to the non-zero anomalous dimension of Z qq . Because of this, below we would only analyse the RG equations (4.12) for the finite remainders defined in eq.(4.10) in 4 dimensions, while switching to the respective MS finite remainders is straightforward as explained in eq. (4.13).

Results for the finite remainders
With all necessary ingredients ready, we present in this section our numerical results for the finite remainders of singlet quark FFs with exact top mass dependence over the full kinematic range from the low-energy (large-m t ) limit s/m 2 t → 0 to the high-energy (smallm t ) limit s/m 2 t → ∞. First, we observe that with our bare results for the axial singlet FFs, all poles in cancel after performing UV renormalization and IR subtraction as explained in the previous section, which itself serves as a welcome check. The parts of the resulting finite remainders featuring explicit logarithms in µ can be revealed via solving the RG equations (4.12) perturbatively in a s . In particular, the µ logarithms in the three-loop coefficients of the finite remainders are entirely determined from the well-tested lower-order results. Our 2-loop results are cross-checked with those in ref. [35] and agreement is found. We therefore show below the numerical results for the 3-loop singlet FFs with µ fixed at s where the purely massless contributions evaluate to constants, serving as convenient reference points. Im -ℱ s,t  [42]. In terms of our notations introduced below in eq.(6.1), it reads C a = Re FĀ ,3 s,b (µ 2 = s, n l = 5) . This plot is made from a sample of these functions at about 2 × 10 5 points all evaluated with very high precision. In the high energy limit x → ∞, one expects the asymptotic relation F A s,t (a s , x) → F A s,b (a s , x), corresponding to the well known result that the total singlet contribution to axial quark FFs vanishes with 6 massless quarks. This is demonstrated in fig. 3 by both the blue curve representing the (normalized) real part, and the dashed brown curve representing the (normalized and sign-flipped) imaginary part, approaching 0 at x → ∞. This serves as a strong check of our result for F A,3 s,t (x), given the correctness of the analytically-known purely massless result. The dotted gray line represents the ratio between −Im FĀ ,3 s,b (x, n l = 5) ≈ 931.771 [42] and C a , which is about 5.318, included for reference purpose. The dashed green curve, standing for −Im F A,3 s,t (x) /C a , does not overlap with this in the high energy limit, and this is simply due to the fact that it has effectively n f = 6 massless quark loops in the gluon self-energy insertion while the reference line has n l = 5 by choice. Because of the same reason, the red curve does not really arrive at -1, deviating by about 1.5%, albeit almost invisible from the plot.
As one lowers the energy down around the top pair threshold at x = 4, one observes the typical behavior due to the Coulomb effect in the real and imaginary part of the top-loop contribution F A,3 s,t (x), of which the former still varies smoothly while the later experiences a (non-smooth) sharp turn. There is no actual divergence observed in this 3-loop result, because we have here at most one virtual gluon exchange between the virtual top pair.
Below the top pair threshold x < 4, one enters the domain where in principle the large top mass expansion can approximate the full result well, given power corrections included to sufficiently high orders. What is special here for the axial FF is that the top quark loop effect does not decouple in the large top mass or low energy limit [46][47][48][49][50]. This is reflected in the plot by both the red and blue curve soaring up as x → 0. 4 This is in contrast to the low energy behavior of the top singlet contribution to the vector counterpart shown in fig. 4, where they are clearly power-suppressed. We note, however, the imaginary part of the F A,3 s,t (x) is power-suppressed in the limit x → 0, while that of F A,3 s,b (x) still features a logarithmic enhancement.
Before we dive into the examination of the quality of the large top mass expansion in this region, let us remark that given the typical size of the normalized axial FF values plotted in fig. 3 over a wide range, without incorporating the coherent top-loop diagrams, the result for the singlet contribution to the axial part of the quark FF would be, in general, completely off. Furthermore, their role in getting the proper µ dependence of the total singlet contribution, as well as stabilizing it in a truncated perturbative result, is evident from the RG equations discussed in the previous section.
We now zoom into the low energy region, and examine the quality of the large top mass expansion. In this region it is more sensible to renormalize the perturbative coupling a s such that the top quark effect in the gluon's self-energy correction is decoupled. This is implemented in the form of a finite decoupling renormalization of a s = ζ αās with ζ α = 1 +ā s 2 3 ln µ 2 Re-expand the finite remainders in powers ofā s , the perturbative coupling in effective QCD with n l = n f − 1 massless quark flavors, and one has where we used symbols with a bar to denote the a s -decoupled counterparts of the quantities appearing in previous equations. Setting µ 2 = s, the perturbative expansion coefficients in eq.(5.1) become univariate functions of x = s/m 2 t , and each can be further expanded as a power-log asymptotic series in x in the limit x → 0. To be specific, theā 3 s coefficient admits the following expansion ansatz where the integer m is bounded within a n-dependent range, and the power n is truncated up to certain value in practical calculation. In particular, for the leading power approximation ofF A,3 s (x) to be plotted below, one keeps only the terms without powers in x. In fig. 5 and fig. 6, we show the comparison of the exact result forF A,3 s (x) with its leading and subleading large-mass approximation as function of x in the range (0, 1.33) (corresponding to √ s ∈ (0, 200) GeV if one sets m t = 173 GeV), respectively for the real and imaginary part. As clearly demonstrated in the lower panels, the accuracy of the leading large-mass approximation within this range is already quite good forF A,3 s (x), deviating from the exact result by at most 3%. Including the subleading power correction O(x 1 ) further reduces the deviation to be below 1% in most of the range covered, hardly visible in the plot. The curves representing our exact result for the a 3 s -coefficient F A,3 s (x) and its leading large-mass approximation are included for reference. In particular, one sees that the logarithmic enhancement in its imaginary part is removed by the a s -decoupling in eq.(5.1). Consequently in the leading power approximation Im F A,3 s (x) is just a constant, given precisely by the purely massless result, leading to the overlap between the solid cyan line overlaps with the dashed gray line.
In the supplemental material associated with this article, we attach the expanded result for the finite remainderF A,3 s (x) truncated to O(x 10 ), which is very compact. 5 For the x-range covered in the plot, the difference between this expanded and the exact result is   Figure 5. Comparison of the real part of the exact result forF A,3 s (x) defined in eq.(5.2) with its leading and subleading large mass approximation, as function of x = s/m 2 t in the range (0, 1.33). The FFs are normalized w.r.t the real part of the 5-flavor massless result C a ≈ 175.218 [42]. The curves representing the exact result for F A,3 s (x) introduced below eq.(4.12) and its leading largemass approximation are included for reference. The lower panel shows the ratios of the leading and subleading approximation to the exact result. A similar plot for the imaginary part is given in fig. 6. below 10 −5 . It is sufficient to approximate the full result at the level of one per-mille up to the point x = 3, corresponding to √ s ≈ 300 GeV with m t = 173 GeV, and below 3% up to x = 3.76. As discussed at the end of the section 4, it is straightforward to transform this finite remainder defined under the convention of eq.(4.10) to others in different IR subtraction schemes that one may use in physical applications.

The Wilson coefficient in the low-energy effective Lagrangian
As mentioned already in the introduction, it is well known [46][47][48][49][50] that the effect of top quark loops in axial FFs does not decouple in the large top mass or low energy limit due to the presence of the axial-anomaly type diagrams. In the large top mass limit, the appearance of these non-decoupled mass logarithms in the (IR-subtracted) finite remainders of axial FFs are accompanied by the dependence on the renormalization scale µ as described by the RG equations of individual singlet contributions discussed in section 4. Therefore, the top contribution F A s,t (a s , m t , µ) can not completely decouple in the naive sense in the  large top mass or low energy limit, because the "massless" contribution F A s,b (a s , m t , µ) still has a non-zero anomalous dimension to be compensated such that the total anomaly-free result has the expected µ dependence [46][47][48][49][50]. Once combined together in the form as in eq.(4.9), the explicit remaining µ dependence is again dictated by just the MS renormalization of the a s . However, the non-power-suppressed m t -logarithms would not drop and remain in the form of ln s m 2 t in the total result. Still, it is quite striking to observe that even the four m t -dependent 3-loop Feynman diagrams contributing toF A, 3 s,b (m t , µ), all with top loop insertion through the gluon selfenergy correction with an example drawn in fig. 7, generate non-power-suppressed m tlogarithms beyond those that would be removed by the usual decoupling renormalization of a s . To be more specific, . Similar non-decoupling terms were determined in refs. [48,50] for the singlet contribution to the inclusive Z boson decay rate in the large top mass limit. We note that only after explicitly decoupling the top loop effect from the a s renormalization, For the axial FF of the flavor-q quark, a coefficient function C w (ā s , µ/m t ) is introduced to encode all the remaining non-decoupling m t -logarithms appearing in the total "physical" (or non-anomalous) combination F A s,b (a s , m t , µ) − F A s,t (a s , m t , µ) in the following way: s,i (ā s , µ) where only the remaining non-decoupling m t -logarithms after performing a s -decoupling are absorbed into C w (ā s , µ/m t ), indicated by the expansion done in powers ofā s . By the virtue of RG invariance of the l.h.s. of eq.(6.2), one can then derive the following RG equation of C w (ā s , µ/m t ): with the aid of whereβ,γ s andγ S = n lγs are the counterparts of β, γ s , and γ S in effective QCD with n l = 5 massless quarks in 4 dimensions. Resummation of the non-decoupling m t -logarithms absorbed in C w (ā s , µ/m t ) can be performed by solving the RG eq. (6.3) with the 4-dimensional γ s (andγ S ). We note that the particular form of eq.(6.3) holds for C w (ā s , µ/m t ) defined with the aforementioned non-decoupling m t -logarithms in F A s,b (a s , m t , µ) included, rather than just from the a t -dependent contribution. For instance, only after this combination, will one observe the number of fermions appearing inγ s (andγ S ) becoming n l , as checked to O(ā 3 s ) explicitly below. Expanded perturbatively up to orderā 3 s considered in present work, eq.(6.3) reduces to a simpler form With our results for the finite remainders of the FFs, presented in the previous section determined in a non-MS renormalization scheme, we can extract the result for the coefficient C w (ā s , µ/m t ) defined in eq.(6.3). It reads We subsequently checked that its µ dependence, hence the structure of the non-decoupling m t -logarithms, satisfies eq.(6.3) or rather its expanded form eq.(6.5) truncated to O(ā 3 s ), which serves as a check of the calculation. Note that the µ-logarithm-independent terms of C w (ā s , µ/m t ) in eq.(6.6) can not be determined by the RG eq.(6.3) alone (with given γ s of course), but have to be extracted from explicit calculations. In general, these constant terms depend on the renormalization scheme in use. In fact, before being combined together, the coefficients of the remaining non-decoupling m t -logarithms from a b -and a tdependent parts in eq.(6.2) are respectively renormalization-scheme dependent starting from 3-loop order, e.g., the coefficient of L µ in eq.(6.1). The O(ā 2 s ) coefficient in eq.(6.6) agrees with refs. [71,72] where it was extracted from explicit calculations of different quantities determined in the same renormalization scheme as used here. The O(ā 3 s ) coefficient in eq.(6.6) is found in agreement with the expression given in a recent publication [73] where it was composed from the MS result in ref. [48] with the aid of ref. [45].
As usual, the large mass limit result parameterized in eq.(6.2) can be generated from, or rather encoded by, an effective Lagrangian with a certain set of relevant local composite operators. To this end, the first point to notice is that the only type of additional operators relevant in the leading large mass approximation for the amplitude in question is the singlet axial current operator [46][47][48] as introduced in eq.(4.4). The part with Z boson couplings in the resulting top-less effective Lagrangian with n l massless quarks, reads, in its renormalized form: a iψ B i γ µ γ 5 ψ B i + a bZsJ µ S,5 + a t C w (ā s , µ/m t ) Z ns + n lZs J µ S,5 (6.7) whereJ µ S,5 = n l i=1ψ B i γ µ γ 5 ψ B i with a non-anticommuting γ 5 , andZ ns ,Z s denote the counterparts of the axial-current renormalization constants (c.f. eq.(4.6)) in this effective QCD. The combination of the non-decoupling m t -logarithms from both a b -and a tdependent singlet contribution is accounted for by the term in the second line of eq.(6.7), appended to the usual n l = 5 purely massless QCD. Only the sum of the second term in the first line and the second line is non-anomalous, as can be checked directly with the aid of RG equations of pieces given previously. 6 We emphasize that the Wilson coefficient appearing in eq.(6.3) and eq.(6.7) is determined in a non-MS renormalization scheme, in particular for the individual axial currents. This is intertwined with the particular renormalized form of axial currents appearing in the low energy effective Lagrangian given in eq.(6.7). The µ-logarithm-independent terms of C w (ā s , µ/m t ) given in eq.(6.6) depend on the renormalization scheme in use. And the C w (ā s , µ/m t ) term so-defined is supposed to be used in calculations made with an effective top-less QCD in the standard way where the chiral Ward identity of theψ B b γ µ γ 5 ψ B b is properly restored. In particular, the top-loop induced contribution to singlet axial FFs determined in this work can be directly combined with the result for the bottom contribution FĀ s,b (ā s , µ) derived in ref. [42].

Conclusion
In this work we determined numerically the finite remainders of the singlet contribution to quark FFs with exact top mass dependence over the full kinematic range, both for the axial and the vector part. We have worked out the renormalization formulae and RG equations for the individual subsets of singlet contributions to the axial FFs, subsequently checked using our explicit results to the perturbative order considered in the present work.
Our numerical investigation shows that without incorporating the coherent top-loop diagrams, the result for the singlet contribution to the axial part of the quark FF would be in general completely off. Furthermore, their role in getting the appropriate scale dependence of the total singlet contribution, as well as stabilizing it in a truncated perturbative result, is evident from the RG equations discussed in this article.
A particular low energy effective Lagrangian is composed to encode the leading large top mass approximation of the full result for the axial quark FFs, with the Wilson coefficient defined and extracted in the non-MS renormalization scheme in use. We note that only with all non-decoupling m t -logarithms from both a t -and a b -dependent singlet diagrams (remaining after a s -decoupling) combined together, will then the Wilson coefficient C w (ā s , µ/m t ) so-defined obey a simple RG equation in effective QCD as presented in this work. The accuracy of the leading large top mass approximation, encoded by such a low energy effective Lagrangian, is examined and is shown to be quite good for the 3-loop coefficient of the finite remainder, deviating from the exact result by ∼ 3% for √ s < 200 GeV.
The result presented in this work provides one of the missing ingredients needed to push the theoretical predictions of Z-mediated Drell-Yan processes to the third order in QCD coupling, especially at the differential level.