Running coupling constant from position-space current-current correlation functions in three-flavor lattice QCD

In this Letter, we provide a determination of the coupling constant in three-flavor quantum chromodynamics (QCD), $\alpha^{\overline{\mathrm{MS}}}_s(\mu)$, for $\overline{\mathrm{MS}}$ renormalization scales $\mu \in (1,\,2)$ GeV. The computation uses gauge field configuration ensembles with $\mathcal{O}(a)$-improved Wilson-clover fermions generated by the Coordinated Lattice Simulations (CLS) consortium. Our approach is based on current-current correlation functions and has never been applied before in this context. We convert the results perturbatively to the QCD $\Lambda$-parameter and obtain $\Lambda_{\overline{\mathrm{MS}}}^{N_f=3} = 342 \pm 16$ MeV, which agrees with the world average published by the Particle Data Group and has competing precision. The latter was made possible by a unique combination of state-of-the-art CLS ensembles with very fine lattice spacings, further reduction of discretization effects from a dedicated numerical stochastic perturbation theory simulation, combining data from vector and axial-vector channels and matching to high-order perturbation theory.


Motivation:
The strength of strong interactions, parametrized by the scale-dependent coupling α s , typically quoted at the Z-boson pole mass, is one of the most important parameters of the Standard Model (SM). It is required in perturbative calculations in collider physics and its uncertainty is one of dominant sources of uncertainty in several SM predictions, as well as in tests of SM extensions [1]. Due to the non-Abelian Yang-Mills nature of quantum chromodynamics (QCD), α s vanishes asymptotically at very high energies [2,3] and experiments are able to follow this energy dependence in various processes over a wide range of energy scales. This allows to determine the value of the strong coupling at several scales by fitting experimental data and matching to a perturbative expansion of an appropriate observable. Equivalently, using renormalization group concepts, one may parametrize the running of α s by a single parameter, Λ, corresponding to the scale where perturbation theory breaks down. Examples of experimental processes for the extraction of the strong coupling or the Λ-parameter are hadronic τ decays, deep inelastic scattering and hadronic final states of e + e − annihilation. For a review of many aspects of such determinations and the obtained values, see the Particle Data Group (PDG) review [4]. However, the strong coupling constant or the Λ-parameter can also be extracted directly from the QCD Lagrangian, using the non-perturbative formulation of QCD on the lattice. This proceeds by calculating appropriately designed short-distance Euclidean observables and, again, matching them to their perturbative expansions. Over the years, several methods how to design such observables have been proposed. Recent investigations employed e.g. step scaling methods [5,6], the static quark-antiquark potential [7][8][9], the vacuum polarization function [10], the heavy-quark current two-point correlation function [11], QCD vertices (e.g. ghost-gluon) [12] or eigenvalues of the lattice Dirac operator [13]. For a discussion and overview of these and older results, see the Flavor Lattice Averaging Group (FLAG) review [14]. The determinations from experiments and from the lattice enter the world average of α s in the PDG review [4], recently with a visibly larger impact of lattice results due to their smaller total uncertainties.
In this Letter, we describe a novel method of estimating the running of the coupling or the Λ-parameter, using numerical simulations of QCD. The proposed method employs large volume simulations, it has a moderate numerical cost and is clean and straightforward from the theoretical point of view. It is based on current-current correlation functions in position space, objects well studied and easily accessible in the lattice QCD framework. Thanks to the combination of very fine lattices generated by the Coordinated Lattice Simulations (CLS) effort [15,16], precise renormalization factors from the chirally rotated Schrödinger functional (χSF) framework [17] and O(a)-improvement coefficients [18,19], subtraction of leading-order and next-to-leading-order discretization effects estimated in the numerical stochastic perturbation theory (NSPT) formulation [20] and various improved analysis techniques, it yields a competitive total uncertainty. Therefore, it may serve as a robust method of estimating the Λ-parameter. The approach presented here is not limited to Λ -a position-space analysis may also be used to reliably estimate other important observables, such as quark and gluon condensates [21], quark masses [22] or operator renormalization functions [23][24][25][26]. Strategy: The strategy proposed in this Letter uses a combination of numerical lattice QCD calculations and high-order perturbative results. We concentrate on correlation functions of flavor non-singlet bilinear quark op-erators of the form where x is the physical distance, m q is the quark mass of degenerate three flavors of quarks, a denotes the lattice spacing, Γ = {V ≡ γ µ , A ≡ γ µ γ 5 }, i = j and Z Γ is the (scale-independent) renormalization factor. For a reliable extraction of α s , we need to work in the regime of distances satisfying a window condition, a x Λ −1 . The former condition guarantees that discretization effects are not enhanced, while the latter establishes that reliable contact to perturbation theory can be made. After extrapolating the correlation functions to the continuum limit and after renormalization, they can be matched to their perturbative expansions in terms of α s , typically in the MS scheme (α MS (µ)), where C Γ (µ) ≡ C Γ (1/x, m q = 0, a = 0). Such an expansion of current-current correlators is presently available up to 4 loops [27]. Knowing C Γ (µ) from numerical simulations and the analytic form of coefficients c (i) Γ , we solve Eq. (2) for α MS (µ). Subsequently, we convert that value to our estimate of the Λ-parameter. Now, we provide details of the different steps needed to reliably obtain C Γ (µ). Crucial elements of the analysis: We start with the bare lattice data for correlation functions C A/V (1/x, m q , a) with the O(a)-improvement of the currents implemented by using improvement coefficients c A from Ref. [18] and c V from Ref. [19]. The ensembles used in this study are summarized in Tab. I. We perform 64 inexact and 2 exact measurements per configuration using the truncated solver method [28] and for every lattice distance x/a, we average correlators evaluated from all sites equivalent with respect to the hypercubic symmetry of the lattice.
At fixed lattice spacing and lattice distance, we extrapolate the correlators to the chiral limit. We use a fitting ansatz linear in the dimensionless combination y = t 0 m 2 π , where t 0 is an intermediate unphysical scale introduced in Ref. [30] and we take the values of t 0 /a 2 from Ref. [29]. The quality of the chiral fit was tested at β = 3.55, where we have four pion masses available. We compared the linear fit in y with the quadratic one for all the relevant distances and we observed that coefficient of the y 2 term was always statistically insignificant. Thus, for other βvalues, we employed the linear ansatz. We denote the massless correlator by C Γ (1/x, a).
The massless correlators are then expressed in the MS scheme, using renormalization factors calculated in Ref. [17], determined using the χSF framework [31].
A significant step to reliably perform the continuum limit extrapolation is to reduce the size of discretization effects present in the data. To this aim, we perturbatively  [15,16]. The gauge action is tree-level Symanzik-improved, while the fermionic one is the Wilson O(a)-improved (clover) action with the improvement coefficient, cSW , determined non-perturbatively. rqcd30, X450, B450, X250, X251 have been generated by the RQCD and Mainz collaborations. For more details, see Refs. [15,16]. The values of t0/a 2 are the reweighted estimates using the symmetric definition of the Yang-Mills action density [29]. The lattice spacings corresponding to different β values are 0.07582(24) fm (β = 3.46), 0.0644 (7) fm (β = 3.55), 0.0499(5) fm (β = 3.7) and 0.0391(15) fm (β = 3.85) [29]. The last column indicates the number of configurations used.
compute O(a ∞ g 2 ) artifacts, i.e. we replace the correlation functions Γ , t ∈ {cont, lat} and the superscript TL/1L denotes the tree-level/1-loop contributions. The massless MS continuum correlators, C cont A/V , are given in Ref. [27]. In turn, C lat A/V are computed in NSPT [20] along the lines of Refs. [32,33] 1 and are expressed in the MS scheme using the renormalization factors for the employed gauge action [35]. Thus, all terms appearing on the RHS of Eq. (3) are correctly normalized correlators in the same scheme. The improved correlator, hence, has leading cutoff effects of O(a 2 g 4 ). We demonstrate the reduction of discretization effects in Fig. 1, depicting the distance dependence of C A (1/x, a) at β = 3.85. We show three data sets: without any correction, with the tree-level correction only and with the full 1-loop NSPT correction. The scatter of data points is clearly reduced, yielding a rather smooth curve, with the remaining outliers being the lattice points that break the rotational symmetry most severely. It is important to emphasize that without the 1-loop subtraction of artifacts (all orders in the lattice spacing), the spread of points at different distances prohibits any meaningful extraction of α s and thus, this step is crucial for the success of the method.
In order to perform the continuum extrapolations, we need to follow the lines of constant physics. In our case, the only relevant scale is the correlator distance x, which we keep fixed in physical units by interpolating to the desired distance at all β values. We use two interpolation ansatzes, linear and quadratic in x 2 , between the two and three closest data points to find the interpolated value at each lattice spacing. We consider three lattice directions, for which hypercubic artifacts are known to be the smallest [24,25]: (0, k, k, k), (k, k, k, k), (0, k, k, 2k) with k ∈ {1, 2, 3} and interpolate independently for each of them. In this way, we keep the discretization effects related to the breaking of rotational symmetry fixed as we change the lattice spacing. We use the difference of the two interpolation models as the systematic uncertainty associated to this step.
If discretization effects are under control, the continuum limits corresponding to the same physical distance should agree for each of the three lattice directions. We checked that this is indeed the case and hence, we performed combined continuum fits of data for all three directions. Depending on the distance (and, thus, the available lattice spacings), we use from 6 to 12 data points and constrain the fit by a single common value in the continuum, C Γ (1/x, a = 0). The fitting ansatz reads (4) and has 4 fit parameters.
The flavor non-singlet axial and vector currents are related by a global SU (2) A transformation. In the chiral limit, this symmetry is spontaneously broken. The difference between the two correlation functions was estimated in various frameworks, for a review see Ref. [36], including lattice QCD [21,37]. Also empirical data exist for this observable [38]. At short distances, the difference between the vector and axial correlators is reliably provided by the operator product expansion [39]. Using estimates from Ref. [40], the relative difference ranges from 0.03% at x = 0.1 fm up to 1.5% at x = 0.3 fm. Hence, within the statistical and systematic precision of our data, the two correlators are indistinguishable in that range of distances, see the inset of Fig. 2. We use this observation in a two-fold way. a) First, we employ it as a test of the reliability of the continuum extrapolation. In further analysis, we consider only the physical distances for which the independently extrapolated axial and vector correlators agree within their uncertainties. On the one hand, this criterion excludes lattice directions and physical distances which are too short and no control over discretization effects is possible, setting the lower limit to 0.1 fm. On the other hand, the two correlators are no longer equivalent at distances larger than around 0.35 fm within our precision, which sets the upper limit on the physical distances where the impact of non-perturbative condensates is negligible. In Fig. 2 know both sides of Eq. (2) and we can determine α MS (µ) for different scales, corresponding to different physical distances 1/µ = x. The results are shown in Fig. 3. At distances above around 0.2 fm (scales below 1 GeV), we observe that the running of the coupling freezes, indicating the breakdown of perturbation theory. We convert our results for the coupling to Λ N f =3 MS [41,42] separately at each distance, see Fig. 4. We show the perturbative running of α s using our final value of the Λ-parameter in Fig. 3 and we discuss it below, after addressing systematic effects in our determination. Final result: We consider several sources of uncertainty in our analysis and we decompose the error of our final result for the Λ-parameter according to these different sources. The raw lattice correlators are, obviously, subject to statistical errors ("lat stat"). The perturbative subtraction of discretization effects via NSPT is also subject to statistical errors ("NSPT stat") and moreover, to a systematic uncertainty of extrapolation of NSPT results to the infinite volume limit ("NSPT infV"). The latter is computed as the difference between a polynomial fit to several volumes ranging from V = 32 4 up to V = 80 4 and the estimates from the largest volume V = 80 4 . The correlator interpolation uncertainty, described above, is denoted by "interpol". Renormalizing the correlators in the MS scheme introduces an uncertainty from the values of Z-factors ("Z A " and "Z V "). The uncertainty of the c V and c A improvement coefficients is completely negligible compared to its other sources. Finally, we estimate the truncation uncertainty of the final Λ-value as the difference between conversions of α s results to Λ using the 4-loop and 5-loop β-functions ("trunc"). These differences are shown in Fig. 4, including also the 2-loop and 3-loop cases. The observed behavior suggests that while 3-loop perturbation theory still shows significant trunca- obtained from αs at different distances. Shown are conversions between αs and Λ [41,42] using from 2-to 5-loop perturbative β-functions [43][44][45]. Our final value (indicated by the green band whose width is the total uncertainty) is obtained by a systematic procedure explained in the text. The inset shows the histogram of results corresponding to different ranges of distances taken into account in the systematic procedure, along with a Gaussian fit.
tion effects in the considered energy range, the 4-loop and 5-loop results evince convergence. We double this uncertainty to cover the truncation of the perturbative series of Eq. (2), where the 5-loop coefficient is not available at present.
To make our final result independent from the choice of the window of physical distances where α s is extracted, we adopt a systematic procedure similar to the one used in Ref. [25]. From all distances smaller than 0.2 fm, above which the coupling freezes, we choose the range 0.13-0.19 fm, where all other systematic uncertainties are under good control. Having 7 determinations of Λ corresponding to these different distances, we calculate all possible weighted averages covering from one to seven subsequent distances. We use the 28 resulting values of Λ to build a weighted histogram, where the weights are taken as the squared inverse error of each individual result. The histogram is approximately Gaussian (see the inset of Fig. 4) and we fit its mean and width to determine the central value, i.e. Λ N f =3 MS , and its uncertainty from the choice of the physical distances ("window"). This central value, along with the total uncertainty, is shown as the green band in Fig. 4.
The final result for the Λ-parameter reads: where we combine the individual uncertainties in quadrature to obtain the total error. Our final value agrees well with earlier lattice determinations, e.g. with the recent one of Ref. [5], Λ N f =3 MS = 341(12) MeV, and with a comparable total error, dominated in our case by the uncertainty from the choice of the physical distance and by the uncertainty from the NSPT correction. Discussion and Conclusions: In this Letter, we presented and tested a novel method to estimate the MS strong coupling constant using numerical simulations of coordinatespace correlators and used it to determine the 3-flavor QCD Λ-parameter. It is based on theoretically clean current-current correlation functions in position space at small distances. Our results suggest that the challenging multiscale problem of evaluating α s at scales which at the same time satisfy the requirement that µ a −1 in order to control discretization effects and µ Λ to match results to an MS perturbative expansion can be addressed using lattices available today. We have shown that using a combination of state-of-the-art simulations and novel analysis techniques, one can find a window of available scales µ and provide an estimate of Λ N f =3 MS with a competitive precision. In particular, the crucial steps are the perturbative subtraction of hypercubic artifacts and a combined continuum extrapolation using four lattice spacings and several lattice directions, which allowed us to control discretization effects at small distances in lattice units. We, furthermore, profited from independent evaluations of axial and vector correlators, which have a common continuum limit at short distances, to design a criterion to characterize the quality of continuum extrapolations and gain confidence in the results.
To conclude, we believe that the techniques described in this Letter provide a robust way of extracting the running of the QCD coupling and the QCD Λ-parameter, with good statistical precision and well-controlled sources of systematic effects. Furthermore, the precision reached in this work can be increased even more in a systematic way. In addition, techniques based on current-current correlators in position space, improved by NSPT reduction of discretization effects, can prove to be useful in determining other quantities, such as the quark and gluon condensates.