Non-perturbative renormalization of tensor currents: strategy and results for $N_f = 0$ and $N_f = 2$ QCD

Tensor currents are the only quark bilinear operators lacking a non-perturbative determination of their renormalisation group (RG) running between hadronic and electroweak scales. We develop the setup to carry out the computation in lattice QCD via standard recursive finite-size scaling techniques, and provide results for the RG running of tensor currents in $N_f = 0$ and $N_f = 2$ QCD in the continuum for various Schr\"odinger Functional schemes. The matching factors between bare and renormalisation group invariant currents are also determined for a range of values of the lattice spacing relevant for large-volume simulations, thus enabling a fully non-perturbative renormalization of physical amplitudes mediated by tensor currents.


Introduction
Hadronic matrix elements of tensor currents play an important rôle in several relevant problems in particle physics. Some prominent examples are rare heavy meson decays that allow to probe the consistency of the Standard Model (SM) flavour sector (see, e.g., [1][2][3] for an overview), or precision measurements of β-decay and limits on the neutron electric dipole moment (see, e.g., [4][5][6] for an up-to-date lattice-QCD perspective).
One of the key ingredients in these computations is the renormalization of the current. Indeed, partial current conservation ensures that non-singlet vector and axial currents require at worst finite normalizations, and fixes the anomalous dimension of scalar and pseudoscalar densities to be minus the quark mass anomalous dimension. They however do not constrain the tensor current, which runs with the only other independent anomalous dimension among quark bilinears. Controlling the current renormalization and running at the non-perturbative level, in the same fashion achieved for quark masses [7][8][9][10], is therefore necessary in order to control systematic uncertainties, and allow for solid conclusions in new physics searches.
The anomalous dimension of tensor currents is known to three-loop order in continuum schemes [11,12], while on the lattice perturbative studies have been carried out to two-loop order [13]. Non-perturbative determinations of renormalization constants in RI/MOM schemes, for the typical few-GeV values of the renormalization scale accessible to the latter, have been obtained for various numbers of dynamical flavours and lattice actions [14][15][16][17][18][19][20]. The purpose of this work is to set up the strategy for the application of finite-size scaling techniques based on the Schrödinger Functional (SF) [21], in order to obtain a fully non-perturbative determination of both current renormalization constants at hadronic energy scales, and the running of renormalized currents to the electroweak scale. This completes the ALPHA Collaboration non-perturbative renormalization programme for non-singlet quark field bilinears [7][8][9][10][22][23][24] and four-quark operators [25][26][27][28][29][30][31].
As part of the strategy, we will set up a family of SF renormalization schemes, and perform a perturbative study with the main purpose of computing the perturbative anomalous dimension up to two loops, in order to make safe contact with perturbative physics at the electroweak scale. Preliminary results of this work have already appeared as proceedings contributions [32]. 1 We will then apply our formalism to the fully non-perturbative renormalization of non-singlet tensor currents in N f = 0 and N f = 2 QCD. Our results for N f = 3 QCD, that build on the nonperturbative determination of the running coupling [34][35][36] and the renormalization of quark masses [9,10,22], will be provided in a separate publication [37].
The layout of this paper is as follows. In section 2 we will introduce our notation and discuss the relevant renormalization group equations. In section 3 we will introduce our SF schemes, generalizing the ones employed for quark mass renormalization. In section 4 we will study these schemes in one-loop perturbation theory, and compute the matching factors that allow to determine the NLO values of anomalous dimensions. In section 5 we will discuss our non-perturbative computations, and provide results for the running of the currents between hadronic and high-energy scales and for the renormalization constants needed to match bare hadronic observables at low energies. Section 6 contains our conclusions. Some technical material, as well as several tables and figures, are gathered in appendices.

Renormalization Group
Theory parameters and operators are renormalized at the renormalization scale µ. The scale dependence of these quantities is given by their Renormalization Group (RG) evolution. The Callan-Symanzik equations satisfied by the gauge coupling and quark masses are of the form where g 0 is the bare coupling, Z O is a renormalization constant, and a is some inverse ultraviolet cutoff -the lattice spacing in this work. We assume a massindependent scheme, such that both the β-function and the anomalous dimensions τ and γ depend only on the coupling and the number of flavours (other than on the number of colours N ); examples of such schemes are the MS scheme of dimensional regularization [38,39], RI schemes [40], or the SF schemes we will use to determine the running non-perturbatively [21,41]. The RG functions then admit asymptotic expansions of the form: −g 2 d 0 + d 1 g 2 + d 2 g 4 + . . . , (2.6) γ(g) ≈ g∼0 −g 2 γ 0 + γ 1 g 2 + γ 2 g 4 + . . . . (2.7) The coefficients b 0 , b 1 and d 0 , γ 0 are independent of the renormalization scheme chosen. In particular [42][43][44][45][46][47][48] b 0 = 1 (4π) 2 and While the value of the Λ QCD parameter depends on the renormalization scheme chosen, M i andÔ are the same for all schemes. In this sense, they can be regarded as meaningful physical quantities, as opposed to their scale-dependent counterparts. The aim of the non-perturbative determination of the RG running of parameters and operators is to connect the RGIs -or, equivalently, the quantity renormalized at a very high energy scale, where perturbation theory applies -to the bare parameters or operator insertions, computed in the hadronic energy regime. In this way the three-orders-of-magnitude leap between the hadronic and weak scales can be bridged without significant uncertainties related to the use of perturbation theory. In order to describe non-perturbatively the scale dependence of the gauge coupling and composite operators, we will use the step-scaling functions (SSFs) σ and σ O , respectively, defined as 15) or, equivalently, is the RG evolution operator for the operator at hand, which connects renormalized operators at different scales as O (µ 2 ) = U (µ 2 , µ 1 )O (µ 1 ). The SSFs are thus completely determined by, and contain the same information as, the RG functions γ and β. In particular, σ O (s, u) corresponds to the RG evolution operator of O between the scales µ/s and µ; from now on, we will set s = 2, and drop the parameter s in the dependence. The SSF can be related to renormalization constants via the identity .

(2.19)
This will be the expression we will employ in practice to determine σ O , and hence operator anomalous dimensions, for a broad range of values of the renormalized coupling u. In this work we will focus on the renormalization of tensor currents. The (flavour non-singlet) tensor bilinear is defined as where σ µν = i 2 [γ µ , γ ν ], and s 1 = s 2 are flavour indices. Since all the Lorentz components have the same anomalous dimension, as far as renormalization is concerned it is enough to consider the "electric" operator T 0k . As already done in the introduction, it is important to observe that the tensor current is the only bilinear operator that evolves under RG transformation in a different way respect to the quark mass -partial conservation of the vector and axial currents protect them from renormalization, and fixes the anomalous dimension of both scalar and pseudoscalar densities to be −τ . The one-loop (universal) coefficient of the tensor anomalous dimension is (2.21)

Schrödinger Functional renormalization schemes
The renormalization schemes we will consider are based on the Schrödinger Functional [21], i.e. on the QCD partition function Z = D[A,ψ, ψ]e −S[A,ψ,ψ] on a finite Euclidean spacetime of dimensions L 3 × T with lattice spacing a, where periodic boundary conditions on space (in the case of fermion fields, up a to a global phase θ) and Dirichlet boundary conditions at times x 0 = 0, T are imposed. A detailed discussion of the implementation and notation that we will follow can be found in [52]. We will always consider L = T and trivial gauge boundary fields (i.e. there is no background field induced by the latter). The main advantage of SF schemes is that they allow to compute the scale evolution via finite-size scaling, based on the identification of the renormalization scale with the inverse box size, i.e. µ = 1/L. To define suitable SF renormalization conditions we can follow the same strategy as in [8,[53][54][55], which has been applied successfully also to several other composite operators both in QCD [23][24][25][26][27][28][29][56][57][58] and other theories. 3 We first introduce the two-point functions and 3 See, e.g., [59] for a recent review. where x,yζ is a source operator built with the x 0 = 0 boundary fields ζ,ζ. A sketch of the correlation function in the SF is provided in Fig.1. The renormalization constant Z T is then defined by where we have already fixed µ = 1/L, m 0 is the bare quark mass, and m cr is the critical mass, needed if Wilson fermions are used in the computation -as will be our case. The factor f 1/2−α 1 k α 1 cancels the renormalization of the boundary fields contained in O[Γ], which holds for any value of the parameter α; we will restrict ourselves to the choices α = 0, 1/2. The only remaining parameter in Eq. (3.5) is the kinematical variable θ entering spatial boundary conditions; once its value is specified alongside the one of α, the scheme is completely fixed. We will consider the values θ = 0, 0.5, 1.0 in the perturbative study discussed in the next section, and in the non-perturbative computation we will set θ = 0.5.
The condition in Eq. (3.5) involves the correlation function k T , which is not O(a) improved. Therefore, the scaling of the renormalized current towards the continuum limit, given by Eq. (2.4), will be affected by O(a) effects. The latter can be removed by subtracting suitable counterterms, following the standard on-shell O(a) improvement strategy for SF correlation functions [52]. On the lattice, and in the chiral limit, the O(a) improvement pattern of the tensor current reads where∂ is the symmetrized lattice derivative and V µ =ψ s 1 γ µ ψ s 2 is the vector current. Focusing again only on the electric part, the above formula reduces to which results in an O(a) improved version of the two-point function k T of the form Note that the contribution involving the spatial derivative vanishes. Inserting k I T in Eq. (3.5), and the resulting Z T in Eq. (2.4) alongside the O(a) improved current, will result in O(a 2 ) residual cutoff effects in the value of the SSF Σ T defined in Eq. (2.19), provided the action and m cr are also O(a) improved.

Perturbative study
We will now study our renormalization conditions in one-loop perturbation theory. The aim is to obtain the next-to-leading (NLO) anomalous dimension of the tensor current in our SF schemes, necessary for a precise connection to RGI currents, or continuum schemes, at high energies; and compute the leading perturbative contribution to cutoff effects, useful to better control continuum limit extrapolations.
We will expand the relevant quantities in powers of the bare coupling g 2 0 as where X can be any of Z T , k T , k V , f 1 , or k 1 . To O(g 2 0 ), Eq. (3.7) can be written as . The renormalization constant for the improved tensor correlator k I T at one-loop is then given by  wherec t is the coefficient of the counterterm that subtracts the O(a) contribution coming from the fermionic action at the boundaries, and am (1) cr is the one-loop value of the critical mass, for which we employ the continuum values of am (1) cr from [26,60]. The one-loop value of the improvement coefficient c T has been obtained using SF techniques in [61]. We have repeated the computation of this latter quantity as a crosscheck of our perturbative setup; a summary is provided in Appendix A.

Perturbative scheme matching
Any two mass-independent renormalization schemes (indicated by primed and unprimed quantites, respectively) can be related by a finite parameter and operator renormalization of the form where we have assumed O to be multiplicatively renormalizable. The scheme change factors χ can be expanded perturbatively as Plugging Eqs. (4.4, 4.5, 4.6) into the Callan-Symanzik equations allows to relate a change in a renormalized quantity to the change in the corresponding RG function, viz. . (4.10) In particular, expanding Eq. (4.10) to order g 2 provides a useful relation between the 2-loop coefficient of the anomalous dimension in the two schemes, viz.
The one-loop matching coefficient χ (1) g for the SF coupling was computed in [62,63], where the logarithm vanishes with our choice µ = 1/L, and for the standard definition of the SF coupling one has The other finite term χ O in Eq. (4.10) will provide the operator matching between the lattice-regulated SF scheme and some reference scheme where the NLO anomalous dimension is known, such as MS or RI, that we will label as "cont". The latter usually are based on variants of the dimensional regularization procedure; our SF schemes will be, on the other hand, regulated by a lattice. The practical application of Eq. (4.11) thus involves a two-step procedure, in which the lattice-regulated SF scheme is first matched to a lattice-regulated continuum scheme, that is in turned matched to the dimensionally-regulated continuum scheme. This yields (4.14) The one-loop matching coefficients [χ (1) O ] cont;lat that we need can be extracted from the literature [13,64,65], while the term [χ (1) O ] SF;lat is obtained from our one-loop calculation of renormalization constants. Indeed, the asymptotic expansion for the one-loop coefficient of a renormalization constant in powers and logarithms of the lattice spacing a has the form  Our results for [χ (1) 0 ] SF;lat have been obtained by computing the one-loop renormalization constants on a series of lattices of sizes ranging from L/a = 4 to L/a = 48, and fitting the results to Eq. (4.15) to extract the expansion coefficients. The computation has been carried out with O(a) improved fermions for three values of θ for each scheme, and without O(a) improvement for θ = 0.5, which allows for a crosscheck of our computation and of the robustness of the continuum limit (see below). The results for the matching factors are provided in Table 1; details about the fitting procedure and the assignment of uncertainties are discussed in Appendix B.
Inserting our results in Eq. (4.11), we computed for the first time the NLO anomalous dimension in our family of SF schemes for the tensor currents, which are given in Table 2. We have crosschecked the computation by performing the matching with and without O(a) improvement, and proceeding through both MS and RI as reference continuum schemes, obtaining the same results in all cases. In this context we observe that the NLO correction to the running is in general fairly large. It is also worth mentioning that the choice of θ = 0.5, which leads to a close-to-minimal value of the NLO mass anomalous dimension in SF schemes analogous to the ones considered here [54], is not the optimal choice for the tensor current. We will still use θ = 0.5 in the non-perturbative computation, since our simulations were set up employing the optimal value for quark mass renormalization.
Finally, as already mentioned in the introduction, parallel to our work Dalla Brida, Sint and Vilaseca have performed a related, fully independent perturbative study as part of the setup of the chirally rotated Schrödinger Functional [33]. Their results for the one-loop matching factors [χ (1) O ] SF;lat are perfectly consistent with ours, providing a very strong crosscheck. θ α r α;θ 0;SF (c sw = 0) r α;θ 0;SF (c sw = 1) Table 1: Finite parts of one-loop renormalization constants in the scheme specified by the parameters θ and α for the unimproved and O(a)-improved fermion actions .

One-loop cutoff effects in the step scaling function
As mentioned above, the RG running is accessed via SSFs, defined in Eq. (2.19). It is thus both interesting and useful to study the scaling of Σ T within perturbation theory. Plugging the one-loop expansion of the renormalization constant in Eq. (2.19), we obtain an expression of the form where In order to extract the cutoff effect which quantifies how fast the continuum limit σ T is approached, we define 19) and the relative cutoff effect δ k The one-loop values of δ k for both the improved and unimproved renormalization conditions are listed in Table 3. The behaviour of δ k as a function of the lattice size is shown in Fig. 4. The figure shows that the bulk of the linear cutoff effect is removed by the improvement of the action, and that the improvement of the current has a comparatively small impact. Note also that θ = 0.5 leads to the smaller perturbative cutoff effects among the values explored, cf. Table 3.

Non-perturbative computations
We will now present non-perturbative results for both N f = 0 and N f = 2 QCD. The simulations underlying each of the two cases are those in [25] (which in turn reproduced and extended the simulations in [7]) and [8], respectively. For N f = 2 simulations are performed with non-perturbatively O(a) improved Wilson fermions, whereas in the quenched case the computation was performed both with and without O(a) improvement, which, along with the finer lattices used, allows for a better control of the continuum limit (cf. below). A gauge plaquette action is always used. In both cases, we rely on the computation of the SF coupling and its nonperturbative running, given in [7,62] for N f = 0 and [66] for N f = 2.

N f = 0
Simulation details for the quenched computation are given in [25]. Simulation parameters have been determined by tuning β such that the value of the renormalized SF coupling is kept constant with changing L/a, and fixing the bare quark mass to the corresponding non-perturbatively tuned value of κ c . As mentioned above, two separate computations have been performed, with and without an O(a) improved fermion action with a non-perturbatively determined c sw coefficient. 4 This allows to improve our control over the continuum limit extrapola- tion for σ T , by imposing a common result for both computations based on universality. It is important to note that the gauge ensembles for the improved and unim-proved computations are different, and therefore the corresponding results are fully uncorrelated. Another important observation is that the c T coefficient for the O(a) improvement counterterm of the tensor current is not known non-perturbatively, but only to leading order in perturbation theory. In our computation of Z T for N f = 0 we have thus never included the improvement counterterm in the renormalization condition, even when the action is improved, and profit only from the above universality constraint to control the continuum limit, as we will discuss in detail below. The resulting numerical values of the renormalization constants and SSFs are reported in Tables 4 and 5.

Continuum extrapolation of SSFs
As discussed above, the continuum limit for Σ T is controlled by studying the scaling of the results obtained with and without an O(a) improved actions. To that respect, we first check that universality holds within our precision, by performing independent continuum extrapolations of both datasets. Given the absence of the c T counterterm, we always assume that the continuum limit is approached linearly in a/L, and parametrize Σ csw=0 We observe that, in general, fits that drop the coarsest lattice, corresponding to the step L/a = 6 → 12, are of better quality; when the Σ T (L/a = 6) datum is dropped, σ csw=0 T (u) and σ csw=NP T (u) always agree within ∼ 1σ. The slopes ρ csw=NP T (u) are systematically smaller than ρ csw=0 T (u), showing that the bulk of the leading cutoff effects in the tensor current is subtracted by including the Sheikholeslami-Wohlert (SW) term in the action.
We thus proceed to obtain our best estimate for σ T (u) from a constrained extrapolation, where we set σ csw=0 T (u) = σ csw=NP T (u) = σ T (u) in Eq. (5.1), and drop the L/a = 6 → 12 step from the fit. The results for both schemes are provided in Table 6, and illustrated in Figs. 11, 12.

Fits to continuum step-scaling functions
In order to compute the RG running of the operator in the continuum limit, we fit the continuum-extrapolated SSFs to a functional form in u. The simplest choice, motivated by the perturbative expression for γ T and β, and assuming that σ T is a smooth function of the renormalized coupling within the covered range of values of the latter, is a polynomial of the form  Table 7. The one-and two-loop perturbative predictions are also shown for comparison.
The perturbative prediction for the first two coefficients of Eq. (5.3) reads Note, in particular, that perturbation theory predicts a dependence on N f only at O(u 2 ). We have considered various fit ansätze, exploring combinations of the order of the polynomial and possible perturbative constraints, imposed by fixing either p 1 or both p 1 and p 2 to the values in Eqs. (5.4,5.5). We always take as input the results from the joint c sw = 0 and c sw = NP extrapolation, discussed above. The results for the various fits are shown in Table 7. All the fits result in a good description of the non-perturbative data, with values of χ 2 /d.o.f. close to unity and little dependence on the ansatz. The coefficients of powers larger than u 3 are consistently compatible with zero within one standard deviation. We quote as our preferred fit the one that fixes p 1 to its perturbative value, and reaches O(u 3 ) (fit B in Table 7). This provides an adequate description of the non-perturbative data, without artificially decreasing the goodness-of-fit by including several coefficients with large relative errors (cf., e.g., fit E). The result for σ T from fit B in our two schemes is illustrated in Fig. 5. It is also worth pointing out that the value for p 2 obtained from fits A and B is compatible with the perturbative prediction within 1 and 1.5 standard deviations, respectively, for the two schemes; this reflects the small observed departure of σ T from its two-loop value until the region u 2 is reached, cf. Fig. 5.

Determination of the non-perturbative running factor
Once a given fit for σ T is chosen, it is possible to compute the running between two well-separated scales through a finite-size recursion. The latter is started from the smallest value of the energy scale µ had = L −1 had , given by the largest value of the coupling for which σ T has been computed, viz.
Using as input the coupling SSF σ(u) determined in [7], we construct recursively the series of coupling values This in turn allows to compute the product where U is the RG evolution operator in Eq. (2.18), here connecting the renormalised operators at scales µ had and 2 n+1 µ had . The number of iterations n is dictated by the smallest value of u at which σ T is computed non-perturbatively, i.e. u = 0.8873. We find u 7 = 0.950(11) and u 8 = 0.865(10), corresponding respectively to 8 and 9 steps of recursion. The latter involves a short extrapolation from the interval in u covered by data, in a region where the SSF is strongly constrained by its perturbative asymptotics. This point is used only to test the robustness of the recursion, but is not considered in the final analysis. The values of u k and the corresponding running factors are given in Tables 8 and 9. Once µ pt = 2 8 µ had has been reached, perturbation theory can be used to make contact with the RGI operator. We thus compute the total running factorĉ(µ) in Eq. (2.13) at µ = µ had asĉ whereĉ(µ pt ) is computed using the highest available orders for γ and β in our schemes (NLO and NNLO, respectively). In order to assess the systematic uncertainty arising from the use of perturbation theory, we have performed two crosschecks: (i) Perform the matching to perturbation theory at all the points in the recursion, and check that the result changes within a small fraction of the error.
(ii) Match to perturbation theory using different combinations of perturbative orders in γ and β: other than our NLO/NNLO preferred choice, labeled "2/3" -after the numbers of loops -in Tables 8 and 9, we have used matchings as a means to have a guesstimate of higher-order truncation uncertainties.
We thus quote as our final numberŝ In Fig. 6 we plot the non-perturbative running of the operator in our two schemes, obtained by running backwards from the perturbative matching point corresponding to the renormalized coupling u 7 = 0.950 (11). with our non-perturbative σ T , and compare it with perturbation theory. In order to set the physical scale corresponding to each value of the coupling, we have used Λ SF /µ had = 0.422 (32), from [7]. The latter work also provides the value of µ had in units of the Sommer scale r 0 [67], viz. (2r 0 µ had ) −1 = 0.718(16) -which, using r 0 = 0.5 fm, translates into µ had = 274(6) MeV. It is important to stress that the results in Eq. (5.10) are given in the continuum, and therefore do not contain any dependence on the regularization procedures employed to obtain them.

Hadronic matching
The final piece required for a full non-perturbative renormalization is to compute renormalization constants at the hadronic scale µ had within the interval of values of the bare gauge coupling covered by non-perturbative simulations in large, hadronic volumes. We have thus proceeded to obtain Z T at four values of the bare coupling, β = {6.0129,6.1628,6.2885,6.4956}, tuned to ensure that L -and hence the renormalized SF coupling -stays constant when L/a = {8, 10, 12, 16}, respectively. The results, both with and without O(a) improvement, are provided in Tables 10 and 11. These numbers can be multiplied by the corresponding value of the running factor in Eq. (5.10) to obtain the quantitŷ which relates bare and RGI operators for a given value of g 2 0 . They are quoted in Table 12; it is important to stress that the results are independent of the scheme within the ∼ 1% precision of our computation -as they should, since the scheme dependence is lost at the level of RGI operators, save for the residual cutoff effects which in this case are not visible within errors. A second-order polynomial fit to the dependence of the results in β for the numbers obtained from the scheme α = 1/2, which turns out to be slightly more precise, yields (5.14) These continuous form can be obtained to renormalize bare matrix elements, computed with the appropriate action, at any convenient value of β.

N f = 2
In this case all our simulations were performed using an O(a) improved Wilson action, with the SW coefficient c sw determined in [68]. , for the computation of the renormalization constant Z T (g 0 , a/L). All simulational details, including those referring to the tuning of β and κ, are provided in [8].
Concerning O(a) improvement, the configurations at the three weaker values of the coupling were produced using the one-loop perturbative estimate of c t [21], while for the three stronger couplings the two-loop value [69] was used. In addition, for L/a = 6, β = 7.5420 and L/a = 8, β = 7.7206 separate simulations were performed with the one-and two-loop value of c t , which results in two different, uncorrelated ensembles, with either value of c t , being available for u = 1.5078. Forc t the oneloop value is used throughout. Finally, since, contrary to the quenched case, we do not have two separate (improved and unimproved) sets of simulations to control the continuum limit, we have included in our analysis the improvement counterterm to the tensor current, with the one-loop value of c T [61].
The resulting values for the renormalization constants Z T and the SSF Σ T are listed in Table 13. The estimate of autocorrelation times has been computed using the "Gamma Method" of [70].

Continuum extrapolation of SSFs
In this case, our continuum limit extrapolations will assume an O(a 2 ) scaling of Σ T . This is based on the fact that we implement O(a) improvement of the action (up to small O(ag 4 0 ) effects inc t and O(ag 4 0 ) or O(ag 6 0 ) in c t , cf. above); and that the residual O(ag 4 0 ) effects associated to the use of the one-loop perturbative value for c T can be expected to be small, based on the findings discussed above for N f = 0. Our ansatz for a linear extrapolation in a 2 is thus of the form Furthermore, in order to ameliorate the scaling we subtract the leading perturbative cutoff effects that have been obtained in Sec. 4, by rescaling our data for Σ T as where the values of the relative cutoff effects δ k (a/L) are taken from Table 3. Continuum extrapolations are performed both taking Σ T and the one-loop improved Σ T as input; the two resulting continuum limits are provided in Tables 14 and 15, respectively. As showed in Fig. 6, the effect of including the perturbative improvement is in general non-negligible only for our coarsest L/a = 6 lattices. The slope of the continuum extrapolation is decreased by subtracting the perturbative cutoff effects at weak coupling, but for u 2 the quality of the extrapolation does not change significantly, and the slope actually flips sign. The u = 1.5078 case is treated separately, and a combined extrapolation to the continuum value is performed using the independent simulations carried out with the two different values of c t . We quote as our best results the extrapolations obtained from Σ T .

Fits to continuum step-scaling functions
Here we follow exactly the same strategy described above for N f = 0, again considering several fit ansätze by varying the combination of the order of the polynomial and the number of coefficients fixed to their perturbative values. The results are listed in Table 16. As in the quenched case, we quote as our preferred result the fit obtained by fixing the first coefficient to its perturbative value and fitting through O(u 3 ) (fit B). The resulting fit, as well as its comparison to perturbative predictions, is illustrated in Fig. 7.

Non-perturbative running
Using as input the continuum SSFs, we follow the same strategy as in the quenched case to recursively compute the running between low and high energy scales. In this case the lowest scale reached in the recursion, following [8], is given by g 2 SF (µ had ) = 4.61. Using the coupling SSF from [66], the smallest value of the coupling that can be reached via the recursion without leaving the interval covered by data is g 2 SF (µ pt ) = 1.017(10), corresponding to n = 7 (i.e. a total factor scale of 2 8 in energy, like in the N f = 0 case). The matching to the RGI at µ pt is again performed using the 2/3-loop values of the γ/β functions, and the same checks to assess the The running is illustrated, and compared with the perturbative prediction, in Fig. 8, where the value of log(Λ SF /µ had ) = −1.298(58) from [8] has been used. Using r 0 Λ SF = 0.30(3) from [66] and r 0 = 0.50 fm, this would correspond to a value of the hadronic matching energy scale µ had ≈ 432(50) MeV.

Hadronic Matching
The computation of the renormalization constants at µ had needed to match bare hadronic quantities proceeds in a somewhat different way to the quenched case.  Table 20. In this case the g 2 0 dependence is barely visible within the quoted errors, and the expected scheme independence holds only up to ∼ 3σ.

Conclusions
In this work we have set up the strategy for a non-perturbative determination of the renormalization constants and anomalous dimension of tensor currents in QCD using SF techniques, and obtained results for N f = 0 and N f = 2. In the former case we employed both O(a) improved and unimproved Wilson fermions, and simulations were performed at four values of the lattice spacing for each of the fourteen different values of the renormalization scale, resulting in an excellent control of the continuum limit. For N f = 2 our simulations were carried out with O(a) improved fermions, at only three values of the lattice for each of the six renormalization scales. The precision of the running factors up to the electroweak scale in the schemes that allow for higher precision is 0.9% and 1.1%, respectively. The somewhat limited quality of our N f = 2 dataset, however, could result in the quoted uncertainty for that case not being fully free of unquantified systematics. We have also provided values of renormalization constants at the lowest energy scales reached by the nonperturbative running, which allows to match bare matrix elements computed with non-perturbatively O(a) improved Wilson fermions and the Wilson plaquette gauge action.
As part of the ALPHA programme, we are currently completing a similar study in N f = 3 QCD [37], that builds upon a high-precision determination of the strong coupling [34][35][36] and mass anomalous dimension [9,10,22]. Preliminary results indicate that a precision ∼ 1% for the running to low-energy scales is possible even for values of the hadronic matching scale well below the one reached for N f = 2. This is an essential ingredient in order to obtain matrix elements of phenomenological interest with fully controlled uncertainties and target precisions in the few percent ballpark.

Appendix A Perturbative improvement
The improvement coefficient c T for the tensor current can, by definition, be determined by requiring an O(a) improved approach to the continuum of the renormalized correlation function at any given order in perturbation theory. As discussed in the main text, the computation of c T to one loop has been carried out in [61]; here we reproduce it, mainly as a crosscheck of our perturbative setup. We introduce the following notation for the renormalized tensor correlator k T;R in the chiral limit evaluated with SF boundary conditions at x 0 = T /2, where the θ as well as the a/L dependence have been made explicit. The one-loop expansion reads ;bi (T /2) + am where Z ξ is the renormalization constant of the boundary fermionic fields, and c T is the coefficient we are interested in, providing the O(a) improvement of the operator. In order to determine c wherek T is a shorthand notation for the correlator including the subtraction of the boundary and mass O(a) terms. The divergent part of Z T , as well as of Z ξ , cancel out in the ratio, since they are independent of θ at one loop. Following [71] T . As a second strategy to determine c T to one loop, one can exploit the tree-level identities obtained in [71], which relate k After some simple algebra we find   The comparison between our determination and the one in [61] is displayed in Fig. 9. In all cases, the continuum extrapolation has been performed using similar techniques to the one employed for the finite part of renormalization constants (see App. B).
T , compared with the result in [61].

Appendix B Continuum extrapolations in perturbation theory
In this appendix we summarize the techniques used to extrapolate our perturbative computations to a/L → 0, a necessary step in order to obtain scheme-matching and improvement coefficients. Our approach is essentially an application to the present context of the techniques discussed in Appendix D of [69], which have been applied in a number of cases, see e.g. [26]. The typical outcome of a perturbative computation is a linear combination of one-loop Feynman diagrams, e.g. the one yielding the one-loop coefficient Z (1) of a renormalization constant, for N values {l 1 , . . . , l N } of the variable l = L/a. We consider the quantity to be a function of l only. It is possible to identify all divergences appearing in the quantity of interest at one-loop, which in general means linear divergences related to the additive renormalization of the quark masses proportional to the one-loop critical mass m (1) cr , and the logarithmic divergences proportional to the (one-loop) anomalous dimension. The latter is particularly relevant for the present analysis, since it allows to check the consistency of the fitting procedure and provides a natural criterion for the choice of the best fitting ansatz. In the following we consider finite quantities, since the leading divergence is subtracted, and the critical mass is appropriately tuned. Considering F (l) as a generic one-loop interesting quantity, following [69] we conservatively assign the error since in this case the computation has been carried out in double precision. As expected, the asymptotic behaviour is (cf. Eq. (4.15)) with a residue R n (l) that decreases faster than any of the terms in the sum as l → ∞.