Controlling quark mass determinations non-perturbatively in three-flavour QCD

The determination of quark masses from lattice QCD simulations requires a non-perturbative renormalization procedure and subsequent scale evolution to high energies, where a conversion to the commonly used MS-bar scheme can be safely established. We present our results for the non-perturbative running of renormalized quark masses in Nf=3 QCD between the electroweak and a hadronic energy scale, where lattice simulations are at our disposal. Recent theoretical advances in combination with well-established techniques allows to follow the scale evolution to very high statistical accuracy, and full control of systematic effects.

Abstract. The determination of quark masses from lattice QCD simulations requires a non-perturbative renormalization procedure and subsequent scale evolution to high energies, where a conversion to the commonly used MS scheme can be safely established. We present our results for the non-perturbative running of renormalized quark masses in N f = 3 QCD between the electroweak and a hadronic energy scale, where lattice simulations are at our disposal. Recent theoretical advances in combination with wellestablished techniques allows to follow the scale evolution to very high statistical accuracy, and full control of systematic effects.

Quantum Chromodynamics at all scales
For numerical simulations one starts from one of many valid discretisations of the (bare) Euclidean action density of Quantum Chromodynamics (QCD), with a fixed number of dynamical quark flavours N f . With N f + 1 input parameters the action is supposed to describe the strong interaction at all energy scales. Lattice QCD is a gauge-invariant regularisation and thus does not require a gauge-fixing and Faddeev-Popov term. It is the only regularisation allowing to determine colorless bound states (m π 0 , m K 0 , m p , . . .) between quarks from first principle calculations, which due to computational and theoretical advances can be systematically improved over time. However, incorporating all physical mass-scales that are relevant for describing (N f = 6) QCD-as it seems to be realised in nature-is an impossible task to present state-of-the-art simulations due to the immense costs to accommodate many orders of magnitude of energy scales in a e-mail: p.fritzsch@csic.es a single simulation, cf. Figure 1. Being mainly interested in non-perturbative, long-distance effects of QCD, the phenomenological applicability of LQCD is limited to a fixed window 1 40 MeV ∼ Λ IR µ had , m π 0 , m p , m K 0 , . . .
Any physically interesting scale µ has to be well separated from the imposed infrared and ultraviolet cutoffs, Λ IR µ Λ UV , to allow for a controlled assessment of the respective cutoff effects, see caption of Figure 1. Although, the lattice community can afford to simulate with three or four dynamical flavours nowadays, one still needs to balance the N f + 1 chosen (bare) input parameters, given in terms of N f + 1 physical scales from the real spectrum, against the intrinsic cutoff values and algorithmic costs. As no practical perfect lattice setup exists, different views on how to best achieve this balance lead to controversies among lattice practitioners. Another of those disputable topics is the application of perturbation theory at energies as low as those given by the low-energy window of lattice QCD as specified in (2), including µ = Λ UV ∼ 3 GeV, see also [1].
In course of our determination of quark masses, we have to fix a renormalization prescription at scale µ had that lies well inside the window [Λ IR , Λ UV ]. Then the natural question arises how to compare or connect our result to other determinations as summarised for instance by the Particle Data Group [2]. Both naturally differ in the renormalization scheme and scale. After removing the UV cutoff dependence, we thus have to evolve our (continuum) result, m i (µ had ), using renormalization group (RG) transformations to the same physical scale and scheme, typically MS.
In a mass-independent renormalization scheme for QCD, the renormalization group equations (RGE) for the running coupling and running quark mass(es) read and can be formally integrated. The resulting, exact solutions are known as renormalization group invariants (RGI), conveniently written as Only the leading, universal (scheme-independent) coefficients b 0 , b 1 , d 0 of the perturbative expansions of the beta-function and mass-anomalous dimension, appear such that the individual integrands are finite when g → 0. Higher order coefficients are schemedependent-indicated by the superscript s-and so are β and τ. Thanks to recently completed efforts these approximations are now known to 5-loop order in the minimal subtraction scheme(s), cf. [6][7][8][9][10][11][12] and [13][14][15][16][17]. However, it is well-known that by using an asymptotic series expansions one is missing non-perturbative contributions 2 which become increasingly more relevant towards low energies, i.e., when the strong coupling becomes large. One has to appreciate that these issues can be overcome by determining β(g) and τ(g) non-perturbatively. In that case, eqs.  Figure 1. Sketch of energy scales relevant to QCD phenomenology. The masses for the six quark flavours found in nature so far (u, d, s, c, b, t) are taken as quoted in [2]. Also the masses of the neutral pseudoscalar meson bound states (m π 0 , m K 0 , m D 0 , m B 0 ) as well as the depicted scale dependence of the strong coupling α s are taken from that reference. Although, α s (µ) can be computed in any sensible renormalization scheme one conveniently quotes it for the 5-flavour theory at the electroweak scale µ = m Z in the MS scheme. The lattice regularisation of QCD provides the only known (non-perturbative and gauge-invariant) framework to compute bound states of the strong interaction. For phenomenological applications one is restricted to a window (Λ IR ≤ µ ≤ Λ UV ) in the low-energy regime of QCD in order to incorporate the important long-distance effects. They are of the order of the Compton wavelength of the lightest particle, λ π 0 = 1/m π 0 . It has to be well separated from the infrared cutoff Λ −1 IR ∼ L λ π 0 in order to not be distorted by finite volume effects. Typical simulations have a physical extent of L 3 fm and a number of lattice points of N = (L/a) 4 ≈ (50 − 100) 4 , thus restricting the lattice spacing to a ∼ Λ −1 UV 0.045 fm. To extract physical quantities one also has to stay away from the ultraviolet cutoff, i.e. µ Λ UV , otherwise control over finite lattice spacing effects is lost, and the continuum limit a → 0 cannot be taken. Although, various advances in LQCD steadily increase the quality and reliability of numerical simulations with different number of flavours N f , the cost of simulations and algorithmic difficulties still prevent us from reaching even smaller lattice spacings. For any non-perturbative renormalization problem posed at fixed renormalization scale µ in the low-energy regime of LQCD the continuum limit has to be taken in a controlled way. Then the remaining question is how to relate the obtained results to experiments at much higher energies, typically two orders of magnitude. For high precision physics the supposedly most convenient renormalization scheme, MS, is of little help at scales much below 10 GeV (or α s 0.2). Its intrinsic perturbative nature, the behaviour of asymptotic series expansions in general, and missing non-perturbative contributions prevent us from quantifying its (non-)reliability at low energies from within the scheme itself. As detailed in the main text the renormalization group running of any operator can be determined purely non-perturbative, thus connecting the low-energy regime of LQCD in a controlled and quantifiable way to scales well above 10 GeV. Recently, it has become possible to non-perturbatively test the accuracy of continuum perturbation theory for α s 0.2 [1], stressing the relevance and necessity of a careful assessment beyond perturbation theory in present high precision physics. Non-perturbatively it is even possible nowadays to follow the renormalization group in QCD down to scales of about 200 MeV, cf. Ref. [3]. Both results have been reported at this conference [4,5]. . Non-perturbative β-function for the SF (α SF 0.2) and GF (α GF 0.17) running coupling scheme; picture from [3]. There also the two schemes have been matched non-perturbatively at fixed physical volume corresponding to about 4 GeV.

An intermediate non-perturbative renormalization scheme
For massless renormalization schemes (RS) it is most natural to first solve the RG equation for the coupling and then for the mass, cf. (3). In LQCD it has become customary to determine the scale evolution of coupling and mass using a finite-volume renormalization scheme in combination with recursive finite-size scaling. For that purpose the physical size L of the simulated volume L 4 is identified with the inverse renormalization scale for the sole purpose of solving the RG equations by rescaling the box size, L → sL. The constant factor s is usually chosen to be 2. In that way one is able to cleanly disentangle the renormalization procedure from the low-energy window of LQCD and to connect the hadronic renormalization scale µ had = 1/L had to the electroweak scale µ = m Z , where for instance a change of renormalization schemes towards the commonly quoted MS scheme can typically be achieved in a controlled way, 3 cf. Figure 1. For that purpose one should ideally cover a range of scales compatible with For example, to connect a scale µ had as low as 200 MeV to the electroweak scale of about 100 GeV using s = 2, at least N s = 9 individual steps are necessary. Each of the N s steps requires a series of lattice simulations at different resolutions L/a but fixed renormalized input parameters, g 2 (L) = u and Lm i (L) = 0. In this way µ = 1/L is kept fixed while only the lattice spacing a is being varied, or equivalently, the bare parameters {g 2 0 , m 0,i } as a function of a. A second set of simulations at the same bare parameters but sL/a number of points in each direction permits to trace the change of renormalization scale L → sL at finite a. With adequate renormalization conditions for the coupling and quark mass, this allows to take the continuum limit a → 0 of their lattice approximants in order to determine the response w.r.t. a change of scale s in the continuum, equivalent to In this way, the relevant information is encoded in the step-scaling functions (SSFs) [18,19], While various ways exist to determine the SSFs numerically via lattice simulations of QCD, the socalled Schrödinger functional (SF) setup [20][21][22] has been develop for exactly that purpose. For instance, in contrast to standard (anti-)periodic boundary conditions in all space time directions, Dirichlet boundary conditions are imposed in time direction. This allows for a natural non-perturbative definition of a strong coupling via non-vanishing boundary gauge fields as well as to simulate at vanishing quark mass. For additional details we point to the review article [23] and references therein. The scale evolution of the strong coupling has been determined in refs. [1,3], and a short account of the whole procedure can be found in [24]. However, a technical complication is inherited from that determination: instead of a single RS, two different schemes are being used and matched nonperturbatively at a scale of approx. 4 GeV. At present, its physical motivation is two-fold: 1) control the perturbative matching at the electroweak scale to high precision [1] using the Schrödinger functional coupling, g SF , and 2) reach good statistical accuracy down to energies of about 200 MeV [3] by applying the Gradient flow coupling scheme, g GF [25,26]. Compared to previous estimates of α s (m Z ) this represents a new quality of rigor from lattice QCD determinations and for the first time allowed to accurately determine the corresponding non-perturbative QCD β-functions, see Figure 2.
For our determination of the scale evolution of renormalized quark masses in three-flavour QCD, this does not pose any additional obstacles but only complicates the overall presentation. In fact, assuming a diagonal quark mass matrix of rank N f in two massless schemes, one can write [27] µ = cµ , c > 0, Invariance of a physical observable P under this change of variables implies where P satisfies the Callan-Symanzik equation in the primed scheme w.r.t.
While this is the general connection between two massless schemes at the level of RG functions, every step in our determination is done non-perturbatively. In practise, when switching between the two schemes, we do not change the imposed renormalization condition-which is relevant for determining σ P (2, u)-but merely the implicit, parametric dependence on the renormalized coupling u ≡ g 2 (L). The connection between SF and GF schemes has been established non-perturbatively at a well-chosen fixed physical box size, L 0 ≈ (4 GeV) −1 , and reads [3] g 2 SF (L 0 ) = 2.012 , g 2 GF (2L 0 ) = 2.6723(64) .  [28,29]. Right: Non-perturbatively determined step-scaling functions σ P (u) in the SF (u 2.0) and GF (u 2.1) running coupling scheme, together with different interpolating fits. We note that for u SF 1.5 the data agrees with NLO perturbation theory such that the matching with SF-PT can be safely established.

The running quark mass
Given the discussions in the previous sections, our strategy of a non-perturbatively controlled determination of quark masses in the three-flavour theory with two renormalization schemes proceeds as follows. We determine the flavour-independent RG factor M/m(µ had ), cf. eq. (5), that connects the RGI mass to the quark mass in the hadronic, low-energy regime Note that we intend to make use of perturbation theory at a scale µ PT close to m Z where truncation errors from using the known 3-/2-loop orders in the perturbative β-/τ-functions in the SF scheme are sufficiently suppressed, 4 see Figure 3. The true challenge is to determine the two non-perturbative factors to a sufficient precision while controlling systematic effects at the sub-percent level for each individual contribution. When this is achieved, the scale-and scheme-independent RGI mass of flavour j can be determined from where m j (µ had ) constitutes a renormalized quark mass computed in the hadronic regime. Ideally, M j is the fundamental parameter to be compared to other determinations. However, it has become customary to quote MS masses. In order to do so one subsequently has to employ perturbation theory and, depending on the individual flavour and renormalization scale, appropriately match at the charm and bottom thresholds. An important fact worth emphasising is that the continuum running factor in eq. (14) can be used by other lattice practitioners which may employ a different lattice discretization in their computations. They just have to recompute m j (µ had ) using the same-but so far  unmentioned-renormalization prescription, such that eq. (15) remains well-defined. For the employed Wilson fermion action, the quark mass m j can be obtained, on the lattice up to a finite renormalization, from the non-singlet axial Ward identity, or current quark mass The bare mass m awi j ≡ m awi j (g 2 0 , {am 0,i }) still depends on the bare parameters of the action, cf. eq. (1). From the previous discussions the connection to the mass-SSF and notation should become evident cf. eqs. (8) and (9). Z P is the multiplicative renormalization factor of the non-singlet pseudoscalar current, which itself is determined from a SF renormalization condition as in [19]. In Figure 3 (left panel) we show preliminary results for some continuum extrapolations as in eq. (17) with s = 2. The corresponding continuum SSF data points are then plotted against the initially fixed target couplings u for both the SF-and GF-coupling schemes (right panel). Note that these values of u can be chosen at will but should cover the range of interest well enough in order to reach the desired precision. Typical fit ansaetze as shown in Fig. 3 are polynomials like σ P (u) = 1 + p LO u + p 2 u 2 + · · · or Padé approximants with the correct, leading (universal) asymptotics fixed. With such a given smooth interpolation formula one is then able to recursively construct the running factor and its uncertainty in steps of s = 2 as summarised in (7), i.e., the two ratios in (14) become Accordingly, we can identify the switching scale L swi appearing in eq. (14) with the scale 2L 0 of (13), which we inherit together with the coupling-SSF σ(u) from the complementary determination of the running coupling [1,[3][4][5]24]. This concludes the computation of (14) and we return to the determination of eq. (15), which we can now rewrite as The disappearance of the hadronic scale µ had reflects the scale-independence of M j , which itself is connected to the current quark mass m awi j only by a scale-independent renormalization Z M . Hence, by determining Z M the renormalization problem is fully solved non-perturbatively using the massless (intermediate) SF scheme. While the costs of simulating the Schrödinger functional setup are good to control in general, towards larger physical volumes the simulations get more and more involved and thus set a natural limit on the values the hadronic scale can take. At the moment we are exploring different ways to fix the value of aµ had in order to determine the Z P (g 2 0 , aµ had ), and thus Z M , with high precision. As the reader may have noticed, forcing s ab in eq. (18) to be a multiple of 2 is a quite stringent condition, especially as we already had to fix the scheme-switching scale in (14). Thanks to the high statistical accuracy and the achieved control of a few systematic effects which are propagated into the final error, we are able to lift this restriction for the first time. By reassessing eqs. (8) and (9), one notices that with the data points at hand and a given parameterisation of the β-function, as shown in Figure 2 and taken from [1,3], only a well motivated ansatz for τ(g) is required for a sensible least square minimisation. We present a preliminary but already very encouraging result in Figure 4, again for both the SF-and GF-running coupling schemes. It should be noted that this result is actually independent of the originally chosen scale-change factor s = 2. As a non-trivial cross-check we can use the resultant to numerically reconstruct σ P (u) for arbitrary u from eq. (8). This is shown in Figure 3 and coincides very well with the original estimate. At the time of this conference we were still accumulating statistics at the two strongest renormalized couplings in use. Now we are finalising our analysis by also employing global fit ansaetze in various ways and include some yet missing correlations between our data, which is going to be published soon. Accordingly, we here refrain from quoting any quantitative numbers.
Compared to the aforementioned simulations, the m awi j estimates in (19) are obtained from largescale lattice simulations that define the accessible low-energy window in the first place. We will employ N f = 2 + 1 ensembles jointly produced by the Coordinated Lattice Simulations (CLS) effort [30], with which we share the lattice discretization in use. Although there have been many advances in lattice QCD, the accessible parameters am 0,i , in general, do not correspond directly to the set of parameters at which m π and m K take their physical values for different values of the cutoffs (even after correcting for isospin splitting and QED effects). Beside setting the physical scale [31], one thus still has to apply chiral perturbation theory. In that respect, eq. (19) is a minor simplification in our discussion, but nevertheless a step that has to be carefully evaluated in the near future.

Conclusions
We have presented a full strategy to determine renormalized quark masses from first principle lattice calculations in three-flavour QCD. Compared to other approaches in the community, we separate and thus disentangle the problem of renormalization from large-scale lattice QCD simulations. In this way the massless renormalization group running for the coupling and quark masses could be mapped out non-perturbatively to high accuracy in the employed Schrödinger functional scheme, and a connection to other schemes can be easily established via renormalization group invariant quantities {Λ, M i } N f . This strategy fully circumvents the unanswerable question regarding the applicability of perturbation theory at parametrically large values of the strong coupling. Furthermore, it provides a quantifiable and systematically improvable uncertainty at each stage of the calculation, something perturbation theory cannot really provide.
For the first time two different non-perturbative schemes, based on the Schrödinger functional coupling and the younger Gradient flow coupling, have been used and matched non-perturbatively. This specific combination has been chosen for the avail of improving on statistical and systematic uncertainties. Also the non-perturbative determination of the continuum quark mass anomalous dimension is a unique achievement for QCD.
The prospects for determining fundamental parameters of QCD valuable for high precision physics are excellent. To obtain an even more realistic approximation of nature this study can be extended to cross the charm flavour threshold (N f = 3 → N f = 4) non-perturbatively at some point in the future. We finally remark that also the renormalization and scale dependence of the tensor current-relevant in rare meson decays (especially B decays), and studies of the neutron electric dipole moment-is being computed along the same lines [32].