Coulomb artifacts and breakdown of perturbative matching in lattice NRQCD

By studying an explicit analytical solution of the Schrödinger equation with the Coulomb potential on the lattice we demonstrate a breakdown of perturbative matching for the description of the Coulomb artifacts in lattice NRQCD, which leads to a large systematic error in the predictions for the heavy quarkonium spectrum. The breakdown is a result of a fine interplay between the short and long distance effects specific to the lattice regularization of NRQCD. We show how the problem can be solved by matching the lattice and continuum results for the solution of the full Schrödinger equation without the expansion in the Coulomb interaction.


Introduction
The lattice simulations within the effective theory of nonrelativistic QCD (NRQCD) [1][2][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. A crucial part of the approach is the matching of lattice NRQCD to the full theory of continuum QCD, which properly takes into account the effect of the hard relativistic modes. For a long time this procedure was thought to be well understood. Recently, however, a problem of the standard perturbative matching in the analysis of the heavy quark-antiquark bound states on the lattice has been pointed out [6]. The problem is related to the description of the lattice Coulomb artifacts resulting from the effect of the space discretization on the Coulomb bound state dynamics. The Coulomb artifacts appear as powers of a dimensionless combination α s am q of the strong coupling constant α s , lattice spacing a, and the heavy quark mass m q in the parameters of the bound states evaluated on the lattice. This dependence on the lattice spacing should be cancelled in the final result for the physical quarkonium spectrum through the matching procedure and an inconsistent treatment of the artifacts may lead to a large systematic error of the lattice NRQCD predictions. In the present paper, by studying an explicit analytical solution of the Schrödinger equation on the lattice, we show that the standard finite-order perturbative matching is insufficient in dealing with the Coulomb artifacts, confirming the result of numerical analysis [6]. This appears counterintuitive since the lattice regularization is usually associated with a momentum cutoff at the scale 1/a much larger than the scale α s m q of Coulomb dynamics and the corresponding short-distance effects are supposed to be systematically described by the standard matching procedure. We find that the failure of perturbative matching is a consequence of a fine interplay between the short and long distance effects specific to the lattice regularization of NRQCD. The paper is organized as follows. In the next JHEP12(2017)007 section we consider the case of the four-quark operators in the NRQCD Lagrangian where the perturbative matching breaks down in one loop, and show how the problem can be related to the properties of the Coulomb wave function on the lattice within the method of "Schrödinger matching". In sections 3 and 4 we describe an exact analytical solution of the Schrödinger equation with Coulomb potential on a spatial lattice, which is used to model the Coulomb artifacts of the real lattice NRQCD simulations without the expansion in α s , and explain the mechanism of the perturbative matching breakdown. In section 5 we discuss the application of our result to the analysis of the bottomonium spectrum in lattice NRQCD.

Perturbative and Schrödinger matching of four-quark operators
The effect under consideration is characteristic to the four-quark operators contribution to the NRQCD Lagrangian ψ (χ c ) are the nonrelativistic Pauli spinors of quark (antiquark) field, Γ i is a matrix in color and spinor space and C i is the corresponding Wilson coefficient. The coefficients C i vanish in the Born approximation and is determined order by order in α s by equating the one-particle irreducible quark-antiquark scattering amplitudes in QCD and NRQCD in a given order in heavy quark velocity v. The Coulomb artifacts are associated with the static Coulomb interaction between the nonrelativistic quark and antiquark. The nonrelativistic Coulomb dynamics is not sensitive to the high momentum region and, to the leading order in v, the result for such a contribution obtained within continuum QCD and NRQCD coincide. Thus the matching is determined by the difference between continuum and lattice NRQCD amplitudes. In the one-loop approximation the relevant NRQCD diagram is given in figure 1, where the gray circle represents the effective O(v 4 ) interaction generated by a tree gluon exchange. 1 It, in particular, includes the contact Fermi and Darwin terms proportional to the delta-function of the distance between the heavy quark and antiquark, which have the structure of eq. (2.2). These terms determine the O(α 2 s ) perturbative corrections to the Coulomb energy levels given, up to an overall factor, by the matrix element O i of the corresponding operators between the quarkonium states, which is proportional to the square of the wave function at the origin |ψ(0)| 2 .
Within the standard perturbative matching the difference of the above diagram computed in the continuum and on the lattice should be included into the one-loop Wilson coefficients C i . The corresponding corrections to the energy level are proportional to C i O i . On the other hand the diagram in figure 1 is quite special since it contains a static Coulomb exchange, which is responsible for the formation of the Coulomb-like heavy quarkonium states. When sandwiched between the Coulomb states such an exchange can be absorbed JHEP12(2017)007 into the Coulomb wave function by the equation of motion, i.e. the diagram is already contained in the Coulomb matrix element of the O(v 4 ) tree gluon exchange. Thus the matching for this diagram as well as for the diagrams with any number of the Coulomb exchanges can be incorporated into the "Schrödinger matching" for the value of the wave function at the origin introduced in ref. [6]. The corresponding correction to the energy level is proportional to where ψ c (ψ l ) stands for the Coulomb wave function computed in the continuum (on the lattice). One may suggest that for α s m q ≪ 1/a the two procedures should give equivalent results and the Coulomb artifact contribution to C i should reproduce eq. (2.3) in a given order of the expansion in α s . However in ref. [6] it has been demonstrated by numerical analysis of the Schrödinger equation that the effects of space discretization on the Coulomb dynamics cannot be accounted for order by order in α s and require the exact solution of the Coulomb problem without the expansion in strong coupling constant. In particular the absence of the linear Coulomb artifacts predicted by the one-loop matching has been observed. At the same time a rather subtle mechanism of the perturbative matching breakdown cannot be easily traced within a numerical approach. In the next sections we study an analytic solution of the Coulomb problem on the lattice and give a detailed explanation of this mechanism.

Analytical solution of the Coulomb problem on the lattice
We consider a stationary Schrödinger equation with the Coulomb potential V (r) = −C F α s /r on a spatial lattice and restrict the analysis to the spherically symmetric Sstates so that only the discretization of the radial coordinate r = an, n = 0, 1, . . . is necessary. The time variable is kept continuous. This simplified model retains the main features of the real lattice NRQCD simulations discussed in the next section. For the central-difference discretization of the radial derivative the Schrödinger equation for the function R(n) = rψ l (n) becomes

JHEP12(2017)007
with the boundary condition R(0) = 0. Throughout this section we use the Coulomb units, with the reduced Bohr radius r B = (C F α s m q /2) −1 being the unit of length. The analytical solution of the eigenstate problem eq. (3.1) is known and the corresponding eigenfunctions for discrete and continuum spectrum are found in terms of a hypergeometric function [7]. For the ground state the solution takes a particulary simple form where E c = − 1 2 and ψ c (r) = 1 √ π e −r are the well known continuum results for the ground state energy and wave function, and the normalization condition reads The wave function in the momentum space is defined by the Fourier transform and readsψ whereψ c (p) = 8 √ π (p 2 +1) 2 is the continuum result in Coulomb units. A peculiar feature of the Fourier transform is that ψ l (0) does not contribute to eq. (3.5) due to the n 2 factor. As a consequence the inverse Fourier transform ψ l (n) = 1 2π 2 π a 0 dp p 2ψ l (p) give the correct value of ψ l (n) only for the sites with n = 0. For n = 0 the corresponding integral 2 ψ p (0) ≡ 1 2π 2 π a 0 dp p 2ψ l (p) (3.8) takes the value while the original solution in the coordinate space is JHEP12 (2017)007 in full agreement with the results of the numerical analysis [6]. Though both values have the same continuum limit, they approach it at a different rate since eq. (3.9) has a nonvanishing O(a) term. In the section 5 we show this to be precisely the origin of the spurious linear artifact in the result of perturbative matching. We should emphasize that eq. (3.9) defines a function, which does not satisfy the Schrödinger equation in the limit a → 0. Indeed, let us consider the (forward) derivative of the wave function at the origin. Then for the solution eq. (3.3) we get which in the limit a → 0 recovers the property of the continuum solution At the same time if eq. (3.9) is used to define the value of the wave function at the origin we get which violates the continuum result at O(1).
Finally we should note that while the linear term in eq. (3.9) is universal for all the S-wave states, the coefficient of the quadratic terms in the expansion of the bound state parameters in a is sensitive to the Coulomb dynamics, e.g. for the first excited state this coefficient in eq. (3.10) changes from −1/4 to 1/16.

Breakdown of perturbative matching
Let us now consider the Coulomb contribution of the diagram figure 1 to a Wilson coefficient C i . As it has been pointed out this contribution is given by the difference of the continuum and lattice NRQCD one-loop expressions and can be computed in the kinematics where the heavy quarks are at rest and on the threshold. In the given order of the nonrelativistic expansion it reads where the integration over the time component of the virtual momentum has been performed by taking the residue of the heavy quark propagator, an auxiliary infrared cutoff λ is introduced to regularize the integrals at small momentum, and we switch back to the standard units. The continuum Coulomb gluon and heavy quark propagators correspond to the leading order NRQCD action. Their lattice counterparts D l (p) = sin(ap) ap (a/2) 2 sin 2 (ap/2) , G l (p) = m q (a/2) 2 sin 2 (ap/2) (4.3)

JHEP12(2017)007
exactly correspond to the discretization of the Schrödinger equation (3.1). After integration we get Numerically, the main effect of the lattice regularization comes from the ultraviolet momentum cutoff at p = π/a. Indeed, if the continuum propagators are used in the second term of eq. (4.1) the coefficient in eq. (4.4) is changed only from 1/2 to 4/π 2 = 0.405 . . .. The lattice contribution to eq. (4.1) can be obtained by the expansion ofψ l (p) in eq. (3.8) to the first order in α s . 3 Hence, as one may expect from the general arguments, eq. We can trace the origin of this phenomenon to the fact that the relation eq. (3.12) violated by eqs. (3.9, 4.4) follows from the cancellation of the singular kinetic and potential energy terms in the Schrödinger equation at r → 0. Indeed in the continuum the radial equation for the S-wave function reads Keeping the most singular terms in the limit r → 0 we get which reproduces eq. (3.12) in Coulomb units. The first (second) term in eq. (4.6) corresponds to the kinetic (potential) energy contribution. Hence at the origin these terms should be considered on an equal footing while the standard matching treats the Coulomb potential as a perturbation. Evidently the above mechanism of the perturbative matching breakdown is specific to the lattice regularization and the contact interaction in the NRQCD Lagrangian. We can regularize the contact interaction by separating the quark and antiquark fields with a small spatial interval r 0 ∼ a and take the limit r 0 → 0 after the matching is done. This introduces an additional sin(r 0 p)/(r 0 p) factor into the integrands of eq. (4.1). The lattice integral then reads m q π a λ dp p 2 sin(r 0 p) r 0 p sin(ap) ap (a/2) sin(ap/2) which removes O(a) contribution from the Wilson coefficient eq. (4.4) and brings it into agreement with the Schrödinger matching result.

JHEP12(2017)007
We should emphasize the difference between matching calculations in lattice NRQCD and in NRQCD with an explicit momentum cutoff Λ U V ∼ 1/a, which is not plagued with the problem discussed above. In the latter theory the properties of the solution of the Coulomb problem in the (continuum) coordinate space are significantly different from the solution of the finite difference equation eq. (3.1). The Schrödinger equation in this case is a differential equation and its regular solution satisfies the conditions of the Fourier inversion theorem. Hence the value of the wave function at the origin is unambiguously determined by the integral of the wave function in momentum space and the problem discussed in the previous section does not exist. The correct behavior of ψ(r) and ψ ′ (r) at r → 0 then follows from the continuity and smoothness of the solution. A comprehensive analysis of the four-fermion operator matching with an explicit momentum cutoff can be found in ref. [8] in a context of the NRQED calculation of the radiative corrections to the orthopositronium decay rate. In [8] a numerical solution of the coordinate-space Schrödinger equation has been obtained for the Hamiltonian defined in a momentum cutoff regularization scheme. It has been found that the dependence of the value of the resulting wave function at the origin on the cutoff includes a linear Coulomb artifact which is cancelled by the O(α s m q /Λ U V ) term in the one-loop Wilson coefficient. Thus in the momentum cutoff scheme the results of the perturbative and Schrödinger matching do agree.
The discretization method and the action of the model discussed above do not exactly reproduce the ones used in the real lattice simulations. The inclusion of the relativistic corrections, a different discretization of the gluon field or the use of a cubic lattice would change the dependence of ψ l (0) and C i on a. However the Wilson coefficient obtained within the perturbative matching always has a linear dependence on a from the momentum cutoff in eq. (4.1). At the same time the solution of the Schrödinger equation with the central difference discretization of the kinetic energy operator has only O(a 2 ) global error (see e.g. [9]) and is free of the linear artifacts. This is confirmed by a numerical analysis of the discretized Schrödinger-Pauli equation at O(v 4 ) on a cubic lattice [10], which with a good precision rules out the linear dependence of the bound state parameters on the lattice spacing. In general the Schrödinger matching can be performed numerically for any given Hamiltonian and the lattice used in the simulations, but even a simple analysis based on eq. (3.10) gives a good estimate of the dependence of the bare lattice NRQCD data for the hyperfine splitting in bottomonium on a [6].

Coulomb artifacts and bottomonium spectrum in lattice NRQCD
Let us now consider how the above analysis affects the determination of the bottomonium spectrum in the radiatively improved lattice NRQCD. The results of nonperturbative lattice NRQCD simulations are typically given for a ∼ 1/(vm b ) [5,11]. The use of relatively large values of the lattice spacing ensures the suppression of the singular ultraviolet cutoff dependence from the higher order 1/(am b ) n terms, which are not removed by the finite order matching and 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. The correct treatment of the artifacts is therefore crucial for the analysis.

JHEP12(2017)007
In actual lattice simulations, the quarkonium bound state parameters are extracted from the asymptotic behavior of the quark-antiquark propagator at large Euclidean time. Neglecting the retardation and long distance nonperturbative effects, which do not affect the Coulomb artifacts under consideration, this method should reproduce the properties of the solution of the Schrödinger equation for a given NRQCD Hamiltonian on the spatial lattice. As we have shown above, the perturbative matching of the four-quark operators does not correctly account for the Coulomb artifacts and, for a ∼ 1/(vm b ) and v ∼ α s , results in O(1) error in the prediction for the spectrum. For example, the one-loop matching of the spin-flip four-quark operator with the spurious linear Coulomb artifact gives the value of the bottomonium hyperfine splitting [12], which overshoots the predictions of perturbative QCD [13] by almost a factor of two, in clear conflict with the general understanding of the heavy quarkonium dynamics.
In practice the effect of the lattice artifacts is reduced by numerical extrapolation of the data to a = 0 [5,11]. The extrapolation below a ∼ 1/m b in this case is justified because for the typical values of lattice spacing the numerical effect of the 1/(am b ) n terms on the data points is small. This extrapolation effectively removes all the lattice artifacts including the (aΛ QCD ) n terms associated with the effect of lattice regularization on the dynamics at the confinement scale Λ QCD . 4 The problem of the perturbative matching breakdown, however, is not fully fixed by this procedure. 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. Since the bare lattice data are free of the linear Coulomb artifacts the one-loop perturbative matching in fact introduces a linear dependence of the radiatively improved result on a, which leads to a systematic error of the fit. For example, in the analysis of the bottomonium hyperfine splitting [11] based on the one-loop perturbative matching [12,14] this error exceeds 10% of the total result and is beyond the estimated uncertainty interval of the radiatively improved O(v 6 ) lattice NRQCD analysis [6]. At the same time the perturbative matching can be used for the self-consistent analysis of the quarkonium spectrum within the above extrapolation scheme if the Coulomb artifacts are removed from the Wilson coefficients by means of the asymptotic expansion [15] or a numerical fit [6].

JHEP12(2017)007
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.