Coulomb Artifacts and Bottomonium Hyperfine Splitting in Lattice NRQCD

We study the role of the lattice artifacts associated with the Coulomb binding effects in the analysis of the heavy quarkonium within lattice NRQCD. We find that a"naive"perturbative matching generates spurious linear Coulomb artifacts, which result in a large systematic error in the lattice predictions for the heavy quarkonium spectrum. This effect is responsible, in particular, for the discrepancy between the recent determinations of the bottomonium hyperfine splitting in the radiatively improved lattice NRQCD [1, 2]. We show that the correct matching procedure which provides full control over discretization errors is based on the asymptotic expansion of the lattice theory about the continuum limit, which gives $M_{\Upsilon(1S)}-M_{\eta_b(1S)}=52.9\pm 5.5~{\rm MeV}$ [1].

The lattice simulations within the effective theory of nonrelativistic QCD (NRQCD) [3,4] has developed into one of the most powerful tools for the theoretical analysis of heavy quarkonium properties [5]. This method is entirely based on first principles, allows for simultaneous treatment of dynamical heavy and light quarks and gives a systematic account of the long distance nonperturbative effects of the strong interaction. The perturbative matching of lattice NRQCD to the full theory of relativistic continuum QCD is thought to be well understood. One of the most interesting applications of the method is the analysis of the bottomonium hyperfine splitting. The latter quantity, defined by the mass difference E hfs = M Υ(1S) − M η b (1S) , has been a subject of much controversy since the first observation of the spinsinglet η b state in radiative decays of the Υ(3S) mesons by the BaBar collaboration [6]. The measured value of the hyperfine splitting overshot the predictions of perturbative QCD [7] by almost a factor of two, well beyond the experimental and theoretical uncertainty bands. Such a discrepancy would indicate a serious failure of perturbative QCD in the description of the bottomonium ground state, in clear conflict with the general concept of the heavy quarkonium dynamics. Further experimental studies [8][9][10] were consistent with the initial measurement, while the Belle collaboration reported a significantly lower value of the splitting with higher experimental precision [11]. The advance of lattice NRQCD is expected to provide an accurate model-independent prediction and solve the problem on the theory side. Surprisingly, the two most recent independent calculations of the hyperfine splitting [1,2] which fully incorporate the one-loop radiative corrections give significantly different values of the splitting, see Table I. The analysis [1] reconciles the theoretical predictions of lattice and continuum QCD, as well as the most accurate experimental data [11]. On the other hand the result of Ref. [2] is close to the PDG average [12] but is not fully consistent with the perturbative QCD estimate [7]. Both calculations are based on the same lattice data and the discrepancy between the results exceeds what one would expect for the perturbative approximations which are formally of the same order in the strong coupling constant α s . At the same time Refs. [1,2] rely on different methods of perturbative matching and the inconsistency of the results indicates that a careful study of the general procedure of the radiative improvement of lattice NRQCD is necessary.
In this paper we study a subtle problem of the lattice NRQCD analysis of the heavy quarkonium spectrum related to the lattice artifacts associated with the Coulomb binding effects. We show that a widely used direct numerical matching procedure [13,14] generates spurious linear Coulomb artifacts and, in particular, leads to a large systematic error in the lattice prediction for the hyperfine splitting [2,5]. The problem is related to the all-order character of the Coulomb binding effects and is naturally solved when the perturbative matching of lattice NRQCD is performed through the asymptotic expansion about the continuum limit [1]. We show that after removing the spurious contribution the result of Ref. [2] is in a good agreement with [1].
The paper is organized as follows. In the next section we outline the general framework and describe different approaches to the fixed order perturbative matching. In Sect. II the structure of the Coulomb lattice artifacts is studied in detail. The result is applied to the analysis of the hyperfine splitting in Sect. III. Sect. IV is our summary and conclusion.

I. RADIATIVE IMPROVEMENT AND MATCHING IN LATTICE NRQCD
Within the NRQCD approach the hard modes, which require a fully relativistic analysis, are separated from the nonrelativistic soft modes. The dynamics of the soft modes is governed by the effective nonrelativistic action given by a series in heavy quark velocity v, while the contribution of the hard modes is encoded in the corresponding Wilson coefficients. The nonrelativistic action can be applied in a systematic perturbative analysis of the heavy quarkonium spectrum [15][16][17]. At the same time the action may be used for lattice simulations of the heavy quarkonium states, which gives full control over nonperturbative long-distance effects [18,19]. In the latter approach the inverse lattice spacing a plays a role of the effective theory cutoff separating the hard scale m q and the soft scale vm q , where m q is the heavy quark mass.
As an example, let us consider the spin-dependent part of the NRQCD Lagrangian, which is responsible for the hyperfine splitting to O(v 4 ). It reads (see e.g. [20,21]) where B is the chromomagnetic field, C F = (N 2 c − 1)/(2N c ) is the SU (N c ) color group factor, ψ (χ c ) are the nonrelativistic Pauli spinors of quark (antiquark) field, and we have projected the four-quark interaction on the color-singlet state. The coefficients c F = 1 + O(α s ) and d σ = O(α s ) parameterize the quark anomalous chromomagnetic moment and the effective local four-quark interaction, respectively. In the given order of the NRQCD expansion in 1/m q they depend logarithmically on the effective theory cutoff 1/a. This dependence can be predicted to all orders of perturbation theory by renormalization group methods (see e.g. [22,23]). The radiative improvement of the action is therefore mandatory for the correct continuum limit.
The effect discussed in this paper is characteristic for the quark-antiquark interaction and we focus on the Wilson coefficient d σ of the four-quark operator. It vanishes in the Born approximation and is determined by matching the one-particle irreducible quark-antiquark scattering amplitudes in QCD and NRQCD. The matching becomes particulary simple when the amplitude is computed at the quark-antiquark threshold and vanishing momentum transfer. In this case the one-loop full QCD amplitude is where C A = N c , T F = 1/2, and we introduced a small auxiliary gluon mass λ to regulate the infrared divergence. The power enhanced 1/λ term corresponds to the Coulomb singularity of the threshold amplitude, while the term proportional to T F is due to the two-gluon annihilation of the quark-antiquark pair.
A. Expansion about the continuum limit This approach has been developed in [1] and relies on the formal asymptotic expansion of the lattice loop integrals about the continuum limit [24] to obtain the NRQCD amplitude as a series in a order by order in the heavy quark mass expansion. To the leading order in 1/m q and a it gives (cf. Eqs. (2,3)) where L = ln(am q ). For the simplest lattice action with no improvement for gluonic and heavy quark fields the method provides the analytical result [1] where the irrational constants b 2 = 0.02401318 . . ., b 3 = 0.00158857 . . . parameterize the lattice tadpole integrals and can be computed with arbitrary precision. For the HPQCD action [5], which is used in real simulations, the nonlogarithmic coefficient has been computed numerically [1]: Note that Eq. (3) has only a logarithmic singularity in a in the formal continuum limit a → 0. In higher orders of the NRQCD expansion in 1/m q the asymptotic expansion includes more singular terms with a negative power of a. Such 1/(am q ) n terms are suppressed with respect to Eq. (3) in the region 1/a ≪ m q , where lattice NRQCD is applied.

B. Direct numerical matching
This approach has been originally used for the radiative improvement of lattice NRQCD. Within this prescription for a given action in a given order in α s the NRQCD amplitude is computed numerically without the expansion in 1/m q and a. The Wilson coefficient is then determined by the difference between the QCD and NRQCD amplitudes in the limit λ → 0. Since no expansion is performed, it has a nontrivial dependence on a dimensionless variable am q and can be written as follows where the logarithmic and annihilation contribution are separated and given in an analytic form. The function ∆(am q ) can formally be expanded in an asymptotic series where the lower summation limit is negative and depends on the approximation used for the NRQCD action. To determine the function ∆(am q ) we use the numerical data of the most recent analysis [2] based on the O(v 6 ) action. 1 In Ref. [2] the numerical values of the Wilson coefficient are given for three different values of the lattice spacing corresponding to am q = 1.95, 2.73, 3.31, where the actual lattice simulations are performed. Numerical simulations [2,14] show that in general the terms with negative n become important for significantly lower values of the lattice spacing corresponding to am q ∼ 1 and can be neglected in the region under consideration. Indeed, the numerical data are well approximated by a linear function with the coefficients The zero-order term of the expansion can be related to the value of the Wilson coefficient obtained through the expansion about the continuum limit, Eq. (4), as follows in a reasonable agreement with Eq. (9). A characteristic feature of the result of the numerical matching is the 1 In Refs. [2,5,13,14] a different basis of the four-quark operators is used and the Wilson coefficient dσ/αs should be identified with the linear combination 9 FIG. 1. One-loop Feynman diagrams with Coulomb singularity contributing to the spin-dependent one-particle irreducible part of the scattering amplitude in QCD (a) and NRQCD (b). The symmetric NRQCD diagram is not shown. In the diagram (b) the double arrow, dashed and wavy lines stand for the nonrelativistic quark, Coulomb and transverse gluon propagators, respectively. The black circles denote the effective spin chromomagnetic interaction proportional to the Wilson coefficient cF in Eq. (1).
linear dependence of the Wilson coefficient on a, which is unusual for the lattice simulations with the improved action. It is related to the Coulomb binding effects in heavy quarkonium discussed in the next section.

II. COULOMB BINDING EFFECTS ON THE LATTICE
In perturbation theory the Coulomb binding effects shows up through the singular (α s /v) n terms in the contribution of the n-loop planar ladder diagrams. Since in an approximately Coulomb bound state v ∼ α s , such terms have to be resummed to all orders. In the perturbative approach [16] this is done by constructing the perturbative expansion about the Coulomb nonrelativistic solution rather than the free quark and antiquark. At the same time the characteristic momentum scale of the Coulomb dynamics is vm q ≪ 1/a and the Coulomb effects are included in the lattice NRQCD simulations along with the nonperturbative effects of strong interactions at the scale Λ QCD . The Coulomb contribution is ultraviolet finite and therefore its effect on the matching coefficients is suppressed by a power of a, i.e. is a lattice artifact. Below we consider the role of such Coulomb artifacts in the calculation of the coefficient d σ .

A. One-loop Coulomb artifacts
The Coulomb singularity is contained in the planar box diagrams of QCD (Fig. 1a) and NRQCD (Fig. 1b), and takes the form α s m q /λ since the matching calculation is performed with v = 0. Let us consider the evaluation of the corresponding contribution to the NRQCD amplitude to O(a). The expansion of the lattice NRQCD Feynman rules in a generates the second or higher order terms so we can use the continuum expressions for the gluon and nonrelativistic heavy quark propagators where k = (k 0 , k). After integrating over the time component of the virtual momentum by taking the residue of the heavy quark propagator, the Coulomb contribution to the scattering amplitudes takes the form where the integration over the spatial virtual momentum is restricted to the first Brillouin zone. Without loss of generality we consider a spherically symmetric lattice with the Brillouin zone defined by |k| < π/a, and after integrating over the angular components obtain The contribution of the first singular term of Eq. .
The correction term in the denominator of Eq. (14) results in an additional contribution to the integral in Eq.
where we neglected the gluon mass since the integral is infrared finite. Thus the O(v 4 ) correction to the nonrelativistic kinetic energy increases the coefficient of the linear term by factor two, which gives For comparison with the direct numerical matching this value should be multiplied by a geometrical factor ν = 0.831 . . ., which converts the result obtained on the spherically symmetric lattice into the one for the standard cubic lattice [1]. This gives ∆ (1) ≈ −1.87, which is slightly above the O(v 6 ) value of Eq. (9), but is in a very good agreement with the value ∆ (1) ≈ −1.82 obtained from the fit of the O(v 4 ) result [13]. Numerically the one-loop linear artifact dominates the series in Eq. (8) for typical values of a and one may argue that its inclusion into the Wilson coefficient is mandatory. However the above analysis takes into account only a single Coulomb gluon exchange while the effect of multiple Coulomb exchanges is not parametrically suppressed and significantly changes the structure of the expansion in a as discussed in the next section.

B. Coulomb artifacts to all orders
Though we consider the properties of the heavy quarkonium bound states, the analysis of the previous sections involved the scattering amplitudes of the free quark and antiquark. This is sufficient if in the matching region the binding effects can be expanded in a regular series in α s . The Coulomb artifacts, however, are related to the dependence of the bound state characteristics on the lattice spacing, which cannot be described within the finite-order perturbation theory. In this case the matching procedure should be applied to the matrix elements of the effective action operators between the quarkonium states with the wave functions computed on the lattice and in the continuum. The relevant nonrelativistic Coulomb wave function in the continuum is well known. On the lattice it can be obtained in a straightforward way by solving the nonrelativistic Schrödinger equation as a difference equation for a given finite a. In the formal limit Λ QCD ≪ v 2 m q one can neglect the nonperturbative dynamics of strong interactions at long distance and the result obtained by numerical solution of the disretized Schrödinger equation provides the same bound state wave function as the real lattice simulations based on the functional integral approach.
Let us apply the above "Schrödinger matching" approach to the analysis of the hyperfine splitting. The relevant four-quark operator is generated by the magnetic gluon exchange and corresponds to the leading order spin-dependent amplitude 2 In coordinate space this local spin-flip operator is proportional to δ(x). The corresponding matrix element, which in fact determines the leading order hyperfine splitting, is proportional to |ψ(0)| 2 , where ψ(x) is the ground state quarkonium wave function. The Coulomb solution for this quantity takes into account the contribution of allorder Coulomb exchange diagrams including Fig. 1b. In the continuum it reads |ψ(0)| 2 = C 3 F α 3 s m 3 q /(8π). The lattice value of the wave function at the origin is obtained by numerical solution of the Schrödinger equation with the Coulomb Hamiltonian. It is performed on a spherically symmetric lattice, which retains the qualitative properties of the solution. To match the setup of real lattice simulations [19] we use the central difference discretization of the kinetic energy operator, which has O(a 4 ) local error. The boundary condition of the eigenstate problem is determind by the value of the exact continuum solution at sufficiently large distance, where the wave function is exponentially suppressed. Though the parameters of the bound state can be obtained for an arbitrary value of lattice spacing, we are interested in their behavior at small a. For the expansion of the ground state energy and the wave function at the origin about their continuum values we get whereā = C F α s am q /2 is the dimensionless lattice spacing in Coulomb units, and the rational coefficients of the expansion are conjectured from the high accuracy numerical result. The expression for the ground state energy is not required for our analysis and is given for completeness. Eq. (19) does not have a linear dependence on a. This may be expected since the integration of a second order difference equation with O(a 4 ) local discretization error gives O(a 2 ) global error of the solution (see e.g. [25]). Eq. (19) determines the difference between the lattice and continuum results for the matrix element of the leading order spin-flip operator Eq. (17). They can be matched by adding to the lattice NRQCD action the four-quark operator with d σ , which has the following coefficients of the expansion in am q Thus when the Coulomb effects are taken into account consistently to all orders in α s , the linear artifact in the four-quark matching coefficient is absent and the first nonvanishing term is quadratic in a. The qualitative difference between the one-loop and all-order expressions can be attributed to the modification of the bound quark propagator at large momentum, which make the all-order result less sensitive to the ultraviolet cutoff. We would like to emphasize that though the coefficient ∆ (2) is proportional to α 2 s , it gets contributions from all-order Coulomb exchange diagrams. This coefficient is changed by the higher order terms in the NRQCD action and is different for the standard cubic lattice, as in the case of the linear artifact discussed in the previous section. In principle, within the same method the Coulomb lattice artifacts can be evaluated for a given NRQCD action on a given lattice. However for practical applications they can be removed along with the nonperturbative artifacts by the extrapolation of the lattice data to a = 0, as it is  [2]. All data points include the statistical error and the uncertainty in the value of the lattice spacing. The error bars of (a) include also the uncertainty due to the higher order perturbative corrections. The difference between (a) and (b) data sets is mainly due to the spurious linear Coulomb artifact contributing to (b). discussed in the next section. The absence of the linear artifact is crucial for this procedure though.

III. DETERMINATION OF THE ENERGY SPECTRUM FROM THE LATTICE DATA
Let us now consider how the Coulomb artifacts affect the determination of the energy spectrum from the lattice data. The results of nonperturbative lattice NRQCD simulations are typically given for a ∼ 1/(vm b ) [2,5]. The use of relatively large values of the lattice spacing ensures the suppression of the unphysical 1/(am b ) n contributions, which become important at a ∼ 1/m b . At the same time it results in sizable Coulomb lattice artifacts proportional to a power of α s am b ∼ 1. In addition the lattice data include the nonperturbative lattice artifacts which scale as (aΛ QCD ) 2 and cannot be removed through the matching procedure discussed above. To minimize these effects the results of the lattice simulations are numerically extrapolated to a = 0. The extrapolation below a ∼ 1/m b in this case is justified because the numerical effect of the 1/(am b ) n terms on the data points is small. Since the radiatively improved lattice result is supposed to be free of linear artifacts, the extrapolation is performed through a constrained fit of the data points by a polynomial in a with vanishing linear term (see e.g. [1,2,5]). The correct treatment of the linear artifacts is therefore crucial for the extrapolation procedure. As it has been shown in the previous section by the analysis of the discretized Schrödinger equation, the linear Coulomb artifacts are absent in the bare lattice data. The contribution of the four-quark interaction to E hfs reads Thus the linear Coulomb artifact in the Wilson coefficient obtained by the direct numerical matching [13,14] results in s spurious linear dependence of the radiatively improved lattice data on a, which leads to a systematic error in the extrapolation procedure based on the fit with the vanishing linear term. At the same time the Wilson coefficient obtained by the asymptotic expansion about the continuum limit is free of the Coulomb artifacts and provides the correct functional dependence of the radiatively improved lattice data on a and therefore can be used for consistent extrapolation procedure. The numerical effect of the spurious linear artifact turns out to be very significant, see Fig. 2. The analysis of the hyperfine splitting [1] based on the O(v 6 ) lattice simulations and the matching coefficient from Eq. (4) gives E hfs = 51.5. At the same time the analysis [2] based on exactly the same lattice data gives E hfs = 60.0. The discrepancy between the results is mainly due to the spurious linear artifact and is well beyond the reported discretization/extrapolation uncertainty, which is below 3 MeV. Thus the analysis of the hyperfine splitting in Refs. [2,5,13,14] contains a systematic error and should be corrected. The result of the direct numerical matching can yet be used for the self-consistent analysis of the quarkonium spectrum through the decomposition of the form of Eqs. (7,8). After separating the logarithmic part, the result for the Wilson coefficient should be fitted by a polynomial in am q and the linear term of the expansion should be subtracted. In the case under consideration only the ∆ (0) term should be retained in d σ . The further analysis follows Ref. [1] with the coefficient ∆ (0) from Eq. (10) substituted by the one from Eq. (9). This gives the central value E hfs = 52.7, which is significantly below the result of Ref. [2] but in a very good agreement with the O(v 6 ) result of Ref. [1] given above.
Though the quadratic Coulomb artifact is eliminated by extrapolation, it is instructive to estimate its contribution to the dependence of the lattice data on a and corresponding uncertainty in the the extracted value of E hfs . The result of the fit for the hyperfine splitting can be represented as follows where Λ is the mass scale characterizing the approach of the lattice approximation to the continuum limit. Numerically one gets Λ ≈ 360 MeV for the O(v 4 ) and Λ ≈ 790 MeV for the O(v 6 ) lattice action [1]. On the other hand the quadratic Coulomb artifact with the coefficient Eq. (20) corresponds to which gives Λ ≈ 530 MeV for the values of the input parameters taken in the middle of a typical interval for the lattice spacing. Though Eq. (23) is obtained in a simplified model with the Coulomb Hamiltonian and on a spherical lattice, we can conclude that the quadratic Coulomb artifact to a large extent determines the dependence of the bare lattice result on a and can be used as a prior for the constrained fit. As we observed in Sect. II A the effect of the lattice artifacts is enhanced by the relativistic corrections since the contribution of the higher dimension operators is more sensitive to the ultraviolet momentum region. This explains a slower approach to the continuum limit and larger discretization errors of the extrapolation based on O(v 6 ) lattice data. The smaller discretization uncertainty balances the larger relativistic corrections in the O(v 4 ) case and both actions provide comparable total errors. The best estimate is obtained as the weighted average of two results [1] E hfs = 52.9 ± 5.5 MeV, which is slightly above the O(v 6 ) result. Hence Eq. (24) can be considered as an unambiguous and the most accurate lattice NRQCD prediction for the bottomonium hyperfine splitting available so far.

IV. SUMMARY AND CONCLUSION
In this paper we critically examined the matching procedure for the radiative improvement of the lattice NRQCD. We have demonstrated that the Wilson coefficients of the effective four-quark interaction obtained by the widely used direct numerical matching suffer from spurious linear Coulomb lattice artifacts, which result in a large systematic error in the predictions for the heavy quarkonium spectrum. This problem is solved by using the matching procedure based on the asymptotic expansion about the continuum limit. We also have shown how the direct numerical matching should be modified for a consistent treatment of the lattice artifacts.
Our analysis resolves the discrepancy between the most recent lattice NRQCD predictions for the bottomonium hyperfine splitting [1,2] in favour of the result of Ref. [1], Eq. (24), which reconciles the lattice predictions, continuum QCD [7], and the most accurate experimental data [11].