Charm quark effects on the strong coupling extracted from the static force

We compute the fermionic contribution to the strong coupling αqq extracted from the static force in Lattice QCD up to order g4 in perturbation theory. This allows us to subtract the leading fermionic lattice artifacts from recent determinations of αqq produced in simulations of two dynamical charm quarks. Moreover, by using a suitable parametrization of the βqq-function, we can evaluate the charm loop effects on αqq in the continuum limit.


Introduction
Many simulations of QCD are carried out with 2 + 1 dynamical quarks, taking into account only the effects of light sea quarks (up, down, strange). So far, this kind of approach has provided really good results and its simulation costs are affordable with modern computer facilities.
However, the discovery of new charm-states in experiments like Belle, CLEO and BABAR has made charm physics really appealing in the last few years and including a dynamical charm quark in Lattice QCD simulations could help us understand better the properties of these new states.
Although nowadays new supercomputers allow to simulate the 2 + 1 + 1 flavor theory, we know [1,2] that the effects of a dynamical charm quark are small on low energy quantities, thus a high statistical precision is needed to disentangle them. Moreover, to resolve the small correlation length associated with a charm quark, really fine lattices are required to control the extrapolation to zero lattice spacing. Therefore it is interesting to understand for which kind of observables it is more worthwhile including a charm quark in simulations of QCD.
In this work, we evaluate the charm loop effects on the strong coupling α qq , whose definition is where C F = 4/3, r is the distance between a static quark-antiquark pair and V (r) ≡ dV dr is the derivative of the static potential V(r) with respect to the distance. The static force F(r) = V (r) can also be used to measure a hadronic scale r 0 defined through r 2 F(r)| r=r 0 = 1.65 [3]. Our first results about these studies can be found in [4], where we saw that charm loop effects on α qq become significant at about 2 GeV. Here we extend our previous results, considering an additional dynamical ensemble, subtracting the leading fermionic lattice artifacts from our non-perturbative data and taking the continuum limit of α qq through a convenient parametrization of the β qq -function.

Numerical setup
To evaluate the dynamical charm effects on α qq , we compare QCD (N f = 2), with two heavy degenerate quarks having the charm mass M c , to quenched QCD, namely QCD (N f = 0).
As a lattice discretization, we consider Wilson's plaquette gauge action [5] and a clover improved doublet of twisted mass Wilson fermions [6,7]. At maximal twist, the inclusion of a clover term is not needed for O(a) improvement of physical observables, but it reduces the O(a 2 ) lattice artifacts [8]. Moreover, we choose open boundary conditions in the time direction and periodic boundary conditions in spatial directions.
To generate the dynamical ensembles, the bare coupling g is chosen such that the lattice spacings cover the range 0.23 fm < a < 0.36 fm, whereas the hopping parameter κ is set to its critical value to achieve maximal twist. Finally, the twisted mass parameter µ is chosen such that the RGI mass in our simulations corresponds to the charm mass M c . The pure gauge theory is simulated at similar or smaller lattice spacings, setting the scale through the hadronic quantity r 0 as described in [9].
In particular, we study 1 three quenched ensembles at β = 6/g 2 = 6.34, 6.672, 6.90 and three dynamical ensembles at β = 5.70, 5.88, 6.00 and M/Λ = 4.87, where M/Λ is the ratio of the RGI mass to the Λ parameter. All the simulation parameters of these ensembles are listed in Table 1 of Ref. [11]. For further details we refer to [4,11].

Computation of the static force
In this section we summarize the main steps that lead us to the computation of the strong coupling α qq in QCD (N f = 2) at M = M c and QCD (N f = 0).

Non-perturbative calculation
From the definition of α qq in the continuum, given in Eq. (1), it is clear that a lattice regularization of the derivative of the static potential V(r) is needed. The most natural choice would be where the static potential V(Ra) can be extracted from the expectation value of a rectangular Wilson loop W(R, T ) in the limit of infinite time separation, T → ∞. However, it has been shown [3,9] that it is better to introduce an improved distance r I such that the static force has no cutoff effects at the tree-level in perturbation theory, namely A table of the improved distances for unsmeared links is provided in [9].
Here, in order to reduce the typical gauge noise which affects the measurements of V(Ra) at large distances we follow Ref. [12]. In particular, before we measure the Wilson loops, all the gauge links are replaced by HYP2-smeared ones [13], which correspond to the following choice of the smearing HYP-parameters: α 1 = 1.0, α 2 = 1.0, α 3 = 0.5. The static potential aV(Ra) is extracted with great accuracy by smearing the initial and final lines of rectangular Wilson Loops up to four levels of HYP-smearing and then solving a generalized eigenvalue problem [14]. Since the improved distances r I depend on the static quark action, we make use of the values of r I for HYP2-smeared links listed in Table 2 of [12]. Wilson loops have been computed using B. Leder's program available at https://github.com/bjoern-leder/wloop/.

Perturbative calculation
One of the goals of this work is to subtract the leading fermionic lattice artifacts from the nonperturbative measurements of α qq realized on our dynamical ensembles.
For this purpose, first we compute the fermionic contribution to the static force up to order g 4 in perturbation theory. Then we extract the leading lattice artifacts comparing our calculations to the one-loop predictions of the continuum theory [15]. Such a computation has been performed in [16,17] for clover improved Wilson fermions using unsmeared links. Here we adopt the same strategy, but we have to consider a clover improved doublet of twisted mass Wilson fermions and HYP2-smeared links. We just summarize the main ideas and we refer to [16,17] for further details.
The perturbative expansion of a Wilson loop W(R, T ) can be written as W 2 (R, T ) involves only gluons, whilst W 4 (R, T ) can be splitted into two terms where W g 4 (R, T ) comes from the pure-gauge theory and W f 4 (R, T ) is a purely fermionic contribution. Once W 2 (R, T ) and W 4 (R, T ) have been calculated, it is possible to access the perturbative expansion of the static potential and the static force up to order g 4 . Finally, converting to the MS scheme we can rewrite F(r I ) as with where N f is the number of flavors and m i are the quark masses. Since the corresponding continuum expressions f 1,g (0) and f 1, f (z, 0) are known [15,17], it follows that the relative lattice artifacts can be written as where for unsmeared (empty markers) and HYP2smeared links (full markers) choosing the same lattice discretization as our dynamical ensembles. In the right panel we subtract the leading lattice artifacts from our non-perturbative data of α qq (full markers) and we compare to the unsubtracted ones (empty markers). α PT is a perturbative 4-loop calculation of α qq in QCD (N f = 0), using the Λ parameter from [19].
In this work we only focus on the fermionic term 4πδ (1, f ) F of Eq. (9), as the calculation of the gluonic term becomes much more intricate when using HYP-smeared links and it would need some extra care 2 . Moreover, as we will see in Section 4, our continuum extrapolation of α qq is already accurate enough and allows to clearly distinguish the charm-loop effects on α qq at high energies. This part of the calculation has been carried out using our computer package written in Mathematica. The fermionic contribution to the static force has been calculated for a sequence of finite lattice sizes (L = [16,64]) and then extrapolated to infinite volume.

One-loop cutoff effects
In this section we summarize our numerical results concerning the extraction of the leading fermionic lattice artifacts. The calculation has been realized for unsmeared and HYP2-smeared links, in order to see what are the main differences between these two cases.
On the left hand side of Figure 1 we show the calculation of 4πδ (1, f ) F using the lattice parameters of our dynamical ensembles. We see clearly that at large distances in lattice units unsmeared and HYP2smeared links produce similar effects, while at small distances using HYP2-smeared links gives rise to larger lattice artifacts. Moreover, as we expect, these unwanted effects become smaller when β → ∞. To obtain the relative size of the lattice artifacts, 4πδ (1, f ) F must be multiplied by g 2 MS (see Eq. (8)) and this means, at our lattice spacings, that we observe around 5% effects.
The second step is to subtract the leading lattice artifacts from our non-perturbative data of α qq . The result of this procedure is depicted in Figure 1 (r.h.s), where also a comparison with the unsubtracted data is shown. The most evident effect of this subtraction is that data corresponding to similar physical distances are closer to each other compared to the unsubtracted ones, as we can see clearly from the measurements of α qq at around 0.38 r/r 0 . This is encouraging because it is a signal that the size of the lattice artifacts has been reduced considerably.
The last part of this work is aimed at extracting the continuum limit of α qq in QCD (N f = 2) at M = M c and QCD (N f = 0), in order to evaluate how large the charm loop effects are on this observable.

Strategy
Studying the step scaling function σ [20] can provide a powerful tool to reach this purpose. If f is a fixed scale factor, σ( f, u) measures how much the coupling changes when the distance scale changes by a factor f Then, from the definition of the β qq -function one arrives at the exact relation Eq. (12) is only true in the continuum, but following [21] we parametrize β qq and the cutoff effects in such a way that the continuum limit of α qq can be easily extracted from our lattice simulations. Let us introduce the following parametrization of β qq , which is not motivated by perturbation theory, but it allows to parametrize our data really well: , P(g qq ) = p 0 + p 1 g 2 qq + p 2 g 4 qq + · · · .
This choice permits to rewrite Eq. (12) as where 2 × (n max + 1) is the degree of the polynomial P(g qq ).
However, on a lattice one can only measure an approximation Σ( f, u, a/r 0 ) of the step scaling function σ( f, u) such that lim This means that if we want to extract the coefficients p 0 , p 1 , . . . , p n max +1 of our parametrization of β qq from lattice simulations, instead of Eq. (14) we have to use where h ≡ h( f, u, a/r 0 ) is a particular function that depends, other than f and u, on the lattice spacing a. Since we expect cutoff effects proportional to a 2 /r 2 0 , we choose to parametrize h( f, u, a/r 0 ) as Thus, we can estimate the coefficients p i of β qq performing a global best-fit to our data at different lattice spacings through the Eqs. (16) and (17). The continuum extrapolation of α qq can be realized solving the ODE which follows from the definition of the β qq -function (Eq. (11)). Since this requires the choice of an initial condition (or the value of the Λ parameter for both theories), we need to know the value of α qq at a reference distance r re f in the continuum limit. This can be achieved by using an interpolation function for the static force [9] F(r re f ) = f 1 + f 2 r −2 re f (18) between the two neighboring points. Setting r re f = 0.75r 0 and taking the continuum limit of the interpolations realized at different lattice spacings we obtain

Results of the best-fits
Before showing our final results, we begin this section with a few remarks. Using improved distances r I , it is not possible to keep the factor f exactly constant. However, we have seen that choosing f ∈ [2.03, 2.10] and f ∈ [1.94, 2.10] for our quenched and dynamical ensembles respectively produces acceptable best-fits for a large number of parametrizations of β qq and ρ(u), see Eqs. (13), (17). The widths of these ranges depend in some manner on the lattice spacings and on the explored distances, therefore some attempts were needed before arriving at the ranges above-mentioned.
We tried different types of (correlated) best-fits, varying the number of parameters both in the parametrization of the β qq -function and in the function ρ(u) that parametrizes the cutoff effects. We only show the parametrizations that provide an acceptable chi-squared using the minimum number of parameters, underlining that we obtain compatible results for different best-fits.  Table 1: Results of the continuum extrapolation in N f = 0 and N f = 2 theories. In the second row we show the results obtained with the original non-perturbative data, while in the third row the ones obtained subtracting the leading fermionic lattice artifacts. Table 1 summarizes the results of our continuum extrapolations. The numbers listed in the table have been produced using 6 parameters for β qq and 1 for ρ(u) in quenched QCD (thus 6 + 1 parameters on the whole), whilst 5 + 1 parameters have been used in QCD (N f = 2). For equal number of parameters, subtracting the leading fermionic lattice artifacts for the 2 flavor theory produces a bit smaller relative errors and a better chi-squared. Moreover, we can see that the coefficient ρ 0 , introduced to parametrize the cutoff effects, is compatible with zero when the one-loop fermionic lattice artifacts are subtracted. Expanding (13) in powers of g qq , we can rewrite β qq in a way similar to the one motivated The blue circles and the red squares denote the measurements of α qq produced with our finest lattices in N f = 0 and N f = 2 theories respectively. The blue and red bands stand for the continuum extrapolations obtained using three different lattice spacings for both theories. The widths of the bands originate from the errors on the data and their correlation is taken into account. The right panel shows instead a comparison between our continuum extrapolation of α qq in quenched QCD (the blue band) and the one obtained in a previous work [9] (the red circles) using a different strategy. The dashed black lines are the predictions in perturbation theory up to four loops. The spread of the lines comes from the uncertainty in the Λ parameter. by perturbation theory. In particular, the first two coefficients are given by This allows us to identify b 0 ≡ 1/p 0 and b 1 ≡ −p 1 /p 2 0 . From our continuum extrapolation of the quenched theory we obtain b 0 = 0.062(3) and b 1 = 0.013(3). These estimates deviate a bit from the perturbative results for N f = 0 (b 0 = 11 (4π) 2 ≈ 0.070 and b 1 = 102 (4π) 4 ≈ 0.004), but this is something that we could expect because the parameters p i have been estimated in a range which is far away from the domain of validity of perturbation theory (α qq < 0.20 in case of a 4-loop calculation). Once we know the parameters of the β qq -function, the continuum limit of α qq can be easily extracted and our final results are depicted in Figure 2. On the left hand side a comparison of α qq in N f = 0 and N f = 2 (subtracting the leading fermionic lattice artifacts) theories is shown, where we use r/r 0 (M c ) on the x-axis for the dynamical points. We see that the continuum limits are accurate enough to distinguish the dynamical charm effects on α qq at distances r/r 0 0.5. In the right panel we compare our continuum extrapolation of α qq in quenched QCD to the one obtained in Ref. [9] and the predictions of perturbation theory up to four loops. We observe a really good agreement with [9] and with perturbation theory at α qq 0.20.

Conclusions and Outlook
In this work we have focused on evaluating the one-loop cutoff effects of a dynamical charm quark on α qq using HYP2-smeared links. We see that these effects are small, but HYP2-smeared links produce bigger lattice artifacts compared to unsmeared links at small distances in lattice units.
We have also tried to extract the continuum limit of α qq studying the step scaling function σ. This strategy allowed us to evaluate the dynamical charm effects on α qq and for quenched QCD we find a really good agreement with perturbation theory at high energies and Ref. [9] at low energies. Our continuum extrapolations also indicate that the dynamical charm effects on α qq are significant at distances r/r 0 0.5.
In future we plan to extend these measurements to other values of quark masses to study the mass-dependence of the strong coupling α qq .