Non-perturbative renormalisation with interpolating momentum schemes

Hadronic matrix elements evaluated on the lattice can be converted to a continuum scheme such as $\MSbar$ using intermediate non-perturbative renormalisation schemes. Discretisation effects on the lattice and convergence of the continuum perturbation theory are both scheme dependent and we explore this dependence in the framework of the Rome-Southampton method for generalised kinematics. In particular, we implement several non-exceptional {\em interpolating} momentum schemes, where the momentum transfer is {\em not} restricted to the symmetric point defined in RI/SMOM. Using flavour non-singlet quark bilinears, we compute the renormalisation factors of the quark mass and wave function for $N_f=3$ flavours of dynamical quarks. We investigate the perturbative and non-perturbative scale-dependencies. Our numerical results are obtained from lattice simulations performed with Domain-Wall fermions, based on ensembles generated by RBC-UKQCD collaborations; we use two different lattice spacings $1/a \sim 1.79 $ and $2.38$ GeV. We also give the numerical values for the relevant anomalous dimensions and matching coefficients at next-to-next-to-leading order.


Introduction
One major goal of lattice QCD is to determine the hadronic matrix elements of hadronic operators, with particular emphasis on flavour non-singlet quark bilinear operators, such as the local scalar, vector, axial vector or pseudoscalar operators. To produce physically meaningful results we need renormalisation factors to relate the bare lattice results to quantities in standard renormalisation schemes, such as the modified minimal subtraction scheme, MS. This conversion can be done in lattice perturbation theory but it is well known that it results in large systematic error because the perturbative series converges poorly. Instead, in the last couple of decades, more efficient methods have been developed, the most popular ones being the Schrödinger Functional [1,2] and Rome-Southampton [3] frameworks. In this work we follow the latter. In both cases, the strategy comprises two distinct steps: First the renormalisation factors are obtained non-perturbatively from lattice simulations in a given intermediate scheme, hence the name Non-Perturbative Renormalisation (NPR). Then these factors are converted to MS providing us with a bridge between low and high energies. This second step is done perturbatively, but in the framework of the Rome-Southampton method, it can be done in standard continuum perturbation theory so that we can take advantage of the multi-loop computation available in the literature. In principle the method is relatively simple: one makes lattice measurements of Green's functions with bilinear operators inserted into quark 2-point functions (for example), and compares the lattice results with continuum perturbative results. This method can also be applied to more complicated Green functions, higher twist bilinears, four-quark operators, etc. It is worth noting that although traditionally we work in massless limit, a massive momentum scheme has been developed in [4].
In early lattice studies the operator insertion was at zero momentum transfer, see for example [3,5] for the early direction of this approach. This choice of kinematics corresponds to an exceptional momentum configuration and is subject to potential infrared issues. The origin of the problem and potential solutions were already discussed in detail in [6]. These infrared problems can become even more severe when the chiral limit is taken, and had lead to significant discrepancies in the case of neutral meson mixing [7,8]. With the improvement in lattice analyses as well as the use of dynamical fermions, the refinement of the non-perturbative structure that is of phenomenological interest became suspect due to the emergence of the latent infrared issues. Therefore in [9][10][11][12], an extension to the zero momentum insertion for this class of Green's functions was developed. In particular the Green's functions were computed with a non-zero momentum flowing through the operator and external quark legs at what is now termed the fully Symmetric MOMentum (SMOM) subtraction point * . This circumvents any potential ambiguity in assessing infrared effects and in particular in the chiral limit. In either scenario of zero or non-zero operator momentum flow lattice measurements were assisted by matching to the continuum limit. By this we mean that one can compute the same Green's function for any of the above operators in the continuum field theory at several loop orders in the MS scheme. In the zero momentum case this has been carried out to three loops, but for the nonzero momentum case this has only been carried out to two loops, [9][10][11] through there has been progress towards the next loop order in [13,14]. Such continuum matching allows one to accurately tune to the ultraviolet before extrapolating to the non-perturbative or low energy region where measurements are made and errors estimated.
The effect of using a momentum configuration which is infrared secure can be seen for example in [8,[15][16][17][18][19] where the renormalisation factors for the various flavour non-singlet quark bilinear or four-quark operators were calculated for both momentum configurations. A quantity which can be used as a benchmark is the ratio of the vector and axial vector operators. On symmetry grounds this has to be unity for Domain-Wall fermions and thus provides a stringent test of determining which momentum configuration was more reliable. It is evident from the analysis of, for example, [15,18,19] that the SMOM ratio is unity over virtually the full momentum range and well into the deep infrared in marked contrast to the zero momentum case. While this augurs well for significant improvement in lattice measurements and errors, one question naturally arises which will form the main focus of our article. Now that it is established that the SMOM non-exceptional momentum configuration gives a clear improvement, it is only one of many possible configurations that can be studied. Indeed, for example in [9,10], a more general off-shell case was developed which is termed IMOM to denote the interpolating momentum subtraction configuration. In the SMOM definition the momentum squared of the three external momenta were all equal. By contrast in the IMOM case the magnitude of the operator momentum flow is allowed to vary and this is governed by an additional kinematic parameter ω. For example the SMOM case corresponds to ω = 1. Varying ω between 0 and 4, which are the respective values giving infrared or collinear singularities, allows one to search for a value which can improve the matching to the continuum, the non-perturbative behaviour or maybe even reduce the lattice artefacts. This therefore leads to our investigation to see if there is an optimal value for ω, possibly different from one, which minimises errors and improves benchmarks such as that studied in [15,18,19]. This would be the appropriate next stage after the success of SMOM lattice studies especially in the chiral limit. For instance the ratio of vector and axial vector renormalisation constants can be examined in the chiral limit to gauge pion mass effects. While ω will be our main parameter the overall mass scale µ at which measurements are made will also be important. Therefore to have a thorough investigation we have performed a lattice computation of various operators on several different lattices for half integer values of ω between 0 and 4 for various values of µ. This has allowed us to identify an optimal area of the (ω, µ 2 )-plane where there is an effective plateau upon which measurements are reliably accurate. Moreover the picture appears to be consistent for the operators we consider, though the actual locations are not necessarily the same point for each case. It suggests that there should be a similar qualitative behaviour for operators with covariant derivatives, which we do not consider here.
Although the article primarily concentrates on the lattice computations and analysis we were indebted throughout to the various two loop MS continuum results accumulated over the last decade, [9-12, 18, 20], to facilitate the continuum matching. Moreover the measurement of quark bilinear currents in the IMOM configuration will give crucial direction to 3-point quark Green's functions such as those relevant to kaon mixing for example, [7,8]. In that scenario the effective 4-point functions can be mimicked by products of the Green's functions we consider here. It is particularly important in this case, as it has been argued that the choice of kinematics was responsible for the discrepancy observed in the literature.
The paper is organised as follows: In Section 2 we outline the kinematics setup, explaining in detail how our incoming and outgoing momenta are chosen to lead to a given value of ω for a given value of µ. In Section 3 we explain how we obtain the Zfactors, and discuss our choice of projectors. In Section 4, we present our results for the Z−factors and their scale-evolutions. We conclude in Section 5. Appendix A shows some perturbative results, such as the values of the conversion factors for matching from IMOM to MS, computed at Next-to-Next-To-Leading Order (NNLO) More details about the lattice simulations can be found in Appendix B.

Kinematics
We consider the interaction in the diagram below, where the incoming momentum is defined as p 1 , the outgoing momentum as p 2 , and the momentum transfer as q. The method is based on the Rome-Southampton regularisation independent momentum scheme: the fermions are off-shell, we work in the chiral limit and the renormalisation scale is given by the choice of external momenta. In the original RI/MOM scheme (or RI /MOM) [3], the incoming and outgoing momenta are equal, p 1 = p 2 , so that there is no momentum transfer, q = 0. The renormalisation scale µ is defined by µ 2 = p 2 1 = p 2 2 . In contrast, in the SMOM schemes, the incoming and outgoing momenta are different: they are chosen such that q 2 = p 2 1 = p 2 2 ; naturally the renormalisation scale is given by µ = q 2 . The advantage is to avoid exceptional kinematics which behave badly in the infrared regime [9].
The above schemes can be generalised further. We retain the condition p 2 1 = p 2 2 = µ 2 but drop the requirement q 2 = µ 2 . We define the variable ω = q 2 /µ 2 to parameterise the kinematics. ω = 0 corresponds to the original RI/MOM and RI'/MOM schemes, ω = 1 corresponds to the SMOM scheme. The possible values of ω range between 0 and 4. Let us denote by α the angle between the incoming and outgoing momenta. We choose our frame such that the momenta take the following form (the time component is represented by the last entry): and where we have defined ω = 2(1 − cos α). Clearly the parameter ω can take any real value between 0 and 4, and the value of ω is directly related to the value of α as shown in Fig 2. We can check that the extreme cases ω = 0 and ω = 4 correspond to p 1 = p 2 and p 1 = −p 2 , respectively, and therefore these situations suffer from infrared and collinear singularities. We also note that ω = 0 is the choice made in the original (exceptional) kinematics of a standard RI/MOM scheme while ω = 1 corresponds to the SMOM kinematic of [9]. This has been known for a while, see for example [9,10]. Different choices of ω are referred to as interpolating momentum schemes [11,20,21]. However, to our knowledge, there is no numerical study or practical implementation of such kinematics.
Let us now turn to the lattice implementation. If L is the physical spatial extent of the lattice (we do not consider the temporal extent here for simplicity), we have where l, m and n are dimensionless (they are the Fourier modes if they are chosen to be integers). In order to obtain any desired value of µ and ω we take advantage of (partially) twisted boundary conditions, which allow l, m and n to be any real numbers. For any pair (µ, ω), the numbers (l, m, n) can be obtained from the kinematic constraints p 2 1 = p 2 2 = µ 2 and q 2 = ωµ 2 by solving (2.11) In order to show a concrete example, let us anticipate and consider one of the ensembles used in this work, with spatial extent of L/a = 32 and inverse lattice spacing of a −1 = 2.383(9) GeV. In that case, we can see that the following choice: Depending on the number of minus signs in p 2 , we can obtain any integer value of ω between 0 and 4. With this choice we also expect the quark wave renormalisation factor so to have similar lattice artefacts, whereas in the setup mentioned above they will have different discretisation effects (except in some special cases). In that case, one can expect the discretisation effects to be ω-independent. We note that the same procedure with fewer than four components can lead to non-integer value of ω: if the last component is fixed to zero, we can achieve the values ω = 0, 4/3, 8/3, 4.
In this work we only consider the setup given in Eqs. (2.1)-(2.2) as we also want to investigate arbitrary non-integer values of ω. We chose to vary ω between 0.5 and 4 with a step of 1/2. In addition, we can use the ω dependence of the discretisation effects to obtain an additional handle to estimate the systematic uncertainties.

Lattice implementation
For the bilinears, omitting flavour and spin-colour indices for clarity (we only consider flavour non-singlet operators) we define O Γ =ψΓψ, with Γ = 1, γ µ , γ µ γ 5 , γ 5 (S,V,A,P). We sketch here the computation of the propagators and of the bilinears needed to renormalise the quark mass and the quark wave function. We drop the volume factors and set a = 1 for simplicity.
Following [22], the momentum source propagators are computed by first solving for a given momentum p, on Landau-gauged fixed configurations. In the previous equation, D denotes a generic Dirac matrix in spin-colour space, regardless of the lattice discretisation. Then we obtain the momentum source propagators by multiplying by the corresponding phase factor So within our conventions, G x (p) is an incoming propagator (with respect to x) with momentum p. After gauge average, the outgoing one is then given by G x (−p) = γ 5 G x (p) † γ 5 . The advantage of the method proposed in [22] compared to a traditional Fourier transom is that there sum over z in Eq.
Defining G(p) = x G x (p), we then amputate the external legs This quantity Π Γ is obviously a matrix in spin-colour space, the renormalisation factors Z Γ are then defined by projecting it to the corresponding tree level value. A possible choice (which defines a so-called γ µ -scheme) is to use the same Γ matrix as in the vertex where the trace is taken in colour space and F is the corresponding tree level value.

Z-factors
Let us now specify our conventions for the renormalisation factors. We denote the renormalised quantities with an index R, whereas the bare quantities do not carry extra indices. The quark wave function renormalisation factor Z q is defined by ψ R = Z q ψ, and G R = Z q G represents the renormalised fermion propagator. For the bilinears defined above, the renormalised quantities are denoted by Note that for the corresponding Green's functions, we have as usual The renormalisation factors for the bilinears are then obtained by imposing where m is the (bare) quark mass (before taking the chiral limit, we work with degenerate quark masses for simplicity). Note that we could take also the chiral limit by sending the renormalised mass to zero, here we assume that there is no additive mass renormalisation ‡ . There are several possible ways to define the quark field renormalisation constant Z q . For example in RI /MOM, one imposes a condition directly on the fermionic propagator (3.23) We now turn to the quark mass; The choice of the field renormalisation fixes the renormalisation factor Z m , which relates the bare and renormalised masses as m R = Z m m. To be specific, we write in Fourier space the inverse bare fermion propagator and similarly, having a generic MOM scheme in mind, we write for the renormalised propagator, Traditionally in a MOM scheme, one defines the renormalisation factors Z m and Z q by imposing that the renormalised parts are finite in the chiral limit at a certain renormalisation point. For example the renormalisation factor Z m can then be extracted from Tr G −1 (p) which is equivalent to In [9], it was suggested to replace this condition by The motivation for this definition of Z m is that it leads to Z m = 1/Z P if Z P is defined from Eq.(3.21) at the same kinematic point. This follows from requiring that the Ward-Takahashi identities (WTI) given below are satisfied both for the bare and renormalised quantities (for degenerate quark masses) in the case of degenerate quark masses. Imposing this condition with the SMOM kinematics leads to maintaining the WTI on the renormalised quantities. In particular the procedure leads to Z A = Z V .

Our definitions
We deviate from the choice made in [9] for several reasons. Firstly, some derivations in [9] rely on the SMOM kinematics, which we want to generalise here. Secondly, rather than imposing the renormalisation condition directly on the quark propagators, we take advantage of the good chiral properties of the Domain-Wall discretisation and extract the renormalisation factors from the bilinears. Thirdly, in our framework, we use the local vector and axial vector current which are related to the conserved (or partially conserved) ones by the finite "renormalisation" factors Z V and Z A respectively. In this paper we are using Domain Wall fermions, so we have Z A = Z V up to small artefacts. However the SMOM scheme can also be used for other lattice fermions, with Z A = Z V so we do not want to restrict ourselves to the case Z A = Z V .In this paper we are using Domain Wall fermions, so we have Z A = Z V up to small artefacts. However the SMOM scheme can also be used for other lattice fermions, with Z A = Z V so we do not want to restrict ourselves to the case Z A = Z V . Our values of Z A = Z V are known from previous studies on the same lattice, and are given in Appendix B for convenience. Additionally, within our framework we define Z m directly as and we also define where the shorthand IMOM stands for µ 2 = p 2 1 = p 2 2 , (p 1 − p 2 ) 2 = ωµ 2 . (These equations are used for the quark wave function and the mass renormalisation factors.) More explicitly, we define where Λ V,S are the bare (and local) vertices defined above. Note that both Z m and Z q acquire an ω-dependence through these definitions. This dependence is discussed in section 4. The main advantage of this approach is that the Green's functions of the bilinears have a much safer infrared behaviour than the quark propagators.
In order to completely fix the definition of Z q and Z m in Eqs. (3.36) and (3.37) we need to specify our the choice of projector P µ for the vector Green's function, We implement the so-called γ µ and / q-projectors (the trace is taken over both Dirac and colour indices) where obviously q = p 1 − p 2 is the momentum transfer of the IMOM kinematic defined above, µ 2 = p 2 1 = p 2 2 , (p 1 − p 2 ) 2 = ωµ 2 . Plugging these definitions in Eqs (3.36) and (3.37), we see that we have defined two kinds of IMOM schemes similarly to was done in the SMOM case. In order to keep track of what projector was used, we introduce the respective notations Z The renormalisation factors for the other bilinears are defined exactly in the same way as Eqs 3.34,3.35. This, together with a lattice formulation which preserves chiral symmetry, leads to Z V = Z A and Z S = Z P , apart from the physical effects of spontaneous chiral symmetry breaking, which may become important at low p 2 and q 2 . This is discussed in Section 4. Naturally for a lattice formulation which break explicitly chiral symmetry, one can follow [9] with the modification q 2 = ωµ 2 .

Nomenclature
To define unambiguously a momentum scheme for a composite operator one needs to specify the kinematics (the choice of momenta) and the choice of projectors, not only for the Green's function corresponding to the operator in question, but also for the quark wave function. § Historically, when the Rome-Southampton method § Strictly speaking we also need to specify the gauge, but it is standard to assume that the computation is performed in the Landau gauge. was introduced in [3], the kinematics was the exceptional case, p 1 = p 2 and Z q was defined through a γ µ -projector, i.e. Eq. (3.39). Although the authors of [3] mention explicitly that other choices were possible, by convention "RI/MOM scheme" refers to this specific choice of kinematics and wave function. Similarly, it also became standard to call "RI /MOM scheme" a scheme in which the kinematics is still the exceptional case but the quark wave function is renormalised through a / q-projector, ie Eq. (3.40). (In [3] the definition of the projector differs from the one given here, but both definitions are equivalent up to lattice artefacts).
We turn now to the SMOM schemes. By definition, the "S" of "SMOM" refers to the situation where p 2 1 = p 2 2 = (p 1 − p 2 ) 2 , ie ω = 1, or "Symmetric" [9]. For the quark wave function, the convention in [9] was that "RI/SMOM" refer to a / q-projector and "RI/SMOM γµ " to a γ µ -projector. Strictly speaking there is also a subtlety in the definition of Z m , see section 3.2.
Here, following the conventions adopted by RBC-UKQCD (see for example [8]), we use the notations Z (X) (µ, ω), where X ∈ (γ µ , / q) indicates the choice of the projector for the wave function. Obviously the choice of kinematics is explicitly given by ω. When needed, we also refer to these schemes as IMOM (γµ) and IMOM ( / q) .

Name
Kinematics

Non-perturbative running
To compute the non-perturbative scale evolution of Z m : we define Σ m as where X can either be γ µ or / q. We can then take the continuum limit a 2 → 0 and define Similarly, for the quark wave function we define σ q , by just substituting Z q for Z m in the above definitions. Since in practice we have only two lattice spacings, we add a systematic error as an estimate of the residual discretisation effects. This error is obtained by adding half the difference beween the extrapolated value and the value on the the finest lattice spacing (see Appendix D). We now have for each value of µ, µ 0 , ω, ω 0 a central value of σ m,q , some statistical error and a systematic error from the continuum extrapolation In the following we will always sum the systematic and statistical uncertainties in quadrature except when stated otherwise. We found that the chiral extrapolations are well under control and that an associated systematic error would be negligible compared to the other sources of error. The main motivation for looking at these quantities rather than the Z-factors themselves is the existence of the continuum limit. We can also compare the non-perturbative running obtained from our lattice simulation with the perturbative prediction.

Results
To determine the relationship between lattice results and MS results we need to compare lattice and perturbation theory in a kinematic region where both perturbation theory and lattice results are reliable. In this region the perturbative and lattice results should run in the same way, so the ratio should be constant. We will see deviations from constancy in regions where perturbation theory converges slowly (for example at small µ or small q 2 = ωµ 2 . We will also see non-constancy if lattice artefacts are important, which is likely to happen at large µ 2 and large ωµ 2 .
We also note that our procedure to estimate the residual discretisation errors could potentially underestimate them in the large µ 2 and large ωµ 2 regions.
In the following sections we will first keep ω constant and compare the µ dependence of the lattice results and the continuum. An advantage of IMOM is that we have two parameters, µ 2 and ω, so we can better test that we are working in a region where it is sensible to trust both perturbative and lattice results. Accordingly, we next will vary both µ and ω and study if a plateau is emerging where we can reliably determine Z.

Quark mass renormalisation factor
We now apply this scheme to determine Z for the mass operator and quark wavefunction. The running of the quark mass is shown in Fig 3 for the γ µ -scheme and Fig. 4 for the q /-scheme. In these plots we show both the non non-perturbative scale evolution σ (γµ) m (µ, µ 0 , ω, ω 0 ) and the perturbative prediction u (γµ) m (µ, µ 0 , ω, ω 0 ). In order to study the µ-evolution for fixed values of ω, we set ω = ω 0 = 0.5, 1.0, 1.5, . . . , 4.0 and let µ vary between 1 and 4 GeV. We find a good agreement for intermediate values of µ and ω, where both perturbation theory and lattice artefacts are expected to be under good control. In fact perturbation theory works surprisingly well even for small values of µ, where significant discrepancies with the lattice results only emerge at µ 1GeV. Out of the two schemes, perturbation theory and lattice results agree best in the γ µ -scheme. The onset of lattice artefacts for large values of µ and ω becomes only relevant for q 2 25GeV 2 . This becomes particularly visible for large values of ω = 4, where perturbation theory also becomes less reliable. The discrepancy of lattice results and perturbation theory occurs already at values of q 2 10GeV 2 in the q /-scheme. We show our results for the scale evolution of the quark mass in the (µ, ω) plane in Fig. 5 for the γ µ -scheme and in Fig. 6 we show the same quantity for the / q-scheme.
Here we allow ω = ω 0 and show the ratio of the non-perturbative running divided by the perturbative prediction as a function of µ and ω at LO, NLO and NNLO, while fixing ω 0 = 2 and µ 0 = 2.5 GeV. The discrepancy from one can then be compared with the combined lattice and perturbative uncertainties that have been added in quadrature.
As expected we observe that increasing the order of the perturbative expansion improves the agreement with the non-perturbative evolution. The corner of the planes are affected by larger systematic errors, in particular where ω = 4 and/or µ > 3.5 GeV where the discretisation effects become more sizeable. However our data agree with NNLO at a few percent level for a large part of the (µ, ω)plane: approximately for 1.5 GeV ≤ µ ≤ 3.5 GeV and 1 ≤ ω ≤ 3.    For completeness, we study the next order in perturbation theory for the case ω = 1, as the three loop matching has been computed in [13,14] and the MS fourloop anomalous dimension in can be found in [23]. We take µ 0 = 2 GeV and µ = 2.5 GeV and compute the running in different schemes, see Table 1. In Table 2, we show the contributions of the higher orders. Although the corrections to the leading order contributions are rather small, we see that the convergence in MS seems to be much better than in the MOM schemes, in the sense the (i + 1)-th correction is much smaller than the i-th order. However, the important observations here are that the leading order gives by far the main contribution and that the perturbative series seems to converge nicely to the non-perturbative values ¶ .   Table 2: Study of the convergence of the perturbative series for running of the quark mass between 2 and 2.5 GeV in MS, SMOM-γ µ and q /. ¶ The attentive reader would have noticed that in Table 1 The NLO for γ µ and q / are identical. This seems to be nothing but a numerical accident: adding more significant figures lead to 0.941112 for γ µ and 0.944132 for q /.
ω/µ =    Tables 3 and 7 we give our results for the nonperturbative running of the quark mass, with µ 0 = 2 GeV and ω = ω 0 = 0.5, 1.0, . . . , 4.0. In Appendix C, we show more numerical results, where we vary our parameters ω, ω 0 , µ and µ 0 in various ways. Our conclusion remains that the non-perturbative results agree very well with the perturbative ones as long as we stay from the corner of the (ω, µ) plane.   Table 6: Study of the convergence of the perturbative series for running of the quark wave function between 2 and 2.5 GeV in MS, SMOM-γ µ and q /.

Quark wave function renormalisation factor
We also study the convergence for the quark wave function, Z (X) q , see Table 5 and Table 6. Here the situation is very different from the quark mass mainly as there is no contribution at leading order (in the Landau gauge). In the q /-scheme, the non-perturbative results differ significantly from the perturbative predictions; however one should keep in mind that the running is very small in magnitude. It is also interesting to have a close look at the q / scheme where there the perturbative prediction is known at N 3 LO. As we can see the perturbative series to converge very poorly in the sense that the relative difference decrease very slowly as we increase the order. The difference between the non-perturbative result and the N 3 LO, namely ∼ 1.0195 − 1.0113 ∼ 0.0082, could be explained by higher order in the perturbative series.
Keeping these facts in mind, we show our results in Fig. 7 in 8 for X = γ µ and X = / q respectively. Once again we divide the non-perturbative running by the perturbative prediction, we fix ω 0 = 2 and µ 0 = 2.5 GeV and let ω and µ vary.

Study of the ω-(in)dependence of σ q
Due to the vector WTI, Eq 3.31, we might expect Z (q /) q to be ω-independent up to lattice artefacts. This is illustrated in Fig. 11 where we show σ (q /) m (µ, µ 0 , ω, ω 0 ) as a function of ω = ω 0 for different values of µ and µ 0 . We can see that this property is rather well satisfied as long as the scales remain moderate, say µ, µ 0 ≤ 3 GeV. However it is worth noting that this invariance is only true in the continuum, as can be seen in Fig. 12. At finite lattice spacing, within our small statistical errors, the lattice artefacts are visible and Z q / q clearly show a ω dependence. These figures also suggest that the lattice artefacts are well under control, with two lattice spacings, as long we do not go too high in energy. Note that in these plots, the systematic errors due to the continuum extrapolations are not included. The non perturbative-scale evolution of Z q , σ q (µ, µ 0 ) with µ 0 = 2 GeV and ω = ω 0 for various values of µ and ω can be found in Tables 7 and 8.
Also it is interesting to notice in Fig 12 how the discretisation effects depend on ω. Clearly for ω = 2.0 and 2.5 the lattice artefacts are much smaller than for ω = 1.0. Of course this effect is quantity-dependent but it this could be useful in future computations.

Study of the chiral symmetry breaking effects
One of the original arguments to motivate the RI/SMOM schemes is a drastic reduction in the effects of spontaneous chiral symmetry. Even though these effects are physical, they can prevent a clean determination of the renormalisation factors because they are absent from the perturbative calculations. This is particularly true for quantities like Π P and Σ S where the presence of pseudo-Goldstone poles can completely dominate the signal. In practice, the vertex functions from which we want to extract the Z-factors are 'polluted' by negative powers of the quark mass or the momentum scale. Of course these chiral symmetry breaking effects are nonperturbative and disappear for high-momentum. However, in practice we want to keep the Rome-Southampton windows as open as possible, so it is always desirable to reduce these infrared contamination or at least keeping them under control.
If chiral symmetry is exactly realised then we should find that Z S = S P and Z V = Z A . Therefore, in order to study these chiral symmetry breaking affects as a function of ω, we study the deviation of Z V from Z A and Z S from Z P . We start by the analysis of the ratio Z V /Z A , see Fig. 13 for ω = 1 using the the γ µ -projector. This ratio collapses quickly to one as the energy scale increases: at µ = 2 GeV, the deviation from one is at most at the per-mille level. We also find the quark mass dependence to be linear within our statistical error -as expected -therefore we extrapolate to the chiral limit using a straight line. Our results for the q /-projector are very similar.
In Fig. 13 we show our results in the chiral limit for the different values of ω and find that Z V /Z A is always very close to one. At our lowest energy point µ = 1 GeV, we can observe a deviation from one at the order of a percent. Our results seem to indicate a trend that Z V /Z A increases toward one when ω increases between 0.5 and 2 (Z V /Z A ∼ 0.98 at ω = 1 and Z V /Z A ∼ 0.99 at ω = 2), but this could well be a statistical effect. All together we find that Z V /Z A = 1 to a very good approximation for all values of ω. This true for both values of the lattice spacing and both projectors X ∈ (γ µ , q /). We also study (Λ S − Λ P )/Λ V , which is proportional to Z V (1/Z S − 1/Z P ). Our results are shown in 15. Here it is clear that the IR effects due to chiral symmetry breaking are affected by the value of ω. We can observe that at µ = 1.5 GeV, (Λ S − Λ P )/Λ V ∼ 0.1. The same quantity reduces to ∼ 0.03 for ω = 2.0. It will be interesting to perform a similar study on four-quark operators, especially those where the IR contamination can be important and a source of disagreement (see for example the section on BSM kaon mixing in FLAG [24] and [8]).

Conclusions and outlook
As a proof of concept we have implemented several IMOM schemes, defined via two different projectors, and for various kinematics. In this framework we have determined the renormalisation factors and non-perturbative scale evolution functions of the quark mass and of the quark wave function. To compare our lattice results with NNLO results in continuum perturbation theory we also present the numerical form of the perturbative scheme conversion factors for these schemes in these general kinematics.
We find that the non-pertubative and perturbative results agree very well as long as we stay from the corner of the ω, µ plane, with one exception, namely Z (q /) q . Where, we argued that the reason for this relatively bad agreement is the poor convergence of the perturbative expansion. Clearly having several values of ω help to have a better handle on the systematic errors coming from the NPR procedure. As an application, we have seen a example where the discretisation effects depend significantly on ω (see σ (q /) q where the discretisation effects for ω = 2, 2.5 are much smaller than for ω = 1) as well as the perturbative convergence. In this proof of concept study only two lattice spacings have been used. Adding a finer lattice could potentially allow us to test the agreement of the perturbative and non-perturbative window even further.
Generally speaking, it is well-known that the SMOM kinematics (ω = 1) leads to much cleaner determinations of the renormalisation factors than the ω = 0 case. In this work, we have shown that increasing the value of ω leads to smaller contributions of the pseudo-Goldstone poles contamination in Λ P − Λ S . It will be interesting to extend this study to the case of four-quark operators where the infrared contaminations due to chiral symmetry breaking are significantly more sizeable. We warmly thank our colleagues of the RBC and UKQCD collaborations. We are particularly indebted to Peter Boyle, Andreas Jüttner, J Tobias Tsang for many interesting discussions. We also thank Peter Boyle for his help with the UKQCD hadron software. We wish to thank Holger Perlt for his early contribution in this area. N.G. thanks his collaborators from the California Lattice (CalLat) collaboration and in particular those working on NPR: David Brantley, Henry Monge-Camacho, Amy Nicholson and André Walker-Loud.    Table 9: Values of the strong running coupling as a function of the dimensional regularisation scheme ν. They are obtained in the MS scheme in the three-flavour theory.

Appendices A Perturbation Theory
In order to make contact to phenomenology, results obtained from lattice simulations need to be converted to a scheme that is used in continuum perturbation theory such as MS. The non-perturbative renormalisation schemes used in this work serve as an intermediate scheme that can be defined on the Lattice and in the continuum. This enables us to compare the non-perturbative with the perturbative change of the renormalisation scheme parameters µ and ω. This appendix provides the ingredients for the perturbative scheme changes at NNLO.
The dependence of the strong coupling constant on the dimensional regularisation scale ν in the MS scheme is governed by the β-function Although it is known at five-loop [25][26][27], we will only use the first three coefficients in the expansion to be consistent with our NNLO analysis of the running. In the following we will need the first two coefficients where we fixed the number of colours N c = 3 in this and the following expressions. Using α S (M Z ) = 0.1179 (10) and M Z = 91.1876(21) GeV [28] we find the required values of the strong coupling constant for N f = 3 flavours given in Table 9. To arrive at three active flavours, we decoupled the bottom and charm quark at ν b = 5 GeV and ν c = 1.5 GeV respectively using central values m b (m b ) = 4.13 GeV and m c (m c ) = 1.27 GeV for the quark masses in the MS scheme [28]. The relevant threshold effects and running of the strong coupling constant are implemented using RunDec [29]. The fields and masses in the IMOM (X) schemes are related to the MS scheme as ψ MS (ν) , where X ∈ (γ µ , / q). The continuum perturbation theory expansion of the conversion factors   are known up to two loops [10,20] for the IMOM scheme and up to three loops [13,14] in the SMOM limit, i.e. where ω = 1. Writing and setting µ = ν we find the numerical conversion factors as a function of ω given in Table. 10. The scheme transformation can then be written as a product of the conversion factor, its inverse, and the and the MS evolution kernel U MS i (ν 1 , ν 0 ). The evolution Kernel fulfils the renormalisation group equation and we expand the anomalous dimensions where X denotes the renormalisation scheme. Hence we can transform for example a light quark masses renormalised in given scheme X at different kinematic points (µ 1 , ω 1 ) and (µ 0 , ω 0 ) via where the logarithms log µ 1 /µ 0 are summed using renormalisation group improved perturbation theory. Explicitly, we find at NNLO X and J (2) X are given by The leading order anomalous dimensions are scheme independent. The expressions for the NLO and NNLO anomalous dimensions are given in Table 11 for different values of ω as a function of the number of flavours N f .

ω B Simulation Details
Our numerical work is based on RBC-UKQCD data, the lattice details can be found in [30]. We compute the propagators using Landau gauge-fixed 2+1, Domain-Wall (Shamir [31])/Iwasaki lattices. The values of the parameters can be found in [7].
In a nutshell, we use two lattice spacings (we refer to them as 24 3 C.1 Study of the ω-dependence for fixed energy scales In order to study the ω-dependence, we first fix µ and µ 0 to some reasonable values, where we expect a rather good control over both the perturbative errors and the lattice artefacts. We choose (µ 0 , µ) = (2.5 GeV, 1.5 GeV). We then vary ω, but first we only consider the 'diagonal' case, ie ω 0 = ω. We show our results in Tables 12  and 13. The middle column shows the running itself obtained from the lattice, while the other three columns on the right show the ratio of the non-perberturbative running over the perturbative prediction at Leading Order (LO), Next-to-Leading Order (NLO) and Next-to-Next-to-Leading Order (NNLO). The errors quoted there combine an estimate of the discretisaiton errors and the statistical one. We observe the non-perturbative running agrees extremely well with NNLO predictions for all values of ω ≤ 3. Looking at the perturbative convergence and stability, our data seem to favour the region ω ∼ 2.
C.3 Study of the running in the non-degenerate ω case.