Short-flow-time expansion of quark bilinears through next-to-next-to-leading order QCD

The gradient-flow formalism proves to be a useful tool in lattice calculations of quantum chromodynamics. For example, it can be used as a scheme to renormalize composite operators by inverting the short-flow-time expansion of the corresponding flowed operators. In this paper, we consider the short-flow-time expansion of five quark bilinear operators, the scalar, pseudoscalar, vector, axialvector, and tensor currents, and compute the matching coefficients through next-to-next-to-leading order QCD. Among other applications, our results constitute one ingredient for calculating bag parameters of mesons within the gradient-flow formalism on the lattice.


Introduction
The gradient-flow formalism (GFF) [1][2][3] extends the fields of QCD in terms of the flow time t and is meanwhile an established tool in lattice gauge theory calculations.Its main application up to now has been in the scale setting procedure, required to determine the lattice spacing in physical units [3,4], as a scheme for defining the strong coupling constant [3,[5][6][7][8][9][10], or simply as a smearing mechanism [1,3].However, it has been shown that the GFF has a much larger potential.One of the key elements for this is the shortflow-time expansion (SFTX), where composite operators of flowed fields are expressed in terms of a regular operator product expansion (OPE).By inversion of a complete basis of operators, this lets one express an effective Lagrangian in regular QCD in terms of flowed operators and corresponding flowed Wilson coefficients [11][12][13][14].Matrix elements of the former do not require renormalization [11,15,16] and are thus ideally suited to be computed on the lattice.The flowed Wilson coefficients, on the other hand, can be obtained from the regular MS results via suitable conversion factors, which can be calculated perturbatively.Obviously, the perturbative order of these conversion factors has to match the one of the regular MS Wilson coefficients.This is why, in many cases, next-to-leading order (NLO) results are not sufficient, but higher orders are required.
The feasibility of the above approach was demonstrated via a flowed formulation of the energy-momentum tensor in QCD [12,13,17] which was subsequently used to extract thermodynamical observables from the lattice [18][19][20][21][22][23][24][25][26][27].In this case, the coefficients of the regular operators are rational numbers.It was shown that the next-to-next-to-leading order (NNLO) corrections to the matching matrix, which determines the conversion to flowed operators, lead to a significant improvement in the extrapolation to the physical limit at t = 0 [23].
As is well known from regular perturbation theory, every additional order leads to an enormous increase in complexity.Fortunately, however, many of the tools and techniques from regular perturbation theory can be adapted to higher orders in the GFF.An outline of this strategy has been described in Ref. [28], where a number of three-loop quantities were evaluated at finite flow time.Using this approach allowed to extend the NLO results for the effective weak |∆F | = 2 Hamiltonian [29] or the magnetic dipole moment operator [30] to the NNLO in Ref. [31] and Ref. [32], respectively.Similarly, the matching matrix between flowed and regular operators and Wilson coefficients was obtained to the same order for the hadronic vacuum polarization [33].
One of the main benefits of the GFF is the exponential suppression of high-momentum modes.As already mentioned above, this implies that composite operators of flowed fields do not require renormalization.Matrix elements of operators which only involve flowed gluons are even finite after renormalization of the regular QCD parameters (strong coupling and masses).Flowed quark fields, on the other hand, still require multiplicative renormalization, typically denoted by Z χ , see Section 2.1.In order to match perturbative and lattice results, one needs to define a suitable renormalization scheme, most conveniently via a Green's function which involves two flowed quark fields.One option is the so-called ringed scheme, originally proposed in Ref. [13], which fixes Z χ via the tree-level vacuum expectation value of the quark kinetic operator.The conversion factor between the MS and the ringed scheme is known through NNLO [28].
However, other options for the scheme of Z χ may be more convenient.For example, one may fix it via the SFTX of some quark bilinear operator, usually called current.A preliminary lattice study of this strategy was recently presented in Ref. [34], where the flowed four-quark operator was normalized to the flowed axialvector current.Which current is most suitable may depend on the specific calculation or observable under consideration.It will therefore be useful to have all the associated results at disposal.Moreover, the simplicity of the quark currents could also be used for systematic studies of the SFTX.First, one can compare perturbative and nonperturbative determinations of the matching coefficients with each other.Some preliminary studies in this direction have already been carried out in Ref. [35] for the CP-violating quark chromoelectric dipole moment operator and for the currents in Ref. [36].Secondly, one can compare results for the renormalized currents obtained through the SFTX with results obtained in more conventional non-perturbative schemes.This may allow one to test non-perturbatively the accuracy of the SFTX and assess the systematics associated with the t → 0 limit.A first study of higher-power terms has been done in the context of the energy-momentum tensor in Ref. [26].Besides these indirect applications, the SFTX of the currents directly contribute to a number of observables in the GFF such as the chiral condensate [19] or semileptonic contributions to the neutron electric dipole moment [37].
Through NLO, the SFTX of the currents has been calculated already several years ago [38,39]. 1 In order to be consistent with the uncertainties expected from the associated lattice calculations, one can expect that higher orders of the matching coefficients will be relevant.In this paper, we will therefore derive the corresponding NNLO results.
The remainder of this paper is structured as follows: In Section 2, we discuss the theoretical basis of our calculation, starting with the GFF in Section 2.1, the definition of the regular and flowed currents in Section 2.2, and the methods to obtain the SFTX in Section 2.3.Our results for the matching coefficients are presented in Section 3. The latter also allow us to evaluate the so-called flowed anomalous dimensions, describing the logarithmic flowtime evolution of the currents.This is presented in Section 4. Section 5 contains our conclusions.

Theoretical framework 2.1 The QCD gradient flow
In this paper, we work in D-dimensional Euclidean space-time with D = 4 − 2ϵ.The GFF continues the gluon and quark fields A µ and ψ of regular 2 QCD to fields B µ (t) and χ(t) through the initial conditions and the flow equations [1,3,15] where the "flow time" t is a parameter of mass dimension [t] = −2, g B is the bare strong coupling, and κ is a gauge parameter which drops out of physical observables.In our calculations, we set κ = 1.
The flowed field-strength tensor is defined as the flowed covariant derivative in the adjoint representation is given by and with the flowed covariant derivative in the fundamental representation, The flow equations can be solved perturbatively, leading to generalized QCD Feynman rules which involve exponential factors for the quark and gluon propagators, plus additional "flow-lines" representing the evolution of the fields in the flow time.The latter couple to the quarks and gluons via "flowed vertices".The general formalism has been worked out in Refs.[11,15], and more details can be found in Ref. [28].
Since the flow time acts as regulator for ultraviolet divergences, the GFF improves the renormalization properties.After renormalization of the fundamental parameters of QCD, the flowed gluon field B µ (t) is finite and does not require field renormalization [3,11,16].The flowed quark fields χ(t), on the other hand, require multiplicative field renormalization [15].Throughout this paper, we will adopt the ringed scheme [13], where Both the MS expression Z MS χ as well as the finite conversion factor to the ringed scheme, ζ χ , are available through NNLO [13,15,17,28].Explicit expressions are collected in Appendix A.
Parameter and field renormalization are sufficient to render physical matrix elements of composite operators finite [11,15,16].Thus, composite flowed operators do not mix under renormalization, which enormously facilitates the lattice evaluation of their matrix elements compared to those of regular operators, in particular if the latter mix with operators of different mass dimension.Connection of the flowed matrix elements to physics at t = 0 can be made through a perturbative calculation, as will be explained in more detail in Section 2.3.

Quark currents
In this paper, we consider n f = n l + n h quark flavors, where n l is the number of massless quarks, while the remaining n h quarks have identical mass m.The bare non-diagonal and diagonal currents are defined as respectively, where p and q are flavor indices3 with p ̸ = q.We furthermore define the bare and renormalized singlet and non-singlet currents as where h a is a traceless flavor generator, and m B is the bare quark mass which is related to the MS renormalized mass m through with Z m ≡ Z MS m given in Eqs.(A.6) and (A.12).Z s , Z ns , and Z 1 are renormalization constants which will be specified below.They depend on the Dirac structure Γ i.e. we consider the parity-even scalar, vector, and tensor currents, and the parity-odd axialvector and pseudoscalar currents.The renormalization constant Z 1 allows for the possibility of the currents to mix with the unit operator.

Short-flow-time expansion
Returning to the non-diagonal and diagonal notation, we define the flowed currents as which means that we adopt the ringed scheme for the flowed quark renormalization throughout this paper, unless stated otherwise.The flowed singlet and non-singlet currents are defined by replacing the diagonal and non-diagonal currents by their flowed versions.
The SFTX [11] for the singlet and non-singlet current can be written as where we have taken into account that we only have n h massive quarks of equal mass m.The terms of order t will be neglected in this paper.
The dependence of the matching coefficients ζ(t) on the flow time t is logarithmic.They are most conveniently computed by defining projectors onto the regular-QCD (i.e.not flowed) diagonal and non-diagonal currents.In our case, we choose 3) [O] = 1 3! and where M pq is a two-quark Green's function defined as the trace is over spinor and color indices, and (2.17) The nullification of the masses and external momenta in Eqs.(2.14) and (2.15) is understood to be taken before any loop integral is evaluated.This means that only tree-level diagrams contribute when the projectors are applied to the r.h.s. of Eq. (2.13), because all higher-order diagrams are scaleless and thus vanish in dimensional regularization.
The matching coefficients are then obtained as where with p ̸ = q, and ζ for a massive quark flavor p.The renormalized matching coefficients follow from this by inserting Eq. (2.9) into Eq.(2.13): Figure 1: Sample diagrams contributing to ζ (1) and ζ (3) .Spiral lines are gluons, straight lines denote quarks; lines with accompanying arrows are the corresponding flow lines (they are always connected to flowed vertices, denoted by small circles).The two fermion lines in diagram (c) can be of different flavor.The vertex with the cross denotes the current.All Feynman diagrams in this paper have been drawn with FeynGame [40]. ζ B and ζ B are just given by the first two terms in an expansion in m 2 B t of the currents' vacuum expectation values.Due to Lorentz and parity invariance, they are non-zero only for the scalar current, corresponding to the so-called quark condensate.Since the one-loop contribution is of order g 0 B , it is required to three-loop order in order to obtain the SFTX up to NNLO QCD.Sample three-loop diagrams are shown in Fig. 1.
In contrast, the projections for the other matching coefficients in Eq. (2.20) require only two-loop calculations.In the diagrams that contribute to ζ ns B , the current is connected to the external states by a single quark line, cf.Fig. 2, as opposed to ζ ∆ B , where the current and the external states belong to different quark lines, cf.Fig. 3.
We will refer to the latter class as triangle diagrams in what follows.The triangle diagrams for the scalar, pseudoscalar, and tensor currents vanish after taking the fermion trace.For the vector current, they only start to contribute from three-loop order due to Furry's theorem.At the perturbative order considered here, we can therefore drop the superscripts "ns" and "s" in these cases and simply write and analogously for the bare matching coefficients.For the axial current, on the other hand, we will find ζ ∆ B ̸ = 0 at the two-loop level.We evaluate all diagrams in D = 4 − 2ϵ space-time dimensions.The occurrence of γ 5 in Eq. (2.11) causes the well-known complications which we take care of by following the strategy outlined in Refs.[41][42][43].This means to replace both in the currents of Eq. (2.12) as well as in the projectors of Eq. (2.15).The resulting products of two (intrinsically four-dimensional) ε tensors are replaced by Figure 2: Examples for contributions to the non-singlet matching coefficients.
The notation is the same as in Fig. 1.
Figure 3: Examples for contributions to the singlet matching coefficients.
The notation is the same as in Fig. 1.
where the square brackets denote the anti-symmetric combination, e.g.
This also affects the normalization factors of Eq. (2.17) via4

Tr( Γµ
This strategy violates the Ward identities, but they can be restored by an additional finite renormalization discussed below.
For the actual calculation, we adopt the framework developed in Ref. [28], which is based on qgraf [44,45] for the generation of the diagrams, q2e/exp [46,47] for inserting the Feynman rules and identifying the momentum topologies, in-house FORM [48][49][50] routines for performing various computer algebraic operations including Dirac and color algebra [51], and Kira⊗FireFly [52][53][54][55] for the reduction to master integrals employing integrationby-parts-like relations [28,56,57] and the Laporta algorithm [58].Up to two-loop level, we find the same master integrals as in Ref. [17].They can be evaluated analytically in terms of the transcendentals5 where Li n (z) = ∞ k=1 z k /k n is the polylogarithm of order n.The three-loop vacuum expectation value contributing to ζ (1) and ζ (3) of Eq. (2.21) leads to 304 diagrams which are reduced to 216 master integrals already contributing to Ref. [33].They have been evaluated numerically following the strategy described in Ref. [59].

Results
Parity-even currents.For the currents which do not involve γ 5 , we adopt the MS scheme, i.e., we define the renormalization constants of Eq. (2.8) as Because of Lorentz and parity invariance, only the scalar current can mix with the identity operator.Thus, we have where Z 0 is the renormalization constant of the vacuum energy which can be found in Eq. (A.14), and µ is the renormalization scale in the MS scheme, and γ E = −Γ ′ (1) = 0.577216 . . . is the Euler-Mascheroni constant, with Euler's gamma function Γ(z).The Ward-Takahashi identities ensure that Z MS V = 1 and Z MS S = Z m , with Z m the quark mass renormalization constant introduced above.For the tensor current, the renormalization constant is given in Eqs.(A.5) and (A.16).
We express our results in terms of the color factors where in QCD, the number of colors is n c = 3.Furthermore, we introduce where we have implicitly defined the t-dependent energy scale µ t , and with g the MS renormalized strong coupling, see Eq. (A.1).We then find the following matching coefficients for the parity-even currents: ) χ is associated with the ringed scheme and would not appear if the fermions were renormalized in the MS scheme.Its explicit form is given in Eqs.(A.24) and (A.25).The results for ζ S (t) and ζ S are already known from Refs.[15,28,30,32]. 6For the sake of brevity, we only quote the three-loop results with four significant digits here.In the ancillary file accompanying this paper, we provide the results with higher numerical accuracy, as well as analytic expressions for the logarithmic terms and the coefficient of S (t) (see Appendix B).Also, while all results in this paper are in the ringed scheme, the ancillary files also contain the matching coefficients in the MS scheme of the flowed fermions.
Parity-odd currents.For the parity-odd currents, one has to introduce a non-minimal renormalization in order to restore the associated Ward identities in regular QCD for the non-singlet cases, and the correct anomaly in the singlet axialvector case, which are broken by adopting Eq. (2.23) combined with minimal subtraction.Therefore, we define the renormalization constants of Eq. (2.8) in these cases as where the MS part is given in Eqs.(A.5) and (A.16).The Z 5,X are finite renormalization constants given in Eq. (A.18).For the non-singlet cases, taking them into account is actually equivalent to working with a naively anti-commuting γ 5 combined with minimal subtraction, which means that We explicitly verified these relations by using Eq.(2.23) and the corresponding renormalization constants of Eq. (3.11).This provides a strong validity check of our calculational setup.
For the renormalized triangle contribution of Eq. (2.21), on the other hand, we find Additional checks.In order to further corroborate the correctness of our results, we performed all calculations in general R ξ gauge and confirmed that the dependence on the gauge parameter drops out in the final result.The only exception to this is ζ (3) , for which the calculation in R ξ gauge exceeds our computing resources.Through NLO, the matching coefficients ζ X (t) were already computed in Refs.[38,39], and we find full agreement after fixing the erroneous finite renormalization7 for ζ P (t) and ζ S (t).For ζ T (t), only the bare NLO result is provided in Ref. [39], and it agrees with our bare result.
Numerical results for ζ V .In order to see the improvement of the impact of the NNLO terms, we display in Fig. 4 the result for ζ V as a function of the unphysical renormalization scale µ at leading order (LO), NLO, and NNLO.The flow-time t has been fixed at µ t = 3 GeV and µ t = 10 GeV, corresponding to t ≈ 0.003 GeV −2 and t = 0.03 GeV −2 , cf.Eq. (3.4).For the plot at µ t = 10 GeV, we set n f = 5 and use α GeV) = 0.1880 as input.For the plot at µ t = 3 GeV, we set n f = 3 and use α We then evaluate α s (µ) through one-, two-, and three-loop running for µ t /3 ≤ µ ≤ 3µ t and insert it into the LO-, NLO-, and NNLO-approximation of ζ V (t) (with the corresponding value for n f ) in order to obtain the three curves in the plots.
Since this quantity is RG invariant, we expect the µ dependence to decrease from NLO to NNLO, and this is indeed what we observe.Taking the variation of µ around µ t by a factor of two as an estimate of the perturbative uncertainty, we find that it decreases from 4.4% to 1.4% at µ t = 3 GeV, and from 1.8% to 0.4% at 10 GeV.This is indicated by the red and blue bands in the plot.Another important observation is that these bands overlap, indicating that µ t as defined in Eq. (3.4) is indeed a reasonable choice for the central renormalization scale.Since the other matching coefficients are not RG invariant, we refrain from showing the analogous plots for them.

Flowed anomalous dimension
As suggested in Ref. [33] for a general set of flowed operators O = (O 1 , . . ., O p ), one may define flowed anomalous dimensions which allow to resum their logarithmic t-dependence.Let us briefly recapitulate the idea behind it.First consider the flow time t and the renormalization scale µ as independent quantities.The regular operators O are then independent of t, the flowed operators Õ are independent of µ, and the elements of the matching matrix ζ are functions of a s (µ) and L µt .Therefore, neglecting terms that vanish as t → 0, and thus where the flowed anomalous dimension matrix γ is a function of a s (µ) and L µt , which is, however, formally independent of µ: Note that the derivative acting on ζ in Eq. (4.2) only affects the logarithmic terms L µt .The latter can be derived at higher orders by noting that flowed operators (in the ringed scheme) are µ-independent, i.e.
where γ is the anomalous dimension matrix of the regular operators O, i.e.
The term in brackets starts at O(a 2 s ), and thus the one-loop term of the flowed anomalous dimension is given by the (negative of the) regular anomalous dimension of the current.Furthermore, knowledge of ζ at order a n s is sufficient to obtain γ through order a n+1 s , given that γ is known to a n+1 s and β to a n s .
It may be interesting to note that Eq. (4.5) can also be derived by tying the flow time and the renormalization scale together from the start, i.e., setting µ = c/ √ t, with some constant c.In this case, also the regular operators become t dependent, while ζ depends on t only through a s (c/ √ t), and thus which again leads to Applied to the current operators considered in this paper, Eq. (4.2) reads and using Eqs.(4.5), (3.8) to (3.10), (A.16), and (A.17), we find   Note that, in these formulas, a s is still renormalized in the MS scheme, and we have set µ = µ t , see Eq. (3.4) (the expression for general µ can be easily reconstructed using Eq.(4.3); it is also given in the ancillary file accompanying this paper, see Appendix B).
In order to eliminate any reference to the MS scheme, one can simply convert a s in these expressions to the gradient-flow scheme according to [3] a s = âs 1 − e 1 âs + â2 s (2e (4.17) The exact expression for the coefficient of C A T R n f can be found in Ref. [28].
Numerical values for the flowed anomalous dimensions are shown in Fig. 5.The input parameters are the same as in Fig. 4. Also here, we observe the expected reduction of the renormalization scale dependence when including higher orders, albeit sometimes less pronounced than for ζ V (t).But also here the NLO and the NNLO uncertainty bands nicely overlap.

Conclusions
In this paper, we have considered the SFTX of the scalar, pseudoscalar, vector, axialvector, and tensor currents and computed the corresponding matching coefficients through NNLO in QCD.Possible applications of these results are the calculation of the chiral condensate on the lattice [19] or the semileptonic contributions to the neutron electric dipole moment [37].
Our results could also serve as alternatives to the ringed renormalization scheme [13], which requires the calculation of the vacuum expectation value of the quark kinetic operator.In certain cases, it may be more efficient to normalize quark matrix elements to one of the currents instead.We believe that this strategy will especially find applications in flavor physics, in particular in combination with the SFTX of the relevant four-quark operators [29,31].A first preliminary study, already employing the result for the axialvector current obtained in the present paper, has been published in Ref. [34].
Finally, the simplicity of the quark currents could also be advantageous for systematic studies of the SFTX, building up on the preliminary studies of Refs.[26,35,36]  The renormalization group equation for the currents is thus given by µ 2 d dµ 2 j(µ) = γ(a s ) + γ fin (a s ) j(µ) , (A.9) where γ fin arises from any finite renormalization as introduced for the parity-odd currents when adopting Eq. (2.23).Specifically, if In this paper, we need the MS quark mass renormalization constant Z m ≡ Z MS m through O(a 3 s ), given by as well as the renormalization constant of the vacuum energy Z 0 through O(a 2 s ).It is related to the corresponding anomalous dimension γ 0 through The first three perturbative coefficients are given by [60,61] γ 0,0 = 1 , γ 0,1 = C F , γ 0,2 = −C (A.17) Recall that the vector current does not require renormalization, and the scalar current renormalizes with Z S ≡ Z m .The finite renormalization constants for the axial and the pseudoscalar current introduced in Eq. (3.11) are given by [42,43]  The flowed-quark field renormalization constant introduced in Eq. (2.7) assumes the same form as Eq.(A.5) with the anomalous dimensions given by [15,17] γ χ,0 = −

Z 5 ,FF
P (a s ) = 1 − 2 a s C F + a T R n f + O(a 3 s ) , Z ns 5,A (a s ) = 1 − a s C F + a T R n f + O(a 3 s ) , Z s 5,A (a s ) = Z ns 5,A + 3 16 a 2 s C F T R n f + O(a 3 s ) .(A.18)Let us remark that quoting the results for Z ns 5,A and Z 5,P is redundant, because they could be derived from Eq. (A.11) andγ S = γ P + γ fin P , γ V = γ ns A + γ ns,fin A .(A.19) The matching coefficient ζ V (t) at two different values of t = e −γ E /(2µ 2 t ) as a function of µ/µ t .

Table 1 :
The expressions of the ancillary file that encode the main results of this paper.

Table 2 :
Notation for the variables in the ancillary file.