Results and techniques for higher order calculations within the gradient-flow formalism

We describe in detail the implementation of a systematic perturbative approach to observables in the QCD gradient-flow formalism. This includes a collection of all relevant Feynman rules of the five-dimensional field theory and the composite operators considered in this paper. Tools from standard perturbative calculations are used to obtain Green's functions at finite flow time $t$ at higher orders in perturbation theory. The three-loop results for the quark condensate at finite $t$ and the conversion factor for the"ringed"quark fields to the $\overline{\mbox{MS}}$ scheme are presented as applications. We also re-evaluate an earlier result for the three-loop gluon condensate, improving on its accuracy.


Introduction
The gradient flow in field theory has proven to be a very useful concept. Originally introduced as a method to regularize divergences and ultraviolet fluctuations in lattice calculations [1,2], it has mainly gained traction through its use for scale-setting in lattice QCD [3][4][5] and has meanwhile become an indispensable tool in many practical calculations in this field (see, e.g. Refs. [6][7][8][9]). While first focused on QCD, generalizations of the gradient-flow formalism (GFF) yield a much wider range of applications, for example the study of dualities in field theory [10][11][12][13][14][15].
In QCD, the formulation as a five-dimensional field theory has been presented in Ref. [16], where the additional dimension is associated with the so-called flow time. The form of the fundamental theory (i.e., actual QCD) serves as the boundary condition at vanishing flow time. It has been proven that, when expressed in terms of renormalized parameters and fields, composite operators at positive flow time are finite [16]. This property and the realization that the gradient flow can be used as a matching scheme between lattice and perturbative calculations shed light on a whole new set of applications that are currently being explored.
One of the first possible cross-fertilizations among perturbative and lattice QCD is given by the definition of a new scheme for the strong coupling, defined by the gluon condensate (QCD action density) at finite flow time [3]. Its perturbative relation to the MS coupling is known through three loops [17]. 1 Another cross-boundary application is the gradient-flow definition of the energy-momentum tensor, which uses the small-flow-time expansion in order to express it in terms of well-defined composite operators at positive flow time and perturbatively accessible coefficient functions [19][20][21][22]. Yet another example is the proposal to relate Euclidean quasi PDFs 2 on the lattice to perturbative light-front PDFs by using the gradient flow [23,24].
All these applications rely on input from the lattice as well as from perturbative calculations, both at finite flow time. It was found that higher order corrections in perturbation theory are crucial for reducing the perturbative truncation error as estimated by the variation of the renormalization scale [17,21]. Such corrections were also found to significantly stabilize the required extrapolations of the corresponding lattice results [22].
Despite the fact that loop integrals in the GFF involve additional exponential factors as well as integrations over flow-time variables, many important techniques for regular perturbative calculations retain their usefulness, albeit in slightly generalized form. It is the goal of this paper to provide the basis for further perturbative calculations within the GFF, and thus to contribute to the further exploration of the capabilities of this approach. Sect. 2 reviews the QCD gradient flow in perturbation theory, recapitulating results of Refs. [3,6,16]. The various stages of automation in a perturbative multi-loop calculation of correlation functions are described in Sect. 3: the generation of Feynman diagrams and insertion and evaluation of Feynman rules, followed by a reduction of the loop integrals to a set of master integrals, and finally the numerical evaluation of the latter. As an application, Sect. 4 presents the three-loop results for the gluon condensate, obtained before in Ref. [17], the quark condensate, as well as the "ringed" quark field renormalization constant, which are new results. These quantities allow one to define a precision gradient-flow coupling and gradient-flow mass scheme and their matching to MS as shown in Sect. 5. Our conclusions are presented in Sect. 6.

The QCD gradient flow in perturbation theory
For completeness, we collect the main steps of the original derivations of Refs. [3,16] in this section.

Definition of the gradient flow
In the following, we work in D-dimensional Euclidean space-time with D = 4 − 2 . The GFF continues the gluon and quark fields A a µ (x) and ψ i α (x) of regular 3 QCD to (D + 1)dimensional fields B a µ (t, x) and χ i α (t, x) through the boundary conditions and the flow equations [3,6] ∂ t B a µ = D ab where the "flow time" t is a parameter of mass dimension minus two, and κ is a gauge parameter which drops out of physical observables (see below).
The (D + 1)-dimensional field-strength tensor is defined as the covariant derivative in the adjoint representation is given by and with the covariant derivative in the fundamental representation, As usual, the color indices of the adjoint representation are denoted by a, b, c, . . ., while µ, ν, ρ, . . . are D-dimensional Lorentz indices. Color indices of the fundamental representation are denoted by i, j, k, . . ., but they are suppressed throughout this paper, unless required by clarity; similarly for spinor indices α, β, γ, . . .. The symmetry generators T a are understood in the fundamental representation. The obey the commutation relation with the structure constants f abc and the trace is normalized to with T R > 0.
Different choices of κ in Eq. (2) correspond to gauge transformations of the form where Of course, all observables are independent of the gauge parameter κ [3]. In perturbative calculations, it is usually most convenient to set κ = 1.
The flow equations (2) can be incorporated in a Lagrangian formalism by defining The first three terms constitute the regular Yang-Mills Lagrangian, with fermions added in the fundamental representation (quarks). Introducing an index f in order to distinguish different quark flavors of mass m f , the classical, gauge-fixing, and Faddeev-Popov ghost part are given by respectively, where is the regular field strength tensor and are the regular covariant derivatives in the fundamental and adjoint representation, respectively, g is the gauge coupling, ξ the QCD gauge parameter, and n F the number of different quark flavors. The flow equations are incorporated by introducing Lagrange multiplier fields of mass dimensions 3 and 5/2 that otherwise carry the same quantum numbers as the flowed gluon and quark/antiquark fields B a µ and χ,χ, respectively. Their Euler-Lagrange equations derived from indeed lead to Eq. (2) [16].
Note that the Lagrangian (11) does not include flowed ghost fields d a (t, x) andd a (t, x). They arise in the same way as the usual Faddeev-Popov ghosts c a (x) andc a (x) due to gauge-fixing, and obey the initial condition Similar to the ghosts of the regular gauge fields, they always form closed loops as long as there are no external d a ord a fields (which can be avoided by considering only physical degrees of freedom in amplitudes with external gluons). As becomes clear later, closed loops of only flowed fields vanish, so that one can omit d a andd a already at the level of the Lagrangian.

Perturbative solution of the flow equations
Let us introduce the short-hand notation where it should be clear from the context whether an integration variable is in position (x, y, z, . . .) or momentum space (p, k, q, . . .).
It is helpful to separate the flow equation (2) of the flowed gauge field B a µ (t, x) into a linear and a non-linear part: The linear equation can be solved by introducing the integration kernel which fulfills lim Taking into account the initial condition (1), the full solution of the flow equation is then given by or, in momentum space, By inserting the solution iteratively into itself, one can express the Fourier transform of the non-linear part of Eq. (19) as where The structure of X 3 is identical to the four-gluon vertex of regular QCD. When formulating the Feynman rules later, X 2 and X 3 describe the three-and four-point vertices of the flowed gluon fields.
The flow equation (2) for the flowed quark fields can be solved by again splitting it into a linear and a non-linear part, The linear equation is solved by the integration kernel with the help of which we can write the full solution as Here and in what follows, we suppress the flavor index f unless required for clarity. The non-linear part of the Fourier-transformed field can be expressed as where These expressions lead to the three-and four-point vertices of the flowed quark fields. For χ one proceeds analogously.

Feynman rules (Flowed) Propagators
Plugging in the solution of the flowed gluon field in momentum space Eq. (23) into the two-point function, one finds [16] B a µ (t, p) B b ν (s, q) where D ab µν (p, t, ξ, κ) = δ ab 1 and we have used the result for the fundamental gluon propagator, The factor g 2 is taken into account in the corresponding vertices further below. Since the flowed gluon propagator coincides with the fundamental gluon propagator at t + s = 0, we can express both of them by the same Feynman rule: We refer to this as the (flowed) gluon propagator. Eq. (35) also applies to the mixed propagator A B , which is obtained by setting one of the two flow-time variables to zero.
The same can be done for flowed quark fields. Inserting their momentum space solution in Eq. (29) into the two-point function results in where and we have used the result for the fundamental quark propagator Since one can express the fundamental propagator through the flowed propagator at vanishing flow times, both can be represented by the same Feynman rule: We refer to this as the (flowed) quark propagator, which again includes the mixed propagators ψχ and χψ .

Flow lines
Since there are no quadratic terms of the Lagrange multiplier fields L a µ , λ, andλ in Eq. (16), there are no propagators for these fields. However, the Lagrangian contains bilinear terms of a flowed field and a Lagrange multiplier field. Using standard methods, one derives the two-point function at leading order as where H µν (t, x) obeys the equation and T R is the usual color trace normalization. Since the fundamental gluon field A a µ (x) = B a µ (0, x) does not couple to the Lagrange multiplier field L a µ , and all flow-time variables are positive, one can impose the initial condition so that the unique solution becomes We refer to the BL bilinear as "gluon flow line". The inverse of the factor 1/(2T R ) appears in the corresponding "flow vertices" further below. We can therefore discard it altogether and write where K µν (t, p) has been defined in Eq. (20), and the adjacent arrow indicates the direction towards increasing flow time as implied by the θ-distribution. As opposed to an actual propagator, the gluon flow line is a regular function for all p.
Similarly, one determines the mixed fermionic two-point function at leading order as where with the condition The unique solution reads The Fourier transformed expression defines the "fermion flow line" Feynman rule, where the θ-distribution again imposes a direction as indicated by the adjacent arrow pointing towards increasing flow time: where K(t, p) has been defined in Eq. (27). As usual, the arrow on the fermion line denotes the "charge flow" of the fermion. In this way, we have a unified Feynman rule for both the χλ and the λχ bilinear, where the latter can be obtained analogously as above and simply corresponds to reversing the direction of the charge flow: Since K(t, p) only depends on p only quadratically, the momentum direction for the fermion flow lines is irrelevant.

Flow vertices
Since L B and L χ of Eq. (16) are proportional to a Lagrange multiplier field, the resulting vertices always involve at least one flow line. Such vertices are always associated with a flow-time parameter which is integrated over. We denote them by "flow vertices" and represent them by empty circles in Feynman diagrams. The corresponding Feynman rules can be derived straightforward. In this paper, we define Feynman rules by assuming all momenta to be outgoing. The three-point gluon flow vertex is governed by X 2 : The interaction terms of L B involve exactly one Lagrange multiplier field L(s), cf. Eq. (16). According to Eqs. (42) and (44) Similarly, the four-point gluon flow vertex is governed by X 3 : Again, there is exactly one outgoing flow line, while the other three lines are either flowed gluons or incoming flow lines. This means that the Feynman rule (52) actually represents four different vertices.
Note that the factor 2T R , arising from the trace in Eq. (16), has been discarded in both vertices in accordance with the normalization of Eq. (44).
Similar considerations applied to L χ lead to vertices involving flowed quark fields. For example, the quark flow vertex with one gluon is described by the following Feynman rule: In this case, the fermion line with outgoing charge flow is also outgoing in the gradient flow, while the other two lines can be either flowed propagators, or ingoing flow lines. Thus, this Feynman rule actually represents four different vertices.
A complete list of the Feynman rules can be found in Appendix A.
Let us summarize the Feynman rules for the GFF and point out a few more features: 1. Propagators and flow lines always carry an exponential factor e ±tp 2 for each vertex, where t is the flow time of the vertex. 3 Automated implementation

Generation of Feynman diagram expressions
We generate the Feynman diagrams including symmetry factors and signs (from closed fermion loops) with the help of the program qgraf [25,26]. In the notation of this program, the propagators for the gluon, the ghost, and a single quark flavor can be defined as The sign in the third entry of each square bracket denotes whether the particle is a boson or a fermion.
The flow lines for the gluon and the quark are implemented by introducing separate fields b, fr, and fs for L, λ, andλ, respectively. They are implemented as where the latter two represent the χλ and λχ bilinears (see Eqs. (49) and (50)). With these fields, we can then define the regular as well as the flow vertices. For example, for the trilinear flow vertex of the pure gauge theory defined in Eq. (51), we define which corresponds to the three combinations of Fig. 1.
Already in regular QCD it is convenient to separate the color structure from the rest of the calculation. This becomes non-trivial in cases which involve the four-gluon vertex. It is thus convenient to introduce an auxiliary "particle" Σ a µν whose "propagator" is given by [27] ρ, σ, b µ, ν, a = δ ab δ µρ δ νσ .
The four-gluon vertex can then be replaced by a trilinear ggΣ vertex whose Feynman rule reads This allows to factorize the color factor off all Feynman diagrams, albeit at the cost of increasing their number. We proceed correspondingly for the quartic gluon flow vertices by introducing an additional Σ particle which is directed in flow time. It forms trilinear vertices with gluons and gluon flow lines, as shown explicitly in Appendix A. Again, they contain exactly one outgoing (Σ or gluon) flow line.
By default, qgraf also generates diagrams with closed flow-line loops. We use a Perl script to parse the qgraf output in order to eliminate such diagrams before the actual calculation. They vanish trivially algebraically though, as discussed earlier. In the same way we multiply diagrams by a factor n F for each closed fermion loop. We then process the diagrams with the help of q2e/exp [28,29] in order to insert the Feynman rules and convert the diagrams to FORM code.
Within FORM [30,31], we contract the Lorentz indices, take the fermion traces, and simplify the resulting expressions to a standard form of scalar integrals, as defined below. The color factor is evaluated separately with the help of the color package [32].
For example at the three-loop level (as needed for the results in this paper) the flow-time integrals can be described in the following form are sets of integers, N ≤ 4, t 0 ≡ t, and the upper limits for the flow-time integrations are linear combinations of the other flow-time variables, t up r = t up r (t 0 , . . . , t r−1 ). Initially, all the c i are zero. It is helpful to introduce these parameters for the subsequent discussion though. The momenta p 4 , p 5 , p 6 are linear combinations of the integration momenta p 1 , p 2 , p 3 . For the three-loop quantities considered in this paper, the flow time t is the only dimensionful external scale. A generalization of this notation to a different number of loops and to additional external scales is straightforward.

IBP reduction of flow-time loop integrals
In the first step, we reduce all occurring integrals to so-called master integrals by employing integration-by-parts (IBP) identities [33,34]. They are based on the observation that integrals over total derivatives vanish in dimensional regularization. Applying the operator On the other hand, explicit evaluation of the derivative at the integrand level either reduces one of the indices c k or one of the indices n i by one. One therefore arrives at linear relations among integrals with different n, c, a, and different number of flow-time integrations N .
Applying the above operations to "seed integrals", i.e., integrals with fixed numerical values for the n, c, and a, allows one to build a system of linear relations among the I(t, n, a, c, D). Defining an ordering (complexity) criterion allows one to solve the system using a Gaussian elimination type reduction for a minimal set of integrals (master integrals) [35]. Through this procedure all integrals of the type defined in Eq. (56) are reduced to the minimal set of master integrals, which are simplest by the ordering criterion introduced. The system of linear relations is process dependent and for the observables in this study we construct it in Mathematica [36]. We use the specialized software Kira [37,38] for the reduction of such linear equation systems to solve it.
While Kira alone is sufficient for the reduction of all our integrals at the two-loop level (see e.g. Ref. [21]), we find that the algebraic solution of the system at the three-loop level would require more than 750 GiB of RAM and thus exceeds our available computing resources. 6 Using finite field and reconstruction techniques, which over the last decade have gained an increased use for higher-order perturbative calculations (see, e.g. Refs. [39][40][41]), one can improve on the required computational resources as follows.
Since t is the only dimensionful scale and can thus be factored out, the coefficients of the integrals in the linear system are simply polynomials in D. The individual steps for the reduction of the linear system only involve elementary arithmetic operations, which means that the coefficients of the master integrals are rational functions in D. One can then solve the linear system numerically over a finite field using a sufficiently large number of different integer values for D ("probes") and reconstruct the rational functions in D exactly. With this approach the size of the coefficients remain simple numbers and the requirements on computational resources can be improved.
While Kira already uses pyRed as a first step to remove linearly dependent equations over a finite field, in our approach pyRed is used to reduce the system itself multiple times over a finite field. The resulting numbers are processed with the library FireFly [42], which provides an efficient implementation of interpolation [43,44] and rational-reconstruction algorithms [45,46]. 7 In our case, we require 201 probes, chosen from three different finite fields, defined as prime fields Z p with p a 63-bit prime number, to reconstruct all coefficients, plus one additional probe in a fourth prime field to verify the reconstruction. For each probe, the solution of the system with Kira now just requires about 70 GiB of RAM and takes about three CPU hours. In total, the reduction took less than three days, using ten threads on two Intel Xeon Gold 6138 processors.
For the observables considered further below in this paper at the three-loop level we start with a total of 3195 integrals and can reduce them to a minimal set of 188 master integrals. 8 Their numerical evaluation is described in the next section.

Numerical computation of flow-time loop integrals
As a first step of the numerical evaluation of the gradient-flow integrals given by Eq. (56), we express the propagators through Schwinger-parameter integrals and map them from x ∈ (0, ∞) to y ∈ (0, 1) using the simple transformation x = y/(1 − y). Momentum integrations are performed as D-dimensional Gaussian integrals after a diagonalization. Similarly, all flow-time integrations are mapped to the unit interval with simple linear transformations. The overall result is an integral over a unit hypercube.
The integrand typically involves a number of (overlapping) singularities, which we factorize using FIESTA [47], an implementation of the sector decomposition algorithm [48]. As opposed to regular Feynman integrals where sector decomposition is typically employed in combination with Feynman parameterization, it appears that we cannot restrict the singularities to the lower integration bound only. 9 If an integral involves singularities both at the lower and the upper bound, we split all integration intervals in the middle, and map the singularities to the lower bound by an appropriate change of variable.
The sector decomposed integrals obtained from FIESTA are then integrated with our implementation of fully symmetric integration rules of order 13 [17,50], which can handle integrable logarithmic-like singularities at the integration boundaries very well. We perform all arithmetic with 256-bit precision using the MPFR library [51], and used a local adaptive bisection in the direction of the largest fourth difference. A high precision arithmetic turns out to be necessary to achieve a relative numerical accuracy of 10 −10 or better. The integration uncertainty is estimated by the difference between the integration results of rules of order 13 and 11.
As a check of our results further below, we apply the numerical integration method to our unreduced set of 3195 integrals as well as to the set of 188 master integrals. Finding full agreement for the results obtained with both sets within numerical uncertainties serves as a check of the numerical integration as well as for the reduction procedure.
Analytical computation. For some master integrals we can find analytical results either through elementary methods with the help of Mathematica [36] or by using HyperInt [52]. The integration with Mathematica sometimes yields hypergeometric functions which can be expanded with the package HypExp [53,54]. For example as shown further below, we find analytical results for all contributions which contain at least one factor of T R in the color structure. A fully analytical calculation of integrals with other color factors requires a more detailed investigation. However, for all practical purposes, the numerical results provided in this paper are (more than) sufficient.

Observables
Among the simplest quantities one can consider within the GFF are vacuum expectation values of gauge-invariant operators at finite flow time. It is one of the remarkable properties of the GFF that these operators do not require any renormalization beyond that of regular QCD, and that of the involved flowed fields. For the gauge coupling and the quark masses, the MS renormalization is given by the replacement in Eq. (12), where µ is the renormalization scale, α s = g 2 /(4π), and γ E = 0.5772 . . . the Euler-Mascheroni constant. Through the perturbative order required in this paper, the renormalization constants are given by with C F and C A are the quadratic Casimir eigenvalues of the fundamental and the adjoint representation of the gauge group, respectively. Furthermore, T F = T R n F , with n F the number of quark flavors, and T R the trace normalization in the fundamental representation. For SU(N c ), it is C F = (N 2 c − 1)/(2N c ), C A = N c , and T R = 1/2. The flowed gauge field B a µ (t, x) does not require renormalization, so that for example matrix elements of the gluon action density, are finite after just the renormalization of g and m f . This allows for a direct comparison of results obtained in different regularization schemes (lattice and perturbation theory, for example).
On the contrary, flowed quark fields require a renormalization factor Z 1/2 χ (α s ) in order to render Green's functions finite. In the MS scheme, it is [3,21] with γ χ,0 = 6 C F , The quantity thus acquires an anomalous dimension, which prevents a direct comparison of results from different regularization schemes. Alternatively, one may work with "ringed quark fields" [20], which amounts to using instead of Z χ in order to renormalize the quark fields, where This corresponds to a "physical" renormalization scheme, which means that the anomalous dimension of the operator vanishes. In Eq. (67), we used the notation for the zero-momentum Fourier transform of the vacuum expectation value of an operator O(t, x), which we adopt in the following.
The Feynman rules for the operators E(t, x), S(t, x), and R(t, x) introduced above are obtained in the same manner as those described in Sect. 2.3. Expressing them in terms of Using the methods described in Sects. 2 and 3, we can calculate the vacuum expectation values of E(t, x), S(t, x), and R(t, x) through three loops. Sample diagrams for these quantities are shown in Fig. 2. The number of diagrams in each case is given in Table 1; recall, however, that this takes into account the re-writing of the four-gluon vertex into trilinear vertices as described in Sect. 3. Not included in this number are diagrams with closed flow-line loops, but integrals with scale-less sub-loops are counted in.
where α s = α s (µ), with µ the renormalization scale. As indicated, quark mass terms are neglected here. Their effect at NLO has been studied in Ref. [17] and it was found to be negligible. In this case, the e i depend only on the product z ≡ µ 2 t. For reasons that will become clear shortly, it is convenient to parameterize this dependence in terms of the variable A natural choice for the renormalization scale is defined by the inverse "smearing radius" q 8 ≡ 1/ √ 8t of the gradient flow. However, from explicit higher order calculations, it was found that the slightly larger value corresponding to L(µ 2 0 t) = 0, improves the stability of the perturbative corrections [21] (see also Ref. [17]). This was corroborated by the behavior of the t → 0 extrapolation for thermodynamical quantities evaluated in Ref. [22] using the two-loop result of Ref. [21].
We thus write the coefficients of Eq. (71) as e 0 (z) = e 0,0 , e 1 (z) = e 1,0 + β 0 L(z) , with β 0 , β 1 from Eq. (62). For the coefficients e i,j , we find e 0,0 = 1 , e 1,0 = 52 9 where ζ(z) is Riemann's ζ function with ζ(3) = 1.20206 . . .. The three dots in the coefficient of T F C A indicate that we were able to obtain an expression in analytical form. It can be found in Appendix B, Eq. (128). Only four decimal places of the numerical result for the C 2 A coefficient are displayed here, while our estimate of the numerical accuracy, obtained by propagating the uncertainty of the individual integrals to the final result, is about six digits beyond that. This estimate matches the observed differences between the numerical and analytical results in the cases where the latter have been computed. The same statements hold for the other observables listed below. The NLO coefficient e 1 was first evaluated in Ref. [3]. Setting µ = q 8 , one finds that the NNLO result agrees with Ref. [17] at the sub-percent level. 10 The logarithmic terms e n,k L k (z) obey the all-order recursion formula ( which follows from renormalization-group invariance, i.e.
This provides a welcome check of our calculation. The residual dependence on µ can be used as probe of the perturbative behavior. We refer to Ref. [17] for such a study.
Since E(t) is also gauge independent, we are free to choose Feynman gauge for the QCD gauge parameter (ξ = 1) and κ = 1 for the gradient-flow gauge parameter for the evaluation of the results presented above, which facilitates the actual calculation significantly. However, we also check gauge invariance explicitly by evaluating the terms with the highest power in the gauge parameter ξ, i.e. ξ 4 , and show that its coefficient vanishes. Indeed, this already happens algebraically after the reduction to master integrals, even before their numerical values are inserted. Keeping the gauge parameter fully general unfortunately leads to intermediate expressions which exceed our available computing resources.

Results for the quark condensate at three loops
To obtain the quark condensate with ringed fields we first consider the conversion factor ζ χ (t, µ) from MS renormalized to ringed quark fields as defined in Eq. (69). Applying the methods described in Sect. 2, we can calculate the diagrams contributing to the Green's function R(t) for massless quarks, which eventually leads to the result The NLO result has been obtained in Ref. [20], and the µ-dependent NNLO terms were derived in Ref. [21]. That paper also used a preliminary result for the coefficient C 2 , derived for the first time in the current paper. Keeping four decimal places for the individual color coefficients, we find where again, as indicated by the three dots, the coefficient of C F T F was obtained in analytical form and is given in Appendix B, Eq. (130).
Let us now turn to the quark condensate S(t) . Since it vanishes in the chiral limit m f = 0, we compute the leading coefficient of the expansion in the quark masses m f , where it is understood that all quark masses are set to zero after taking the derivative. This expansion can be applied before the loop and flow-time integrations are carried out, which means that, in the quark propagators, we write and keep only the leading linear term in m of the integrand. The resulting integrals are then again of the form shown in Eq. (56). We thus write Renormalizing the quark masses in the MS scheme according to Eq. (60), we find s 0 (z) = s 0,0 , s 1 (z) = s 1,0 + γ m,0 2 L(z) , with z = µ 2 t and the γ m,i and β i from Eq. (62). The non-logarithmic terms are given by where again the T F C F term is known analytically, as indicated by the three dots. It can be found in Appendix B, Eq. (131). The NLO term has been obtained in Ref. [55], the NNLO term is new. Again, only four decimal places are displayed and our uncertainty estimate for the numerical integration is several digits beyond that.
The logarithmic terms s n,k L k (z) obey the all-order recursion formula (1 ≤ k ≤ n) which follows from renormalization group invariance, i.e.
and the renormalization group equation for α s given in Eq. (77). This again provides a welcome check of our calculation. The residual dependence on µ can be used as probe of the perturbative behavior. It is displayed in Fig. 3 for three values of the central energy scale µ 0 = (2te γ E ) −1/2 , see Eq. (73). We observe a qualitatively similar behavior as was found for t 2 E(t) in Ref. [17]. While for large µ 0 , the reduction of the renormalization scale with increasing perturbative order is quite significant, this improvement reduces as µ 0 gets closer to the non-perturbative regime. Taking the variation of µ by a factor of two around the central scale as an estimate of the uncertainty due to missing higher orders, one finds that they overlap between successive perturbative orders for all three values of µ 0 . However, the trend of the curves does not explicitly support µ 0 as the central scale in this case, but seems to favor a slightly smaller scale.
Another important consistency check is that S(t) and R(t) can be renormalized with the same flowed-quark field renormalization constant Z χ as obtained in Ref. [21]. S(t) and R(t) are also gauge independent. For the calculation of the full NNLO result, we again chose Feynman gauge (ξ = 1) and κ = 1. We check gauge invariance explicitly by evaluating the terms with the highest power in the gauge parameter ξ, which is ξ 3 for S(t) and R(t) , and show that their coefficients vanish. This happens here also immediately after the reduction to master integrals.     The renormalization group invariance of E(t) allows us to fix a scale µ = ρ/t with constant ρ ∈ R, and its proportionality to α s suggests to define a "gradient-flow coupling" For the usual definition of the gradient-flow coupling [3], ρ is set to 1/8 = 0.125, while ρ = 1/(2e γ E ) = 0.281 . . . appears to be more suitable from a perturbative point of view (cf. Eq. (73)).
Perturbatively solving Eq. (71) for α s results in with the coefficients e n (ρ) given in Eq. (74). Such a transformation between mass-independent renormalization schemes leaves the first two coefficients of the β function invariant, i.e.
with β 0 and β 1 from Eq. (62). The third coefficient is given bŷ with the MS coefficient β 2 which we quote here for the SU(3) gauge group:

Conclusions
A framework for the calculation of correlation functions in the perturbative gradientflow formalism using tools and techniques from standard perturbative QCD calculations has been presented. It employs a Feynman-diagrammatic approach to the gradient-flow formalism, based on the five-dimensional field-theoretical formulation of Ref. [16]. We implemented the corresponding QCD Feynman rules within the gradient-flow formalism as well as those of the composite operators relevant for the quantities considered in this paper, and automated the reduction and calculation of algebraic expressions and loop integrals to a large extent. An important ingredient of our setup is the reduction to master integrals with the help of integration-by-parts identities that also involve the flowtime variables. Except for a few cases where an analytic result was derived, the master integrals were evaluated numerically.
To demonstrate the viability of this setup, we reproduced and improved the three-loop result for the gluon condensate of an earlier calculation, which was based on explicit field-theoretic calculations at the level of Wick contractions [17].
Using this setup, we then evaluated the three-loop approximations to the quark condensate and the conversion factor from MS renormalized to ringed quark fields, both of which are new results. Checks on gauge parameter independence, renormalization group invariance, as well as alternative calculational approaches (with and without reduction to master integrals) convinced us of the correctness of these results.
The gluon and the quark condensate lend themselves for straightforward definitions of a gradient-flow gauge coupling and a gradient-flow quark mass, for which we derive the three-loop matching relations to the MS scheme. Provided that sufficiently accurate lattice results for these quantities become available, the conversion factors evaluated in this paper should allow for precise first-principle determinations of the gauge coupling, 11 and possibly also quark masses.