Non-perturbative renormalization of quark bilinear operators with Nf=2 (tmQCD) Wilson fermions and the tree-level improved gauge action

We present results for the renormalization constants of bilinear quark operators obtained by using the tree-level Symanzik improved gauge action and the Nf=2 twisted mass fermion action at maximal twist, which guarantees automatic O(a)-improvement. Our results are also relevant for the corresponding standard (un-twisted) Wilson fermionic action since the two actions only differ, in the massless limit, by a chiral rotation of the quark fields. The scale-independent renormalization constants ZV, ZA and the ratio ZP/ZS have been computed using the RI-MOM approach, as well as other alternative methods. For ZA and ZP/ZS, the latter are based on both standard twisted mass and Osterwalder-Seiler fermions, while for ZV a Ward Identity has been used. The quark field renormalization constant Zq and the scale dependent renormalization constants ZS, ZP and ZT are determined in the RI-MOM scheme. Leading discretization effects of O(g^2 a^2), evaluated in one-loop perturbation theory, are explicitly subtracted from the RI-MOM estimates.


Introduction
Non-perturbative renormalization is an essential ingredient of lattice QCD calculations which aim at a percent level accuracy. In this paper, we present calculation methods and results for the renormalization constants (RCs) of bilinear quark operators for the lattice action used by the ETM Collaboration. The simulations with N f = 2 dynamical flavours employ the tree-level Symanzik improved gauge action and the twisted mass fermionic action at maximal twist. In the chiral limit, the latter is related to the standard Wilson fermionic action by a chiral transformation, under which quark composite operators behave in a definite and simple way. Therefore the RCs we compute here can be also employed to renormalize bilinear quark operator matrix elements computed with the same glue action but standard (un-twisted) Wilson quarks.
For the computation of the bilinear operator RCs we have used the RI-MOM method [1]. For the calculation of the scale independent RCs, Z V , Z A and the ratio Z P /Z S , a set of alternative methods have also been used. In particular, we have determined Z V by imposing the validity of the axial Ward identity whereas for Z A and Z P /Z S we have employed a new method which makes a combined use of standard twisted mass (tm) and Osterwalder-Seiler (OS) formalism. A comparison of the results obtained with different approaches provides an estimate of the size of the systematic errors affecting the calculation of the RCs at fixed lattice bare coupling. Needless to say, in renormalized quantities the discretization errors coming from both RCs and bare quantities will be removed by continuum extrapolation, as illustrated with an example in sect. 3.3.
Throughout the paper we follow the standard practice of labelling the RCs according to the notation adopted in the un-twisted Wilson case. Thus, for instance, Z V and Z A denote the RCs of the local vector and axial-vector currents of the standard Wilson action even though, in the maximal twist case, they renormalize in the physical basis the local (charged) axial-vector and polar-vector currents respectively. We refer to Table 1 for the explicit presentation of the renormalization pattern.
We have computed the RCs at three values of the gauge coupling, namely β = 3.80, 3.90 and 4.05, which correspond to inverse lattice spacing a −1 ≃ 2.0, 2.3 and 2.9 GeV. Preliminary results of this work have been presented in ref. [2]. Further details concerning the ETMC simulations can be found in refs. [3]- [6].
The paper's contents are the following. In Section 2 we present the novel approach of computing the scale independent RCs Z A and Z P /Z S and the Ward identity determination of Z V . In Section 3 we discuss the RI-MOM approach applied to tm quarks and provide the details of the numerical analysis. A collection of our final results for the RCs can be found in Tables 4, 5 and 6. In Tables 4 and 5 the predictions of one loop boosted perturbation theory [7,8] for the various RCs are also given. Finally in the Appendix we present the proof of the automatic O(a)-improvement of the RI-MOM RCs calculated with maximally twisted quarks.
An RI-MOM calculation of the RCs of bilinear quark operators similar to the one described in the present paper but for the standard Wilson's plaquette gauge action and nonperturbatively improved Wilson (clover) fermions has been recently presented in ref. [9].

A novel approach to the calculation of the scale independent RCs
In this section we present a calculation of the scale independent RCs, namely Z V , Z A and Z P /Z S . The evaluation of Z V is based on the PCAC Ward identity, which in the quenched approximation has led to very precise results (see e.g. refs. [10,11]). On the other hand, Z A and Z P /Z S are computed from two-point correlators only with a new method, which is based on the simultaneous use of two regularizations of the valence quark action. One is the standard tm action at maximal twist, while the other is the OS variant [12]. In the so called physical basis these actions can be written in the form: x q=u,dq with W cr = − a 2 µ ∇ * µ ∇ µ +M cr . The Wilson parameters are r u = −r d = 1 for the standard tm case and r u = r d = 1 for OS quarks. In the sea sector we have two degenerate quarks, regularized in the standard tm framework.
Let us illustrate the general idea of the calculation of the scale independent RCs, Z A and Z P /Z S . It is based on the observation that the axial transformations of the quark fields from the physical basis to the so called twisted basis (primed fields), q = exp (iγ 5 r q π/4) q ′ andq =q ′ exp (iγ 5 r q π/4) , transform the tm and OS actions of Eq. (1) into an action with the Wilson term in the standard form (i.e. having no γ 5 and r u = r d = 1) and a mass term which takes the form iµ uū ′ γ 5 u ′ − iµ dd ′ γ 5 d ′ in the tm case and iµ uū ′ γ 5 u ′ + iµ dd ′ γ 5 d ′ in the OS case. This also implies, in turn, that the RCs for the operators in the twisted basis, being defined in the chiral limit, are the same for the Wilson, tm and OS cases. Consider, now, a non-singlet quark bilinear operator, O Γ =ūΓd, defined in terms of the fields of the physical basis. Under the axial transformations of Eq. (2) this operator transforms into an operator in the twisted basis, denoted as OΓ and OΓ for the tm and OS case respectively. In general, the two operators are not of the same form. However their renormalized matrix elements between given physical states, say Z OΓ α|OΓ|β tm and Z OΓ α|OΓ|β OS are estimates of the same physical matrix element α|(O Γ ) R |β up to O(a 2 ) errors: RCs are named, as anticipated, after the twisted basis, in which the Wilson term has its standard form. The operator renormalization pattern in the physical and twisted bases is given in Table 1 for both OS and tm formulations at maximal twist. The primed operators refer to the twisted basis while the unprimed ones to the physical basis. Notice that the condition of maximal twist ensures that cut-off effects are of order O(a 2 ) [13]. The main point of Eq. (3) is that the tm and OS determinations of the same continuum matrix element are equal up to discretization effects. Our proposal for the determination of Z P /Z S and Z A , described in the following subsections, is based on this observation. Table 1: Renormalization pattern of the bilinear quark operators for the OS and tm case at maximal twist. The unprimed operators refer to the physical basis while the primed ones to the twisted basis. The symbols P ud andT µν,ud indicate the operatorsūγ 5 d andūσ µν γ 5 d.

Calculation of Z P /Z S
The method is based on comparing the amplitude g π = 0|P |π , computed both in tm and OS regularizations. We start by considering, in the physical basis, the correlation function which at large time separations behaves like From Table 1 we see that in the twisted basis this corresponds to C S ′ S ′ (t) (for OS regularization) and C P ′ P ′ (t) (for tm regularization): which implies The ratio Z P /Z S is extracted from the above equation in the chiral limit. As our simulations are performed at finite quark masses, an extrapolation of our Z P /Z S -estimators to zero quark mass is necessary.

Calculation of Z A
In order to compute Z A , we consider the calculation of the charged pseudoscalar meson decay constant f π , in both OS and tm regularizations. We start with the tm regularization. From the axial Ward identity in the physical basis, we have the well-known result for the pseudoscalar meson decay constant [14]: Note that in this case no RC is needed [14]. Thus the pion decay constant can be extracted from the large time asymptotic behaviour of [C P ′ P ′ (t)] tm .
In the OS formulation we consider, besides the correlator C P P of Eq. (4), also C AP , which is defined by Its large time asymptotic behaviour is: with The last equation may be solved for f π , which is thus obtained from the correlation functions of Eqs. (5) and (10). With the aid of Table 1, this solution is expressed in terms of OS-regularized quantities in the twisted basis (note ξ AP = iξ A ′ S ′ ) as follows: Combining Eqs. (8) and (12), we obtain from which an estimate of Z A can be extracted. Again, results obtained at finite quark masses need to be extrapolated to the chiral limit. In ref. [15] the discretization effects affecting the quenched pseudoscalar decay constant in the tm and OS regularizations have been studied in detail. The tm and OS results were compatible within one standard deviation for three values of the lattice spacing, indicating the smallness of O(a 2 ) cut-off effects.

Calculation of Z V
The determination of the renormalization constant Z V can be done using solely tm quarks. It is based on the comparison of the point-like axial current A µ,ud , which renormalizes with Z V (see Table 1), with the exactly conserved one-point split current ud . The four-divergence of the latter is exactly equal to the sum of the valence quark mass values, (µ 1 + µ 2 ), times the corresponding pseudoscalar density. Therefore, we extract Z V by solving the equation: where∂ 0 is the symmetric lattice time derivative, and extrapolating results to the chiral limit.

Results
In   highest sea quark mass is around half the strange quark mass. For the inversions in the valence sector we have made use of the stochastic method (one-end trick of ref. [16]) in order to increase the signal to noise ratio. Propagator sources have been placed at randomly located timeslices. This turned out to be an optimal way to further reduce the autocorrelation time. Our correlation functions are computed for all combinations of valence quark masses appearing in Table 2. In order to deal with effectively independent measurements of RC-estimators, we have selected gauge field configurations separated by 20 HMC (length 1/2) trajectories. Within each ensemble of such gauge configurations statistical errors are evaluated using a jackknife procedure. For results involving an extrapolation in the sea quark mass and/or combined fits at several lattice couplings, statistical errors have been estimated with a bootstrap procedure for a 1000 bootstrap events. Typical plots on the data quality of Z V , Z A and Z P /Z S are shown in Fig. 1; in this example we show the case aµ min sea = aµ 1 = aµ 2 for the three values of the gauge coupling. Various methods have been implemented in taking the chiral limit (in the sea and valence quark sector). A first method consists in calculating the RCs at fixed sea quark mass for three choices of valence quark pairs, namely with aµ 1 = aµ 2 , or for all pairs of (aµ 1 , aµ 2 ) or for (aµ 1 , aµ 2 ) ≥ aµ sea and taking the "valence chiral limit" using a linear fit in terms of the sum of the valence quark masses. Subsequently, the RCs were quadratically extrapolated to the sea quark chiral limit. Note that a quadratic dependence on aµ sea is expected from the form of the sea quark determinant, assuming that lattice artefacts on the RCs are not sensitive to spontaneous chiral symmetry breaking. However we have verified that a linear fit in µ sea leads to compatible results, albeit with larger final errors. A second method consists in inverting the order of the two chiral limits. A third method is simply (c) Figure 1: The plateaux quality for Z V (a), Z A (b) and Z P /Z S (c) for three gauge couplings and masses aµ min sea = aµ 1 = aµ 2 . From Table 2 we have aµ min sea (β = 3.80) = 0.0080, aµ min sea (β = 3.90) = 0.0040 and aµ min sea (β = 4.05) = 0.0030. Note that in Fig.(b) the data have been shifted for clarity.
the extraction of the RCs from the subset of data satisfying µ valence = µ sea . This allows reaching the chiral limit with a single, linear extrapolation in the quark mass. In most cases the quality of the fits is very good and the final results, obtained from these different extrapolations, are compatible within one standard deviation. We have opted to quote as final results those produced by the first method, for all pairs of (aµ 1 , aµ 2 ), followed by a linear fit in terms of a 2 µ 2 sea . In Fig. 2 (left panel) we show the "valence chiral limit" extrapolation of the three scale independent RCs, computed at the lightest sea quark mass for each gauge coupling. In Fig. 2 (right panel) we present the chiral limit in the sea sector for all the RCs for each gauge coupling.
Our final results are collected in Table 4 (see rows labelled as "Alt. methods"). For comparison with the RI-MOM results and conclusions on the precision of the methods used, see below. Here we point out that the Z V value at β = 3.90 compares nicely with the result of Ref. [17], produced from 3-and 2-point correlation functions, in the context of a calculation of the pion form factor.

Renormalization constants in the RI-MOM scheme
The non-perturbative determination of the RCs in the RI-MOM scheme [1] is based on the numerical evaluation, in momentum space, of correlation functions of the relevant operators between external quark and/or gluon states. In this work we are interested in the case of the bilinear quark operators O Γ =ūΓd, where Γ = S, P, V, A, T stands for I, γ 5 , γ µ , γ µ γ 5 , σ µν respectively. The fields u and d of the operator O Γ belong to a twisted doublet of quarks, regularized by the tm action of Eq. (1), with r u = −r d = 1.
The relevant Green functions are the quark propagator the forward vertex function and the amputated Green function 1 The RI-MOM renormalization scheme consists in imposing that the forward amputated Green function, computed in the chiral limit in the Landau gauge and at a given (large Euclidean) scale p 2 = µ 2 , is equal to its tree-level value. In practice, the renormalization condition is implemented by requiring in the chiral limit that where PΓ is a Dirac projector satisfying Tr [Γ PΓ] = 1. 2 The quark field RC Z q , which enters Eq. (18), is obtained by imposing, in the chiral limit, the condition 3 In order to minimize discretization effects, we select momenta with components p ν = (2π/L ν ) n ν in the intervals [2,3], [2,3], [4,7]) , for L = 24 (β = 3.8 and 3.9) 1 The calculation of G Γ (p) in Eq. (16) involves the propagator S d (0, y) which, in the tm approach, equals The appearance ofΓ in Eq. (18) is a trivial consequence of the fact that RCs, as we said, are named after the (here unphysical) twisted quark basis while the operator O Γ =ūΓd is expressed in the physical quark basis. 3 Strictly speaking, the renormalization condition of Eq. (19) defines the so called RI ′ scheme. In the original RI-MOM scheme the quark field renormalization condition reads The two schemes differ in the Landau gauge at the N 2 LO. In this work, when perturbative predictions are used, the difference between the two schemes has been properly taken into account.
and [2,5], [2,5], [4,9]) , for L = 32 (β = 4.05) (21) (L ν is the lattice size in the direction ν). The time component of the four-momentum is shifted by the constant ∆p 4 = π/L 4 , in order to account for the use of anti-periodic boundary conditions on the quark fields in the time direction. Furthermore, we have only considered for the final RI-MOM analysis the momenta satisfying which helps in minimizing the contribution of Lorentz non-invariant discretization effects (cfr Eq. (31)).
In order to improve the statistical accuracy of the RCs, we have averaged the results obtained from the correlation functions of the operators O Γ =ūΓd and O ′ Γ =dΓu. Similarly, in the case of the quark field RC, we have computed Z q by averaging the results obtained for the up and down quark propagators.
The RCs, calculated in the chiral limit in the way described above, are automatically improved at O(a). Actually in the maximal twist situation O(a) improvement holds for all external momenta p and non-zero (valence and sea) quark mass at the level of the basic quantities Tr[Λ Γ (p) P Γ ] and Tr [(/ p S q (p) −1 )/p 2 ] entering in Eqs. (18)- (19). A proof of this statement, which follows from an analysis based on the symmetries of the tm action at maximal twist and the O(4) symmetry of the underling continuum theory, is given in the appendix. 4 Details of the RI-MOM simulation are collected in Table 2. The lattice parameters are those used also in the determination of the scale independent constants with the methods discussed in the previous section. However, a different ensemble of independent gauge configurations has been analyzed in each case. Moreover, the inversions in the valence sector for the RI-MOM study have been performed without using the stochastic approach, but with point-like sources randomly located on the lattice for each gauge configuration.

Analysis of the twisted mass quark propagator
The necessary Green functions for the RI-MOM determination of RCs of bilinear quark operators are the quark propagator S q (p) and the amputated vertex Λ Γ (p), evaluated in momentum space in the Landau gauge. In this section, we first illustrate the results on the tm quark propagator.
At O(a), the spin-flavour structure of the lattice artefacts of the quark propagator differs from that of the standard Wilson case, due to the twisted Wilson term in the action. With tm fermions, the explicit breaking of parity at finite lattice spacing allows for the presence of a parity violating contribution in the quark propagator. By neglecting O(4) symmetry violating effects, which only appear at O(a 2 ), the inverse quark propagator can be expressed in terms of three scalar form factors, as follows: At large p 2 , Σ 1 and Σ 2 are related to the quark field RC Z q (see Eq. (19)), and to the renormalized quark massμ q respectively [19]. The parity violating term proportional to Σ 3 represents an O(a) discretization effect, induced by the twisted Wilson term in the action. 5 At maximal twist, the form factors Σ 1,2,3 are given at tree-level by In Fig. 3, the non-perturbatively determined form factors Σ 1,2,3 are plotted against a 2p2 , for β = 3.90 and aµ sea = 0.0040. The following observations are in order: • The form factor Σ 1 exhibits a nice plateau at large p 2 , with a tiny slope which is partly due to a NLO renormalization scale dependence (the quark field anomalous dimension vanishes at LO in the Landau gauge) and partly to O(a 2 ) discretization effects. The quark field RC Z q is obtained from the Σ 1 results at large p 2 , extrapolated to the chiral limit (see Eq. (19)).
• The form factor Σ 2 also exhibits a plateau at large p 2 , which is expected to be proportional, up to the quark field RC, to the renormalized valence quark massμ q in the RI-MOM scheme. Indeed, at variance with the cases of Σ 1 and Σ 3 , a clear separation among the results obtained for Σ 2 at different bare valence quark masses is visible in the plot.
• As expected, the form factor Σ 3 , which represents a pure O(a) discretization effect, has opposite sign for the up and down quark propagators (r u = −r d = 1). Its behaviour is quite close to the tree level estimate of Eq. (25). We find Σ 3 ≃ ± c q ap 2 /2, with c q ≃ 0.8 ÷ 0.9.

Renormalization constants
The RCs in the RI-MOM scheme are determined by closely following the procedure summarised in section 3 above and described in detail in ref. [18]. This procedure involves two main steps: the chiral limit extrapolation of the RCs, at fixed coupling and renormalization scale, and the study of the renormalization scale dependence.

Chiral extrapolations
RI-MOM is a mass independent renormalization scheme. Since in practice the RCs are obtained at non vanishing values of the quark masses, an extrapolation of the results to the chiral limit, both in the valence and the sea quark masses, must be performed. The validity of the RI-MOM approach relies on the fact that, at large p 2 , Green functions are expected to depend smoothly on the quark masses, since their non-perturbative contributions vanish asymptotically in that limit. However, specific care must be taken in the study of the pseudoscalar Green function V P , since the leading 1/p 2 -suppressed contribution is divergent in the chiral limit, due to the coupling with the Goldstone boson [1,20,21]. Moreover, for tm fermions, the explicit breaking of parity at finite lattice spacing induces a coupling of the Goldstone boson also with the scalar operator. Though the latter is suppressed at O(a 2 ), its effect may be not negligible at the couplings considered in the present simulations.
In order to subtract the pseudoscalar mass dependence of the amputated vertex functions V P and V S , at each given p 2 we fit their values at different quark masses to the Ansatz where M P S (µ 1 , µ 2 ) is the mass of the pseudoscalar meson composed by valence quarks of mass µ 1 and µ 2 (with r 1 = −r 2 ) 6 . The first term in the fit, i.e. the function A P (S) (p 2 ), represents the Green function V P (S) in the (valence) chiral limit, from which the RC Z P (S) is eventually extracted. The last term in Eq. (26) accounts for the presence of the Goldstone pole. Using the results of the fit, we can define the subtracted correlators The effect of the subtractions is illustrated in Fig. 4. We find that, while in the case of V P the contribution of the Goldstone pole is clearly visible, the subtraction is practically irrelevant in the case of V S , indicating a strong O(a 2 ) suppression of the parity violating coupling. We have also verified that, the alternative procedure [22] for the Goldstone pole subtraction, based on the definition .
(28) leads to consistent results in the chiral limit, within our statistical errors.
The dependence of Z V , Z A and Z T , which are not affected by the Goldstone pole contamination, on the valence quark masses, is shown in Fig. 5, for all three couplings. For illustration purposes, a specific value of the sea quark mass has been chosen in each case. We see that the valence quark mass dependence of the RCs is rather weak and well consistent with a linear behaviour.
In the tm formulation of lattice QCD with N f = 2 flavours of degenerate sea quarks, the fermionic determinant is a function of the sea quark mass squared. Therefore, we expect that the RCs, which by definition are insensitive to the effects of spontaneous chiral symmetry breaking, will also exhibit the same dependence. This is shown in Fig. 6. We find that all RCs depend very weakly on the sea quark mass squared. We thus obtain our mass independent RCs by performing the sea quark mass extrapolation to the chiral limit linearly in a 2 µ 2 sea . We also checked that chiral extrapolations based on a first or second order polynomial fit in aµ sea lead within errors to practically equivalent results.

Renormalization scale dependence and subtraction of the O(g 2 a 2 ) effects
Once the RCs have been extrapolated to the chiral limit, we investigate their dependence on the renormalization scale by evolving, at fixed coupling, the RCs to a reference scale µ 0 = 1/a. This is done using where the evolution function C Γ is expressed in terms of the beta function β(α) and of the anomalous dimension of the relevant operator γ Γ (α) by This function is known, in the RI-MOM scheme, at the N 2 LO for Z T [23] and at the N 3 LO for Z S and Z P [24]. Since Z V , Z A and the ratio Z P /Z S are scale independent, they have vanishing anomalous dimensions; thus for these quantities we have C Γ (µ 0 , µ) = 1.
In the non-perturbative calculation, the RCs evolved to the reference scale µ 0 = 1/a maintain a dependence on the renormalization scale p 2 at which they have been initially computed. We will keep track of this dependence, which (at large enough p 2 ) is mostly due to discretization effects, by denoting these RCs as Z Γ (1/a; a 2 p 2 ), where the first variable indicates the renormalization scale, µ 0 = 1/a, and the second the dependence on the initial momentum.
In Fig. 7 we show the results for all Z Γ (1/a; a 2 p 2 ) at β = 3.9 as a function of a 2p2 (empty symbols). The residual a 2 p 2 -dependence which is observed in these results in the large momentum region is practically linear, suggesting that leading discretization effects are O(a 2 p 2 ). As illustrated in the plots, this dependence is particularly pronounced in the case of the pseudoscalar RC Z P .
In order to reduce the size of discretization errors, we analytically subtract from the quark propagator and the amputated vertex functions the O(g 2 a 2 ) contributions, recently computed in lattice perturbation theory [8]. Up to 1-loop and up to O(a 2 ), the amputated projected Green functions V Γ (p), defined in Eq. (18), have the simple and general expression The values of the coefficients b  Table 3 for the lattice action used in the present study, namely the tree-level Symanzik improved gluon action in the Landau gauge and the tm fermionic action at maximal twist 7 . A result similar to Eq. (31) also holds for the form factor Σ 1 (p) of the inverse quark propagator, from which the quark field RC Z q is evaluated: Our definition of this form factor on the lattice, which is equivalent to Eq. (19) up to terms of O(a 2 ), is where the sum ′ ρ only runs over the Lorentz indices for which p ρ is different from zero and N(p) = ′ ρ 1. With this definition, the coefficients of Eq. (32) take the values Using Eqs. (31) and (32), we then define the corrected amputated Green functions and the corrected form factor Σ 1 as which are free of O(g 2 a 2 ) effects. Note that the O(a 2 ) terms depend not only on the magnitude, ρ p 2 ρ , but also on the direction of the momentum, p ρ , as manifested by the presence of ρ p 4 ρ . As a consequence, different renormalization prescriptions, involving different directions of the renormalization scale p ρ , treat lattice artifacts differently. In the numerical evaluation of the perturbative correction in Eq. (35), we replaced g 2 with a  simple boosted coupling [25], defined asg 2 = g 2 0 / P , where the average plaquette P is computed non-perturbatively.
The effect of the subtraction (35) is also illustrated in Fig. 7, which shows that the a 2 p 2 -dependence of the RCs is significantly reduced by the perturbative correction. By fitting the RCs as in the large momentum region a 2p2 > ∼ 1 (with λ Γ just constrained to depend smoothly on g 2 -see Eq. (38)), we find that the slopes λ Γ are reduced, after the perturbative subtraction, at the level of 10 −2 or smaller for Z V , Z A , Z T and Z q . For Z S we find that the slope, which is also about 10 −2 before implementing the correction, increases slightly after the subtraction. For Z P , on the other hand, the slope is rather large before the subtraction; i.e. λ P ≃ 5 · 10 −2 . As shown in Fig. 7, the effect of the subtraction is beneficial also in this case, but inadequate for correcting the bulk of the observed a 2 p 2 -dependence. This is not completely unexpected: when discretization effects are large, the subtraction of only the leading O(g 2 a 2 ) terms may be not sufficient to reduce them to a negligible level. Similar results, with approximately equal values of the slopes λ Γ , are obtained at all three values of the lattice spacing.
A useful way to investigate the size of lattice artifacts, which are left after the perturbative subtraction, consists in comparing the renormalization scale dependence of the RCs as observed at different lattice spacings. Specifically, at each value of the lattice spacing we fix two common values of the renormalization scale, µ and µ ′ , and compute the step scaling functions If not for discretization effects, the step scaling functions should be independent of the cut-off and equal to the evolution functions C Γ (µ, µ ′ ) of Eq. (29).
In fig. 8 we show the results obtained for the step scaling functions of the quark bilinear operators and of the quark field RCs, as a function of (a/r 0 ) 2 , for two common pairs renormalization scales, (µ 2 , µ ′ 2 ) = (8.5, 6.0) GeV 2 and (µ 2 , µ ′ 2 ) = (8.5, 7.25) GeV 2 , which both lie in the interval of momenta accessible at all values of the lattice spacing. We see from these plots that the dependence of the step scaling functions on the lattice cut-off is tiny in most of the cases, being clearly visible only in the case of Σ P . Moreover, for the scale independent RCs Z V and Z A , a linear extrapolation of the results to the continuum limit leads to values which are perfectly consistent with unity, within the statistical errors. For the other RCs, the continuum limit of the step scaling functions leads to non-perturbative results which are in good agreement with the perturbative determinations of the evolution functions C Γ (µ, µ ′ ), computed at the N 2 LO for Z T and at the N 3 LO for Z S , Z P and Z q . These results provide evidence that higher order perturbative contributions to the anomalous dimensions of the various operator are not relevant for describing the renormalization scale dependence of the RCs in the region of momenta explored in the present calculation.
In order to account for the residual discretization effects in the calculation of the RCs, we follow two different approaches: -Extrapolation method (M1): after subtraction of the O(g 2 a 2 ) terms according to Eq. (35), we extrapolate the RCs linearly to a 2 p 2 → 0. Specifically, we fit Z Γ (1/a; a 2 p 2 ) with Eq. (36) in the large momentum region, 1.0 < ∼ a 2p2 < ∼ 2.2. The slopes λ Γ exhibit  only a very mild dependence on the coupling. We parameterize this dependence by performing a simultaneous extrapolation at the 3 values of the coupling and writing the slopes as where g 0 is the coupling corresponding to the intermediate value β = 3.9. The intercept of the extrapolation as determined from the fit, and indicated with Z Γ (1/a) in Eq. (36), provides our final estimate of the RC with the extrapolation method. The fit is illustrated, for all RCs, in fig. 9.
-p 2 -window method (M2): in this approach we do not attempt any additional subtraction of discretization effects, besides the perturbative one. The final estimates of RCs are obtained by taking directly the results of Z Γ (1/a; a 2 p 2 ) at a fixed value of p 2 (in physical units). In practice, this is done by fitting the RCs to a constant in the small momentum intervalp 2 ∈ [8.0, 9.5] GeV 2 . The idea behind this approach is that, once RCs are combined with bare quantities, so as to construct the physical observables of interest, the residual O(a 2 ) effects, which are present in both RCs and bare quantities, will be extrapolated away in the continuum limit.
The two methods are compared below, in a specific example. This consists in the scaling of the pseudoscalar meson mass, composed by two degenerate quarks of fixed (renormalized) mass. As expected, the two methods give different results at fixed lattice cut-off, but the difference disappears in the continuum limit.

Analysis of systematic errors
Before presenting our final results for the RCs, we list and discuss in some detail the main systematic uncertainties.
-Finite size effects: the RCs, determined in the RI-MOM scheme at large momenta, are short distance quantities and, as such, should not be affected by significant finite size effects. In order to verify this expectation, an additional RI-MOM calculation at β = 3.90 has been performed, on a 32 3 × 64 lattice. A comparison of these results with those obtained on the smaller 24 3 × 48 lattice is illustrated in Fig. 10. We find that, for the momentum range in which results on both volumes are available, the differences are smaller than statistical errors.
-Subtraction of the Goldstone pole and chiral extrapolations: the chiral extrapolation of the RCs is rather delicate for Z P (and to a lesser extent for Z S ), due to the presence of the Goldstone pole which has to be subtracted. The uncertainty introduced by the subtraction of the Goldstone pole contribution to Z S and Z P can be estimated by comparing the results obtained following two different subtraction approaches: the fit of Eq. (26) or, alternatively, the procedure based on Eq. (28). We find that differences are always negligible for Z S , whereas for Z P they are at the level of our statistical errors.
As illustrated in Figs. 5 and 6, the chiral extrapolations in the valence quark mass for the other RCs and in the sea quark mass for all RCs are quite smooth. They are also well consistent with the expected leading linear and purely quadratic dependence on the valence and sea quark masses respectively. Therefore, we estimate that the uncertainties associated with these chiral extrapolations are negligible with respect to the statistical errors.
-Discretization effects: the analysis of these effects has been presented in section 3.2.2. The leading O(g 2 a 2 ) discretization errors have been accounted for, using the analytical results of ref. [8]. With the method denoted as M1 in section 3.2.2, after the perturbative subtraction, residual O(a 2 p 2 ) effects are extrapolated away, with the Ansatz (36). These contributions are not subtracted, instead, in the approach denoted as M2. Discretization effects affecting both RCs and the bare matrix elements will be corrected for, by extrapolating the physical (renormalized) quantity to the continuum limit.
As an example of the results obtained by implementing the two approaches, we show in Fig. 11 the continuum extrapolation of the pseudoscalar mass squared, computed at a fixed value of the renormalized quark mass µ R . The renormalization of the latter has been achieved using the estimates M1 and M2 for Z P given in Tables 4 and 5 8 . Although the two Z P determinations differ by visible discretization effects at finite lattice spacing, particularly at the two coarsest lattice, the continuum limit results of the pseudoscalar mass squared are consistent within the two determinations. This is evidence that the discretization errors as well as other possible systematic effects in the analysis of the RI-MOM RCs are well under control, at least within our current statistical errors (say at the level of one standard deviation and possibly even less).

Final results and comparison with perturbation theory
The final results for the bilinear quark operators and the quark field RCs, obtained with the RI-MOM method, are collected in Tables 4 and 5 (labelled as RI-MOM (M1) and RI-MOM (M2)). These results are compared with those obtained in Sec. 2 for the scale independent Z V , Z A and Z P /Z S (Alt. methods). In addition, we give in Tables 4 and 5 the predictions of one loop boosted perturbation theory (BPT) [7,8], obtained with two definitions of the boosted coupling: the first isg 2 = g 2 0 / P , also used in the previous section, and the second is based on the one-loop tadpole-improved expansion of log P [25].
From Tables 4 and 5, we see that the largest deviations between the central values of the RI-MOM (M1) determinations and the alternative determinations of Section 2 are at the level of 3-4% or smaller for Z A and Z V . These differences are not always accounted for by the errors quoted for each result. When statistically significant, these differences provide an estimate of the systematic uncertainties affecting the calculations. As discussed in the previous section, we expect these uncertainties to be dominated by O(a 2 ) discretization effects. Nevertheless, we do not include these differences in the final errors. The reason is that such estimates of discretization effects, though meaningful at finite lattice spacing, are of no practical significance when the continuum limit of a physical quantity is eventually computed.
The comparison of the Z P /Z S results, obtained with the RI-MOM and the alternative method, deserves a further comment. The difference between the two results is about 15% for the coarsest lattice and 4% for the finest one. The rather big difference noticed in the former case may be attributed not only to discretization effects, but also to an imprecise tuning of the twist angle to maximal twist at this coupling. This is expected to influence substantially the OS pseudoscalar density.
As shown in Tables 4 and 5, the non-perturbative results also lie in the range of oneloop BPT estimates, with the exception of Z S and, more notably, of Z P . We emphasise that the accuracy of Z P is crucial in the determination of quark masses, since the bare quark mass, computed with maximally twisted fermions, renormalizes with of Z −1 P ; see for example ref. [26].   In most phenomenological applications, the renormalization scheme of choice is MS and the preferred reference scale is 2 GeV. For this reason, we provide in Table 6 the values of the scheme dependent RCs Z S , Z P , Z T and Z q in the MS scheme at the scale of 2 GeV. The conversion from the RI-MOM scheme at µ = 1/a to MS at µ =2 GeV has been performed using renormalization group-improved perturbation theory at the N 2 LO for Z T and at the N 3 LO for Z S , Z P and Z q . For completeness, we also report in the table the results for the scale independent Z V and Z A . Note, however, that the Ward identity determination of Z V , given in Tables 4 and 5, has a much better statistical accuracy.
It should be remarked that in Table 6 we only quote the results obtained from RI-MOM renormalization conditions with the method M1, since they are in general affected by smaller discretization effects and should hence be preferably used in computations at fixed gauge coupling without continuum extrapolation 9 . The case of Fig. 11 is an exception in this sense, because there an accidental cancellation (enhancement) of the cutoff effects coming from the bare pseudoscalar meson mass M ps and from Z P computed with the method M2 (M1) leads to a similar scaling behaviour of the two sets of data corresponding to the M2 and M1 determinations of Z P . Moreover, with respect to the method M1, the method M2 tends to give results that are affected by smaller statistical uncertainties, since the latter procedure does not involve the extrapolation to a 2 p 2 = 0 of Eq. (36).

Conclusions
In this paper we have presented a non-perturbative determination of the RCs of bilinear quark operators for the tree-level Symanzik improved gauge action and the N f = 2 twisted mass fermion action at maximal twist, which guarantees automatic O(a)-improvement. The values of these constants do not depend on whether tm, OS or Wilson fermions are used. The RCs of the five bilinear quark operators and the quark field have been determined with the RI-MOM method. We have also obtained an independent estimate of Z V , based on a Ward identity, and have introduced a new method for Z A and Z S /Z P , based on a combined use of standard tm and OS fermions. The main results for the RCs are collected in Tables 4, 5 and 6, where they are also compared with the predictions of 9 Anyway, the values of the RCs in the MS scheme at the scale of 2 GeV obtained with the method M2 can easily be obtained for each a = a(β) from their M1-method counterparts and the RI-MOM results in Tables 4 and 5 with the help of the relation Z MS one-loop perturbation theory. We find that the differences between perturbative and nonperturbative determinations are large in the cases of Z S and Z P . This result is of particular relevance to the lattice determinations of quark masses and the chiral condensate.
quantities which are invariant under both O(4) and parity. (ii) Z q carries no flavour index f since, owing to its definition (41), it depends quadratically on r q . The proof that our estimators of the RI-MOM RCs are O(a)-improved at all quark masses and external momenta in maximally twisted QCD follows closely the argument of ref. [13] (subsequently refined in ref. [27]) for automatic improvement of correlation functions. It is based on the exact invariance of the maximally twisted QCD action under the transformation P × D d × (µ q → −µ q ), where µ q is the twisted quark mass. For completeness we recall the definition of the discrete transformations P (continuum parity) and D d in the physical quark basis P : where x P = (−x, x 0 ) andμ is the unit vector in the µ-direction. Actually the proof is also valid in the more general "mixed action" setting, where Osterwalder-Seiler fermions are used as valence quarks. This is because the corresponding action is similarly invariant under P × D d × (M q → −M q ), where M q denotes the set of all quark (sea and valence) twisted mass parameters.

The proof
The proof can be given in three steps: • Due to H(4) symmetry, the lattice form factors V q (p) and V f h Γ (p) have no terms of the form p k 0 + p k 1 + p k 2 + p k 3 and (p 0 p 1 p 2 p 3 ) k , with odd k. Thus they are invariant under parity transformation; i.e. they satisfy the relations V q (p) = V q (p P ) , V f h Γ (p) = V f h Γ (p P ) , p P = (−p, p 0 ) .
• As has been proved in [27], the invariance of maximally twisted QCD under P ×D d × (µ q → −µ q ) implies that O(a j ) cut-off artefacts of lattice correlators arise, in their Symanzik expansion, from the insertion of Symanzik operators of parity (−1) j .
• Subsequently, O(a j ) terms with odd j must then be absent in the Symanzik expansion of parity-even form factors, like V q (p) and V f h Γ (p). The Symanzik expansion of the form factors above is straightforwardly obtained from that of the correlators in terms of which they are defined.

Improvement at a generic value of the twist angle
Since RCs are UV quantities, they should be independent of the QCD infrared structure and properties; in particular they should be unaffected by spontaneous chiral symmetry breaking. Thus, if we compute RCs at zero quark mass and at high momentum scale p 2 such that chiral breaking effects are negligible, we can fully exploit the R 5 -symmetry of the continuum theory. This will be useful because in turn the lattice action is invariant under R 5 × D d , where

The proof
The proof of the O(a) improvement of RCs in lattice QCD at generic values of the twist angle, e.g. of the bare twist angle ω 0 = arctan(µ q /(m 0 − M cr )), runs on lines similar to those of the previous section, with R 5 replacing P. One can, in fact, argue as follows.
• RCs constants can be extracted from the large-p 2 behaviour of the V q (p) and V f h Γ (p) form factors in the chiral limit. In this regime chiral breaking effects can be ignored and V q (p) and V f h Γ (p) are invariant under R 5 . • The invariance of the theory under R 5 × D d implies that the R 5 -parity of an operator equals the parity of its naive dimension. Thus, O(a j ) lattice artefacts of correlators or form factors arise in the Symanzik expansion from the insertion of operators of R 5 -parity equal to (−1) j .
• Therefore O(a j ) terms with odd j cannot be present in the Symanzik expansion of R 5even form factors, like V q (p) and V f h Γ (p). In fact, the continuum target theory is invariant under the chiral group (including R 5 ), while O(a j ) terms would be odd under R 5 .

A final comment
The O(a) improvement of RI-MOM RCs at any value (including zero) of the twist angle ω 0 hardly comes as a surprise because the RC-estimators computed at generic ω 0 become ω 0independent in the chiral limit, provided they are analytic in the complex mass parameter (m 0 − M cr ) + iµ q , which is certainly the case at sufficiently large p 2 . As a consequence, a unique set of RI-MOM RCs, all free from O(a j ) (j odd) artifacts, exist for lattice QCD with Wilson quarks and a given pure gauge lattice action. These RCs can in principle be extracted from simulations at arbitrary (and not necessarily all equal) values of the twist angle. Computing the RI-MOM RCs at maximal twist (|ω 0 | = π/2), as we do in the present study, is just a simple and technically convenient choice. In our case this choice was a quite natural one, as we could then evaluate the necessary correlators on the N f = 2 gauge ensembles which were produced for studying large volume physics.