Nf=2+1 QCD thermodynamics with gradient flow using two-loop matching coefficients

We study thermodynamic properties of Nf=2+1 QCD on the lattice adopting O(a)-improved Wilson quark action and Iwasaki gauge action. To cope with the problems due to explicit violation of the Poincare and chiral symmetries, we apply the Small Flow-time eXpansion (SFtX) method based on the gradient flow, which is a general method to correctly calculate any renormalized observables on the lattice. In this method, the matching coefficients in front of operators in the small flow-time expansion are calculated by perturbation theory. In a previous study using one-loop matching coefficients, we found that the SFtX method works well for the equation of state, chiral condensates and susceptibilities. In this paper, we study the effect of two-loop matching coefficients by Harlander et al. We also test the influence of the renormalization scale in the SFtX method. We find that, by adopting the mu_0 renormalization scale of Harlander et al. instead of the conventional mu_d=1/sqrt{8t} scale, the linear behavior at large t is improved so that we can perform the t ->0 extrapolation of the SFtX method more confidently. In the calculation of the two-loop matching coefficients by Harlander et al., the equation of motion for quark fields was used. For the entropy density in which the equation of motion has no effects, we find that the results using the two-loop coefficients agree well with those using one-loop coefficients. On the other hand, for the trace anomaly which is affected by the equation of motion, we find discrepancies between the one- and two-loop results at high temperatures. By comparing the results of one-loop coefficients with and without using the equation of motion, the main origin of the discrepancies is suggested to be attributed to O((aT)^2)=O(1/N_t^2) discretization errors in the equation of motion at N_t =<10.

The GF modifies the fields according to a flow equation, which is given by the gradient At small flow-times, operators in the GF-scheme ("flowed operators") can be expanded in terms of operators at t = 0 in a conventional renormalization scheme, say the MSscheme [4]. By inverting the relation, we can also expand correctly renormalized physical observables in conventional schemes in terms of flowed operators at small t. Thanks to the asymptotic freedom of QCD around the small flow-time limit, we can calculate the matching coefficients relating the operators in two schemes by perturbation theory. Basic idea of the SFtX method is that, because the flowed operators are finite, we can evaluate their expectation values directly on the lattice without further renormalization [9]. We can thus extract the values of correctly renormalized physical observables by extrapolating proper combinations of expectation values in the GF-scheme to the small flow-time limit t → 0. In these extrapolations, the matching coefficients act not only to match the renormalization schemes but also to make the t-dependence milder by removing calculable part of operator mixings and t-dependences around the small flow-time limit. Note that the method is applicable also to observables whose founding symmetry is violated explicitly on the lattice, provided that the lattice theory has the correct continuum limit in which the symmetry is restored.
The SFtX method was first tested in quenched QCD to calculate the energy-momentum tensor (EMT) [13][14][15][16], which has not been easy to evaluate on the lattice due to the explicit violation of the Poincaré invariance by the lattice regularization. 1 It was found that the equation of state (EoS) calculated from the diagonal components of the EMT correctly reproduces previous estimation using the conventional integral methods [19][20][21][22]. The SFtX method was tested successfully also in solvable models [23,24].
We note that the SFtX method is applicable also to chiral observables [11,12]. We thus apply the method to QCD with dynamical Wilson-type quarks, with which the correct continuum limit is guaranteed, to cope with the problems due to explicit violation of the chiral symmetry on the lattice [25][26][27][28]. To reduce the finite lattice spacing effects, we adopt the renormalization-group improved Iwasaki gauge action [29,30]

and the O(a)-improved
Wilson quark action [31] with a non-perturbatively estimated clover coefficient using the Schrödinger functional method [32].
In our previous study of (2+1)-flavor QCD with slightly heavy u and d quarks (m π /m ρ 0.63) and approximately physical s quark (m ηss /m φ 0.74), we calculated the EMT as well as chiral condensates and disconnected chiral susceptibilities in the temperature range of T 174-697 MeV (N t = 16-4, where N t is the lattice size in the temporal direction) [25].
The lattices are relatively fine with the lattice spacing a 0.07 fm. Adopting one-loop matching coefficients calculated in Refs. [10,12], we found that the EoS extracted from diagonal components of EMT by the SFtX method is well consistent with that estimated with the conventional T -integration method at T < ∼ 280 MeV (N t > ∼ 10) [33]. At the same time, the two estimates of EoS deviate at T > ∼ 348 MeV, suggesting contamination of aindependent lattice artifacts of O((aT ) 2 = 1/N 2 t ) at N t < ∼ 8 in the EMT. We also found that the chiral condensates bend sharply and the disconnected chiral susceptibilities show peak at T 190 MeV which was suggested as the pseudo-critical temperature from other observables [33]. We have further studied topological properties of QCD by the SFtX method on these lattices [26]. We found that the topological susceptibilities estimated with the gluonic and fermionic definitions agree well with each other at T < ∼ 279 MeV even at finite lattice spacing of a = 0.07 fm. This is in clear contrast to their conventional lattice estimations: For example, a study with improved staggered quarks reports more than hundred times larger gluonic susceptibility than fermionic one at similar lattice spacings [34]. These suggest that the SFtX method is powerful in calculating observables from lattice simulations.
Recently, Harlander, Kluth, and Lange have completed the calculation of the matching 1 For a recent development in lattice determination of the EMT, see Ref. [17,18]. coefficients for EMT up to the two-loop order [35]. Some details of their calculation are given in [36]. Removing more known small-t behaviors, we may perhaps expect milder t-dependence in the t → 0 extrapolation. The two-loop coefficients were first tested in quenched QCD [15]. It was found that the results of EoS with one-and two-loop coefficients are well consistent with each other. It was also noted that the two-loop coefficients lead to a milder t-dependence such that systematic errors from the t → 0 extrapolation are reduced.
In this paper, we extend the test of two-loop coefficients to QCD with (2 + 1)-flavors of dynamical quarks. The lattice setup is the same as in Ref. [25]. A point to be noted here is that, unlike the one-loop coefficients of Ref. [10], in the calculation of two-loop coefficients of Ref. [35], the equation of motion (EoM) in the continuum, is used for quark operators, where In Appendix A, we define the group factors appearing in perturbative expressions of the matching coefficients. In Appendix B, we confirm that the one-loop coefficients of Ref. [10] are consistent with the one-loop part of Ref. [35].

II. THE SFtX METHOD
A. Gradient flow Our flow equations are identical to those given in Refs. [3] and [5]. That is, for the gauge field, we set 2 where the field strength and the covariant derivative of the flowed gauge field are and respectively. For the quark fields, we set where f = u, d, s, denotes the flavor index, and Note that our flow equations are independent of the flavor.

B. Energy-momentum tensor with two-loop matching coefficients
In terms of unflowed operators, the EMT under the dimensional regularization is given by 2 In what follows, summations are always understood over repeated Lorentz indices, µ, ν, . . . , and adjoint indices, a, b, . . . . On the other hand, without stated otherwise, summation over repeated flavor indices, with F a µν (x) the field strength of unflowed gauge field A a µ (x).
Here and in what follows, we assume for notational simplicity that the vacuum expectation value (VEV), i.e., the expectation value at zero-temperature, is subtracted in each operator.
In terms of flowed operators, the EMT is expressed as [10,25] whereχ are "ringed" quark fields introduced in Ref. [10] to carry out a wave function renormalization of quark fields non-perturbatively, where the operators in the denominator are not VEVsubtracted. The O(t) term in the right-hand side of Eq. (20) can be removed by a t → 0 extrapolation. Explicit forms of the matching coefficients c 1 , · · · , c 5 to the one-loop order (NLO) are given in Ref. [10].
Recently, two-loop (NNLO) calculation of the matching coefficients has been completed by Harlander, Kluth, and Lange [35]. Unlike the calculation of Ref. [10], the EoM (1) was used in Ref. [35] assuming that the EMT operators are spatially separated from other composite operators. See also Sec. 5 of Ref. [10]. In terms of unflowed operators, Eq. (1) reads as Note that Eq. (23) implies that the last two terms of Eq. (9) cancel with each other: In terms of the flowed operators, the EoM is where, in the one-loop level, the coefficients d i (t) are given by (cf. Eqs. (5.3)-(5.5) of Ref. [10]) where g is the running coupling in the MS-scheme at the renormalization scale µ, and with γ E the Euler-Mascheroni constant. The definitions of the group factors (T F etc.) are summarized in Appendix A.
Using Eq. (25), we can eliminateÕ 5,µν (t, x) order by order in perturbation theory (up to The matching coefficientsč 1 (t) andč 2 (t) to the two-loop order are given in Ref. [35] aš where Li 2 (z) is the dilogarithm function with Li 2 (1/4) = 0.26765263908 · · · , and β 0 and β 1 are the first two coefficients of the beta function, The matching coefficientsc 3 (t) andc 4 (t) are given bẙ usinǧ c 3 (t) andč given in Ref. [35], where Here, ζ χ (t) adjusts the normalization of quark fields to that of the ringed variables in Eqs. (21) and (22), and is given by with Its inverse reads as We then obtaiň wherec In Appendix B, we confirm that the results of Ref. [10] are consistent with the one-loop parts of Eqs. (31)-(44).
From the diagonal components of the EMT, we compute the pressure and the energy density as The entropy density and trace anomaly are then computed as ( + p)/T and − 3p, respectively.
C. Extrapolation to t → 0 To extract physical results of EMT in Eqs. (20), (30), etc., we remove contamination of O(t) terms in the right-hand side of these equations by extrapolating them to t → 0. On finite lattices, lattice artifacts contaminate additionally. With O(a)-improved Wilson quarks we adopt, the lattice artifacts start with O(a 2 ) and we expect the EMT, for example, to be where T µν (x) is the physical EMT, S µν and S µν are contaminations of dimension-six operators with the same quantum number, and A µν , B f µν , C µν , and D µν are those from dimension-four operators. Though both a → 0 and t → 0 extrapolations are needed to extract the physical EMT, it is often attractive to reserve the a → 0 extrapolation for a late stage of numerical analyses. To perform a t → 0 extrapolation on finite lattices a = 0, the influence of singular terms such as a 2 /t must be suppressed. This will be possible when we have a window in t in which the linear terms dominate ("linear window") [25].
We found in the study of Ref. [25] that, depending on the observable and simulation parameters, we do have ranges of t in which the data show well linear behavior. We performed linear t → 0 extrapolation of observables when a linear window is available and obtained reasonable results, as introduced in Sec. I. We think that the success of the SFtX method in Ref. [25] is largely due to the fineness of the lattices studied. In this paper, we adopt the same strategy. We also discuss a method which may be used to improve linear behaviors in Sec. IV.
We identify linear windows as follows: First of all, we require the flow-time to satisfy √ 2a ≤ √ 8t ≤ min(N t a/2, N s a/2), i.e., the smearing range √ 8t by the gradient flow covers the minimal lattice separations to make the smearing effective, and, simultaneously, is smaller than the half of the smallest lattice extent to avoid finite-size effects due to overlapped smearing. In terms of the dimensionless flow-time t/a 2 , these conditions read 3 For each observable O, we then look for a range of t in which terms linear in t look dominating, and try linear extrapolations of the form with various choices of the fitting range. We then select a temporally best linear fit whose fitting range is the widest within the range (47), under the condition that χ 2 /N dof is smaller than a cutoff value. In this study, due to limitation of the statistics, we disregard correlations among data at different t. Thus the absolute value of χ 2 /N dof does not have a strong sense.
We thus repeat the test varying the cutoff value widely. In Fig. 1, we show some typical results of this test. One-loop matching coefficients of Ref. [10] and the µ 0 -scale discussed in Sec. IV are adopted. The cutoffs for the linear fits 1, 2, 3, · · · , and 8 are 10 −5 , 10 −4 , 3 In practice, t/a 2 is bounded also by the maximum value of t/a 2 in the calculation of flowed fields. In this study, we calculate them up to t/a 2 = 2.0. 10 −3 , · · · , and 10 +2 , respectively. The selected fitting range for each fit is shown by a line at the bottom of each plot with the same color. From Fig. 1, we note that the linear fits are stable when the cutoff value is large. On the other hand, when we require χ 2 /N dof < 10 −3 or smaller, our selection procedure for the temporal linear window becomes sometimes unstable and fails to give a window. Consulting these plots and also requiring that the resulting linear windows are common among similar observables, we decide to choose the fits 5 which require χ 2 /N dof < 0.1 to select optimum linear windows, for all observables we study in this paper.
To confirm the validity of the linear window and to estimate a systematic error due to the fit ansatz, we also make additional fits using the data within the same window: One is a "nonlinear fit" inspired from Eq. (46), Another is a "linear+log fit" to estimate the effects of neglected higher-order loop corrections in the matching coefficients. For the case of one-loop matching coefficients, possible correc- of O(g 4 ) terms may be taken by fits of the form For the case of two-loop matching coefficients, we instead try to estimate the O(g 6 ) contaminations. We take the difference between the linear fit and the nonlinear or linear+log fits as an estimate of the systematic error due to the choice of the fit ansatz.

III. SIMULATION PARAMETERS AND OUR NUMERICAL METHODS
The numerical setup for this study is the same as that of Ref. [25]. We study 2 + 1 flavor QCD with slightly heavy u and d quarks (m π /m ρ 0.63) and approximately physical s quark (m ηss /m φ 0.74) on a relatively fine lattice with the lattice spacing a 0.07 fm (a −1 2.79 GeV) [33,37]. To reduce the finite lattice spacing effects, we adopt the renormalization-group improved Iwasaki gauge action [29,30] and the O(a)-improved Wilson quark action [31]. For the clover coefficients of the improved Wilson quark action, we adopt non-perturbatively evaluated values [32]. The bare gauge coupling parameter and the hopping parameters are set to β = 2.05, κ ud = 0.1356, and κ s = 0.1351.
Finite temperature configurations in the range of T 174-348 MeV (N t = 16-8) have been generated adopting the fixed-scale approach [21,33]. The values of temperature at each N t are given in Table I. The spatial box size is 32 3 for T > 0 and 28 3 for T = 0. Although we have configurations also at T 464 MeV (N t = 6) and T 697 MeV (N t = 4) [25], limitations by t 1/2 = 1.125 and 0.5, respectively, are too severe to obtain a stable linear window. It was also noted that the EoS on N t < ∼ 8 lattices has large O((aT ) 2 = 1/N 2 t ) lattice artifacts [25]. We thus just omit these configurations in this study.  Our numerical algorithm for gradient flow is given in Ref. [25]. We adopt the third order Runge-Kutta method with the step size of = 0.02. For the quadratic terms of the field strength tensor G µν (x), we adopt clover operator with four plaquette Wilson loops and that with eight 1 × 2 rectangle Wilson loops such that the tree-level improved field strength squared is obtained [39]. To calculate the MS running coupling and masses in the matching coefficients, we adopt five-loop beta and gamma functions [38], instead of the four-loop functions adopted in our previous study [25].
We evaluate fermionic observables by the noisy estimator method. The number of noise vectors we adopt is 20 for each color and spinor component. To measure fermionic bilinear observables at t > 0, instead of the adjoint Runge-Kutta integration adopted in Ref. [25] (see Appendix B of Ref. [25]), we apply an alternative method using usual forward Runge-Kutta integrations only by locating the noise vectors at t = 0. See Appendix C for details. We confirm that the both methods give consistent results within statistical errors of the noise method. This alternative method is applicable to fermionic bi-linear observables and reduces the computational cost. On the other hand, because the data at all t's are estimated with the same noise vector at t = 0, the correlation among different t's is stronger than the study of Ref. [25] in which independent noise vectors were generated at each t. A conventional choice of µ is which is a natural scale of flowed operators because √ 8t is the physical smearing extent of flowed fields. On the other hand, the authors of Ref. [35] argued that the choice We note that, since the µ 0 -scale is more perturbative than the µ d -scale in asymptotically free theories. Thus, the µ 0 -scale may improve the quality of perturbative expressions, in particular at large t. 4

A. Test of renormalization scales using one-loop matching coefficients
In Fig. 2, we examine effects of the renormalization scale by comparing the results for the entropy density, at finite t with µ 0 -and µ d -scales. Corresponding results for the trace anomaly, − 3p are shown in Fig. 3. We adopt one-loop matching coefficients of Ref. [10] in this test. In Table II, Table III. We find that the results with the µ d -scale are consistent within the errors with both those given in Table II and also in Ref. [25]. a too big value for µ will make L(µ, t) large and thus may invalidate the perturbative expansion, i.e., µ should be O(µ d ).
where connected quark loop contribution is dropped from the scalar density two-point function and N Γ is the lattice volume. Here, because the VEV-subtraction is not required for χ disc. f f , we also study the case of T = 0. We find that, though a linear window is sometimes not clear with the conventional µ d -scale, the linear behavior is much improved by adopting the µ 0 -scale in particular at large t/a 2 . The µ 0 -scale extends the applicability of the t → 0 extrapolation method based on the linear window. In the followings, we perform t → 0 extrapolations with the µ 0 -scale TABLE IV. Results for chiral condensates and disconnected chiral susceptibilities by the SFtX method with the µ 0 -scale using the one-loop matching coefficients of Ref. [10]. The chiral condensates are in unit of GeV 3 , and the disconnected chiral susceptibilities are in unit of GeV 6 . The first parenthesis is for the statistical error, and the second for the systematic error due to the fit ansatz.
only. 5 In Figs. 4, 5, 6, and 7, the arrow at the bottom of each plot is the linear window we adopt, and the symbols at t ∼ 0 shows the results of t → 0 extrapolations, from the right to the left, using linear, nonlinear, and linear+log fits, respectively. In Table IV, we summarize the final results for the chiral condensates and disconnected chiral susceptibilities, obtained in the t → 0 limit using the one-loop matching coefficients of Ref. [10] with the µ 0 -scale. 5 Similar but more drastic improvement with the µ 0 -scale was observed in our preliminary study of 2 + 1 flavor QCD at the physical point on a less fine lattice [40]. On the other hand, no apparent improvement with the µ 0 -scale was reported in the study of quenched QCD [15]. This may be understood as follows: Because T c in quenched QCD is higher than T pc in full QCD, the effective coupling in quenched QCD is smaller at similar T /T c/pc , and thus the relevant range of t is well perturbative already with the µ d -scale.
We also note that T pc in full QCD decreases as we decrease the quark mass towards the physical point.  Black dots in Fig. 8 show the results of EoS obtained previously by the conventional T -integration method using the same configurations [33]. Our conclusions are the same as Ref. [25]: At T < ∼ 279 MeV (N t ≥ 10), the SFtX method leads to EoS which is well consistent with the result of the conventional method, while a-independent lattice artifacts of O((aT ) 2 = 1/N 2 t ) are suggested for N t < ∼ 8. Because the continuum extrapolation is not taken yet, the good agreement of different estimations at N t ≥ 10 suggests that the remaining O(a 2 T 2 , a 2 m 2 , a 2 Λ 2 QCD ) lattice artifacts are small with our lattice action at a 0.07 fm. Results for the VEV-subtracted chiral condensates are shown in the left panel of Fig. 9.
In the right panel of Fig. 9, we show the results for disconnected chiral susceptibilities as function of temperature. Because the VEV-subtraction has no effects in this quantity, we also show the results at T = 0. We find a clear peak at T 199 MeV, which may be indicating the pseudo-critical point around T pc ∼ 190 MeV previously suggested using the Polyakov loop etc. [33]. We also note that, although the errors are large yet, the height of the peak looks increasing as we decrease the valence quark mass from s quark to u (or d) quark.

V. TEST OF TWO-LOOP MATCHING COEFFICIENTS
We now test the effects of two-loop matching coefficients for EMT by Harlander et al. [35] in QCD with 2+1 flavors of dynamical quarks. Following the discussion in Sec. IV, we adopt the µ 0 -scale in this test. As mentioned in Sec. I, unlike the one-loop coefficients of Ref. [10], the EoM is used in the two-loop coefficients of Refs. [35].

A. Entropy density
In Fig. 10, we compare the results of the entropy density at finite t using the one-loop matching coefficients of Ref. [10] (blue squares) and those using the two-loop coefficients of Ref. [35] (red circles). Note that, because the contribution of the EoM to the EMT given by Eq. (23) is proportional to δ µν , only the trace part of the EMT is affected by the use of the EoM. Thus, the EoM has no effects in the entropy density which is a trace-less combination of the EMT. We find that the entropy density with the two-loop coefficients is larger than its one-loop value at finite t, but the difference becomes smaller in the t → 0 limit.
Physical results for the entropy density with one-and two-loop matching coefficients are shown in Fig. 11. The errors include the systematic error due to the t → 0 extrapolation. Here, it should be recalled that EMT data at T > ∼ 348 MeV are contaminated with O((aT ) 2 = 1/N 2 t ) lattice artifacts at N t < ∼ 8 [25]. We find that one-and two-loop results agree well within the errors at T > ∼ 279 MeV (N t < ∼ 10).

B. Trace anomaly
In Fig. 12, we show the results for the trace anomaly ( − 3p)/T 4 as function of the flow-time. The trace anomaly is just the trace part of the EMT and thus will be sensitively affected by the EoM on finite lattices. In order to identify the effects of EoM clearly, we also compute the trace anomaly using the one-loop part of the matching coefficients of Ref. [35] in which the EoM is used. We find that one-and two-loop results both using the EoM are close with each other, while the one-loop results without using the EoM deviates from the results using the EoM. We thus conclude that the deviation is mainly due to the use of the EoM. The deviation increases with increasing T (decreasing N t ) and becomes sizable at high temperatures, T > ∼ 279 MeV (N t < ∼ 10). the EoM are compared with results using two-loop matching coefficients of Ref. [35] in which the EoM is used. Black dots are the results of the T -integration method [33]. Errors include systematic error due to the fit ansatz for the t → 0 extrapolation. The symbols are slightly shifted horizontally to avoid overlapping.
Physical results of the trace anomaly extracted by the t → 0 extrapolation are summarized in Fig. 13. We find that though the analyses with and without using the EoM lead to consistent results within errors at low temperatures, they show visible discrepancy at high temperatures, T > ∼ 279 MeV (N t < ∼ 10). Even with disregarding the data at T > ∼ 348 MeV lattice artifacts is suggested in EMT, we see discrepancy at T 279 MeV (N t = 10). The O((aT ) 2 ) lattice artifacts will contaminate the EoM too. Our results suggest that the EoM suffers from larger O((aT ) 2 = 1/N 2 t ) discretization errors than the EMT, and has visible effect at N t < ∼ 10. the EoM are compared with results using two-loop matching coefficients of Ref. [35] in which the EoM is used. Also shown are the results using one-loop part of the matching coefficients of Ref. [35] in which the EoM is used. Black dots are the results of the T -integration method [33].
Errors include systematic error due to the fit ansatz for the t → 0 extrapolation. The symbols are slightly shifted horizontally to avoid overlapping.

VI. CONCLUSIONS
We presented the results of our tests of the µ 0 = e −γ E /2 / √ 2t renormalization scale and two-loop matching coefficients recently calculated by Harlander et al. [35]. For this test, we revisited the case of QCD with heavy u and d quarks [25].
We find that, comparing with the results using the conventional µ d = 1/ √ 8t scale, the µ 0 -scale improves the quality of perturbative expressions, in particular at large t, and thus leads to clearer and wider linear windows so that we can carry out t → 0 extrapolations much confidently. We also find that, for observables for which the linear window is clear with the conventional µ d -scale, the results using µ 0 -and µ d -scales are consistent with each other, i.e., the results extrapolated to the t → 0 limit are insensitive to the choice of the renormalization scale, as expected.
Concerning the test of two-loop matching coefficients, unlike the case of the one-loop matching coefficients of Ref. [10], the equation of motion for quark fields in the continuum limit is used by Harlander et al. in their calculation of the two-loop matching coefficients [35]. We are attempting to extend applications of the SFtX method in various directions: thermodynamics of 2+ 1 flavor QCD at the physical point [27,40], shear and bulk viscosities from two-point correlation functions of the energy-momentum tensor [28], end-point of firstorder deconfining transition region in QCD near the quenched limit [41], PCAC quark masses [42], etc. The choice of the µ 0 -scale as well as higher order coefficients may help improving these calculations too. We normalize the gauge group generators as The structure constant defined by [T a , T b ] = f abc T c has quadratic Casimirs as In the perturbative expressions in Sec. II B, group factors are defined by For G = SU (N ) and R = N , In particular, for the N f = 2 + 1 QCD, and for the quenched QCD (the SU (3) pure Yang-Mills) these coefficientsc i (t) in Eq. (B2) to the coefficients in Eq. (30), we have to eliminate the operatorÕ 5,µν (t, x) from Eq. (B2) in favor ofÕ 4,µν (t, x) andÕ 2,µν (t, x) by using the relation (25). In the present order of approximation, this is easy because we can use the relation 1 2Õ 4,µν (t, x) +Õ 5.µν (t, x) = 0 (B9) that holds in the tree-level in Eq. (B2) becausec 4 (t) andc 5 (t) are already one-loop quantities.

Two-loop confirmation
As pointed out in Ref. [9], in the case of the pure Yang-Mills theory, a certain part of the two-loop order coefficients can be extracted by using the trace anomaly without any higher order calculations. The expression in our present notation is (see Eq. (4.69) of Ref. [10]), where the second term in the right-hand side contains the information in the two-loop order.
This coincides with the result obtained from Eqs. (31) and (32) for the pure Yang-Mills theory.
We can similarly deduce three-loop c 2 (t) for the pure Yang-Mills theory from the two-loop coefficients. A concrete form is given in Ref. [15].
Appendix C: Alternative method for flowed fermionic bi-linear observables As discussed in Appendix A of Ref. [25], fermionic parts of the EMT are given in terms and where the covariant derivatives in Eq. (C1) refer to the flowed gauge field B µ (t, x), N Γ = x is the number of lattice points, α and i denote the spinor and color indices, respectively, and c fl is an improvement coefficient associated with the flowed quark field [5]. Here, we note that Eqs. (C1) and (C2) can be equally written as where we have used the relation S r (v, w) = γ 5 S r (w, v) † γ 5 , and s r (t) = − 1 N Γ x,y,v,w α,i K(t, x; 0, y) † K(t, x; 0, v) [S r (v, w) − c fl δ v,w ] αi,αi δ w,y . (C5) We thus introduce a new noise field η αi (x) η = 0, η αi (x)η βj (y) * η = δ αβ δ ij δ x,y , and define Ξ(t, x) ≡ y K(t, x; 0, y) η(y), (C7) where contraction of spinor and color indices is understood. We then obtain compact expressions as where the gauge field in the covariant derivative D ν is the flowed B µ (t, x), and The building blocks obey the forward flow equations as which can be solved by the time-forward Runge-Kutta method as explained in Appendix D.2 of Ref. [5]. That is, setting ∂ t χ t = ∆(V t )χ t , the Runge-Kutta proceeds as where ∆ i = ∆(W i ), i = 0, 1, 2, and These new representations are advantageous in the sense that they do not require the Runge-Kutta steps proceeding backward in the flow time (see Appendix B of Ref. [25]) which is numerically demanding.