Power corrections from decoupling of the charm quark

Decoupling of heavy quarks at low energies can be described by means of an effective theory as shown by S. Weinberg in Ref. [1]. We study the decoupling of the charm quark by lattice simulations. We simulate a model, QCD with two degenerate charm quarks. In this case the leading order term in the effective theory is a pure gauge theory. The higher order terms are proportional to inverse powers of the charm quark mass $M$ starting at $M^{-2}$. Ratios of hadronic scales are equal to their value in the pure gauge theory up to power corrections. We show, by precise measurements of ratios of scales defined from the Wilson flow, that these corrections are very small and that they can be described by a term proportional to $M^{-2}$ down to masses in the region of the charm quark mass.


Introduction
In a field theory which contains light (mass-less) fields and fields of a heavy mass M , the functional integral over the latter can be performed resulting in an effective theory for the light fields which was formulated by Weinberg [1]. The action of the effective theory contains the action of the light fields (without the heavy fields) and an infinite number of non-renormalizable terms. The latter are suppressed by powers of E/M at low energies E M . Moreover, the non-renormalizable couplings do not contribute to the renormalization group equations of the renormalizable couplings of the light fields. This property holds for mass-independent renormalization schemes like the MS scheme as shown in [1]. The heavy fields still affect the value of the renormalized couplings of the light fields through the decoupling relations, which result from the matching of the effective and the fundamental theory at low energies.
Assuming the validity of perturbation theory at the matching scale, the decoupling relations can be computed perturbatively. In the case of QCD and one heavy quark, such as the charm or the bottom quark, the decoupling relation for the renormalized strong coupling is known to four loops [1,2,3,4]. The strong coupling of the five-flavor theory can be extracted in this way from the coupling computed non-perturbatively in the three-flavor theory using lattice simulations [5]. We remark that the decoupling relation for the strong coupling can be equivalently expressed as a relation between the Λ parameters of the effective and the fundamental theory [6].
Simulations of QCD on the lattice are often carried out with three light sea quarks [7,8,9,10,11,12]. The inclusion of a charm sea quark increases significantly the computational cost and introduces additional tuning to set the bare quark masses on a line of constant physics. Moreover, in the case of simulations with Wilson fermions, Symanzik O(a) improvement requires the computation of coefficients which multiplies terms proportional to the bare quark masses in lattice units am [13,14]. The contribution of these terms is significant for the charm quark am c > 0.3 at the affordable lattice spacings a > 0.05fm. Some of these coefficients, like the one of the gluon action, are difficult to extract non-perturbatively. Relying on decoupling of the charm quark at low energies allows to simulate the cheaper and simpler effective theory with three flavors only.
The applicability of decoupling for the charm quark has to be justified. In [6] this was studied in a model, QCD with two heavy mass-degenerate quarks and no light quarks. The decoupling of the heavy quarks leave a pure gauge theory 1 up to power corrections (which are due to the non-renormalizable interactions) at low energies. The latter were extracted by computing low energy quantities related to the Wilson flow [16,17,18,19]. Ratios of two such quantities are insensitive to the matching of the gauge couplings, and after taking the continuum limit, can be compared to their counterparts in the pure gauge theory. The differences are due to the power corrections. By interpolating data obtained from simulations at quark masses ranging from 1/8th up to one half of the charm quark mass with data from simulations of the pure gauge theory, the size of the power corrections due to one sea charm quark was estimated to be at the sub-percent level [6].
In [6] it was noted that the simulated masses were not large enough to see the leading behavior of the power corrections which start at 1/M 2 in the effective theory.
Instead a behavior more like 1/M was observed. In this article we study the same model as in [6] but extend the simulated quark masses to the charm quark mass and slightly above. Thus we can directly compute the size of the power corrections from decoupling of the charm quark. Furthermore we perform a non-perturbative test of the validity of the effective theory of decoupling for the charm quark. Our goal is to determine whether the leading power corrections in the inverse heavy quark mass behave in the charm region as 1/M 2 .
The article is organized as follows. In Sect. 2 we briefly review the theoretical framework of the effective theory of decoupling for QCD with two heavy mass-

Decoupling
To avoid a multi-scale problem, we consider a simplified version of QCD, namely an SU ( L YM is the Lagrangian of the SU (3) Yang-Mills (pure gauge) theory. Due to gauge invariance there are no fields of mass dimension equal to five. A complete set of fields of where F µν is the SU (3) field strength tensor and D µ F νρ its covariant derivative.
At leading order the effective theory, eq. (1), is a Yang-Mills theory. It has only one free parameter, the renormalized gauge coupling. This coupling is fixed by matching the effective theory to the fundamental theory. Equivalently one can fix the Λ parameter of the Yang-Mills theory, Λ YM , which becomes a function Λ YM = Λ dec (M, Λ), see [6,21]. Matching requires that low energy physical observables are the same in the two theories up to power corrections. Let us denote a low energy observable by m had where, for example, it can represent a hadronic scale such as 1/ √ t 0 [22] or 1/r 0 [23].
After matching where m had (M ) is the hadronic scale in QCD with N f = 2 heavy quarks of mass M and m had YM is the hadronic scale in the Yang-Mills theory. Note that m had YM depends on M through the matching, in particular m had YM /Λ YM is a pure number independent of M . Therefore we consider two hadronic scales, m had,1 (M ) and m had,2 (M ), whose values in the Yang-Mills theory are m had,1 YM and m had,2 YM respectively. An immediate consequence of eq. (2) is The matching of the coupling is irrelevant for the ratios and we have direct access to the power corrections [24]. The effective theory of decoupling predicts that the ratios like in eq. (3) are equal to their value R(M = ∞) in the Yang-Mills theory with a leading power correction in the inverse heavy quark mass given by where k is a parameter which depends on the hadronic scales which are taken to form the ratio. The goal of this work is to verify the behavior in eq. (4) and to establish whether it applies for masses around the charm quark mass.

Monte Carlo simulations
We simulate QCD with two mass-degenerate flavors of quarks (N f = 2). Wilson's plaquette gauge action [25] is employed in the the Yang-Mills sector and a doublet of quarks is realized either as standard or as twisted mass [26] Wilson quarks. In both cases a clover term [27,13] with non-perturbatively determined improvement coefficient c sw [28] is added. It is not needed for the O(a) improvement of the twisted mass action at maximal twist, but was found to reduce the O(a 2 ) lattice artifacts, see e.g. [29].
The bare coupling β of the gauge action was chosen such that the lattice spacings cover the range 0.023 fm a 0.066 fm. The lattice spacing is determined from the hadronic scale L 1 [30,31]. The scale L 1 /a is defined at vanishing quark mass, where the standard and twisted mass Wilson quark formulations are equivalent. Therefore, the lattice spacing for a given bare coupling β is the same for both formulations. In order to obtain the scale L 1 in lattice units at a particular value of β, we fitted the data in  [6,21]. For the present work we also simulated twisted mass Wilson quarks at these quark masses. In addition, we also simulated directly at the charm quark mass and at one mass larger than that of the charm quark. In the simulations of twisted mass Wilson quarks, the hopping parameter κ was set to its critical value in order to achieve maximal twist. The critical values were obtained from an interpolation of published data [31,30,35]. The twisted mass parameter aµ was chosen to correspond to certain values of M/Λ listed in Table 1. More precisely, at a given value of the bare coupling the twisted mass parameter is set by The pseudo-scalar renormalization constant at renormalization scale L −1 1 computed in the Schrödinger Functional scheme Z P (L 1 ) = 0.5184(53) (valid for 5.2 ≤ β ≤ 6.0), the relation between the running and the RGI mass M/m(L 1 ) = 1.308 (16) and Λ = 310 (20) MeV are known from [31,36]. We take ΛL 1 = 0.649(45) from [37]. Some of the quantities entering eq. (5) have rather large errors. The dominant error comes from ΛL 1 . Note that it is common to all our simulation points and amounts to a change in the target values of M/Λ. For the charm quark we set M c /Λ = 4.8700, where we use the preliminary value M c = 1510 MeV of [38] which agrees with [39].
In order to determine the power corrections in eq. (3), we also simulate the pure gauge theory at values of the scales t 0 /a 2 [22], r 0 /a [23] which are similar or larger.

Results
On the ensembles generated with two flavors of O(a) improved Wilson fermions reported in [6,21] and those generated in the twisted mass simulations at maximal twist and the pure gauge simulations which are listed in Table 1 At flow times t > 0, E(x, t) is a renormalized quantity [22,40]. The reference scale t 0 is defined as in Ref. [22] by Similarly, the scale t c is defined by The numerical solutions to eq. (7) and eq. (8) are found by quadratic interpolation of the data of E(t). We use the clover (symmetric) definition of the action density E on the lattice, cf. [22]. The scale w 0 is defined as in Ref. [41] by where E (t) = d dt E(t). The numerical solution to eq. (9) is found by first computing the symmetric finite differences of t 2 E(t) on each configuration and then by quadratic interpolation of the data. The error analysis is based on the method of Ref. [42] and takes into account the coupling to slow modes following Ref. [43]. We estimate the exponential autocorrelation time τ exp from the tail of the autocorrelation function of t 0 . We use the ensembles at our largest masses, for which we perform simulations at our finest lattice spacings. As expected with open boundary conditions [32], the observed critical slowing down is compatible with a τ exp ∝ a −2 behavior and from a least squares fit we obtain τ exp = −32(23) + 17.4(2.8) t 0 /a 2 in molecular dynamics units (MDUs). The errors of the fit coefficients are given in parantheses. With periodic boundary conditions the autocorrelation times would be much larger, cf. [43] where an effective scaling proportional to a −5 was observed.
We compute the ratios of hadronic scales In such ratios the bare coupling (or equivalently the lattice spacing) drops out. After taking the continuum limit we can directly compare the ratios in the N f = 2 theory to their value in the pure gauge theory and so determine the size of the 1/M 2 effects in eq. (4).
To determine the values of the ratios in the continuum limit, we fit our data to R(a, M/Λ, A) = R cont (M/Λ) + a 2 t0 c(M/Λ, A), where A is the action, "W" for Wilson, "tm" for twisted mass and "q" for quenched (pure gauge, M = ∞). The functional form is motivated by Symanzik's effective theory for our actions. For a given mass M/Λ we have two fit parameters (continuum value and slope) for cases where the calculation was performed with one action and three parameters (continuum value and two slopes) for cases with two lattice actions. We apply a cut, a 2 /t 0 (M ) < 0.32, to the data to be fitted. The data and the fits for R = t c /t 0 and R = √ t 0 /w 0 are shown in Fig. 1. The pentagrams represent the twisted mass data, the squares are the  The continuum values R cont (M/Λ) obtained from the fit are summarized in Table 2 Figure 2 and on the right plot in Figure 3. In Fig. 2 we also show the fits to the data down to masses corresponding to M c /4.
They are represented by the dashed lines, blue for the quadratic and green for the linear fit. While the linear fit is almost unchanged with respect to the linear fit down to M c /2, the quadratic fit is clearly excluded. Still the quadratic fit to the points down to M c /2 has much better χ 2 values per degree of freedom than the linear fit down to M c /4, see Table 3 . These results demonstrate that the linear behavior is disfavored only when the data for masses smaller than M c /2 are excluded from the fits. The data then begin to show the leading quadratic behavior which is expected from the effective theory for large enough masses. This conclusion is further supported by the dashed-dotted lines shown in Fig. 2. They represent a fit down to M c /4 when a next-to-leading correction term Λ 4 /M 4 (whose coefficient is fitted) is added to the leading behavior of eq. (4).
This fit has very good χ 2 values per degree of freedom: 0.49/2 for R = t c /t 0 and  and subsequently measuring the Wilson loops where the initial and final line of gauge links are smeared using up to four levels of spatial HYP smearing [45]. This allows us to extract the static-quark potential a V (r) very reliably by solving a generalized eigenvalue problem [46]. The static force F (r) = V (r) can then be used to measure the hadronic scale r 0 defined implicitly through [23] r 2 0 F (r 0 ) = 1.65 .
Even a careful state of the art determination of r 0 does not yield a precision high enough to resolve the power corrections we are interested in, as can be seen in Fig. 3. We note that the coefficient k of the (Λ/M ) 2 term in eq. (4) depends on the observables used to form the ratio.

Conclusions
In this work we simulated a model, QCD with two heavy mass-degenerate quarks.
At low energies this theory is described by an effective theory which is a pure gauge theory up to power corrections in the inverse heavy quark mass. By comparing ratios of low energy physical quantities computed in both theories and extrapolated to the continuum limit, we could determine the size of the power corrections. They have been found to be very small. We are now confident that the effects of neglecting the charm quark in N f = 2 + 1 simulations is far below a percent in dimensionless low-energy quantities. The power corrections are expected for sufficiently large heavy quark masses to be proportional to the square of the inverse quark mass, see eq. (4).
Our data shown in Fig. 2 follow very well this expectation down to masses equal to half of the charm quark mass. A behavior linear in the inverse quark mass, which is possible for masses outside the range of validity of the effective theory, is strongly disfavored by the data.
In order to achieve a stronger conclusion and find the range of quark masses where the linear behavior can be completely excluded, the statistics of our simulations should be increased and larger quark masses up to the bottom quark mass be simulated. The resources needed to carry this out are beyond our computational budget. We emphasize that the computational resources used to produce the data for this article amount to a large scale project already. Simulating with a yet larger statistical precision and heavier quark masses would require a computational effort which is comparable to simulations with light sea quarks and are therefore beyond the scope of this model calculation.