Hadronic vacuum polarization using gradient flow

The gradient-flow operator product expansion for QCD current correlators including operators up to mass dimension four is calculated through NNLO. This paves an alternative way for efficient lattice evaluations of hadronic vacuum polarization functions. In addition, flow-time evolution equations for flowed composite operators are derived. Their explicit form for the non-trivial dimension-four operators of QCD is given through order $\alpha_s^3$.


Introduction
The vacuum polarization functions (VPFs) for (axial-)vector and (pseudo-)scalar particles are among the most important objects when studying QCD. On the one hand, this is because their imaginary part is directly related to physical observables such as the decay rates of the Z-or the Higgs boson, or the hadronic R-ratio. If the characteristic energy scale is far above the QCD scale, a perturbative evaluation of the polarization functions is sufficient in these cases to arrive at high-precision results (see, e.g., Ref. [1]).
But VPFs also contribute indirectly to physical observables such as anomalous magnetic moments [2,3], the definition of short-distance quark masses [4], or hadronic contributions to the QED coupling [5,6]. These applications involve an integration of the VPFs over the non-perturbative regime, which is typically achieved with the help of experimental data and dispersion relations. Only very recently, first-principle lattice calculations have become competitive with these dispersive approaches. In the case of the hadronic vacuum polarization contribution to the muon's anomalous magnetic moment, the two approaches turn out to lead to incompatible results [7] 1 . It would therefore be highly desirable to have additional independent first-principle calculations of the VPF. About ten years ago, the gradient-flow formalism (GFF) was suggested as a mechanism to improve the efficiency of lattice calculations [10][11][12] (see also Refs. [13,14]). Since then, it has become a standard for the scale-setting procedure [15,16]. However, also other applications of the GFF have been studied, among them a new way to determine the energy-momentum tensor on the lattice. The underlying idea in this case is the smallflow-time expansion of composite operators [11], leading to the flowed Operator Product Expansion (OPE) (also named smeared OPE in Ref. [17,18]), where the regular operators are replaced by operators taken at finite flow time. Its main advantages with respect to (w.r.t.) the regular OPE is the absence of operator mixing, and the improved efficiency of the evaluation of operator matrix elements. The translation of the regular to the flowed operators can be done perturbatively. For the energy-momentum tensor, it is available through next-to-next-to-leading order (NNLO) [19][20][21]. Quite recently, the small-flow-time expansion was applied at next-to-leading order (NLO) to CP-violating operators [22], and to four-quark operators [23].
In this paper, we present the flowed OPE for the time-ordered product of two currents through NNLO QCD. Taking the vacuum expectation value (VEV) leads to the VPF. This should thus allow for an alternative first-principle evaluation of VPFs on the lattice. In addition, we derive a general logarithmic flow-time evolution equation for flowed operators which resembles the renormalization group (RG) equation of regular operators.
The remainder of this paper is organized as follows. Section 2 introduces the regular OPE of current correlators with operators up to mass dimension four. This includes the renormalization of these operators as well as an overview of the literature which provides the corresponding perturbative Wilson coefficients. (The coefficients for the dimensionfour operators for various currents are reproduced in Appendix B.) The transition to the flowed OPE is presented in Sect. 3. Section 4 describes the calculation of the mixing matrix between regular and flowed operators in the small-flow-time limit. While a large part of this mixing matrix is already known [21] and recollected in Appendix C, the missing components require higher order mass terms of the VEV of the flowed dimension-four operators and their renormalization with the help of the vacuum-energy renormalization constant. These results complete the ingredients required for the flowed OPE of the VPF through NNLO. In Sect. 5, we derive a logarithmic evolution equation from the generic flow-time dependence of the mixing matrix. Section 6 presents our conclusions and gives a short outlook on possible extensions of this work.

Current-current correlators
Our results are presented for a general non-Abelian gauge theory based on a simple compact Lie group with n f quark fields ψ 1 , . . . , ψ n f in the fundamental representation, of which the first n h are degenerate with mass m, while the remaining n l are massless. The generators T a of the fundamental representation are normalized as Tr(T a T b ) = −T R δ ab , and the structure constants f abc are defined through the Lie algebra [T a , T b ] = f abc T c . The dimensions of the fundamental and the adjoint representation are n c and n A , respectively, and their quadratic Casimir eigenvalues are denoted by C F and C A . For SU(N ), it is and QCD is recovered for T R = 1/2 and N = 3, i.e. C F = 4/3 and C A = 3. For brevity, we often use "QCD" also to refer to the more general gauge theory in the following.

Operator product expansion
The role of the perturbative and the non-perturbative regime of VPFs can be made most explicit through the OPE (see, e.g., Ref. [24]): where j(x) generically stands for a scalar, pseudo-scalar, vector, axial-vector, or tensor current, and k labels the mass dimension. Fig. 1 shows sample Feynman diagrams which arise from the perturbative evaluation of the current correlator in Eq. (2.2). In the following, we only consider the so-called non-singlet diagrams, where the currents are connected by a common quark line. An example for a singlet-diagram, on the other hand, is shown in Fig. 1 (e).
The coefficients C (k),B n on the right-hand side of Eq. (2.2) depend on the quantum numbers of the current and may thus carry Lorentz indices. Apart from the explicit results for specific currents in Appendix B, we suppress these indices throughout the paper. We furthermore assume that, upon transition from the left-to the right-hand side, possible global divergences are subtracted off of T (Q). Figure 1: Sample diagrams contributing to the perturbative calculation of the VPF, i.e., the left-hand side of Eq. (2.2). The currents are symbolized by wavy lines, gluons by spirals. We consider the case where n l quarks are massless (thin straight lines), and n h quarks are degenerate with mass m (thick straight lines).
(a) One-loop contribution for non-diagonal currents; (b-e) sample three-loop diagrams for diagonal currents. In (d), the currents couple to massless quarks. (e) is a "singlet" diagram. All diagrams in this paper were produced with the help of FeynGame [25].
Up to mass dimension two, the only operators of QCD which contribute to physical matrix elements are proportional to unity, i.e., where m B is the bare mass of the n h degenerate massive quarks. This means that are ultra-violet (UV)-finite, where Z m is the renormalization constant of the quark mass defined in Appendix A.
At mass dimension four, we choose the following basis of operators (the space-time argument is suppressed in most of what follows): where with the regular (as opposed to "flowed", see Sect. 3) quark and gluon fields ψ q (x) and A a µ (x), respectively, and the bare coupling g B . We employ Euclidean space-time, but translation of the intermediate formulas and the final results to Minkowski space is possible without difficulty. Working in space-time dimensions, the mass dimensions of O 1 and O 2 are actually equal to d, while that of O 3 is equal to 4. Higher dimensional operators are neglected in the following.
The set in Eq. (2.5) does not contain gauge dependent operators, or operators that vanish due to the equations of motion when acting on physical states, since they are irrelevant for the scope of this paper. In fact, in this respect the upper limit of the sum over q in O 2 could be replaced by n h , because the terms with massless quarks vanish on-shell. For the same reason, one could use O 2 ≡ −2m B n h q=1ψ q ψ q instead of O 2 in the definition of the operator basis (2.5). Other choices are possible as well, but the basis in Eq. (2.5) is particularly suitable for our purposes, because it is most directly related to the operator basis used in Refs. [20,21,26].
Matrix elements of the dimension-four operators are divergent in general. However, one may define "renormalized operators" O R n as linear combinations among them, for which physical matrix elements become finite: Analogously, one defines renormalized coefficient functions through the condition where C B n ≡ C (4),B n , cf. Eqs. (2.2) and (2.5). It is well known that, since the operators of Eq. (2.5) are part of the QCD Lagrangian, the renormalization matrix Z can be expressed in terms of the anomalous dimensions of QCD [27,28]: (2.10) The 't Hooft mass µ ensures that each renormalized operator O R n in Eq. (2.8) has the same mass dimension as the corresponding bare one, andμ appears because we will adopt the MS scheme by default (γ E = −Γ (1) = 0.577216 . . .). We also introduced the quantity a s = α s /π = g 2 /(4π 2 ) here, where g is the renormalized strong coupling in the MS scheme. Z 0 is the MS renormalization constant for the vacuum energy [28]. It is given in Appendix A, together with the anomalous quark mass dimension γ m and the d-dimensional beta function β .
Since the dimension-zero and -two operators in Eq. (2.2) are proportional to unity, the coefficients C (0) and C (2) are immediately determined from the small-mass expansion of these perturbative results for the VPF. They are thus known up to the four-loop level at the moment [48][49][50][51][52][53][54][55][56]. 2 The Wilson coefficients C n of the dimension-four operators, on the other hand, require a dedicated calculation which keeps track of the contributions from the individual operators. This has been done through O(a 3 s ) for C 1 and C 2 , and through O(a 2 s ) for C 3 in Refs. [61][62][63][64]. For the purpose of this paper, only the O(a 2 s ) results are required. For completeness, we include them in Appendix B.

Flowed operator product expansion
Having introduced the setup in the "regular" theory, we now translate this to the flowed OPE for the current correlators.

Flowed operators
We introduce the flowed operators as The flowed gauge and quark fields B a µ = B a µ (x, t) and χ q = χ q (x, t) obey the equations [10,12] Here we used the flowed covariant derivative in the adjoint representation, 5) and the flowed Laplace operators Z χ is the non-minimal renormalization constant for the flowed quark fields χ q , defined by the all-order condition [20] where · denotes the VEV. It readsZ where Z χ is the MS part, and (3.10) The short-hand notation reflects our choice of the "central" renormalization scale µ t [21].
Eq. (3.14) displays only the first six leading digits in numerical results. Results with higher accuracy are provided in an ancillary file, which also includes the L µt terms, see Appendix D. We expect that these floating-point numbers can be considered equivalent to their exact results for all practical purposes. This is why we often use the numerical values for the coefficients in what follows, even if the exact result is available.
Similar to the regular operators in Eq. (2.5), one could trade the flowed operatorÕ 2 (t) forÕ 2 (t) = −2mZ χ n f q=1χ q (t)χ q (t). However, in this case the final results to be derived below would be different, because the equations of motion for the flowed operators relatẽ O 2 (t) to bothÕ 2 (t) andÕ 1 (t) (see Refs. [20,21]). A transformation of the results in this paper toÕ 2 (t) is straightforward though.

Small-flow-time expansion
The small-flow-time expansion allows us to relate the QCD operators and coefficients with the flowed quantities as follows: where the ellipsis denotes terms that vanish as t → 0, and are the renormalized, finite mixing coefficients. Inversion of Eq. (3.15) gives where ζ −1 nk is the nk-element of the inverse of the mixing matrix ζ. This lets one define the "flowed OPE" for the current correlator: where the corresponding coefficient functions are related to the regular Wilson coefficients throughC The regular QCD coefficients C (0) and C (2) are given by the first two terms in m 2 /Q 2 of the large-Q 2 expansion of the VPFs. Through the required order, they can be found in Ref. [52] for vector-, in Ref. [53] for axial-vector-, and in Ref. [54] for scalar-and pseudoscalar currents, for example. The dimension-four coefficients can be found in Refs. [61,64]. For convenience of the reader, they are also collected in Appendix B.

Calculation of the mixing matrix
We now determine the mixing matrix ζ in a perturbative calculation through NNLO. By using the known results for the regular Wilson coefficients given in Appendix B, one can determine the flowed coefficients of Eq. (3.19) to the same order. Together with an evaluation of the flowed operator matrix elements on the lattice, the VPFs can be extracted and used in the determination of various physical quantities.
The bare mixing matrix ζ B can be determined with the help of the method of projectors: where the action of P (n) is to take suitable derivatives of a specific Green's function of the operator such that for n, m ∈ {0, 2} and k, l ∈ {1, 2, 3}. For details, see Refs. [21,66,67].
Specifically, the projectors onto 1, m 2 B 1, and O B 3 are given by derivatives of vacuum matrix elements w.r.t. m B : .

(4.3)
A crucial point of the method of projectors is that the derivatives and limits are taken before loop integration [66,67]. This guarantees that all loop contributions of projections on the r.h.s. of Eq. (3.15) lead to scaleless integrals and thus vanish in dimensional regularization. On the other hand, it implies that the projections on the l.h.s. of Eq. (3.15) can be divergent, even though physical matrix elements ofÕ n (t) are finite. This is why we need to carefully account for a possible non-integer mass dimension of these operators, see Eq. (3.1).
We directly obtain where the third set of equations follows fromÕ 3 (t) = m 4 = O 3 and the projector property P . The bare and renormalized mixing matrices for the dimension-four operators thus take the form where 0 T = (0, 0), and (4.6) ζ 2×2 can be obtained from the mixing matrix of the operators occurring in the energymomentum tensor and is accordingly known through NNLO [21]. Explicit results are given in Appendix C.
The coefficient ζ whereā s = a s (µ t ). Due to the RG invariance of a sÕ1 (t) [10,11], the result for general values of the 't Hooft mass µ can be obtained by multiplying the result given in Eq. (4.7) byā s /a s , replacinḡ Therefore, the only coefficients that are not yet known through NNLO are 2 , ζ 13 and ζ 23 . (4.10) According to Eq. (4.3), they require the calculation of derivatives w.r.t. m B in Õ 1 (t) and Õ 2 (t) up to three loops (see Fig. 2 for a set of associated diagrams). We achieved this using the setup described in Ref. [65], which employs qgraf [69,70] and q2e/exp [71,72] for the generation and subsequent categorization of the Feynman diagrams, FORM [73,74] for the manipulation of the ensuing algebraic expressions, the color package [75] of FORM for the calculation of the gauge group factors, and Kira⊕FireFly [76][77][78][79][80] for the Feynman integral reduction using integration-by-parts identities and the Laporta algorithm [44][45][46] over finite fields [81][82][83][84]. For the evaluation of the master integrals, we adopt the method described in Ref. [68], which performs sector decomposition [85] with the help of FIESTA [86,87] in order to extract the UV poles, along with a fully symmetric integration rule of order 13 for the numerical evaluation of their coefficients [88], implemented with high precision arithmetics by using the MPFR library [89]. Some intermediate steps of the calculation are done within Mathematica [90].
Multiplying the result by Z 2 m suffices to obtain the renormalized expressions for ζ (2) n , and we find, setting µ = µ t , ζ (m ∈ {0, 2}) and ζ n3 , for n = 1 (a-d) and n = 2 (e-h). The notation is the same as in Fig. 1; in addition, white circles denote flow-vertices, and lines with arrows next to them denote flow-lines (see Ref. [65] for details). Diagrams (a) and (b) only contribute to ζ where againā s = a s (µ t ). Since quark loops appear in Õ 1 (t) only at the two-loop level, ζ 1 (t) starts at O(a s ). In the case of ζ 1 (t) one needs to multiply by Eq. (4.12), and in addition byā s /a s . ζ 13 and ζ 23 require the more sophisticated renormalization given in Eq. (4.6). It is important here to work consistently in d space-time dimensions. Since Z 3 contains a 1/ pole already at O(a 0 s ), we need to keep the O( 2 ) terms of ζ 2×2 at NLO, and the O( ) terms at NNLO. They were not required in the calculation of Ref. [21], so we recalculated ζ 2×2 , keeping these higher terms in . Using the identity n A T R = n c C F [75], our final result for ζ 3 reads: The logarithmic terms at O(a n s ) are determined by the RG equation derived in Sect. 5. Nevertheless, for the convenience of the reader, we provide the result for general µ in this case.
This, together with the results for ζ 2×2 obtained in Refs. [21] and explicitly given in Eq. (C.5), completes the result for the small-flow-time coefficients of the OPE up to dimension four of Eq. (3.18).

Flow-time evolution
In the final section of this paper we derive a general flow-time evolution equation for flowed operators. It resembles the RG equation for regular operators but with a "flowed anomalous dimension matrix". While studies of the relation between the RG and the flowtime evolution have also been performed elsewhere in the literature (see, e.g., Refs. [26,[91][92][93][94]), to our knowledge the treatment described here has not been discussed before.
Let us return to the small-flow-time expansion of the operatorsŌ defined in Eqs. (3.15), (3.17), employing a matrix rather than component-wise notation for the sake of clarity: Since we work in the small-flow-time limit, the dependence of ζ(t) on t can be only through L µt , defined in Eq.
So far the discussion is general and holds for any flowed OPE. Specializing to our case of the QCD dimension-four operators, we can write the "flowed anomalous dimension" matrix as

(5.4)
Through O(a 2 s ), the result can be directly evaluated from Eqs. (C.5) and (4.13). A consistency check is obtained by noting that ζ(t) depends on t only through L µt : On the other hand, we know that a sÕ1 (t) andÕ 2 (t) are RG invariant [10,11,20] and therefore, with Eq. (3.15), where Since operators of different mass dimensions do not mix under RG evolution and ζ (0,2) 3 (t) = 0, we can drop the first two terms in the brackets on the r.h.s. of Eq. (5.6). 5 We thus arrive at where γ O is the anomalous dimension of the operators O R , defined through It can be written as   [21]. Thus, for γ f 2×2 , Eq. (5.5) is not just a consistency check, but a means to derive higher order terms. In our case, we can obtain the result through O(a 3 s ): We verified that this agrees through O(a 2 s ) with the result which is obtained by directly inserting Eq. (C.5) into Eq. (5.4). Due to the factor of 1/g 2 inÕ 1 (see Eq.
It may be useful to note that, by subtracting the VEVs off ofÕ 1 andÕ 2 ,Õ We checked, of course, that Eq. (5.8) is consistent with the results for ζ 13 and ζ 23 of Eq. (4.13).

Conclusions
We presented the flowed OPE for general current correlators and its matching to regular QCD through NNLO in the strong coupling α s and through mass dimension four by using the small-flow-time expansion. Our calculation is based on the renormalization procedure for the regular QCD dimension-four operators worked out in Ref. [27,28], the mixing matrix between flowed and regular operators derived in Ref. [21], the method of projectors [66], and the tools and results for perturbative calculations in the GFF presented in Ref. [65].
Overall, our results allow to combine the known perturbative results for the regular QCD current correlators from the literature to gradient-flow lattice calculations. This lays out the path for an alternative determination of hadronic contributions to observables such as the anomalous magnetic moment of the muon. In addition, we derived a general logarithmic flow-time evolution equation for flowed operators and presented its explicit form for the dimension-four operators considered in this paper.
Our methods are sufficiently general to be applied to similar problems at higher orders in perturbation theory, such as CP violating operators [22] relevant for the electric dipole moment of the neutron, or four-quark operators occurring in flavor physics [23].

A Renormalization group functions
The d-dimensional beta function is defined as where a s ≡ α s /π ≡ g 2 /(4π 2 ). The renormalized coupling g = g(µ) is related to the bare one through g B =μ Z Through Sect. 5, we only need the first two perturbative coefficients, while β 2 is required in order to derive the O(a 3 s ) terms of γ f in Eq. (5.12): The anomalous dimension of the quark mass is defined through with the first three perturbative coefficients given by Similarly to β , the third coefficient γ m,2 is needed only in Sect. 5.
The renormalization constant of the vacuum energy Z 0 is related to the corresponding anomalous dimension γ 0 through (A.8) The first three perturbative coefficients are given by [28,62]

B Perturbative coefficient functions
This appendix cites the results for the coefficient functions C n of the regular dimension-four operators appearing in the OPE of the current correlators defined in Eq. (2.2). We consider scalar, pseudo-scalar, vector-and axial-vector currents, both diagonal and non-diagonal, i.e. the currents assume the form This means that ψ k and ψ l can be either both massive with mass m (k, l ∈ M ), or both massless (k, l ∈ N ), or one of them is massless, the other massive (e.g. k ∈ M , l ∈ N ). While C 1 is independent of k and l, the coefficient C 2 of the quark operator takes the form where δ kM = 1 for k ∈ M , and 0 otherwise. Also the results for C 3 depend on whether the quarks k and l are massive or not. This dependence will be indicated explicitly below, using the δ kM symbol defined above.
For convenience, we introduce the short-hand notation and the dimensionless quantitieŝ The extra factor (−2) between C 2 andĈ 2 arises from using O 2 in Ref. [64] rather than O 2 from Eq. (2.5). For the sake of brevity, we insert the SU(3) color factors.

B.1 Vector and axial-vector currents
In this case, the current correlator can be decomposed into a transversal and a longitudinal part, each associated with a set of coefficient functions: The upper sign refers to the vector, the lower sign to the axial-vector case. The results have been taken from Ref. [64]: For the reader's convenience, we provide the relation between the mixing matrix ζ used in this paper and its definition in Ref. [21], referred to as "HKL" in what follows. This is easily derived from the relation between the operators O n defined in Eq. (2.5) and the O n,µν of HKL: with H 2×2 defined in Eq. (5.7). The mixing matrix between regular and flowed operators in HKL, restricted to the two operators which are relevant for this paper, is defined through   with ζ χ (t) from Eq. (3.10), gives the renormalized mixing matrix used in the current paper in terms of the bare mixing matrix of HKL:  where c (2) χ is given in Eq. (3.13). The results including higher orders in are provided in the ancillary file to this paper, see Appendix D.

D Ancillary File
The main results of this paper are provided in computer readable format (e.g. with Mathematica [90]). The notation is described in Table 1. All coefficients are represented by floating-point numbers in this file. The relative uncertainty for our numerically evaluated coefficients is estimated to be 10 −10 or better.