Nonminimal gradient flows in QCD-like theories

The Yang-Mills gradient flow for QCD-like theories is generalized by including a fermionic matter term in the gauge field flow equation. We combine this with two different flow equations for the fermionic degrees of freedom. The solutions for the different gradient flow setups are used in the perturbative computations of the vacuum expectation value of the Yang-Mills Lagrangian density and the field renormalization factor of the evolved fermions up to next-to-leading order in the coupling. We find a one-parameter family of flow systems for which there exists a renormalization scheme in which the evolved fermion anomalous dimension vanishes to all orders in perturbation theory. The fermion number dependence of different flows is studied and applications to lattice studies are anticipated.


Introduction
The past decade has seen a lot of interest in the Yang-Mills gradient flow. Its classical version was formulated in a pure mathematical context in [1,2], and the applications to quantum field theory followed in [3,4]. See also [5] for a nice review on 'smearing', as the effect of the gradient flow on the fields involved is often referred to.
The gradient flow evolved fields have very particular renormalization properties, first demonstrated in [6], see also [7]. Correlators consisting of evolved gauge fields B µ -and possibly composite operators thereof -do not require any renormalization beyond that of the usual renormalization of the fields and parameters of the nonevolved theory over which the path integral average is being taken, and correlators containing evolved fermion fields require a universal 1 field renormalization factor [8].
A much studied quantity is the observable with L Y M the Euclidean pure Yang-Mills Lagrangian density, g 0 the bare gauge coupling, and t the gradient flow 'time' which is of dimension length 2 , see appendix A for conventions used throughout this paper. Since E(t, x) solely consists of evolved gauge fields 2 , it acquires no multiplicative renormalization factor. Its expectation value E(t) is used in lattice studies for e.g. scale setting [9], and defining finite volume running couping schemes 3 [12]. It is customary to apply the Yang-Mills gradient flow, defined by to QCD-like theories, where also (nonevolved) fermions are present. We consider this equation more in depth in section 2, but for now it is sufficient to note that the gradient flow (1.2) drives the gauge field B µ (t, x) towards the stationary points of the Yang-Mills action for increasing flow time t. Since some properties of the expectation value of E(t, x) in pure Yang-Mills theory are attributed to this fact [4], one might wonder what happens if one adds the fermionic term to the action in (1.2), i.e. S Y M → S Y M + S F (see appendix B for the explicit expressions). The consequences of this change are the subject of this paper.
There exists a close relation between the gradient flow and the Langevin equation used in stochastic quantization [13]. Since it is customary to include the fermionic matter term in the Langevin equation when stochastically quantizing QCD-like theories, we will make this relation more explicit. 1 With universal we mean that any correlator consisting of n evolved fermion fields χ, n fieldsχ and any number of Wilsonian normalized evolved gauge fields Bµ is rendered finite by multiplicative renormalization with Z n χ , also if evolved composite operators are present. We will elaborate on this in section 2.4. 2 This is only true for Wilsonian normalized gauge fields. We choose to work with canonical instead of Wilsonian normalization from the outset. This will make for a neater comparison with stochastic quantization, and we are mainly concerned with perturbative analyses. All statements about Wilsonian normalized Bµ holds equally true for canonically normalized Bµ multiplied by g0. 3 There is some degree of disagreement in the lattice community on the validity of certain schemes defined in this way, especially in relatively large N f studies, see e.g. [10,11].

Motivation from stochastic quantization
In stochastic quantization the basic idea is to introduce an extra scale t -we will shortly see how it is related to the flow time t of (1.2) -and a Gaussian noise field η(t, x), and describe the t evolution of some field φ(t, x) by the generalized Langevin equation where M (x, y) is some kernel 4 , S is the Euclidean action for the theory under study, and δS δφ(t,y) is short for δS δφ(y) φ(y)=φ(t,y) . The Gaussian noise satisfies η(t, x)η(s, y) η = 2M (x, y)δ(t − s), η(t, x) η = 0 (1.4) where the η subscript on the angle brakets indicates that the average over η is being taken, see e.g. [14] for details. Note that the presence of the noise in (1.3) makes the value of φ at large t completely independent of the boundary condition at t = 0, and a common choice is φ(0, x) = 0. We now imagine that the fields are coupled to some heat resevoir. It can be shown that in the large t limit we reach equilibrium [14], and all stochastic correlators converge to their Euclidean quantum field theory counterparts: Now, the gradient flow evolution can be understood as follows 5 : at some large time t = t 0 we abruptly switch off the heat resevoir and the noise field. At some later time t 1 > t 0 , the fields will have evolved via the flow equation which is exactly the gradient flow equation.
When stochastically quantizing pure Yang-Mills theory, the following Langevin equation is used [14] -here we neglect the possible presence of a 'gauge fixing term', usually called a Zwanziger term in this context: the field strenght. However, when considering QCD-like theories, a different Langevin equation is being used for the gauge field in the presence of fermionic matter: with T a the SU (N ) generators, where for the Langevin evolution of the fermions one can take [15] where θ,θ are Gaussian noise fields. 4 The standard Langevin equation is the special case M (x, y) = δ (d) (x − y). 5 This reasoning is based on [5].

Nonminimal gradient flows
In gradient flow studies it is customary instead to use the Yang-Mills gradient flow of (1.2), i.e. (1.7) with 6 η = 0, also when fermions are present. We set out to investigate what will change if we instead use (1.8) with η = 0 as the basis for the gradient flow, i.e. the case where we include the fermion bilinear color non-singlet vector current in the right hand side of the flow equation. We will refer to this as the nonminimal gradient flow equation.
For the evolution of the fermions we are considering two different flow equations. Firstly we use the one that is standard in gradient flow studies [8] using the operator D 2 , and secondly we use the operator / D 2 instead, i.e. (1.9) with θ =θ = 0.
All flows must respect the symmetries of the nonevolved theory, and the Yang-Mills gradient flow is indeed the simplest possibility. However, the nonminimal cases we will investigate have the conceptual advantage that they are actually driving the gauge field towards the stationary points of the full QCD-like action, not just of the Yang-Mills action.
That said, we also carry out the exercise of writing the most general flow equations for the gauge field and the fermions, respecting the symmetries of the nonevolved theory, for completeness.
We use the nonminimal flows to calculate the expectation value of the operator E(t, x) from (1.1) to next-to-leading order in the coupling, generalizing the result first presented in [4]. We also generalize the demonstration of the renormalization properties of the evolved fields of [6,8] (see also [7]), and calculate the evolved fermions' field renormalization factor Z χ for the different nonminimal flows. We find that there exists a one-parameter family of flow equation systems for which there is a renormalization scheme in which the evolved fermion anomalous dimension vanishes to all orders in perturbation theory.
We solely consider massless fermions in this work. The quantity E(t) will have a mass dependence at next-to-leading order in the gauge coupling, see e.g. [16] for the Yang-Mills gradient flow case. However, the MS expressions for Z χ and the other renormalization properties discussed in this paper are independent of the presence of fermion masses.

Organization of the paper
In section 2 we review the basics of the Yang-Mills gradient flow (YMGF), the next-toleading order in the gauge coupling result for E with E in (1.1), and the renormalization properties of the evolved fields. In section 3 we introduce the nonminimal gradient flow (NMGF), calculate E , present the generalization of the arguments for the renormalization properties for the nonminimally evolved fields, and calculate the field renormalization factor of the evolved fermions. Subsequently, in section 4 we include the effects of replacing the operator D 2 in the fermion flow equation with / D 2 , coining this case the 'slashed' nonminimal gradient flow (sNMGF). In section 5 we investigate the N f dependence of E with N = 3, calculated with the three different flows. In section 6 we generalize the flow equations by putting them in the most general form compliant with all the symmetries of the nonevolved theory. We present our conclusions and outlook in section 7. Appendices A to G contain relevant conventions and calculations.

Yang-Mills gradient flow
The Yang-Mills gradient flow equation is given by For further conventions and the explicit expression for the action we refer the reader to appendices A and B, respectively. When performing perturbative calculations, this flow equation is usually conveniently modified by adding a term which in the context of stochastic quantization is known as a Zwanziger term [14]: and subsequently setting α 0 = 1 [4]. The solution to this modified flow equation is related to the solution of (2.1) via a t-dependent gauge transformation (see appendix C), and thus observables will be independent of α 0 . The solution to (2.2) with α 0 = 1 is given by where the scalar kernel reads We will often use the 'exponential-of-Laplacian' notation, with the Laplacian ∆ and

Fermion flow
The conventional flow equations for the fermions introduced in [8] read In section 4 we investigate the consequences of using / D 2 instead of D 2 , establishing a more direct connection to the stochasic quantization interpretation outlined in section 1.1.
Similarly to the addition of the Zwanziger term in (2.2), these flow equations are conveniently modified by adding an α 0 -dependent term: (2.10) and the solutions to (2.10) are related to those of (2.8) by the same t-dependent gauge transformation as in the pure gauge case from section 2.1 (see appendix C). The solutions to (2.10) with α 0 = 1 are given by

Perturbative calculation of E
For later comparison, we briefly review the perturbative calculation of E up to O(g 4 0 ) of [4]. At this order in the coupling it only involves the evolved gauge fields given in (2.7); the fermionic contribution will solely be coming from the nonevolved path integral averaging.
It is useful to write the evolved gauge field as a series in powers of the bare coupling: with the first few orders iteratively given by These different orders can be conveniently expressed in diagrammatic form, see figure 1.
The observable E from (1.1) written in terms of the B µ field reads: (2.15) and the expectation value E can be computed order-by-order in the coupling [4,16], using (2.14) combined with the contributions stemming from the nonevolved path integral (a) (b) (c) Figure 1: Diagrammatic representation of the first two nontrivial orders of the evolved gauge field's expansion in the coupling g 0 from (2.13). (a) represents B µ,1 from (2.14), (b) collectively represents the first two terms in B µ,2 from (2.14), and (c) represents the last term in B µ,2 from (2.14). In this notation, double lines stand for kernels ('exponentials-of-Laplacian') together with a flow time integration, white blobs represent GF vertices, black blobs are the usual QCD vertices, and the crossed blobs are the external points.
averaging. The result up to next-to-leading order in the coupling is given by (we use dimensional regularization with d = 4 − 2ε) This observable is rendered finite by the coupling renormalization of the nonevolved theory.
In the MS scheme one has: with g(µ) the renormalized coupling, and After renormalization 7 we express the renormalized coupling α(µ) = g 2 (µ) 4π in terms of the renormalization group invariant coupling α(q) by inverting (G.6), and subsequently setting q = (8t) −1/2 and N = 3 we arrive at For the renormalized result for general N we refer the reader to appendix G. 7 We would like to emphasize that the fact that E is rendered finite by the coupling renormalization at this order does not constitute an example of the 'nonrenormalization' of the evolved gauge field Bµ, to be discussed in section 2.4. In order to achieve this, one has to use the 2-loop universality of Zg from (2.19) and calculate E to next-to-next-to-leading order. This has been done in [16], where the 'nonrenormalization' is indeed confirmed.

Renormalization of evolved fields
It has been shown in [6] (see also [7]) that in Wilsonian normalization all correlators consisting of only B µ fields are finite, i.e. they do not require renormalization beyond the usual renormalization of the parameters of the four dimensional theory over which the path integral average is being taken. When switching to canonical normalization, i.e. taking B µ → g 0 B µ , any correlator consisting of m number of B µ fields now automatically contains a factor of g m 0 , which will renormalize as given in (2.18) in the case of MS. The fermion fields χ andχ do require a renormalization factor [8], namely and the contributing self-energy diagrams together with their values are presented in appendix E. Any bare evolved correlator consisting of an arbitrary number of g 0 B µ 's, and n number of χ's andχ's is now rendered finite by multiplicative renormalization with Z n χ , supplemented with the usual renormalization of the parameters 8 of the nonevolved theory. An important feature is that this also holds true when some of the evolved fields have coinciding spacetime positions at stricktly positive flow time; i.e. composite operators of evolved fields do not aquire any additional renormalization factor, in contrast to the nonevolved case.
The all order demonstration of these statements is found in [6] and [7]. In the remainder of this section we present the general reasoning for this demonstration, without any intention or pretention of being complete. We nevertheless choose to include it here, since these will be the arguments that we generalize later on in sections 3.3 and 4.2.
The starting point of the demonstration from [6] is to describe the theory in d + 1 dimensions, the extra dimension being the flow time direction t ∈ [0, ∞). The evolved fields are being considered as independent from the nonevolved elementary fields, apart from their implicit dependence through the boundary conditions. The d + 1 dimensional 'bulk' action is given by: where we have introduced the lie algebra valued lagrange multiplier field L µ (t, x) = L a µ (t, x)T a with purely imaginary components, the two fermionic lagrange multipliers λ(t, x) and λ(t, x) which carry the same indices as the quark fields, and the bulk ghost fields d(t, x) andd(t, x). The ghost d has boundary condition while the other new fields do not obey any. Note that the lagrange multipliers L µ , λ and λ are there to enforce the flow equations upon variation. Thed field acts as a lagrange multiplier enforcing the correct diffusion equation on the d + 1 dimension ghost, such that the flow equations are invariant under an infinitesimal gauge transformation . This generalizes the BRST symmetry to the d + 1 dimensional bulk action: with δ = δ BRST , which can be checked using the BRST identities in appendix D. Together with the BRST invariance of the d dimensional QCD action, and the fact that the path integral measure is invariant, we have with O any combination of evolved and nonevolved elementary fields. Equation (2.29) has some important consequences. For our purposes the most important relation that follows, for reasons becoming obvious below, is The next step in the demonstration is to exclude the possibilty of counterterms. These are either localized in the d + 1 dimensional R d × (0, ∞) bulk or at the d dimensional R d boundary [17]. We first focus on the former.
The impossibility of bulk divergences is shown by first eliminating the elementary nonevolved fields from the quadratic part of the bulk action, by plugging in where b µ , Ψ,Ψ obey homogenous boundary conditions. The terms involving the nonevolved elementary fields all drop out, since they solve the linearized flow equations. Therefore, the only nonzero 2-point bulk correlators are bL , Ψλ , λΨ and dd . By subsequently analyzing the interaction part of the bulk action it becomes clear that the only nonzero bulk correlators of these fields are trees, i.e. they can never form loops, and thus will be finite. Therefore there will be no bulk counterterms.
The remaining possibility is the presence of boundary counterterms. Dimensional analysis together with symmetry requirements and ghost number conservation tells us that the only allowed new counterterms are of the form At this point, the importance of (2.30) becomes apparent; the tree level computation of its renormalized version 9 forces z 1 = z 2 = 0, i.e. there is no field renormalization for L µ andd, and thus no field renormalization for B µ and d (since there can be no bulk counterterms). However, there is no restriction on z 3 , and consequentlyλ and λ do get a field renormalization factor which implies the reciprocal field renormalization factors for χ andχ in order not to have a bulk counterterm. This completes the reasoning for the case of correlators consisting of the evolved fields at distinct points in spacetime. In order to include composite operators it is sufficient to notice the following: the small -or indeed vanishing -distance limit of any number of evolved fields inside a correlator will always remain finite for t > 0 due to the exponential suppression of the high frequency modes.
Missing details of the demonstration can be found in [6] and [7].

Defining the nonminimal gradient flow
We now include fermionic matter in the flow equation for the gauge field, giving rise to the nonminimal gradient flow (NMGF) equation: is the gauge covariant matter current. The quark fields χ andχ are the same as in section 2.2, apart from the fact that the B µ field in their flow equations is now the solution of the NMGF equation (3.1). Importantly, the NMGF equation is gauge invariant, as shown in appendix C.
We can modify the flow equation (3.1) by the same Zwanziger term as in (2.2) by utilizing a particular t-dependent gauge transformation, see appendix C. This gives rise to the modified NMGF equation For α 0 = 1 the solution to the modified nonminimal gradient flow equation reads: where R µ is still the same as in (2.5), though now the B µ 's are the solution to (3.3) instead of (2.2), and we have explicitly put the superscript 'YM' on the solution (2.3) of the Yang-Mills gradient flow equation (2.2). Again, we write the fields as a series in powers of the bare coupling: and thus we have for the first three orders of B µ These quantities are diagrammatically represented in figure 2.

Perturbative calculation of E
We present the calculation of the new contributions -collectively called E χ -to E stemming from the fermionic term F µ in B µ , see (3.4). The full observable has the same form as in (2.15), though now the B µ 's represent solutions to the NMGF, so we have 14) The new contributions start at O(g 4 0 ), and are given by and the integrals I 1 -I 9 are listed in appendix F. These contributions are diagrammatically represented in figure 3. The total fermionic contribution is subsequently given by Figure 3: Diagrammatic representation of the contributions to E due to the NMGF, whose total value is given in (3.21). (a) represents E χ,1 from (3.15), (b) represents E χ,2 from (3.16), and the combination of (c) and (d) represents E χ,3 from (3.17).
The integrals I 1 , I 5 -I 8 , are solved using Schwinger parametrization: and Gaussian integration The expressions obtained in this way contain gamma functions, incomplete beta functions and hypergeometric functions, which depend on d only; the t dependence factors out.
The ε-expansions with d = 4 − 2ε of I 1 , I 5 -I 8 are provided in appendix F. Adding all contributions, we find for (3.18): which is manifestly finite, as is to be expected anticipating the results of section 3.3. We now have for the full observable (3.14): Following the same steps as in section 2.3, i.e. renormalizing the coupling as in (2.18), expressing the renormalized coupling α(µ) in terms of the RG invariant coupling α(q) by inverting (G.6) and setting q = (8t) −1/2 and N = 3, we now have: Note that the N f contribution in k 1 is negative, in contrast to the result from (2.20) that is obtained using the Yang-Mills gradient flow. See appendix G for the renormalized result for general N . We further comment on our result (3.24) in section 5, where we will make a more elaborate comparison with the YMGF result from (2.20).

Renormalization of evolved fields
The result for E from the previous section seems to indicate that the nonminimally evolved gauge field remains devoid of any renormalization. However, one might suspect that at the next order in g 0 , i.e. O(g 6 0 ), divergences arise from the fermionic contribution; at that order the self-energy contribution to the evolved fermionic propagator might have to be canceled with Z χ factors. Our aim here is to show that this is not the case, and any correlator consisting only of evolved gauge fields B µ -g 0 B µ when canonically normalized -is still finite, without any renormalization beyond that of the nonevolved d dimensional theory that is being used in the path integral averaging.
The demonstration is a straighforward generalization of the one from [6,8], see also [7], that we have reviewed in section 2.4. We will now show where the changes are with respect to the YMGF case, and why the same reasoning still holds.
First of all, the gauge bulk action S G,f l from (2.24) changes in order to comply with the NMGF equation (3.3): Using the BRST variations provided in appendix D, we establish that the complete bulk action is again invariant: Also, the path integral measure has not been changed and is still BRST invariant, therefore δO = 0 still holds, where O is any combination of evolved and nonevolved elementary fields. Therefore, also (2.30) still holds.
The argument for the absence of bulk counterterms translates one-to-one from the YMGF to the NMGF case, with one addition. There now is an additional interaction term proportional toχL µ χ, which upon the substitution of (2.31) yields the bulk interaction ΨL µ Ψ. However, recalling that the only nonzero propagators are bL , Ψλ , λΨ and dd , in order to form a loop in the bulk the presence of the interaction termλB µ λ is essential. It is clear from the fermionic bulk action (2.25) that this term is absent, and thus no loops can be formed. Therefore, all bulk correlators still solely consist of trees and no divergences can arise, i.e. there will be no bulk counterterms.
The allowed boundary counterterms remain to be of the form given in (2.32). Also, since we still have δO = 0, and thus in particular (2.30) still holds, we have z 1 = z 2 = 0, i.e. there is again no field renormalization factor for L µ , and consequently also no field renormalization factor for B µ . Now, returning to the issue raised at the beginning of this section: might it be that the χ andχ field in j a µ in the solution (3.4) of the NMGF equation (3.3) have to be renormalized in order to render correlators of B µ finite? It is now clear that the answer is no. A Z −1 χ factor in the flow equation (3.3) would imply a Z −1 χ factor multiplying j a µ in the gauge bulk action (3.25). This in turn would imply a bulk counterterm, unless L µ = Z χ (L R ) µ i.e. that L µ renormalizes. The bulk counterterm is excluded due to the absence of loops in the bulk, and the renormalization of L µ is excluded due to BRST symmetry. NM1 NM2 Figure 4: New next-to-leading order diagrams for the NMGF evolved fermion 2-point function. Combined with the diagrams from figure 10 they represent the full NMGF evolved fermion 2-point function, which is rendered finite by multiplying with the renormalization factor Z χ from (3.28). Table 1: UV divergence of the diagrams from figure 4 in units of C 2 (R) The discussion presented above does not make any statement about the value of Z χ , and there is no reason why it should remain the same as when using the YMGF, (2.22). The new vertex at O(g 0 ) due to the NMGF -see (3.9) and corresponding figure 2.c -gives rise to a new contribution to Z χ . In figure 4 we show the diagrams that contribute to the evolved fermion 2-point function, and their values are presented in table 1 in units of C 2 (R) is the evolved fermion propagator: Combining these new contributions with the self-energy diagrams presented in appendix E we find NMGF : 4 Using / D 2 in the fermion flow equations Next we present the fermion flow equation using / D 2 instead of D 2 , and calculate its contribution to E . The well-known relation between these operators is: The modified fermion flow equations from (2.8) now read with The solutions to these flow equations are where ∆ and ← − ∆ are given in (2.12). The combination of these fermion flow equations and the nonminimal gauge field flow from section 3 will be referred to as the 'slashed' nonminimal gradient flow (sNMGF).
It is convenient to define the extra contribution in χ(t, x) stemming from the σ µν term in (4.1): where χ / ∆ is representing the solution from (4.3), and χ ∆ the solution from (2.11). The extra contributions χ * ,χ * are given by Writing these contributions as a series in powers of g 0 (a) (b) Figure 6: Diagrammatic representation of the contribution to E due to the sNMGF, whose total value is given in (4.14). (a) and (b) collectively represent E * χ from (4.12).
we have for the first two orders: The first -and at this order in g 0 only -nonzero new contribution to B µ from (3.4) is then given by These new contributions are diagrammatically represented in figure 5. Note that the new interaction vertices due to the sNMGF are represented by white blocks instead of blobs.

Perturbative calculation of E
The calculation of E now consists of calculating one additional contribution with respect to the calculation presented in section 3.2. Comparing to (3.14), we now have: with E χ as presented in section 3.2, and which is diagrammatically represented in figure 6. Employing the expressions presented in the previous section, we evaluate this term as  The integrals I 10 , I 11 and their ε-expansions are listed in appendix F. The new contribution is then given by and, similarly to E χ from (3.21), this is manifestly finite. Again this is to be expected, anticipating the arguments provided in section 4.2.
We obtain for the total fermionic contribution: We therefore now have for the full observable from (4.11): (4.17) Repeating the procedure from sections 2.3 and 3.2, i.e. renormalizing the coupling as in (2.18), expressing the renormalized coupling α(µ) in terms of the RG invariant coupling α(q) by inverting (G.6) and setting q = (8t) −1/2 and N = 3, we now obtain: Again, for the renormalized result for general N we refer the reader to appendix G. We will compare the sNMGF result (4.18) with the NMGF and YMGF results (3.24) and (2.20), respectively, in section 5.

Renormalization of evolved fields
The arguments presented in section 3.3 are fully applicable to the case presented here. The added gamma matrix structure in (4.    figure 7 in units of C 2 (R) The full sNMGF evolved fermion 2-point function is given by the combined contributions from figures 10, 4 and 7, and it is rendered finite by multiplying with the sNMGF evolved fermion renormalization factor sNMGF :  It is convenient to express α(q) in terms of the fundamental scale Λ QCD 10 using the universal 2-loop renormalization group improved UV asymptotic expression [18], see (G.7). We expand the coupling in powers of 1/l, with l = − log 8tΛ 2 QCD , and we obtain for (5.1): The values for t 2 E(t) i are plotted for N f = 4, 8 in figure 8.
In [4] it was found, using lattice simulations, that for N = 3, N f = 0, beyond the perturbative regime, t 2 E(t) grows roughly linearly with t. There, as a possible explanation 10 In general the renormalization group invariant scale ΛQCD depends on N , N f , the value of the coupling at some given scale, and the renormalization scheme used. We will solely use it as a measure for the flow time t and will not worry about its value or how it changes with N f ; when comparing theories with different N f we will keep the dimensionless product tΛ 2 QCD fixed. for this behaviour the author alluded to the fact that the gradient flow drives the gauge fields towards stationary points of the YM action. For these configurations the right hand side of the flow equation (2.1) is small, implying that E will change relatively little with t. When considering cases with N f = 0, this same reasoning should result in a more pronounced slowing down of t 2 E(t) in the (s)NM cases compared to the YMGF, since the latter now no longer drives the gauge field towards the stationary points of the full action of the theory.
The graphs in figure 8 seem to confirm this reasoning. Moreover, as N f is increased, the difference between the use of the different flow equations becomes more clear. It would be interesting to see this behaviour being reproduced in lattice simulations when the different flows are implemented.

Sensitivity to change in N f
We can study the sensitivity to changes in N f of the different E i 's from (5.1), (5.2) by considering the relative change: The results are plotted in figure 9. It shows that the sNMGF provides a much more stable value for E under changes of N f than the YMGF. For example, for N f = 4 → 5, at t = 0.008/Λ 2 QCD , the sNMGF result changes by 2%, while for the YMGF we have a change of around 13%.
We would again like to stress that although Λ QCD will also change with N f , the dimensionless product tΛ 2 QCD of the flow time, or 'smearing radius' √ t, and the fundamental physical scale Λ QCD is kept fixed.

Generalized flows
The extensions of the Yang-Mills gradient flow presented in sections 3 and 4 are motivated by the similarities between the Langevin equations used in stochastic quantization and the gradient flow, and by the observation that the YMGF does not drive the gauge field towards the stationary points of the full action of the QCD-like theory when fermionic matter is present.
However, we also noted that in principle any gradient flow equation respecting all the symmetries of the nonevolved theory is equally valid from a purely mathematical point of view. Therefore we also present the most general flow equations: where we include the modification terms proportional to α 0 , and a 1 , a 2 are real constants parametrizing the amount of 'nonminimality'. The YMGF is the minimal case: Using these solutions we calculate the expectation value of the observable E(t, x) from (1.1). For general N we obtain The evolved fermion renormalization factor is also calculated using the generalized flow equations, and reads Z χ (a 1 , a 2 ; g(µ), ε) Thus the anomalous dimension associated with the evolved fermion fields is given by Note that when a 1 (1 + 6a 2 ) = 3 the first coefficient of the anomalous dimension vanishes. This implies that there exists a renormalization scheme such that γ χ vanishes to all orders in the gauge coupling, see e.g. [19].

Conclusion and outlook
In this paper we have investigated some of the properties and consequences of nonminimal gradient flows in non-Abelian gauge theories containing fermions, i.e. flow equations which are defined using the full action, instead of only the Yang-Mills action. The most important new feature is that the flow equation for the gauge field B µ acquires a purely fermionic term proportional to the gauge covariant vector current. In addition to the conventional fermion flow equations based on the operator D 2 , we also studied those based on / D 2 .
We found that the same qualitative renormalization properties hold for the nonminimally evolved fields as in the Yang-Mills gradient flow case. That is, any correlator -possibly containing evolved composite operators -consisting of n (s)NMGF evolved fermion fields χ and n (s)NMGF evolved fieldsχ and any number of (s)NMGF evolved gauge fields B µ with Wilsonian normalization is rendered finite by multiplicative renormalization with Z n χ , supplemented with the usual gauge coupling renormalization, where the different Z χ factors for the different flows are given by Furthermore, we have extended the nonminimal flows to the most general form compliant with the symmetries of the nonevolved theory. We have shown that in general Z χ depends on two constants a 1 , a 2 , and that there exists a particular combination of these constants for which the O(g 2 ) term vanishes. This implies that there exists a one-parameter family of flow systems for which there is a renormalization scheme in which the anomalous dimension of the evolved fermion field vanishes to all orders in the gauge coupling.
Additionally, we calculated E(t) up to next-to-leading order in the coupling using the generalized nonminimal gradient flow, of which the YMGF, NMGF and sNMGF are special cases. We found that out of these three different cases E(t) calculated with the sNMGF has the smallest dependence on N f . Also, for fixed N f = 0, t 2 E(t) sN M grows slower with t than t 2 E(t) N M , which in turn grows slower with t than the conventionally employed t 2 E(t) Y M . This behaviour might be attributed to the fact that the gauge fields are now driven towards the stationary points of the full action, instead of to the stationary points of the Yang-Mills action only.

Possible lattice applications
It may be interesting to see if the above mentioned properties of E(t) will persist beyond the pertubative regime. If so, the nonminimal gradient flows might have useful applications in lattice studies involving fermions.
First of all, the gradient flow derived scales t 0 and w 0 (see e.g. [9]) might have a smaller N f dependence when derived using the nonminimal flows. A much used definition for t 0 is and figure 8 seems to imply that when using the YMGF at N f = 8 this definition leads to a relatively small value for the dimensionless quantity t 0 Λ 2 QCD . The nonminimal flows will lead to a significantly larger value for t 0 Λ 2 QCD when using the same definition (7.2)more comparable to the pure Yang-Mills case [4] -possibly making them practically more useful in studies including fermions.
Secondly, the application of the generalized nonminimal gradient flow to the lattice will give the scales t 0 , w 0 a continuous dependence on the parameters a 1 , a 2 which are defining the flow. This might provide additional tests on lattice implementations of the gradient flow.

Acknowledgments
The author would like to thank Elisabetta Pallante for useful discussions and suggestions. This work has been supported by the foundation of Dutch research institutes (NWO-I).

A Notation and conventions
We work with anti-Hermitian generators T a , a = 1, ..., N 2 − 1, with real structure constants f abc . For the normalization of the trace we use where T (R) and d(R) are the index resp. dimension of repesentation R. We usually omit the subscript R, since we have (unless stated otherwise): with F denoting the (anti)fundamental representation. The quadratic casimir of representation R is related to the generators via and thus We work with Euclidean metric δ µν , δ µµ = d. The gamma matrices satisfy For the Euclidean space-and momentum integrals we use the abbreviations:

B Euclidean QCD-like action
The Euclidean gauge-fixed d-dimensional action of QCD-like theories with massless fermions and canonically normalized gauge fields is given by and

C Gauge transformations
We parametrize a gauge transformation by Λ ∈ SU (N ): Nonevolved case: t = 0 This is the well-known case, repeated here for convenience with our conventions: and the matter current j a µ T a =ψγ µ T a ψT a transforms as: Evolved case: t > 0 The modified flow equations are given in (2.2) and Λ(t, x) must satisfy Note that the boudary condition at t = 0 is unconstrained, and thus these transformations generalize the gauge symmetry of the SU (N ) theory to all flow times.
and it is then straightforward to verify that The same holds when ∆ and ← − ∆ are replaced by / ∆ and ← − / ∆, as in (4.2).

D BRST transformations
The BRST variations of the d-dimensional unrenormalized fields are given by where δ can be viewed as a nilpotent operator, which anticommutes with grassman variables, i.e. the (anti-)ghosts and (anti-)quarks.
The BRST variations of the bulk fields are given by δd = D µ L µ − g 0 {d,d} + g 0λ T a χT a − g 0χ T a λT a (D.13) The nontrivial variation ofd is chosen such that δS bulk = 0, and δ 2d = 0. Useful identities in checking δS bulk = 0 are obtained by defining The same identities hold when ∆ and ← − ∆ are replaced by / ∆ and ← − / ∆, as in section 4. In the Yang-Mills gradient flow case the g 0 j a µ T a term is absent in (D.14), but (D.18) will still hold.

YM1
YM2 YM3 YM4 YM5 YM8 YM6 YM7 Figure 10: Next-to-leading order diagrams for the YMGF evolved fermion 2-point function. White blobs represent GF vertices, black blobs the usual QCD vertices, and double lines indicate a kernel including integration over the associated flow time. Table 3: UV divergences of the diagrams from figure 10 in units of C 2 (R) E Z χ using the Yang-Mills gradient flow We present the results of the calculation of the evolved fermion renormalization factor Z χ using the evolved fermion 2-point function χ(t, x)χ(t, y) . We find Z χ by imposing: χ R (t, x)χ R (t, y) = Z χ χ(t, x)χ(t, y) = O(ε 0 ) (E.1) The diagrams contributing at O(g 2 0 ) are shown in figure 10. We collect the results in table 3, where S(x t −ȳ t ) is the evolved fermion propagator given in (3.27). Adding the results from table 3, we have at O(g 2 0 ) in d = 4 − 2ε: and therefore the renormalization factor for fermions evolved with the flow equation (2.8), while the gauge fields are evolved with the Yang-Mills gradient flow (2.1), is given by as is given in (2.22). This calculation is generalized to the case where the gauge fields are evolved with the nonminimal gradient flow from (3.1) in section 3.3, and the case where the fermions are evolved by the 'slashed' flow equation (4.2) together with the gauge fields evolved by the nonminimal gradient flow in section 4.2. The generalized case is presented in (6.5).