Non-perturbative running of renormalization constants from correlators in coordinate space using step scaling

Working in a quenched setup with Wilson twisted mass valence fermions, we explore the possibility to compute non-perturbatively the step scaling function using the coordinate (X-space) renormalization scheme. This scheme has the advantage of being on-shell and gauge invariant. The step scaling method allows us to calculate the running of the renormalization constants of quark bilinear operators. We describe here the details of this calculation. The aim of this exploratory study is to identify the feasibility of the X-space scheme when used in small volume simulations required by the step scaling technique. Eventually, we translate our final results to the continuum MSbar scheme and compare against four-loop analytic formulae finding satisfactory agreement.


Introduction
In lattice Monte Carlo simulations of QCD, connecting the high energy regime where perturbation theory can be safely applied, to the low energy regime of large volume simulations where hadron matrix elements are evaluated, usually cannot be done using a single ensemble. The problem to safely treat the two very seperated scales within a single lattice extent and lattice spacing is usually called the 'window' problem and is summarized by the following equation in momentum space or equivalently as in position space. The latter says that although the relevant physical distances at which we evaluate our observables should be sufficiently small, we are limited by the discretization errors. At the same time, when we try to study larger distances where cut-off effects are expected to be smaller, we may approach the intrinsic scale of QCD too closely and perturbation theory will not be applicable. A method to connect those two regimes while keeping all systematic errors under control was developed in Refs. [1,2,3,4]. It is called step scaling and consists in performing several simulations at a smaller and smaller lattice spacing as well as smaller and smaller volumes starting at the coarse, large volume level and repeating the procedure until the mentioned high energy regime of the theory can be reached.
The most popular implementation of the step scaling method is used to compute the running of the strong coupling constant and uses Schrödinger functional (SF) boundary conditions [5]. Although it necessitates the implementation of boundary conditions which differ from the ones usually used in large volume simulations, such framework offers many advantages. By choosing apropriate field values at the time boundaries t = 0 and t = T , one removes problematic zero modes of the gauge field, hence significantly simplifying perturbative calculations. By the same token, the Dirac operator develops an energy gap which allows to simulate dynamical fermions at vanishing quark mass. The SF framework offers also additional freedom in constructing correlation functions, which can be used to implement renormalization or improvement conditions that have significantly reduced cut-off effects. Within the SF scheme, also the running of renormalization factors for the first moment of the quark non-singlet parton distribution has been computed [6,7].
An alternative implementation of the step scaling technique was described in Ref. [8]. It uses the RI-MOM renormalization scheme [9], with twisted boundary conditions in order to enable an easy implementation of renormalization conditions at a fixed physical momentum on the set of ensembles. Hence, instead of tuning the parameters of the entire simulated volume to stay on the line of constant physics, one chooses the twist angles such that the particular momentum used to define the renormalization condition remains fixed. The results of this implementation were used for phenomenology among others in Ref. [10].
In the current work, we investigate an alternative to the RI-MOM implementation of the step scaling technique. We test renormalization conditions imposed on correlation functions in coordinate space (X-space) following Refs. [9,11,12,13,14] at sufficiently small distances. The same renormalization scheme can be used in a large volume simulation, hence we avoid implementing different boundary conditions. The aim of our work is to test the precision and the cost of this implementation of the step scaling technique.
Due to the known difficulties of perturbation theory for QCD in a small periodic box (see for example Ref. [15] or the recent discussion in Ref. [16] in the context of gradient flow), we look at correlation functions at very short distances (of the order of 0.01-0.1 fm) compared to the physical extent of the simulated box (which is an order of magnitude larger). This allows us to use infinite volume perturbation theory to translate our results to the MS scheme. This solution comes with a certain price, as we have to simulate effectively large lattices in lattice units.
The remainder of the present paper is organized as follows. In Section 2, we recall the most important details of the X-space renormalization scheme and the characteristics of the twisted mass formulation of lattice QCD which is used as the discretization of valence quarks. We describe the generation and tuning of the small volume ensembles. We pay particular attention to the matching of physical volumes, which allows us to perform continuum extrapolations along the lines of constant physics. Next, we consider possible sources of systematic errors in Section 3. We discuss, in particular, uncertainties related to the mismatch of the tuned ensembles, discretization effects and finite volume effects. Finally, we present numerical results for the running of renormalization constants between energy scales of around 1.5 GeV and 17 GeV in all four channels (pseudoscalar, scalar, vector and axial vector) in Section 4. Discussion of results and an outlook is presented in Section 5.

Lattice setup
This feasibility study is performed in the quenched approximation. This means that sea quarks are infinitely massive, i.e. effectively there are N f = 0 active quark flavours. Therefore, the obtained results, extrapolated to the continuum limit, can be compared with continuum perturbation theory setting N f = 0.
To generate the configurations, we use the CHROMA software [17], in its scalar or parallel version, depending on the lattice size. The chosen simulation algorithm is the standard heatbath with four overrelaxation steps.

Valence quarks
The computation of X-space correlators needs an introduction of valence quarks and for this, we use the twisted mass (TM) formulation of lattice QCD [18,19,20,21], which is given in the so-called twisted basis by where τ 3 is the third Pauli matrix acting in flavour space and χ l = (χ u , χ d ) is a two-component vector in flavour space, related to the one in the physical basis by a chiral rotation. m 0 and µ v are the bare untwisted and twisted quark masses. The massless Wilson-Dirac operator D W reads: where ∇ µ and ∇ * µ are the forward/backward covariant derivatives. One of the major advantages of TM fermions is that they are automatically O(a)-improved by tuning just one parameter, the twist angle ω, to maximal twist (ω = π/2). This can be achieved by setting the hopping parameter κ = (8 + 2am 0 ) −1 to its critical value, such that the PCAC quark mass vanishes [19,22,23,24,25,26].

Coordinate space (X-space) renormalization scheme
The X-space renormalization scheme was initially suggested in Ref. [9] and first applied in Refs. [11,12,13] in the quenched approximation. The first application in the dynamical case was reported by us in Ref. [14], where important improvements were implemented, and recently the method has been used, with further improvements, in Refs. [27,28].
Here, we shortly summarize the main ideas of this renormalization scheme. For details, we refer to the above mentioned references. We consider flavour non-singlet correlation functions of two operators of the form where We impose the following coordinate space conditions in the chiral limit, where we denote four-vectors with capital letters, e.g. X = (x, y, z, t) and X 0 is the renormalization point, which we choose to satisfy in order to keep discretization effects under control and to ease contact to perturbation theory. Here, Λ denotes a low-energy scale of the considered theory and is of the order of a few hundred MeV. It was shown in Ref. [14] that the method works best for the so-called "democratic" points, which are defined as points with a direction close to the diagonal of the hypercube (with an angle of at maximum 30 degrees between a given vector X and (1, 1, 1, 1)). However, further studies of this issue in the free theory suggest that certain points with an angle of 45 degrees are also expected to yield relatively small cut-off effects.
To alleviate the impact of cut-off effects, we use a tree-level correction, as defined in Ref. [14]. This boils down to computing the ratio of the tree-level lattice and continuum correlators, where c = 3 for (pseudo)scalar correlators (later referred to also as PP/SS correlators), c = 6 for (axial) vector ones (AA/VV correlators). The corrected correlation function, C ′ Γ (X), reads The (tree-level corrected) renormalization constants, at the scale µ = 1/ X 2 0 , are then given in the X-space scheme by In a certain case, we will also consider non-corrected renormalization constants, and we will indicate it explicitly when this is the case (else we use the corrected ones). One now needs to decide which points X 0 are best suited to extract the renormalization constants. We discuss this issue in the next subsection.

Step scaling
In the step scaling method, one aims at estimating the step scaling function describing the running of a scale dependent quantity as the renormalization scale is changed by some factor. In this work, we fix this factor to 2. The finite lattice spacing version of the step scaling function can be constructed as a ratio of renormalization constants in the X-space scheme at two renormalization scales which differ by a given factor. Then, such ratio, which we denote by Σ X Γ (µ, 2µ), has a well-defined continuum limit, with µ = 1/ X 2 0 and where we have added explicitly the lattice spacing a as an argument of Z X Γ to indicate that it was regularized on the lattice. The step scaling function enables one to non-perturbatively compute the renormalization scale running of scale-dependent quantities.

Main characteristics of the setup
We decide to perform three repetitions of the step scaling procedure, i.e. we start at a renormalization scale µ 0 and then compute succesively Σ X Γ (µ 0 , µ 1 = 2µ 0 ), Σ X Γ (µ 1 , µ 2 = 2µ 1 ) and Σ X Γ (µ 2 , µ 3 = 2µ 2 ). As it will turn out, these three values of the step scaling function will allow us to link nonperturbatively the scales ranging from around 1.5 GeV up to 17 GeV. The lower energy is accessible in a large volume simulation where the hadron matrix elements are evaluated, whereas the upper energy allows for a safe translation to the infinite volume continuum MS scheme.
In order to stay on the line of constant physics, we fix the renormalization scale to be a fraction of the total simulated spatial extent, Since we always extrapolate our results to the chiral limit, µ, L and the the lattice spacing are the only relevant scales in our problem. Hence, the estimation of Σ X Γ (µ, 2µ) necessitates two sets of lattices: one with spatial extent L and the second with 2L.
Our starting point is an ensemble withL ≡ L/a = 24 at β = 9.0 (with a corresponding, but unknown before the step scaling procedure, lattice spacing) which fixes our spatial extent L in the first step of our procedure. We  denote it by L 1 , and the corresponding one of theL = 48 lattice by L 1 = 2L 1 . We fixL to 8, 16, 24 and 32 and find the corresponding lattice spacing (or inverse bare coupling β) such that the spatial extent is always equal to L 1 (we describe the details of this procedure in the following section). We end up with a set of lattices as described by the first line of Tab. 1. We have at our disposal four pairs of ensembles with equal physical volumes (and their doubles) and with different lattice spacings. This defines our line of constant physics and allows to perfom the continuum extrapolation of the step scaling function.
Having four estimates of the step scaling function Σ X Γ (µ, 2µ, a) at different lattice spacings, allows us to extrapolate to the continuum and compute the step scaling function Σ X Γ (µ, 2µ), Eq. (15). In the second step of the step scaling procedure, we work in a larger physical spatial extent L 2 = 2L 1 and L 2 = 2L 1 . Hence, we can reuse some of the ensembles from the previous step. For example, theL = 16 ensemble at β = 7.90 used in the first step can be used again in the second step. Similarly, in the third step of the procedure we can reuse two ensembles generated during the second step. Tab. 1 summarizes all the β andL values needed for our study. For the coarsest ensemble, with the physical spatial extent L 3 = 4L 1 and β = 6.61, we can make contact with large volume simulations where the dimensionful lattice spacing can be fixed in terms of the r 0 distance (see Sec. 2.6).
In the discussion above, we demonstrated the concept of the step scaling method using a particular renormalization scale given by µ 0 = 1/ √ X 2 with X/a = (1, 1, 1, 1). However, one has also the possibility of computing the step scaling function for other renormalization scales as well. In particular, we repeated our analysis for the renormalization scales given by µ 0 = 1/ √ X 2 with X/a = (1, 1, 1, 0) (four permutations averaged over 1 ) and µ 0 = 1/ √ X 2 with X/a = (1, 1, 0, 0) (six permutations with averaging). We excluded the data for the type of points X/a = (1, 0, 0, 0), since they are known to be affected by very large cut-off effects. We label the three above possibilities by points of type IV, III and II, respectively.

Translation to MS
Since we want to compare finally with the running obtained in the MS scheme, after checking that our volumes are large enough, we convert Σ X Γ (µ, 2µ) to the latter, denoted by Σ MS Γ (μ, 2μ), using 4-loop conversion formulae [29]. Note thatμ = 2e −γ E µ ≈ 1.12µ, due to using a redefined scale in this reference. For more details on this step, we refer to our earlier work [14], Section 2.4. We note that in the following, we skip the superscript indicating the renormalization scheme, since it is always clear from the context which one is meant at a given stage.

Matching of ensembles for step scaling
As we discussed in the previous section, in order to perform three steps of the step scaling procedure, we need a set of gauge ensembles with a broad range of bare coupling values, β ∈ [6.61, 9.50] and lattice sizes, L/a = 8, 12, 16, 24, 32, 48, 64. Moreover, ensembles with L/a = 8, 16, 24, 32 of the same step scaling step need to have matched physical volumes, such that one can relate the energy scale yielded by a given type of points of one ensemble to the energy scales for the other ones (i.e. that the scale µ = 1/ √ X 2 corresponds to the same value in physical units for e.g. pointŝ X = (1, 1, 1, 1) (L = 8),X = (2, 2, 2, 2) (L = 16),X = (3, 3, 3, 3) (L = 24) andX = (4, 4, 4, 4) (L = 32)). Such matching for lattices of size L 1 /a 1 and L 2 /a 2 can be done in terms of an effective renormalized coupling 2 . In this work, we define it using a certain ratio of Wilson loops proposed by Creutz [30], where W (l, t) is a Wilson-loop with spatial size l and temporal size t. Plugging W (l, t) ∝ e −αt/l (for t ≫ 1), we get Hence, the effective coupling can be defined as For l, we always use l = L/4. If the thus defined renormalized couplings are the same (up to the precision of the calculation) on lattices with sizeL 1 = L 1 /a 1 andL 2 = L 2 /a 2 , this implies that L 1 ≈ L 2 (physical volumes are the same) and, equivalently,  In Fig. 1, we show an example extraction of the effective coupling, for one of our ensembles, at β = 8.62 on a 16 3 × 64 lattice 3 . The effective coupling reaches a plateau at large t/a, as expected. We adopt a systematic procedure to extract a value from the plateau that takes into account the possibility of choosing different fitting ranges. The procedure is analogous to the one described in Refs. [31,32]. In short, we consider all possible fit ranges [t min /a, t max /a] with t min /a = 25, 26, . . . and t max /a = 39, 40, . . ., i.e. each fit contains at least 15 consecutive data points. Then, we build a weighted histogram of all fit results and define the systematic error as the average difference of the 16th/84th percentile with respect to the median (such that 68% of the fit results are within one systematic error from the median). The weights for each fit are given as exp(−χ 2 /dof), i.e. fits with larger values of χ 2 /dof are suppressed. The smallest value of t min /a is chosen in such a way that even the worst fits in terms of χ 2 /dof still contribute to the weighted histogram (i.e. have χ 2 /dof only slightly larger than 1, with best fits having χ 2 /dof ≈ O(0.1 − 0.2)). However, it is important to observe that the procedure is almost independent of the chosen smallest t min /a -fits that start before the actual pleateau has been reached are almost completely suppressed by the exponential weighting factor and hence only fits with χ 2 /dof 1 have a significant contribution to the finally extracted value of the effective coupling. Note that the number of fits that enter the histograms is of the order of a hundred or a few hundred (depending on the temporal extent of the lattice) and hence the histograms are to a good approximation Gaussian. For each case, we repeat the whole procedure on 1000 bootstrap samples (with blocking) to obtain also the statistical error of the median and confirm that autocorrelations are under control (in all cases, we find τ int = O(0.4 −0.8) for the Wilson loop observables). In the example shown in Fig. 1, the dominant error is the statistical one. However, in several cases, also the systematic error is important and thus our error estimation procedure is essential to obtain reliable results.
Having established the procedure to extract the effective coupling, we can perform the matching of lattices. To show the reliability of our matching procedure, we show an example (Fig. 2) using an auxiliary ensemble (not used for evaluating the step scaling function) with inverse bare coupling β = 6.92 and lattice size 24 3 ×96, for which we can cross-check the result by comparing to the prediction from Eq. (2.6) in Ref. [33]. This formula expresses the Sommer parameter r 0 /a as a function of β in the interval β ∈ [5.7, 6.92] and yields that β = 6.92 with L/a = 24 has the same volume as β = 6.613 and L/a = 16. In our procedure, we want to find values of β that yield the same effective coupling (and thus the same physical volume) on lattice sizes L/a = 16. We have generated several 16 3 × 64 ensembles with values of β in the interval β ∈ [6.4, 6.75] and extracted the effective coupling for each ensemble. Then, we have performed a linear fit (red solid line) to all data points and found the matching point as the intersection of the line from the fit with the line given by the effective coupling for the β = 6.92, 24 3 × 96 ensemble. The corresponding inverse bare coupling is β = 6.618. Considering the error of the effective coupling for β = 6.92, we also determined the error of the matching value of β to be 0.017, by observing intersections of the fit line with lines representing α eff ± 1σ. In Sec. 3, we discuss how to propagate this matching error to the level of renormalization constants, to account for effects of non-ideal matching. Comparing the value 6.613 from Eq. (2.6) in Ref. [33] with our estimate from the Creutz ratio effective coupling, β = 6.618 (17), we conclude that our matching procedure is consistent with Ref. [33], well within the uncertainty that we find. This gives us additional confidence in the adopted matching procedure. Note also that the inverse bare coupling value β = 7.90 appears in all three steps of the step scaling procedure, hence making the matching between different steps more robust.
Our next step is tuning the ensembles to maximal twist. We discuss it in the next subsection.

Tuning to maximal twist
Twisted mass lattice QCD enjoys the property of automatic O(a)-improvement when tuned to maximal twist [19,22,23,24,25,26]. The commonly used criterion to achieve maximal twist is the vanishing of the PCAC quark mass, m PCAC . In practice, it is enough if the criterion m PCAC ≤ 0.1µ v holds [34].
For all our values of β resulting from the matching procedure, we adopt the following strategy. We choose the smallest available volume such that L/a ≥ 16, we choose our intermediate valence quark mass, aµ v = 0.0032 (in addition, we also use aµ v = 0.0016 and 0.0048, see below), and we guess a value of κ 1 for which the PCAC mass is small. Then, we compute the PCAC mass at another value of κ = κ 2 , closer to its critical value, κ c (we know that the sign of the derivative ∂am PCAC ∂κ is negative). This allows us to estimate the derivative ∂am PCAC ∂κ numerically, assuming linear dependence of m PCAC on κ. This, in turn, gives us an estimate of κ c . We illustrate this prescription graphically in Fig. 3 (left). The PCAC masses for the two inital guesses,  Fig. 3 (left)), a value compatible with the criterion m PCAC ≤ 0.1µ v . For the other two valence quark masses, aµ v = 0.0016 and aµ v = 0.0048, we assume that the value of κ c can be taken as the one found for aµ v = 0.0032. The right panel of Fig. 3 illustrates that this is a reasonable procedure -the dependence of am PCAC on aµ v is rather small (and approximately linear) for the considered changes in aµ v .
We take the value of κ c found in the above way also for other volumes at the same value of β. We monitor the value of am PCAC for all ensembles and we observe that the condition m PCAC ≤ 0.1µ v is always satisfied within statistical uncertainty. Hence, we conclude that tuning to maximal twist is very robust for our ensembles and thus we can safely take our continuum limits assuming O(a 2 ) leading cut-off effects.

Summary of generated ensembles
Finally, we list in Tab. 2 all our generated ensembles for which computation of correlation functions is performed, i.e. we don't list the ones only used for matching. As we mentioned above, at β = 6.61 we can express the value of the lattice spacing in terms of the r 0 value. Using the parametrization of Ref. [33], we find r 0 /a = 12.76 for this β, which together with our chosen value of r 0 in physical units (0.48(2) fm), yields 0.0376(16) fm for the lattice spacing. This, together with our results for the matching of physical volumes, sets the scale for all our ensembles (see Tab. 2).

Estimating systematic uncertainties of the step scaling function
Lattice simulations are necessarily influenced by different kinds of systematic effects. To compare to continuum, infinite-volume, massless perturbation theory, we need to reliably extrapolate to the continuum, infinite-volume and chiral limits, respectively. L/a T /a κ nr of confs step 9.50 0.00235 (10) Table 2: Complete list of ensembles used for step scaling. We give the inverse bare coupling β, the lattice spacing in fm (with its uncertainty given by our choice of r 0 = 0.48(2) fm, the lattice size, the value of κ that yields maximal twist, the number of generated (saved) configurations and the number of heatbath updates between saved configurations.

Chiral limit
As already mentioned above, the chiral limit extrapolation is performed using three valence quark masses. This is illustrated with a detailed example in Sec. 4.1. In all cases, we observe behaviour compatible with a linear dependence on the valence quark mass. In several cases, actually, the dependence is mild enough such that a constant fit would also be appropriate. In total, quark mass effects are well under control and they are taken into account by propagating the error from the linear extrapolation to the massless limit.  : Continuum limit scaling of the massless PP correlator (X 2 ) 3 C P (X) in the free theory. We look at points (x/a, x/a, x/a, x/a) satisfying x/L = 1/4, 1/8, 1/12, 1/16, 1/24 and change L/a, up to 336, to approach the finite-volume continuum limit. The case of x/L = 1/4 is not shown as it yields a continuum limit more than 50% higher than the infinite-volume one (3/π 4 ). For x/L = 1/8, this deviation is around 0.5% and for x/L ≤ 1/12 it is below 0.025%.

Finite volume effects
Another source of systematic uncertainties are finite volume effects (FVE). To minimize them, we consider coordinate space points (x/a, x/a, x/a, x/a), (x/a, x/a, x/a, 0) and (x/a, x/a, 0, 0) (where 0 is at different positions) that always satisfy the condition x/L = 1/8 (for the temporal coordinate, we have t/T = 1/16 or t/T = 1/32). This is the smallest value of x/L for which FVE are expected to be small in the free theory, as shown in Fig. 4. In case of larger values of x/L, we observe that the continuum limit in the free theory is different from the one of the infinite-volume free theory. With x/L = 1/4, x/L = 1/8 and x/L = 1/12, the deviations of the continuum-extrapolated result from lattice data with respect to the free theory infinite-volume and continuum result are around 50%, 0.5% and 0.022%, respectively. Ideally, one should then choose x/L < 1/8 in the interacting case. However, the  computational cost of taking e.g. x/L = 1/12 is prohibitive, i.e. it would require generating lattices with L/a = 96 for our finest lattice spacing of each step scaling step or using only three coarser lattice spacings and going to L/a = 72. Hence, it is essential to confirm that FVE are relatively small for x/L = 1/8 that we choose for the interacting case. We perform this check by investigating the dependence of the step scaling function Σ Γ (µ, 2µ, a) at β = 7.90, where we have four volumes available, with L/a = 8, 16, 32, 64. We look at Σ Γ (µ, 2µ, a) calculated for points of type II, III and IV: • with x/a = 1 and x/a = 2 and pairs of lattices with L/a = 8, 16, L/a = 16, 32 and L/a = 32, 64; this implies that x/L changes from 1/8 to 1/16 and 1/32, respectively -see Fig. 5, • with x/a = 2 and x/a = 4 and pairs of lattices with L/a = 16, 32 and L/a = 32, 64; thus x/L changes from 1/8 to 1/16 -see Fig. 6  Looking at these two figures, we arrive at very important conclusions. We confirm that FVE are smaller than statistical errors even with x/L = 1/8, however only when the pair of lattice extents L/a = 8, 16 is excluded. For this pair, FVE can be statistically significant and up to 3% (red squares and circles in Fig. 5 for points of type IV). The pair L/a = 16, 32 is still fine from the point of view of FVE, even with x/L = 1/8. Hence, the ratio x/L is not the only factor that determines the size of FVE -in addition, the value of (X/a) 2 is also important and apparently FVE can be large if (X/a) 2 = 4, as happens on L/a = 8 lattices. The outcome of this test of FVE leads us to using the results from the pairs of lattices with L/a = 8, 16 in our continuum limit extrapolations only for checks of systematic effects from the fitting ansatz, while the central values are taken excluding such lattices. Following this strategy, we can obtain good estimates of the infinite volume step scaling functions and hence apply infinite volume perturbation theory.

Continuum limit
We now move on to discuss the continuum limit extrapolations, which are the next source of systematic effects. We consider renormalization constants extracted from two types of correlation functions, ones to which tree-level correction has or has not been applied, see Sec. 2. We use four types of continuum extrapolation fitting functions: i.e. a linear fit in a 2 , we take into account the three finest lattice spacings (for the fourth lattice spacing, we observe that deviations from the linear behaviour are usually significant) and only values of Σ Γ obtained with the tree-level correction (the non-corrected ones are not well described by a linear fitting ansatz due to larger discretization effects), in total three data points and two fitting parameters, Σ Γ (µ, 2µ) cont,A and c A , which is a quadratic fit in a 2 to all four lattice spacings, again with only tree-level corrected ratios of renormalization constants, four data points and three fitting coefficients, Σ Γ (µ, 2µ) cont,B , c B and d B , i.e. a combined fit linear in a 2 , using only three finest lattice spacings, in total six data points and three fitting parameters, Σ Γ (µ, 2µ) cont,C , c C1 and c C2 , which is a combined fit quadratic in a 2 , using all four lattice spacings, in total eight data points and five fitting parameters, Σ Γ (µ, 2µ) cont,D , c D1 , c D2 , d D1 and d D2 .

Non-ideal matching
We also consider systematic effects related to non-ideal matching of physical volumes. For given lattice volumes L and 2L and a given value of β, we obtain the estimates of the lattice step scaling function. Moreover, we know the matching uncertainty of β, denoted by ∆β and the values of the step scaling function obtained from ensembles at neighbouring value(s) of β, a smaller one, β − (for the first and second step scaling steps) and a larger one, β + (for the second and third step). This allows us to estimate the derivative of the lattice step scaling function with respect to β: ∂Σ Γ (µ, 2µ, a(β)) ∂β Our estimate of the matching uncertainty is then obtained as the product of

Detailed example
We now illustrate our procedure, including our systematic error estimates, with a detailed example. We choose the first step of step scaling and points of type IV, (x/a, x/a, x/a, x/a). We consider the scalar case, i.e. Γ = ½ ≡ S.
We aim to construct the step scaling function for the scale change from approx. 5.911 GeV to 11.822 GeV 4 , using the β values: (31) To construct the step scaling function at each lattice spacing, we first extrapolate the appropriate correlation functions to the chiral limit. This is illustrated in Fig. 7. In all cases, the behaviour is consistent with our fitting ansatz, linear in the valence quark mass, aµ v . We also note that the correlators at our lowest mass, aµ v = 0.0016, are always compatible with the chiral limit extracted value for this case. This confirms that quark mass effects are small.
The values of the correlation functions for appropriately chosen points are then used to construct the step scaling function Σ S (µ, 2µ, a) at each of the four lattice spacings. Our next step is to extrapolate to the continuum limit, using the four fitting ansätze defined in Eqs. (23)- (28). The extrapolations are shown in Fig. 8 and the continuum limit values, given in the key of the plot, agree well with one another within their errors. As our preferred value, we take the one for fit C, since it yields the smallest statistical error. The systematic uncertainty from the fitting ansatz is taken as the difference with respect to the continuum limit value for fit D, which includes one higher power of the lattice spacing and the data points from our coarsest lattice spacings.
Since we want to compare with continuum perturbation theory in the MS scheme, we now convert our continuum limit value from fit C to this scheme. The conversion factor is around 1.002 for these scales (for our third step of step scaling, involving scale change from around 1.5 GeV to 3 GeV, this factor can slightly exceed 1.01). Hence, our final estimate of the step scaling function in the MS scheme is: where the errors are, respectively, statistical, one from the fitting ansatz uncertainty, one from the matching and the propagated uncertainty of Λ   from the X-space scheme to the MS scheme and is hence very small in the first step of step scaling (large energy scales). The r 0 uncertainty enters via the conversion formulae and also via the absolute scale setting, thus influencing the estimate of the scale µ.
The obtained value can be compared to the prediction of perturbation theory (PT) [35,36]: Σ S (5.911 GeV, 11.822GeV) PT = 1.0713 (18), (33) where the uncertainty comes from the uncertainty of Λ MS . We thus obtain good agreement between our computation and the prediction of continuum perturbation theory (with N f = 0).

Full results for the step scaling function
In all other cases that we have computed, the analysis proceeds in a similar manner. We summarize the procedure: 1. Compute the relevant correlation functions in X-space, at three values of the valence quark mass.
2. Extrapolate to the chiral limit.
3. Apply the tree-level correction. 4. Use the chirally extrapolated and tree-level corrected values of correlators to compute the step scaling function. 5. Extrapolate to the continuum, using fit C or fit A (see below). 6. Convert from the X-space to the MS renormalization scheme. 7. Calculate systematic uncertainties from the fitting ansatz (by comparing fits C/D or A/B, see below) and non-ideal matching and the uncertainties of Λ MS and r 0 . Concerning points 5 and 7, we remark that combined fits (C and D) to treelevel corrected and non-corrected data points are possible only when cut-off effects are under control for the non-corrected case. We find that this happens for the PP/SS correlators for points of type IV and III. For these correlators with points of type II, as well as for VV/AA correlators (all types of points), we note that the χ 2 /d.o.f. of the fits including non-corrected points indicate bad fits, most probably due to higher-order discretization effects in the noncorrected correlators. Thus, in such cases, we only consider fits A and B and take the central value from the former, while the latter serve to estimate the systematic uncertainty.
Our results for the continuum-extrapolated step scaling function, determined from both the scalar and the pseudoscalar correlator, are shown in Tab. 3 and the comparison with 4-loop N f = 0 continuum perturbation theory [35,36] is given in Tab. 5. In most cases, we obtain good agreement between our lattice-extracted results extrapolated to the continuum limit and PT. This concerns in particular the results from the scalar correlator, which are typically (0 − 1.5)-σ away from PT. However, we observe some regularities depending on the type of points that we consider -points of type II(IV) tend to lie above(below) the PT result and points of type III tend to agree best with PT. This suggests that cut-off effects are very different for these kinds of points, which is indeed plausible, taking into account that very different cut-off effects are observed even in the free theory. A systematic treatment of these discretization effects (hypercubic artefacts) is highly desirable in order to obtain reliable results from the X-space scheme.
The results obtained from the pseudoscalar correlator are slightly worse in terms of agreement with PT, i.e. they tend to lie systematically below PT. In most cases, they are still within 1.5-σ, with the exception of the lowest scale µ (more than 2-σ too low). We note, however, that much better   agreement with PT is observed for points of type III and IV when not using the tree-level correction for the step scaling function (for points of type II this makes the agreement even worse). This strongly suggests that the cut-off effects for the PP case are much worse under control and, again, their better understanding is very important for the reliability of the X-space scheme.  We also computed the step scaling function Σ V /A (µ, 2µ) for the vector and axial vector case (Tabs. 4,5). Obviously, its continuum value is always 1.0, but its lattice computation provides a cross-check of the method and of the cut-off effects for different types of points in coordinate space. We observe that the expected value of 1.0 is well reproduced by our approach when using points of type III and IV. The deviations from unity are at most 1.3-σ and are of both signs, hence they are plausibly only statistical fluctuations. However, for points of type II, there are more significant deviations, as large as 3-σ. This suggests that cut-off effects in the vector and axial vector correlators are significant for the points of this type.

Comparison of the running with perturbation theory
As a summary of our results, we show two plots that illustrate the running of the renormalization constants, see Figs. 9 and 10.
The former shows the scalar/pseudoscalar case, i.e. the quark mass evolution with the energy scale. We show all our step scaling steps and all types of points. As already indicated, the running obtained from the scalar correlator is well reproduced, particularly for points of type IV and III (the left and the middle points in each group of three points corresponding to the same step). Points of type II (rightmost points in the group of three) give a result that is approximately 1-σ above the continuum curve. When the running is continuum pertubation theory lattice from X-space, PP correlator lattice from X-space, SS correlator To the left of these, there are three groups of three points, corresponding to the three steps of step scaling and to the three types of points. Each rightmost point from a given group corresponds to points of type II, the middle one to type III and the leftmost to type IV. extracted from the pseudoscalar correlator, we observe a possible tendency that the step scaling function is below its continuum value. This is clearly visible for all types of points, although the result is always within 1-σ of the continuum result even in the last step. It is desirable to investigate more the source of this observation, which can still only be a statistical fluctuation. We plan to address this issue in a forthcoming publication, where we will attempt to reduce hypercubic artefacts using strategies similar to ones applied already, with success, in momentum space, see e.g. Ref. [37].
In Fig. 10, we show the analogous "running" for the vector/axial vector case. Obviously, the corresponding renormalization constants are scaleindependent, hence the continuum value is exactly 1.0. We obtain good agreement with this result for points of type III and IV. However, for points continuum lattice from X-space, AA correlator lattice from X-space, VV correlator The three rightmost points are the starting points for our step scaling procedure and hence have no errors. To the left of these, there are three groups of three points, corresponding to the three steps of step scaling and to the three types of points. Each rightmost point from a given group corresponds to points of type II, the middle one to type III and the leftmost to type IV. of type II, deviations from unity are increasing when going to smaller energy scales, with 2-2.5-σ discrepancy in the final step scaling step. Interestingly, the behaviour is opposite for the vector and the axial vector case -the former systematically overestimates the ratio and the latter underestimates it. As for the (pseudo)scalar case, we plan to investigate this issue in more depth in our upcoming work. Using techniques from Ref. [37] could help to reduce discretizations effects and bring also points of type II closer to unity.

Discussion and outlook
In this work, we investigated, for the first time, the step scaling technique using the coordinate space renormalization scheme. This scheme has certain appealing properties with respect to e.g. the RI-MOM scheme, in particular it is gauge invariant and hence no gauge fixing is needed. A feasibility analysis was performed in the quenched approximation, to reduce the computational cost to a tractable level. We used the Creutz ratio to match lattices of different size such that they have the same physical volume and we performed three steps of step scaling to evaluate the running of the renormalization constants of the pseudoscalar and scalar densities, as well as of the vector and axial vector currents. In principle, the computed running of the former, i.e. of the quark mass, allows for an evaluation of the renormalization group invariant renormalization constants [38,39,40,41], which we have not attempted here, since our aim was a feasibility study of the method.
We considered three types of points in position space and showed that they lead to somewhat different levels of consistency between the latticeextracted and continuum-extrapolated results and continuum quenched perturbation theory. The overall consistency with perturbation theory is satisfactory. Only for the (axial) vector case with points of type II, we have encountered some difficulties, which motivates to explore other ways to reduce cut-off effects, which are different for different kinds of points in coordinate space.
We demonstrated that the X-space method can provide reliable results. It will be interesting to investigate whether the remaining sources of systematic uncertainties, especially the hypercubic artefacts, can be controlled with even better precision than the one obtained here. It would also be desirable to perform the analysis with points that are less likely to be influenced by residual finite volume effects, i.e. with a smaller ratio of the x/y/z coordinates to the spatial lattice size L, e.g. of 1/12 instead of 1/8 that we used for this study. This will, however, lead to significant increase of the computational cost, unless one uses only three lattice spacings, instead of four as in our current work. Finally, it would be of interest to perform this study in the dynamical case. Of course, the cost of such a step scaling procedure in the X-space scheme is substantial with dynamical fermions, but it provides a well-controlled computation of renormalization constants, allowing for a reliable non-perturbative renormalization.