Coulomb Artifacts and Breakdown of Perturbative Matching in Lattice NRQCD
aa r X i v : . [ h e p - l a t ] D ec ALBERTA-THY-11-17
Coulomb Artifacts and Breakdown of PerturbativeMatching in Lattice NRQCD
A.A. Penin and A. Rayyan
Department of Physics, University of Alberta,Edmonton, Alberta T6G 2J1, Canada
Abstract
By studying an explicit analytical solution of the Schr¨odinger equation with theCoulomb potential on the lattice we demonstrate a breakdown of perturbative matchingfor the description of the Coulomb artifacts in lattice NRQCD, which leads to a largesystematic error in the predictions for the heavy quarkonium spectrum. The breakdownis a result of a fine interplay between the short and long distance effects specific tothe lattice regularization of NRQCD. We show how the problem can be solved bymatching the lattice and continuum results for the solution of the full Schr¨odingerequation without the expansion in the Coulomb interaction.
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 theoret-ical analysis of heavy quarkonium properties [5]. This method is entirely based on firstprinciples, allows for simultaneous treatment of dynamical heavy and light quarks, and givesa 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 contin-uum QCD, which properly takes into account the effect of the hard relativistic modes. For along time this procedure was thought to be well understood. Recently, however, a problemof the standard perturbative matching in the analysis of the heavy quark-antiquark boundstates on the lattice has been pointed out [6]. The problem is related to the description of the1attice Coulomb artifacts resulting from the effect of the space discretization on the Coulombbound state dynamics. The Coulomb artifacts appear as powers of a dimensionless combina-tion α 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 thelattice spacing should be cancelled in the final result for the physical quarkonium spectrumthrough the matching procedure and an inconsistent treatment of the artifacts may lead to alarge systematic error of the lattice NRQCD predictions. In the present paper, by studyingan explicit analytical solution of the Schr¨odinger equation on the lattice, we show that thestandard finite-order perturbative matching is insufficient in dealing with the Coulomb arti-facts, confirming the result of numerical analysis [6]. This appears counterintuitive since thelattice regularization is usually associated with a momentum cutoff at the scale 1 /a muchlarger than the scale α s m q of Coulomb dynamics and the corresponding short-distance effectsare supposed to be systematically described by the standard matching procedure. We findthat the failure of perturbative matching is a consequence of a fine interplay between theshort and long distance effects specific to the lattice regularization of NRQCD. The paperis organized as follows. In the next section we consider the case of the four-quark operatorsin the NRQCD Lagrangian where the perturbative matching breaks down in one loop, andshow how the problem can be related to the properties of the Coulomb wave function on thelattice within the method of “Schr¨odinger matching”. In Sects. 3 and 4 we describe an exactanalytical solution of the Schr¨odinger equation with Coulomb potential on a spatial lattice,which is used to model the Coulomb artifacts of the real lattice NRQCD simulations withoutthe expansion in α s , and explain the mechanism of the perturbative matching breakdown. InSect. 5 we discuss the application of our result to the analysis of the bottomonium spectrumin lattice NRQCD. The effect under consideration is characteristic to the four-quark operators contributionto the NRQCD Lagrangian δ L NRQCD = X i C F α s m q C i O i , (1)where O i = ψ † Γ i ψχ † c Γ i χ c , (2)2 ( χ c ) are the nonrelativistic Pauli spinors of quark (antiquark) field, Γ i is a matrix incolor 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 theone-particle irreducible quark-antiquark scattering amplitudes in QCD and NRQCD in agiven order in heavy quark velocity v . The Coulomb artifacts are associated with the staticCoulomb interaction between the nonrelativistic quark and antiquark. The nonrelativisticCoulomb 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 NRQCDamplitudes. In the one-loop approximation the relevant NRQCD diagram is given in Fig. 1,where the gray circle represents the effective O ( v ) interaction generated by a tree gluonexchange. It, in particular, includes the contact Fermi and Darwin terms proportionalto the delta-function of the distance between the heavy quark and antiquark, which havethe structure of Eq. (2). These terms determine the O ( α s ) perturbative corrections to theCoulomb energy levels given, up to an overall factor, by the matrix element h O i i of thecorresponding operators between the quarkonium states, which is proportional to the squareof the wave function at the origin | ψ (0) | .Within the standard perturbative matching the difference of the above diagram computedin 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 h O i i . On theother hand the diagram in Fig. 1 is quite special since it contains a static Coulomb exchange,which is responsible for the formation of the Coulomb-like heavy quarkonium states. Whensandwiched between the Coulomb states such an exchange can be absorbed into the Coulombwave function by the equation of motion, i.e. the diagram is already contained in theCoulomb matrix element of the O ( v ) tree gluon exchange. Thus the matching for thisdiagram as well as for the diagrams with any number of the Coulomb exchanges can beincorporated into the “Schr¨odinger matching” for the value of the wave function at the originintroduced in Ref. [6]. The corresponding correction to the energy level is proportional to (cid:18) − | ψ l (0) | | ψ c (0) | (cid:19) h O i i , (3)where ψ c ( ψ l ) stands for the Coulomb wave function computed in the continuum (on thelattice). One may suggest that for α s m q ≪ /a the two procedures should give equivalentresults and the Coulomb artifact contribution to C i should reproduce Eq. (3) in a given orderof the expansion in α s . However in Ref. [6] it has been demonstrated by numerical analysis The leading order NRQCD Lagrangian is O ( v ). O ( v ) one-particle irreduciblescattering amplitude in NRQCD. The double arrow and dashed lines stand for the nonrel-ativistic quark and Coulomb gluon propagators, respectively. The gray circle denotes theeffective O ( v ) interaction generated by a single gluon exchange. The symmetric diagram isnot shown.of the Schr¨odinger 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 Coulombproblem without the expansion in strong coupling constant. In particular the absence of thelinear Coulomb artifacts predicted by the one-loop matching has been observed. At the sametime a rather subtle mechanism of the perturbative matching breakdown cannot be easilytraced within a numerical approach. In the next sections we study an analytic solution ofthe Coulomb problem on the lattice and give a detailed explanation of this mechanism. We consider a stationary Schr¨odinger equation with the Coulomb potential V ( r ) = − C F α s /r on a spatial lattice and restrict the analysis to the spherically symmetric S -statesso that only the discretization of the radial coordinate r = an , n = 0 , , . . . is necessary.The time variable is kept continuous. This simplified model retains the main features of thereal lattice NRQCD simulations discussed in the next section. For the central-difference dis-cretization of the radial derivative the Schr¨odinger equation for the function R ( n ) = rψ l ( n )becomes R ( n + 1) − R ( n ) + R ( n −
1) + 2 a (cid:18) E + 1 na (cid:19) R ( n ) = 0 , (4)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 / − being the unit of length. The analyticalsolution of the eigenstate problem Eq. (4) is known and the corresponding eigenfunctions fordiscrete and continuum spectrum are found in terms of a hypergeometric function [7]. For4he ground state the solution takes a particulary simple form E l = − a (cid:0) (1 + a ) / − (cid:1) = E c (cid:18) − a O ( a ) (cid:19) , (5) ψ l ( n ) = 1 π / (1 + a ) / e − n sinh − a = ψ c ( r ) + O ( a ) , (6)where E c = − and ψ c ( r ) = √ π e − r are the well known continuum results for the groundstate energy and wave function, and the normalization condition reads4 πa ∞ X n =0 n | ψ l ( n ) | = 1 . (7)The wave function in the momentum space is defined by the Fourier transform˜ ψ l ( p ) = 4 πa ∞ X n =0 n ψ l ( n ) sin( nap ) nap (8)and reads ˜ ψ l ( p ) = 8 √ π (1 + a ) / ( a/ (cid:0) sin ( pa/ − a E l / (cid:1) sin( ap ) ap = ˜ ψ c ( p ) + O ( a ) , (9)where ˜ ψ c ( p ) = √ π ( p +1) is the continuum result in Coulomb units. A peculiar feature of theFourier transform is that ψ l (0) does not contribute to Eq. (8) due to the n factor. As aconsequence the inverse Fourier transform ψ l ( n ) = 12 π Z πa d p p ˜ ψ l ( p ) sin( nap ) nap (10)give the correct value of ψ l ( n ) only for the sites with n = 0. For n = 0 the correspondingintegral ψ p (0) ≡ π Z πa d p p ˜ ψ l ( p ) (11)takes the value ψ p (0) = 1 √ π − √ a − a ! = ψ c (0) (cid:16) − a O ( a ) (cid:17) , (12)while the original solution in the coordinate space is ψ l (0) = 1 √ π (1 + a ) / = ψ c (0) (cid:18) − a O ( a ) (cid:19) , (13) The subscript p in the Eq. (11) is used to distinguish the formal result of the inverse Fourier transformfrom the actual value of the solution Eq. (6) at n = 0.
5n full agreement with the results of the numerical analysis [6]. Though both values have thesame continuum limit, they approach it at a different rate since Eq. (12) has a nonvanishing O ( a ) term. In the Sect. 5 we show this to be precisely the origin of the spurious linearartifact in the result of perturbative matching. We should emphasize that Eq. (12) defines afunction, which does not satisfy the Schr¨odinger equation in the limit a →
0. Indeed, let usconsider the (forward) derivative of the wave function at the origin. Then for the solutionEq. (6) we get ψ ′ l (0) = ψ l (1) − ψ l (0) a = − ψ c (0) + O ( a ) , (14)which in the limit a → ψ ′ c (0) = − ψ c (0) . (15)At the same time if Eq. (12) is used to define the value of the wave function at the originwe get ψ ′ p (0) = ψ l (1) − ψ p (0) a = − ψ c (0) + O ( a ) , (16)which violates the continuum result at O (1).Finally we should note that while the linear term in Eq. (12) is universal for all the S -wave states, the coefficient of the quadratic terms in the expansion of the bound stateparameters in a is sensitive to the Coulomb dynamics, e.g. for the first excited state thiscoefficient in Eq. (13) changes from − / / Let us now consider the Coulomb contribution of the diagram Fig. 1 to a Wilson coefficient C i . As it has been pointed out this contribution is given by the difference of the continuumand lattice NRQCD one-loop expressions and can be computed in the kinematics where theheavy quarks are at rest and on the threshold. In the given order of the nonrelativisticexpansion it reads δC i = 4 C F α s π Z ∞ λ d p p D c ( p ) G c ( p ) − Z πa λ d p p D l ( p ) G l ( p ) ! , (17)where the integration over the time component of the virtual momentum has been performedby taking the residue of the heavy quark propagator, an auxiliary infrared cutoff λ is intro-duced to regularize the integrals at small momentum, and we switch back to the standardunits. The continuum Coulomb gluon and heavy quark propagators D c ( p ) = 1 p , G c ( p ) = m q p , (18)6orrespond to the leading order NRQCD action. Their lattice counterparts D l ( p ) = sin( ap ) ap ( a/ sin ( ap/ , G l ( p ) = m q ( a/ sin ( ap/
2) (19)exactly correspond to the discretization of the Schr¨odinger equation (4). After integrationwe get δC i = 12 C F α s am q + O (cid:0) ( λa ) (cid:1) . (20)Numerically, the main effect of the lattice regularization comes from the ultraviolet momen-tum cutoff at p = π/a . Indeed, if the continuum propagators are used in the second term ofEq. (17) the coefficient in Eq. (20) is changed only from 1 / /π = 0 . . . . .The lattice contribution to Eq. (17) can be obtained by the expansion of ˜ ψ l ( p ) in Eq. (11)to the first order in α s . Hence, as one may expect from the general arguments, Eq. (20)agrees with the expansion of the first factor in Eq. (3) to the first order in α s if the momentumspace result Eq. (12) is used to define the value of the wave function at the origin. HoweverEq. (12) does not agree with the actual solution Eq. (13) and results in a pathological wavefunction which does not satisfy the Schr¨odinger equation in the continuum limit at r = 0.Thus the Wilson coefficient Eq. (20) does not cancel the dependence of the O ( v ) tree gluonexchange matrix element on a and we observe a breakdown of perturbative matching in theanalysis of the Coulomb artifacts already in one loop. By contrast the Schr¨odinger matchingwith the exact solution of the Coulomb problem on the lattice gives the correct result forthe Coulomb artifacts to all orders in α s .We can trace the origin of this phenomenon to the fact that the relation Eq. (15) violatedby Eqs. (12, 20) follows from the cancellation of the singular kinetic and potential energyterms in the Schr¨odinger equation at r →
0. Indeed in the continuum the radial equationfor the S -wave function reads (cid:18) d dr + 2 r ddr − C F α s m q r + m q E (cid:19) ψ ( r ) = 0 . (21)Keeping the most singular terms in the limit r → (cid:18) ddr − r B (cid:19) ψ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) r → = 0 (22)which reproduces Eq. (15) in Coulomb units. The first (second) term in Eq. (22) correspondsto the kinetic (potential) energy contribution. Hence at the origin these terms should beconsidered on an equal footing while the standard matching treats the Coulomb potential asa perturbation. Evidently the above mechanism of the perturbative matching breakdown is In the Coulomb units this is equivalent to the expansion in a with fixed ap ∼ r ∼ a and take the limit r → r p ) / ( r p ) factor into the integrands of Eq. (17). The latticeintegral then reads m q Z πa λ d p p sin( r p ) r p sin( ap ) ap (cid:18) ( a/ ap/ (cid:19) = m q (cid:18) λ − π r (cid:19) + O (cid:0) ( λa ) (cid:1) , (23)which removes O ( a ) contribution from the Wilson coefficient Eq. (20) and brings it intoagreement with the Schr¨odinger matching result.We should emphasize the difference between matching calculations in lattice NRQCD andin NRQCD with an explicit momentum cutoff Λ UV ∼ /a , which is not plagued with theproblem discussed above. In the latter theory the properties of the solution of the Coulombproblem in the (continuum) coordinate space are significantly different from the solution ofthe finite difference equation Eq. (4). The Schr¨odinger equation in this case is a differentialequation 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 theintegral of the wave function in momentum space and the problem discussed in the previoussection does not exist. The correct behavior of ψ ( r ) and ψ ′ ( r ) at r → O ( α s m q / Λ UV ) term in the one-loop Wilson coefficient. Thus in the momentum cutoff scheme the results of the perturbativeand Schr¨odinger matching do agree.The discretization method and the action of the model discussed above do not exactlyreproduce the ones used in the real lattice simulations. The inclusion of the relativisticcorrections, a different discretization of the gluon field or the use of a cubic lattice wouldchange the dependence of ψ l (0) and C i on a . However the Wilson coefficient obtainedwithin the perturbative matching always has a linear dependence on a from the momentumcutoff in Eq. (17). At the same time the solution of the Schr¨odinger equation with thecentral difference discretization of the kinetic energy operator has only O ( a ) global error(see e.g. [9]) and is free of the linear artifacts. This is confirmed by a numerical analysisof the discretized Schr¨odinger-Pauli equation at O ( v ) on a cubic lattice [10], which with a8ood precision rules out the linear dependence of the bound state parameters on the latticespacing. In general the Schr¨odinger matching can be performed numerically for any givenHamiltonian and the lattice used in the simulations, but even a simple analysis based onEq. (13) gives a good estimate of the dependence of the bare lattice NRQCD data for thehyperfine splitting in bottomonium on a [6]. Let us now consider how the above analysis affects the determination of the bottomoniumspectrum in the radiatively improved lattice NRQCD. The results of nonperturbative latticeNRQCD simulations are typically given for a ∼ / ( vm b ) [5, 11]. The use of relativelylarge values of the lattice spacing ensures the suppression of the singular ultraviolet cutoffdependence from the higher order 1 / ( am b ) n terms, which are not removed by the finite ordermatching and become important at a ∼ /m b . At the same time it results in sizable Coulomblattice artifacts proportional to a power of α s am b ∼
1. The correct treatment of the artifactsis therefore crucial for the analysis.In actual lattice simulations, the quarkonium bound state parameters are extracted fromthe asymptotic behavior of the quark-antiquark propagator at large Euclidean time. Ne-glecting the retardation and long distance nonperturbative effects, which do not affect theCoulomb artifacts under consideration, this method should reproduce the properties of thesolution of the Schr¨odinger equation for a given NRQCD Hamiltonian on the spatial lattice.As we have shown above, the perturbative matching of the four-quark operators does notcorrectly account for the Coulomb artifacts and, for a ∼ / ( vm b ) and v ∼ α s , results in O (1) error in the prediction for the spectrum. For example, the one-loop matching of thespin-flip four-quark operator with the spurious linear Coulomb artifact gives the value of thebottomonium 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 heavyquarkonium dynamics.In practice the effect of the lattice artifacts is reduced by numerical extrapolation of thedata to a = 0 [5, 11]. The extrapolation below a ∼ /m b in this case is justified because forthe typical values of lattice spacing the numerical effect of the 1 / ( am b ) n terms on the datapoints 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 the9onfinement scale Λ QCD . The problem of the perturbative matching breakdown, however,is not fully fixed by this procedure. Since the radiatively improved lattice result is supposedto be free of linear artifacts, the extrapolation is performed through a constrained fit of thedata points by a polynomial in a with vanishing linear term. Since the bare lattice data arefree of the linear Coulomb artifacts the one-loop perturbative matching in fact introduces alinear dependence of the radiatively improved result on a , which leads to a systematic errorof the fit. For example, in the analysis of the bottomonium hyperfine splitting [11] based onthe one-loop perturbative matching [12, 14] this error exceeds 10% of the total result and isbeyond the estimated uncertainty interval of the radiatively improved O ( v ) lattice NRQCDanalysis [6]. At the same time the perturbative matching can be used for the self-consistentanalysis of the quarkonium spectrum within the above extrapolation scheme if the Coulombartifacts are removed from the Wilson coefficients by means of the asymptotic expansion [15]or a numerical fit [6]. Acknowledgements
We would like to thank Tao Liu for useful discussions and comments on the manuscript.This work is supported in part by NSERC. The work of A.P. is supported in part by thePerimeter Institute for Theoretical Physics.
References [1] W. E. Caswell and G. P. Lepage, Phys. Lett. B , 437 (1986).[2] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D , 1125 (1995).[3] B. A. Thacker and G. P. Lepage, Phys. Rev. D , 196 (1991).[4] G. P. Lepage, L. Magnea, C. Nakhleh, U. Magnea and K. Hornbostel, Phys. Rev. D ,4052 (1992).[5] R. J. Dowdall et al. [HPQCD Collaboration], Phys. Rev. D , 054509 (2012).[6] T. Liu, A. A. Penin and A. Rayyan, JHEP , 084 (2017).[7] A. A. Kvitsinsky, J. Phys. A: Math. Gen. , 65 (1992). For the bottomonium ground state this contribution is numerically suppressed with respect to theCoulomb artifacts since Λ
QCD ≪ α s m q . , 111301 (2000).[9] E. Hairer, S.P. Nørsett and G. Wanner, “ Solving Ordinary Differential Equations I:Nonstiff Problems ”, Springer Berlin Heidelberg (2008).[10] G. S. Bali and P. Boyle, Phys. Rev. D , 114504 (1999).[11] R. J. Dowdall et al. [HPQCD Collaboration], Phys. Rev. D , 031502 (2014) [Erratum-ibid. , 039904 (2015)].[12] T. C. Hammant, A. G. Hart, G. M. von Hippel, R. R. Horgan and C. J. Monahan,Phys. Rev. Lett. , 112002 (2011); [Erratum-ibid. 115 (2015) 039901].[13] B. A. Kniehl, A. A. Penin, A. Pineda, V. A. Smirnov and M. Steinhauser, Phys. Rev.Lett. , 242001 (2004). [Erratum-ibid. , 199901 (2010)].[14] T. C. Hammant, A. G. Hart, G. M. von Hippel, R. R. Horgan and C. J. Monahan,Phys. Rev. D , 014505 (2013); [Erratum-ibid. , 119904 (2015)].[15] M. Baker, A. A. Penin, D. Seidel and N. Zerf, Phys. Rev. D92