Massive Fermions in Lattice Gauge Theory

TL;DR

This paper develops a unified lattice fermion formulation valid for all quark masses, ensuring smooth connections between light and heavy quark regimes, and analyzes cutoff artifacts and improvement criteria.

Contribution

It introduces a comprehensive lattice fermion approach applicable to all quark masses, incorporating interactions and external currents without mass limit expansions.

Findings

Cutoff artifacts factorize into mass- and momentum-dependent terms.

Order $(m_q a)^s$ artifacts are avoided even for large $m_q a$.

Provides interpretation of Wilson and Sheikholeslami-Wohlert actions for nonrelativistic fermions.

Abstract

This paper presents a formulation of lattice fermions applicable to all quark masses, large and small. We incorporate interactions from previous light-fermion and heavy-fermion methods, and thus ensure a smooth connection to these limiting cases. The couplings in improved actions are evaluated for arbitrary fermion mass~$m_q$, without expansions around small- or large-mass limits. We treat both the action and external currents. By interpreting on-shell improvement criteria through the lattice theory's Hamiltonian, one finds that cutoff artifacts factorize into the form $b_n(m_qa)[\vek{p}a]^{s_n}$, where $\vek{p}$ is a momentum characteristic of the system under study, $s_n$ is related to the dimension of the $n$th interaction, and $b_n(m_qa)$ is a bounded function, numerically always~$\order(1)$ or less. In heavy-quark systems $\vek{p}$ is typically rather smaller than the fermion mass~$m_q$. Therefore, artifacts of order $(m_qa)^s$ do not arise, even when $m_qa\gsim1$. An important by-product of our analysis is an interpretation of the Wilson and Sheikholeslami-Wohlert actions applied to nonrelativistic fermions.

Tables

Table 1. Table 1: Interactions that could appear in the action, with and without axis-interchange symmetry (a.i.s.).

dim	w/ a.i.s.	w/o a.i.s.
3	$\overline{\psi }\psi$	$\overline{\psi }\psi$
4	$\overline{\psi }/D\psi$	$\overline{\psi }{\gamma }_0{D}_0\psi$	$\overline{\psi }𝜸\cdot 𝑫\psi$
5	$\overline{\psi }/D^2\psi$	$\overline{\psi }{D}_0^2\psi$	$\overline{\psi }{\left(𝜸\cdot 𝑫\right)}^2\psi$
	$i\overline{\psi }{\sigma }_{\mu \nu }{F}_{\mu \nu }\psi$	$i\overline{\psi }𝚺\cdot 𝑩\psi$	$\overline{\psi }𝜶\cdot 𝑬\psi$
		$\overline{\psi }[{\gamma }_0{D}_0,𝜸\cdot 𝑫]\psi$

Table 2. Table 2: Estimates of the size of dimension-three, -four, and -five interactions in systems with only light quarks, with one heavy quark, and in quarkonia. The latter two columns use m Q subscript 𝑚 𝑄 m_{Q} to emphasize the heavy -quark mass.

$H_n$	only light	heavy-light	quarkonia
$E_0$	${\Lambda }_{\mathrm{QCD}}$	${\Lambda }_{\mathrm{QCD}}$	$m_q{v}^2$
$\overline{\Psi }\Psi$	$1$	$1$	$1$
$\overline{\Psi }𝜸\cdot 𝑫\Psi$	${\Lambda }_{\mathrm{QCD}}$	${\Lambda }_{\mathrm{QCD}}^2/m_Q$	$m_Q{v}^2$
$\overline{\Psi }{𝑫}^2\Psi$	${\Lambda }_{\mathrm{QCD}}^2$	${\Lambda }_{\mathrm{QCD}}^2$	$m_Q^2{v}^2$
$i\overline{\Psi }𝚺\cdot 𝑩\Psi$	${\Lambda }_{\mathrm{QCD}}^2$	${\Lambda }_{\mathrm{QCD}}^2$	$m_Q^2{v}^4$
$\overline{\Psi }𝜶\cdot 𝑬\Psi$	${\Lambda }_{\mathrm{QCD}}^2$	${\Lambda }_{\mathrm{QCD}}^3/m_Q$	$m_Q^2{v}^4$

Table 3. Table 3: Estimates of the size of dimension-six bilinear interactions in systems with only light quarks, with one heavy quark, and in quarkonia. For the actions considered here, the interactions below the line do not arise at tree level.

$H_n$	only light	heavy-light	quarkonia
$\overline{\Psi }{\gamma }_0[𝜸\cdot 𝑫,𝜸\cdot 𝑬]\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^3$	$m_Q^3{v}^4$
$\overline{\Psi }{\left(𝜸\cdot 𝑫\right)}^3\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^4/m_Q$	$m_Q^3{v}^4$
$\overline{\Psi }{\gamma }_i{D}_i^3\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^4/m_Q$	$m_Q^3{v}^4$
$\overline{\Psi }\{𝜸\cdot 𝑫,i𝚺\cdot 𝑩\}\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^4/m_Q$	$m_Q^3{v}^6$
$i\overline{\Psi }{\gamma }_0𝚺\cdot \left(D_0𝑩\right)\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^4/m_Q$	$m_Q^3{v}^6$
$\overline{\Psi }𝜸\cdot \left(D_0𝑬\right)\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^5/m_Q^2$	$m_Q^3{v}^6$
$\overline{\Psi }{\gamma }_0\left(𝑫\cdot 𝑬-𝑬\cdot 𝑫\right)\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^3$	$m_Q^3{v}^4$
$i\overline{\Psi }𝜸\cdot \left(𝑫\times 𝑩+𝑫\times 𝑩\right)\Psi$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^4/m_Q$	$m_Q^3{v}^6$

Table 4. Table 4: Estimates of the size of dimension-six four-fermion interactions in systems with only light quarks, with one heavy quark, and in quarkonia. The role of these interactions can only be treated consistently in concert with the gauge-field action. The interactions above the line can then be considered redundant, while those below the line do not arise at tree level.

$H_n$	only light	heavy-light	quarkonia
${\left(\overline{\psi }{\gamma }_0{t}^a\psi \right)}^2$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^3$	$m_Q^3{v}^6$
${\left(\overline{\psi }{\gamma }_i{t}^a\psi \right)}^2$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^3$	$m_Q^3{v}^6$
${\left(\overline{\psi }\Gamma t^a\psi \right)}^2$, $\Gamma \ne {\gamma }_{\mu }$	${\Lambda }_{\mathrm{QCD}}^3$	${\Lambda }_{\mathrm{QCD}}^3$	$m_Q^3{v}^6$

Table 5. Table 5: Estimates of the size of the dimension-seven interactions that arise in designing an action for quarkonia with an accuracy of m Q ​ v 4 subscript 𝑚 𝑄 superscript 𝑣 4 m_{Q}v^{4} . (There are many other interactions needed to ensure O ​ ( ( Λ QCD ​ a ) 4 ) O superscript subscript Λ QCD 𝑎 4 {\rm O}((\Lambda_{\rm QCD}a)^{4}) accuracy in all-light and heavy-light systems.)

$H_n$	only light	heavy-light	quarkonia
$\overline{\Psi }{\left({𝑫}^2\right)}^2\Psi$	${\Lambda }_{\mathrm{QCD}}^4$	${\Lambda }_{\mathrm{QCD}}^4$	$m_Q^4{v}^4$
$\overline{\Psi }{D}_i^4\Psi$	${\Lambda }_{\mathrm{QCD}}^4$	${\Lambda }_{\mathrm{QCD}}^4$	$m_Q^4{v}^4$

Equations

$$S=\sum _n{c}_n{S}_n,$$

$$S={\mathcal{S}}_{\mathrm{RT}}-\delta S,$$

$${\mathcal{S}}_{\mathrm{RT}}(m_q/{\Lambda }_{\mathrm{QCD}},{\Lambda }_{\mathrm{QCD}}a)=\sum _n{c}_n(m_{0}a,g_0^2)S_n,$$

$$c_n(m_{0}a,g_0^2)=\sum _{l=0}^{\infty }{g}_0^{2l}{c}_n^{[l]}(m_{0}a).$$

$$\delta S=\sum _n\sum _{l=L_n+1}^{\infty }{g}_0^{2l}{c}_n^{[l]}(m_{0}a)S_n.$$

$$\begin{array}[]{l}S_{0}={\displaystyle\sum_{n}}\bar{\psi}_{n}\psi_{n}-\kappa_{t}{\displaystyle\sum_{n}}\left[\rule{0.0pt}{9.70001pt}\bar{\psi}_{n}(1-\gamma_{0})U_{n,0}\psi_{n+\hat{0}}+\bar{\psi}_{n+\hat{0}}(1+\gamma_{0})U^{\dagger}_{n,0}\psi_{n}\right]\\[10.00002pt] \hskip 63.00012pt-\,\kappa_{s}{\displaystyle\sum_{n,i}}\left[\rule{0.0pt}{9.70001pt}\bar{\psi}_{n}(r_{s}-\gamma_{i})U_{n,i}\psi_{n+\hat{\imath}}+\bar{\psi}_{n+\hat{\imath}}(r_{s}+\gamma_{i})U_{n,i}^{\dagger}\psi_{n}\right].\end{array}$$

$$S_{B}=\mbox{\small$\frac{i}{2}$}c_{B}\kappa_{s}\sum_{n;i,j,k}\varepsilon_{ijk}\bar{\psi}_{n}\sigma_{ij}B_{n;k}\psi_{n},$$

$$S_E=i{c}_E{\kappa }_s\sum _{n;i}{\bar{\psi }}_n{\sigma }_{0i}{E}_{n;i}{\psi }_n,$$

$${\psi }_n=\frac{a^{3/2}}{\sqrt{2{\kappa }_t}}\psi (x)$$

$$m_{0}a=\frac{1}{2{\kappa }_t}-[1+r_s\zeta (d-1)],$$

$$\begin{array}[]{c}S_{0}=m_{0}{\displaystyle\int}\bar{\psi}(x)\psi(x)+{\displaystyle\int}\left[\bar{\psi}(x)\mbox{\small$\frac{1}{2}$}(1+\gamma_{0})D_{0}^{-}\psi(x)-\bar{\psi}(x)\mbox{\small$\frac{1}{2}$}(1-\gamma_{0})D_{0}^{+}\psi(x)\right]\\[10.00002pt] +\zeta{\displaystyle\int}\bar{\psi}(x)\mbox{\boldmath$\gamma$}\!\cdot\!\mbox{\boldmath$D$}\psi(x)-\,\mbox{\small$\frac{1}{2}$}ar_{s}\zeta{\displaystyle\int}\bar{\psi}(x)\triangle^{(3)}\psi(x),\end{array}$$

$$T_{\pm\mu}\psi(x)=U_{\pm\mu}(x)\psi(x\pm a\hat{\mu}),\hskip 14.22636pt\bar{\psi}(x)\vbox{\halign{#\cr$\scriptstyle\leftarrow$\crcr\nointerlineskip\cr$\hfil\displaystyle{T}\hfil$\crcr}}_{\pm\mu}=\bar{\psi}(x\mp a\hat{\mu})U_{\pm\mu}(x\mp a\hat{\mu}),$$

$$\begin{array}[]{l}D_{0}^{+}\psi=a^{-1}(T_{0}-1)\psi,\\[11.99998pt] D_{0}^{-}\psi=a^{-1}(1-T_{-0})\psi,\\[11.99998pt] D_{i}\psi=(2a)^{-1}(T_{i}-T_{-i})\psi,\\[11.99998pt] \triangle^{(3)}\psi=a^{-2}{\displaystyle\sum_{i=1}^{3}}(T_{i}+T_{-i}-2)\psi,\end{array}$$

$$\gamma_{0}=\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right),\hskip 28.45274pt\mbox{\boldmath$\gamma$}=\left(\begin{array}[]{cc}0&\mbox{\boldmath$\sigma$}\\ \mbox{\boldmath$\sigma$}&0\end{array}\right),$$

$$\mbox{\boldmath$\alpha$}=\gamma_{0}\mbox{\boldmath$\gamma$}=\left(\begin{array}[]{cc}0&\mbox{\boldmath$\sigma$}\\ -\mbox{\boldmath$\sigma$}&0\end{array}\right),\hskip 28.45274pt\mbox{\boldmath$\Sigma$}=\left(\begin{array}[]{cc}\mbox{\boldmath$\sigma$}&0\\ 0&\mbox{\boldmath$\sigma$}\end{array}\right).$$

$$\psi=\left(\begin{array}[]{c}\phi\\ \chi^{*}\end{array}\right),\hskip 20.00003pt\bar{\psi}=(\phi^{\dagger}\hskip 10.00002pt-\chi^{\rm T})$$

$$\begin{array}[]{r@{\,}l}S_{0}=m_{0}&{\displaystyle\int}\left[\phi^{\dagger}(x)\phi(x)+\chi^{\dagger}(x)\chi(x)\right]\\[10.00002pt] +&{\displaystyle\int}\left[\phi^{\dagger}(x)D_{0}^{-}\phi(x)+\chi^{\dagger}(x)D_{0}^{-}\chi(x)\right]\\[10.00002pt] -\,\mbox{\small$\frac{1}{2}$}ar_{s}\zeta&{\displaystyle\int}\left[\phi^{\dagger}(x)\triangle^{(3)}\phi(x)+\chi^{\dagger}(x)\triangle^{(3)}\chi(x)\right]\\[10.00002pt] +\,\zeta&{\displaystyle\int}\left[\phi^{\dagger}(x)\mbox{\boldmath$\sigma$}\!\cdot\!\mbox{\boldmath$D$}\chi^{*}(x)-\chi^{\rm T}(x)\mbox{\boldmath$\sigma$}\!\cdot\!\mbox{\boldmath$D$}\phi(x)\right]\end{array}$$

$$\begin{array}[]{r@{\;}l}S_{B}&=-\mbox{\small$\frac{i}{2}$}ac_{B}\zeta{\displaystyle\int}\bar{\psi}(x)\mbox{\boldmath$\Sigma$}\!\cdot\!\mbox{\boldmath$B$}(x)\psi(x)\\[10.00002pt] &=-\mbox{\small$\frac{i}{2}$}ac_{B}\zeta{\displaystyle\int}\phi^{\dagger}(x)\mbox{\boldmath$\sigma$}\!\cdot\!\mbox{\boldmath$B$}(x)\phi(x)-\chi^{\dagger}(x)\mbox{\boldmath$\sigma$}^{*}\!\cdot\!\mbox{\boldmath$B$}(x)\chi(x),\end{array}$$

$$\begin{array}[]{r@{\;}l}S_{E}&=-\mbox{\small$\frac{1}{2}$}ac_{E}\zeta{\displaystyle\int}\bar{\psi}(x)\mbox{\boldmath$\alpha$}\!\cdot\!\mbox{\boldmath$E$}(x)\psi(x)\\[10.00002pt] &=-\mbox{\small$\frac{1}{2}$}ac_{E}\zeta{\displaystyle\int}\phi^{\dagger}(x)\mbox{\boldmath$\sigma$}\!\cdot\!\mbox{\boldmath$E$}(x)\chi^{*}(x)+\chi^{\rm T}(x)\mbox{\boldmath$\sigma$}\!\cdot\!\mbox{\boldmath$E$}(x)\phi(x),\end{array}$$

$$F_{\mu\nu}(x)=\frac{1}{8a^{2}}\sum_{\begin{array}[]{c}\scriptstyle\bar{\mu}=\pm\mu\\ \scriptstyle\bar{\nu}=\pm\nu\end{array}}\mathop{\rm sign}(\bar{\mu}\bar{\nu})\,T_{\bar{\mu}}T_{\bar{\nu}}T_{-\bar{\mu}}T_{-\bar{\nu}}-\mbox{h.c.},$$

$$\psi\mapsto e^{\varepsilon a(\hbox to0.0pt{$\,/$\hss}D+m)}\psi,\quad\bar{\psi}\mapsto\bar{\psi}e^{\bar{\varepsilon}a(\hbox to0.0pt{$\,/$\hss}D+m)},$$

$$\begin{array}[]{r@{\,\mapsto\,}l}\psi&\exp\Big{(}\varepsilon a(\hbox to0.0pt{$\,/$\hss}D+m)+\delta a\mbox{\boldmath$\gamma$}\!\cdot\!\mbox{\boldmath$D$}\Big{)}\psi,\\[6.99997pt] \bar{\psi}&\bar{\psi}\exp\Big{(}\bar{\varepsilon}a(\hbox to0.0pt{$\,/$\hss}D+m)+\bar{\delta}a\mbox{\boldmath$\gamma$}\!\cdot\!\mbox{\boldmath$D$}\Big{)}.\end{array}$$

$$\left\langle\psi(t^{\prime},\mbox{\boldmath$p$}^{\prime})\bar{\psi}(t,\mbox{\boldmath$p$})\right\rangle=(2\pi)^{3}\delta(\mbox{\boldmath$p$}^{\prime}-\mbox{\boldmath$p$})C(t^{\prime}-t,\mbox{\boldmath$p$}),$$

$$C(t,\mbox{\boldmath$p$})=\int^{\pi}_{-\pi}\frac{dp_{0}}{2\pi}\;\frac{e^{ip_{0}t}}{i\gamma_{0}\sin p_{0}+i\zeta\mbox{\boldmath$\gamma$}\!\cdot\!\mbox{\boldmath$S$}+m_{0}+1-\cos p_{0}+\mbox{\small$\frac{1}{2}$}r_{s}\zeta\hat{\mbox{\boldmath$p$}}^{2}},$$

$$C(t,\mbox{\boldmath$p$})={\cal Z}_{2}e^{-E|t|}\frac{\gamma_{0}\mathop{\rm sign}t\sinh E-i\zeta\mbox{\boldmath$\gamma$}\!\cdot\!\mbox{\boldmath$S$}+m_{0}+1-\cosh E+\mbox{\small$\frac{1}{2}$}r_{s}\zeta\hat{\mbox{\boldmath$p$}}^{2}}{2\sinh E}$$

$$\cosh Ea=1+\frac{(m_{0}a+\mbox{\small$\frac{1}{2}$}r_{s}\zeta\hat{\mbox{\boldmath$p$}}^{2}a^{2})^{2}+\zeta^{2}\mbox{\boldmath$S$}^{2}a^{2}}{2(1+m_{0}a+\mbox{\small$\frac{1}{2}$}r_{s}\zeta\hat{\mbox{\boldmath$p$}}^{2}a^{2})}$$

$$E^{2}=M_{1}^{2}+\frac{M_{1}}{M_{2}}\mbox{\boldmath$p$}^{2}+\ldots,$$

$$M_1=E(0),$$

$$M_{2}^{-1}=(\partial^{2}E/\partial p_{i}^{2})_{\hbox to0.0pt{$\hskip 0.35002pt\scriptstyle p$\hss}\scriptstyle p\hskip 0.35002pt={\bf 0}}.$$

$$M_1=a^{-1}\log (1+m_{0}a),$$

Full text

1 Introduction

The most promising avenue for a quantitative understanding of nonperturbative quantum chromodynamics—and other field theories—is via numerical (Monte Carlo) integration of functional integrals defined on a lattice [1]. Like any numerical technique this method has uncertainties that must be understood and controlled before the results are useful. In particular, although the continuum theory is defined by the limit of a sequence of lattice theories, the numerical calculations are never carried out at the limit. Because the Monte Carlo introduces statistical errors, the extrapolation to the continuum limit is imperfect. The results for physical quantities are consequently contaminated by lattice artifacts. For a practical result, this uncertainty must be smaller than, say, relevant experimental uncertainties.

The way to reduce lattice artifacts is based on the renormalization group [2]. One starts with a general action

[TABLE]

where the $S_{n}$ include all interactions with the desired field content and the appropriate symmetries. One approach to the continuum limit, which might be called brute force, is to choose the $c_{n}$ in any way that drives the lattice spacing to zero. An ideal approach would be to choose the $c_{n}$ to lie on a renormalized trajectory [2], where there are no lattice artifacts even though the lattice spacing ${a\neq 0}$. In the space of all possible actions specified by eq. (1.1), the renormalized trajectories lie in a subspace, whose dimension equals the number of relevant parameters. Once the relevant parameters have been fixed by physics, they and the renormalization scheme determine all the $c_{n}$.

Unfortunately, a renormalized trajectory is mostly of abstract value, because on one infinitely many $c_{n}$ are nonzero. All practical schemes, such as blocking fields [3] or Symanzik improvement [4] use criteria such as locality [3] or the scaling dimension [4] to truncate the space of actions. (For an asymptotically free theory, such as QCD, these two criteria are not very different.) Furthermore, the calculations of the $c_{n}$ are, in practice, only approximate. For these reasons an improved action is only partially renormalized. Nevertheless, any practical action can be written

[TABLE]

where ${\cal S}_{\rm RT}$ denotes (an action on) the renormalized trajectory. Usually the truncations and/or approximations used to generate $S$ will also yield estimates for the remaining cutoff effects $\delta S$.

This paper treats massive fermions coupled to non-Abelian gauge fields. The relevant couplings are the fermion masses and the gauge coupling. So the renormalized trajectory takes the form

[TABLE]

where $m_{q}$ denotes the fermion mass,111We use $m_{q}$ for a quark mass defined by a physical condition and $m_{0}a$ for the coupling appearing in the action. $\Lambda_{\rm QCD}$ is the scale characteristic of the gauge theory, and the argument $\Lambda_{\rm QCD}a$ labels the renormalization point. The $S_{n}$ are gauge-invariant combinations of four-component fermion and anti-fermion fields ($\psi$ and $\bar{\psi}$) and the lattice gauge field ($U_{\mu}$). For later calculational convenience we choose the bare, rather than some physical, fermion mass $m_{0}a$ and gauge coupling $g_{0}^{2}$ to parameterize the couplings $c_{n}$.

As $a\to 0$ the fermion mass $m_{q}a$ is formally smaller than $g_{0}^{2}$. (By asymptotic freedom $g_{0}^{2}=[\beta_{0}\log(\Lambda_{\rm QCD}^{2}a^{2})]^{-1}$ as $a\to 0$.) It is therefore tempting to expand the couplings $c_{n}(m_{0}a,g_{0}^{2})$ in $m_{0}a$, as in previous analyses [5, 6, 7]. But there may be fermions satisfying $m_{q}/\Lambda_{\rm QCD}\gg 1$; the charm, bottom, and top quarks are examples in nature. If, in practice, $m_{q}a$ is not small, perturbation theory in $m_{0}a$ need not be useful, even though perturbation theory in $g_{0}^{2}$ might be. Indeed, this regime includes the charm and bottom quarks at currently accessible lattice spacings.

The static [8, 9] and nonrelativistic [10, 11, 12] effective theories address the problems of heavy fermions. Their restriction to $m_{q}\gg\Lambda_{\rm QCD}$ implies that couplings of interactions between particle and antiparticle states may be chosen to vanish, and the remaining interactions in eq. (1.1) are organized according to a \mbox{\boldmathp}/m_{q} expansion. But for some $m_{q}\sim 2\Lambda_{\rm QCD}$ the expansion is no longer useful. Furthermore, radiative corrections induce power-law terms, e.g. $g_{0}^{2}/(m_{q}a)$, which must be canceled by adjusting the $c_{n}$. These terms, which diverge as $a\to 0$, are a reminder that the effective theories are to be used at scales below (large) $m_{q}$. Their presence implies that cutoff effects in the effective theories should be removed not by brute force, but by keeping $a\sim m_{q}^{-1}$ and constructing actions systematically closer to the renormalized trajectory [11, 12].

This paper presents a way to encompass both the small and large mass formulations. The doubling problem is handled by Wilson’s method [13]. In addition to treating the $m_{0}a$ dependence of the couplings exactly, the central idea is to enlarge the class of interactions considered in eq. (1.3) to include those from both the small $m_{0}a$ and the large $m_{q}/\Lambda_{\rm QCD}$ limits. In particular, we do not impose a symmetry between couplings of interactions that would be related by interchanging the time axis with a spatial axis. For any $m_{q}a$, the four-momentum of a quark in most interesting physical systems satisfies [E(\mbox{\boldmathp})-E({\bf 0})]a\ll 1 whenever \mbox{\boldmathp}a\ll 1. But when $m_{q}a>1$ the characteristic four-momentum of the physics does not respect time-space axis interchange. Under such circumstances it is inconvenient and unnecessary to choose an axis-interchange symmetric action.222Relinquishing axis-interchange symmetry is common in treatments of heavy quarks with momentum-space and dimensional regulators. It is possible to derive deviations from heavy-quark symmetry from the Dirac action [14, 15], while maintaining time-space axis interchange invariance as a corollary to Euclidean invariance. But usually the derivations are easier with a nonrelativistic action [10, 16], the so-called the heavy-quark effective theory [17].

In the appropriate limits, our formulation of lattice fermions shares properties of previous ones. On one hand, at dimension five or less, couplings related by axis interchange become identical in the limit $m_{0}a\to 0$: the Wilson action and the improved action of ref. [6] are recovered. But when $m_{0}a\neq 0$ the mass-dependent renormalization leaves lattice artifacts that are proportional to \mbox{\boldmathp}a and $\Lambda_{\rm QCD}a$, not $m_{q}a$. At higher dimension, however, we retain Wilson’s time derivative and incorporate “spatial-only” interactions into eq. (1.1).

On the other hand, for ${m_{q}\gg\Lambda_{\rm QCD}}$ one can interpret the lattice theory in a nonrelativistic light. Indeed, all members of our class of actions approach a universal static limit as $m_{0}a\to\infty$. For $m_{0}a$ large but finite, the $1/m_{q}^{s}$ corrections to the static limit can be recovered systematically, provided the fermion mass is defined through the kinetic energy, and provided the general action, eq. (1.3), is truncated only at dimension $4+s$ (or higher). Unlike previous implementations of nonrelativistic fermions, however, our approach crosses smoothly over into the regime of tiny lattice spacings, where $m_{0}a\ll 1$ even for a heavy quark. Thus, after several $c_{n}$ have been tuned close to a renormalized trajectory, thereby removing the worst lattice artifacts, a little brute force can remove the rest.

Because we make no assumptions about the ratio of fermion masses to other scales, our formulation is especially well suited to fermions too heavy for small $m_{0}a$ methods yet too light for large $m_{q}/\Lambda_{\rm QCD}$ methods. With the actions given below one can test whether a given fermion is heavy enough to be treated nonrelativistically, without resorting to brute-force simulations. A practical example might be the charm quark, which has a mass only a few times $\Lambda_{\rm QCD}$, yet even on the finest lattices available today $m_{\rm ch}a$ is largish, at least $\frac{1}{3}$.

For a concrete determination of the $c_{n}$, one must choose a renormalization group, a criterion for truncating the sum in eq. (1.3), and a strategy for determining the $c_{n}$. For illustration we adopt here a Symanzik-like procedure [4], organizing the interactions by dimension. Carrying out this program to arbitrarily high dimension would produce a renormalized trajectory of a renormalization group generated by infinitesimal changes in $a$. For simplicity, however, most of this paper treats interactions only up to dimension five. Although a nonperturbative determination of the couplings is possible in principle, this paper makes the further approximation of perturbation theory in $g_{0}^{2}$, that is

[TABLE]

Except for sect. 8 we work to tree level, so we often abbreviate $c_{n}^{[0]}(m_{0}a)$ by $c_{n}(m_{0}a)$. (Explicit one-loop calculations are in progress [18].) Within these approximations we determine the $c_{n}$ by insisting that on-shell quantities take their desired values, as first suggested by Lüscher and Weisz [19].

One calculational procedure is to work out $n$-point on-shell Green functions via Feynman diagrams and tune them to the continuum limit, to the appropriate order in $p$. An example is in sect. 4. Because this strategy is limited to a finite number of quantities, it is nontrivial to assume that other quantities are improved too [4]. An alternative is developed in sect. 5. Starting from the transfer matrix we derive an expression for the fermion Hamiltonian, valid (at tree level) for states with \mbox{\boldmathp}a\ll 1. Because the Hamiltonian is an operator, improving it to some accuracy guarantees the improvement to the same accuracy of infinitely many $c$-numbers. We show that by adjusting the couplings $c_{n}(m_{0}a)$ correctly, and allowing physically unobservable redefinitions of the fermion field, one can tune the Hamiltonian to the continuum limit, i.e. to the Dirac Hamiltonian, to the appropriate order in $p$.

In addition to equations for the $c_{n}$, the analysis of the Hamiltonian yields two important results. One is that the Wilson time derivative needs no improvement. The other is a canonically normalized fermion field that, to the accuracy of the improvement, obeys Dirac’s (continuum) equation of motion. This field is a potent ingredient in calculations of matrix elements involving heavy-quark systems; see sect. 7.

Owing to the approximations introduced—the truncation of interactions and perturbation theory—some cutoff effects remain. If these errors are small, they may be estimated by insertions of $\delta S$ in correlation functions. If the series for $c_{n}(m_{0}a,g_{0}^{2})$ has been developed to $L_{n}$-th order

[TABLE]

A typical term in $\delta S$ distorts an on-shell correlation function by an amount of order g_{0}^{2l}c_{n}^{[l]}(m_{0}a)(\mbox{\boldmathp}a)^{s_{n}}, where ${s_{n}}=\dim S_{n}-4$, and $l\geq L_{n}+1$. (If $S_{n}$ is omitted from the action altogether, then here $l=0$.) The analysis presented below shows that the tree-level, lower-dimension $c_{n}^{[0]}(m_{0}a)$ are well-behaved for all masses. We also show that loop diagrams have the same or more benign behavior at large mass. In particular, as $m_{0}a\to\infty$ the $c_{n}^{[l]}(m_{0}a)$ either approach a constant or fall as $(m_{q}a)^{-s}$, for some $s\leq s_{n}$. These results provide evidence that the higher-dimension $c_{n}(m_{0}a)$ are well-behaved too.

In Monte Carlo programs it is customary to parameterize the action by the hopping parameter instead of the mass. In this notation the dimension-three and -four interactions are written

[TABLE]

Some terms of dimension five—to solve the doubling problem [13]—are included here too. The relation between the bare mass $m_{0}a$ and the hopping parameters $\kappa_{t}$ and $\kappa_{s}$ is given below in sect. 2. Eq. (1.6) illustrates how our program subsumes properties of the familiar small-mass and large-mass formulations. Imposing axis-interchange invariance would set $\kappa_{s}=\kappa_{t}$ and $r_{s}=1$, and then $S_{0}$ reduces to the Wilson action [13]. Rewriting $r_{s}\kappa_{s}=c_{s}$ and setting $\kappa_{s}$ to zero with $c_{s}\neq 0$ produces the simplest nonrelativistic action [11].

This pattern continues for dimension-five interactions. Aside from the Wilson terms in $S_{0}$, the other dimension-five interactions are the chromomagnetic interaction

[TABLE]

and the chromoelectric interaction

[TABLE]

where $B$ and $E$ are suitable functions of the lattice gauge field $U$, as in sect. 2. The light-quark formalism of refs. [6, 7] considers the special case $c_{B}=c_{E}$, whereas the heavy-quark formalism of refs. [11, 12] sets $c_{E}=0$.

The couplings $r_{s}$, $\zeta=\kappa_{s}/\kappa_{t}$, $c_{B}$, and $c_{E}$ are specific examples of the couplings $c_{n}$ in eq. (1.3). On the renormalized trajectory they are, therefore, all functions of $m_{0}a$. Sect. 4 shows how to adjust $r_{s}$ and $\zeta$ so that the relativistic energy-momentum relation E=\sqrt{m_{q}^{2}+\mbox{\boldmathp}^{2}}+\delta E_{\rm lat} is obtained for all $m_{q}a$. With the correct choice, for which $\zeta\neq 1$ when $m_{0}a\neq 0$, the (tree-level) lattice artifact $\delta E_{\rm lat}$ is proportional to \mbox{\boldmathp}^{4}a^{3} for $m_{q}a\ll 1$ and \mbox{\boldmathp}^{4}a/m_{q}^{2} for $m_{q}a\gg 1$. Similarly, results in sect. 5 include functions $c_{B}(m_{0}a)$ and $c_{E}(m_{0}a)$ that reduce lattice artifacts in the quark-gluon vertex functions to {\rm O}(\mbox{\boldmathp}^{2}a^{2}) for $m_{q}a\ll 1$, or (yet smaller) {\rm O}(\mbox{\boldmathp}^{2}a/m_{q}) for $m_{q}a\gg 1$.

In their on-shell improvement program refs. [19, 6] introduced changes of variables, or isospectral transformations, to expose redundant interactions. Since the coupling of a redundant interaction does not influence physical quantities, one can choose it according to theoretical or computational criteria distinct from improvement. Sect. 3 examines the isospectral transformations when time-space axis-interchange symmetry is not imposed. In our formulation many isospectral transformations are exploited to keep the time discretization, and hence the transfer matrix, as simple as possible.

The remaining redundant directions can be classified in the Hamiltonian approach developed in sect. 5. In a Euclidean version of the standard Dirac-matrix basis, given in sect. 2, matrices are either block diagonal or block off-diagonal. Block-diagonal transformations are absorbed into a generalized field normalization. On the other hand, block–off-diagonal transformations (called Foldy-Wouthuysen-Tani transformations [20] in atomic physics) generate leeway in choosing the mass dependence of associated couplings.333The off-diagonal interactions are precisely the ones usually omitted from nonrelativistic formulations, yet their presence in our formulation permits a smooth transition from large to small mass. For example, in the action $S_{0}+S_{B}+S_{E}$ one may freely choose $r_{s}$, as long as the choice circumvents the doubling problem.

Our approach breaks down, just as any lattice theory does, when \mbox{\boldmathp}a is too large. Fortunately, the typical momenta and mass splittings of hadronic systems usually are bounded; the energy scale around the fermion mass is dynamically unimportant. In quarkonia the typical energy-momentum scales are $m_{q}v$ and $m_{q}v^{2}/2$, i.e. 200–800 MeV for charmonium and 200–1400 MeV for bottomonium. Similarly, in light and heavy-light hadrons the typical momentum scale is $\Lambda_{\rm QCD}$, i.e. 100–300 MeV. In some processes, such as a decaying heavy-quark system that transfers all its energy into light hadrons, a large three-momentum |\mbox{\boldmathp}|\approx m_{q} does arise. Then our formulation and its predecessors all require further extensions. One should appreciate, however, that the breakdown arises not from the large fermion mass per se, but from the large momentum of the decay products.

A by-product of our formalism applies to existing numerical calculations, done with axis-interchange invariant actions. For m_{q}a\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{>}}1 (and, furthermore, ${m_{q}/\Lambda_{\rm QCD}\gg 1}$) we derive in sect. 9 a nonrelativistic interpretation of such actions. One then sees that, with a proper definition of the fermion mass, any action described by $S_{0}+S_{B}+S_{E}$, including the Wilson and Sheikholeslami-Wohlert fermion actions, has the lattice-spacing and/or relativistic inaccuracies of a typical nonrelativistic action. A practical bonus of the nonrelativistic regime is that it is no longer necessary to adjust $\kappa_{s}$ differently from $\kappa_{t}$. In heavy-light systems, one may also set $c_{B}=c_{E}$.

This paper is organized as follows: Sect. 2 introduces some notation, including a form of the action better suited to perturbation theory, and the Dirac-matrix convention used in later sections. The isospectral transformation of ref. [6] is reviewed and generalized in sect. 3, to determine which couplings are redundant. Then, to derive improvement conditions, Feynman-diagram methods are discussed in sect. 4, and the Hamiltonian method in sect. 5. With a Hamiltonian description of the lattice theory in hand, sect. 6 estimates cutoff effects in various hadronic systems. Sect. 7 considers perturbations from the electroweak interactions, needed for the phenomenology of the Standard Model [21]. Some of the issues beyond tree level are outlined in sect. 8. The relationship of our work, for $m_{q}/\Lambda_{\rm QCD}\gg 1$, to nonrelativistic QCD is pursued in sect. 9. We discuss a few phenomenologically relevant applications more thoroughly in sect. 10. Finally, sect. 11 contains selected concluding remarks, and the appendices contain various technical details.

2 Notation

We shall call the form of the action in eq. (1.6) the “hopping-parameter form.” For studying the continuum limit and developing perturbation theory it is useful to present a different form. Let us introduce some notation. The lattice spacing is $a$ and the site labels are $n=x/a$. Rescale the fields:

[TABLE]

and similarly for $\bar{\psi}_{n}$. The bare mass is

[TABLE]

where $d\;(=4)$ is the spacetime dimension, and $\zeta=\kappa_{s}/\kappa_{t}$. With these substitutions the action reads

[TABLE]

where the integral sign abbreviates $a^{4}\sum_{x}$. The covariant difference operators are conveniently defined via covariant translation operators

[TABLE]

where $U_{-\mu}(x)=U_{\mu}^{\dagger}(x-a\hat{\mu})$. Then

[TABLE]

define various covariant difference operators and the three-dimensional discrete Laplacian. We shall call the form of the action in eq. (2.3) the “mass form.”

The temporal kinetic energy in eq. (2.3) is written in a way that does not make the temporal Wilson term explicit. Eq. (2.3) is more convenient, however, for constructing the transfer matrix, and for comparing with nonrelativistic QCD. The spacelike Wilson term, the one proportional to $r_{s}$, is needed to prevent doubling. A convenient choice in computer programs is $r_{s}=1$, but we keep it arbitrary, because other choices may have other advantages.

For constructing the transfer matrix and for examining the nonrelativistic limit, a useful representation of the Euclidean gamma matrices is

[TABLE]

satisfying $\{\gamma_{\mu},\gamma_{\nu}\}=2\delta_{\mu\nu}$. Another convention that we use is $\sigma_{\mu\nu}=\frac{i}{2}[\gamma_{\mu},\gamma_{\nu}]$ so that $\sigma_{\mu\nu}^{\dagger}=+\sigma_{\mu\nu}$. Using eq. (2.6), $\sigma_{0i}=i\alpha_{i}$ and $\sigma_{ij}=-\varepsilon_{ijk}\Sigma_{k}$, where

[TABLE]

The following split into upper and lower two-component spinors

[TABLE]

follows from eq. (2.6). This convention is chosen so that (the operators corresponding to) $\phi$ and $\chi$ annihilate particle and anti-particle states, respectively. With these two-component fields the mass form of the action is

[TABLE]

This form of the action exhibits explicitly that particles and anti-particles are treated on the same footing. (Anti-particles transform under the complex-conjugate representation of the gauge group, however, so $U^{*}$ appears instead of $U$ in the rules (eqs. (2.4)) for covariant translations acting on $\chi$.)

Writing \mbox{\boldmathB}_{n}=a^{2}\mbox{\boldmathB}(x) and \mbox{\boldmathE}_{n}=a^{2}\mbox{\boldmathE}(x), the four- and two-component mass forms of the chromomagnetic and chromoelectric interactions are

[TABLE]

and

[TABLE]

respectively. Except in a technical step in sect. 5 we take the “four-leaf clover” lattice approximant to the field strength

[TABLE]

as introduced in ref. [22]. In eqs. (2.10) and (2.11), B_{i}=\mbox{\small\frac{1}{2}}\varepsilon_{ijk}F_{jk} and $E_{i}=F_{0i}$. As defined here, $F_{\mu\nu}$, $B_{i}$, and $E_{i}$ are anti-Hermitian; similarly we take anti-Hermitian gauge-group generators ${t^{a}}^{\dagger}=-t^{a}$.

3 Redundant Couplings

Before trying to determine the mass dependence of the couplings $\zeta$, $r_{s}$, $c_{B}$, and $c_{E}$, one should establish which combinations are physical. The fields in functional integrals are integration variables, and a change of variables cannot change the integrals. Interactions that are induced by changes of variables are redundant; their couplings can be chosen with some leeway, dictated by calculational or technical convenience, rather than by physical criteria.

A subtle example of a redundancy in the space of interactions is wavefunction (re)normalization, which multiplies the field by a constant. For fermions, for example, it is sometimes convenient to fix the kinetic energy \bar{\psi}\hbox to0.0pt{,/\hss}D\psi to have coefficient unity, which is the mass form of the action, and sometimes to fix the local term $\bar{\psi}\psi$ to have coefficient unity, which is the hopping-parameter form. But neither interaction is redundant, even though the normalization convention drops out of physical quantities.

Otherwise redundant directions are exposed by redefinitions of the field. In the analysis of ref. [6], with axis-interchange symmetry, one considers the transformation

[TABLE]

where $a$ is chosen so that a(\hbox to0.0pt{,/\hss}D+m) is “small.” After carrying out the transformation on the target action \int\bar{\psi}(\hbox to0.0pt{,/\hss}D+m)\psi, one expands the transformed action to ${\rm O}(a)$. One finds changes in the normalization of the lower-dimension terms and the additional interaction a(\varepsilon+\bar{\varepsilon})\bar{\psi}\hbox to0.0pt{,/\hss}D^{2}\psi: from the two independent transformation parameters, only one combination survives. Hence, of the dimension-five interactions listed in Table 1,

one (i.e. \bar{\psi}\hbox to0.0pt{,/\hss}D^{2}\psi) is redundant, and the other is not.

On the lattice the nearest-neighbor discretization of \hbox to0.0pt{,/\hss}D suffers from the doubling problem. Wilson’s prescription adds a nearest-neighbor discretization of $D^{2}$ to eliminate the unwanted states. By the preceding analysis [6], using instead \hbox to0.0pt{,/\hss}D^{2}=D^{2}-i\sigma_{\mu\nu}F_{\mu\nu} would not change the spectrum at ${\rm O}(a)$. When the discretization is chosen to solve the doubling problem, however, the $D^{2}$ interaction does change the spectrum of high-momentum states. Since they communicate with the low-momentum states through virtual processes, lattice artifacts proportional to $g_{0}^{2}$ remain. They can be eliminated with the other dimension-five interaction, $i\bar{\psi}\sigma_{\mu\nu}F_{\mu\nu}\psi$, with a coupling proportional to $g_{0}^{2}$.

Thus, with axis-interchange symmetry there are four interactions up to dimension five, one of which goes with wavefunction normalization (e.g. \bar{\psi}\hbox to0.0pt{,/\hss}D\psi). One coupling is redundant, and it can be chosen to solve species doubling (\bar{\psi}\hbox to0.0pt{,/\hss}D^{2}\psi). The other couplings are fixed by the fermion mass ($\bar{\psi}\psi$) and a physical improvement condition ($i\bar{\psi}\sigma_{\mu\nu}F_{\mu\nu}\psi$).

When axis-interchange symmetry is given up, the transformation in eq. (3.1) should be generalized to

[TABLE]

Applying this transformation to the target action induces the dimension-five interactions listed in Table 1. From the four independent transformation parameters, only three combinations survive: ${\varepsilon+\bar{\varepsilon}}$, ${\delta+\bar{\delta}}$, and ${\delta-\bar{\delta}}$. Therefore, the coefficients of $\bar{\psi}D_{0}^{2}\psi$, \bar{\psi}(\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD})^{2}\psi, and \bar{\psi}[\gamma_{0}D_{0},\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}]\psi can be chosen arbitrarily. The last of these has no redeeming features, so ${\delta-\bar{\delta}}$ should be chosen so that it never appears.

The other two redundant interactions are again used to solve the doubling problem. The $D_{0}^{2}$ term is used to eliminate states that would make contributions to the fermion propagator proportional to $(-1)^{t}$; the factors ${1\pm\gamma_{0}}$ in eqs. (1.6) or (2.3) provide the unique choice. Low-energy states with $p_{i}\sim\pi/a$ are lifted by adding the interaction proportional to $r_{s}$ in eqs. (1.6) or (2.3). When the mass is nonzero, it may prove convenient choose $r_{s}$ to be a function of $m_{0}a$, so we leave it arbitrary.

As with axis-interchange symmetry, the chromomagnetic and chromoelectric interactions are not redundant. Their couplings can be used to remove cutoff effects once the doubling problem has been eliminated.

Thus, without axis-interchange symmetry there are eight interactions up to dimension five, one of which goes with wavefunction normalization (e.g. $\bar{\psi}\gamma_{0}D_{0}\psi$). Three couplings are redundant; two can be used to solve species doubling ($\bar{\psi}D_{0}^{2}\psi$ and \bar{\psi}(\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD})^{2}\psi), and the other to eliminate \bar{\psi}[\gamma_{0}D_{0},\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}]\psi. The other couplings are fixed by the fermion mass ($m_{0}\bar{\psi}\psi$) and three physical improvement conditions (\zeta\bar{\psi}\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}\psi, {ic_{B}\bar{\psi}\mbox{\boldmath\Sigma}\!\cdot\!\mbox{\boldmathB}\psi}, and c_{E}\bar{\psi}\mbox{\boldmath\alpha}\!\cdot\!\mbox{\boldmathE}\psi).

Redundant combinations of higher-dimension interactions can be exposed by generalizing the transformation of eq. (3.2). In particular, after dispensing with axis-interchange symmetry, it is possible to transform away interactions with higher time derivatives of $\psi$ and $\bar{\psi}$, in favor of spatial derivatives of $\psi$ and $\bar{\psi}$, $B$ and $E$, and time derivatives of the latter. Indeed, any action with the Wilson time difference—the first line of eq. (2.3)—has an easy-to-construct transfer matrix. (This is reviewed in sect. 5.) Consequently, it is possible to implement eq. (1.3) by adding “spatial-only” interactions to $S_{0}$.

4 On-shell Correlation Functions

We now turn to the mass dependence of the tree-level couplings, generically denoted $c_{n}^{[0]}(m_{0}a)$ in eq. (1.4), needed to bring the action closer to the renormalized trajectory. This section uses the fermion propagator to obtain the relation between the physical mass and the coupling $m_{0}a$, the correct tuning of the coupling $\zeta(m_{0}a)$, and the normalization of the field $\psi(x)$. Since we are interested in the full mass dependence, we do not expand in the fermion mass. Sect. 5 uses the Hamiltonian of the lattice theory to clarify and extend the analysis to $c_{B}$ and $c_{E}$.

A well-known procedure for determining the couplings [4] is to calculate $n$-point correlation functions and expand in momentum $p$. In gauge theories, however, it is not known whether lattice artifacts can be removed systematically from Green functions off mass shell. Hence, one expands “on-shell” quantities instead [19]. The (lattice) mass shell specifies the energy at given spatial momentum $p$, so on-shell improvement amounts to an expansion in \mbox{\boldmathp}a. Previous analyses [5, 6, 7] also expanded in the coupling $m_{0}a$. We simply avoid the latter expansion, and thus obtain the full mass dependence.

The simplest on-shell correlation function is the fermion propagator as a function of time and spatial momentum. It is used to relate the bare mass to a physical mass and to derive the mass dependence of $\zeta$. In the language of sect. 3, it probes the interactions $\bar{\psi}\psi$ and \bar{\psi}\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}\psi, relative to $\bar{\psi}\gamma_{0}D_{0}\psi$.

Define C(t,\mbox{\boldmathp}) through

[TABLE]

where \psi(t,\mbox{\boldmathp})=a^{3}\sum_{\hbox to0.0pt{\hskip 0.35002pt\scriptstyle x\hss}\scriptstyle x\hskip 0.35002pt}e^{-i\hbox to0.0pt{\hskip 0.35002pt\scriptstyle p\hss}\scriptstyle p\hskip 0.35002pt\cdot\hbox to0.0pt{\hskip 0.35002pt\scriptstyle x\hss}\scriptstyle x\hskip 0.35002pt}\psi(t,\mbox{\boldmathx}) and similarly for \bar{\psi}(t,\mbox{\boldmathp}). Then from eq. (2.3)

[TABLE]

where $S_{i}=a^{-1}\sin p_{i}a$ and $\hat{p}_{i}=2a^{-1}\sin(p_{i}a/2)$, but for brevity eq. (4.2) is given in lattice units. To integrate over $p_{0}$, proceed as follows: rationalize the denominator; for $t\geq 0$ let $z=e^{ip_{0}}$, and for $t<0$ let $z=e^{-ip_{0}}$, yielding a contour integral over the circle $|z|=1$; apply the residue theorem to obtain

[TABLE]

for $t\neq 0$,444To obtain $C(0)$ from eq. (4.3), replace $\gamma_{0}\mathop{\rm sign}t$ by 1 on the right-hand side. where (restoring $a$)

[TABLE]

implicitly defines the energy of a state with momentum $p$. The residue {\cal Z}_{2}(\mbox{\boldmathp}) is given below in eq. (4.12).

Expanding the energy-momentum relation in powers of \mbox{\boldmathp}a yields

[TABLE]

where the “rest mass”

[TABLE]

and the “kinetic mass”

[TABLE]

(Any axis $i$ will do to define $M_{2}$, by spatial axis-interchange symmetry.) The relativistic mass shell has $m_{q}=M_{1}=M_{2}$, and it terminates at \mbox{\boldmathp}^{2}. From the tree-level eq. (4.4)

[TABLE]

and

[TABLE]

Eq. (4.8) shows how to adjust $m_{0}a$ so that $m_{q}=M_{1}$. Similarly, eq. (4.9) shows how to adjust $\zeta$ and $r_{s}$ so that $m_{q}=M_{2}$.555In the $\kappa_{t}$-$\kappa_{s}$ parametrization (eliminate $m_{0}$ with eq. (2.2)) this condition is an implicit transcendental equation. In the $m_{0}$-$\zeta$ parametrization one can solve for $\zeta$ explicitly. Setting $M_{1}=M_{2}$ and solving for $\zeta$ yields the (tree-level) condition (setting $a=1$ again)

[TABLE]

The dimension-five coupling $r_{s}$ is treated here as redundant; it is determined not by physics, but to solve the doubling problem. To alleviate doubling any ${r_{s}(0)>0}$ will do, and the most natural choice is ${r_{s}(0)=1}$.

For small mass the Taylor expansion of eq. (4.10) is

[TABLE]

At $m_{0}=0$ the redundant coupling $r_{s}$ drops out, leaving $\zeta(0)=1$ unambiguously. On the other hand, the full mass dependence of $\zeta$ can only be prescribed hand-in-hand with $r_{s}$. The origin of the link between the two couplings is that both the kinetic (\bar{\psi}\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}\psi) and Wilson ($\bar{\psi}\triangle^{(3)}\psi$) terms contribute to $E^{2}$ at {\rm O}(\mbox{\boldmathp}^{2}). This and analogous links between couplings’ mass dependence are examined further in sect. 5 and Appendix A.

Beyond tree level (in perturbation theory or in Monte Carlo calculations) one would tune $\zeta$ according to the same physical principle that led to eq. (4.10): determine the momentum dependence of the energy of a suitable state and demand that $M_{1}=E({\bf 0})$ and M_{2}=(\partial^{2}E/\partial p_{i}^{2})^{-1}_{\hbox to0.0pt{\hskip 0.35002pt\scriptstyle p\hss}\scriptstyle p\hskip 0.35002pt={\bf 0}} be equal.

When $\zeta$ and $m_{0}a$ have been adjusted so that $M_{1}=M_{2}=m_{q}$, one can rewrite eq. (4.5) as E=\sqrt{m_{q}^{2}+\mbox{\boldmathp}^{2}}+\delta E_{\rm lat}. Expanding eq. (4.4) to \mbox{\boldmathp}^{4}, one finds the lattice artifact \delta E_{\rm lat}\sim\mbox{\boldmathp}^{4}a^{3} at small mass and \delta E_{\rm lat}\sim\mbox{\boldmathp}^{4}a/M_{2}^{2} at large mass. To reduce $\delta E_{\rm lat}$ further, one must incorporate higher-dimension interactions into the analysis.

Finally, let us return to the residue ${\cal Z}_{2}$ in eq. (4.3). In general, the residue is a scalar function of four-momentum $p$, evaluated on shell. With a Euclidean invariant cutoff, scalar functions can depend only on $p^{2}$; on shell, with $p^{2}=-m^{2}$, the spatial momentum $p$ drops out. With the lattice cutoff, however, the mass shell is distorted, cf. eq. (4.4), so three-momentum $p$ dependence can remain. Indeed, after integrating eq. (4.2) over $p_{0}$ one finds

[TABLE]

Normally one identifies the residue with a (re)normalization of the fermion field. Now, however, it is appropriate to expand {\cal Z}_{2}(\mbox{\boldmathp})=Z_{2}+{\rm O}(\mbox{\boldmathp}^{2}), where

[TABLE]

Then $Z_{2}^{-1/2}\psi(x)=e^{M_{1}a/2}\psi(x)$ has the canonical normalization. In the hopping-parameter notation the canonically normalized field is $\sqrt{1-6r_{s}\kappa_{s}}\,\psi_{n}$. This notation shows clearly that the approach to the static limit, $\kappa_{s}\ll 1$, smooth. Indeed, eq. (4.13) captures the dominant mass dependence of the field normalization to all orders in perturbation theory, cf. sect. 8 and ref. [18].

One might ask what to make of the momentum dependence of ${\cal Z}_{2}$, when the action is improved to higher dimensions. The residue itself is not observable; physical quantities are given by ratios of $n$-point functions and the propagator. With the correct on-shell improvement, the \mbox{\boldmathp}^{2} dependence of untruncated $n$-point functions combines with that of ${\cal Z}_{2}$ to yield the desired results (to the order considered).

5 The Hamiltonian

This section introduces another method for deriving conditions on the couplings in the action. The strategy is to obtain an expansion in the lattice spacing for the Hamiltonian. For concreteness, we focus on the action $S=S_{0}+S_{B}+S_{E}$. The couplings are then adjusted so that the Hamiltonian of the lattice theory is equivalent to the Dirac Hamiltonian. The idea is conceptually the same as on-shell improvement, because the “spectral quantities” of ref. [19] are just eigenvalues of the Hamiltonian. But since the Hamiltonian is an operator, it contains the information of infinitely many quantities, rather than the finite number accessible when one computes correlation functions.

This approach reproduces the $c_{n}(m_{0}a)$ derived with on-shell correlation functions. But the analysis is explicitly relativistic, if noncovariant, so one sees clearly that the results are general. On the other hand, we have not attempted to extend the method to four-fermion operators, or to higher orders in $g_{0}^{2}$. The calculations required by those extensions seem simpler with Feynman diagrams.

There is a further conceptual advantage to the Hamiltonian. Lattice field theories are almost always formulated in imaginary time. The interpretation of the results in real time hinges on a good Hamiltonian fixing the dynamics of the Hilbert space of states [2]. Hence the implicit, but seldom stated, goal of improvement is an improved Hamiltonian; this section merely takes direct aim on that goal. Moreover, once one accepts the central role of the Hamiltonian, one appreciates why a satisfactory Hamiltonian $\hat{H}$ implies a satisfactory time evolution $e^{-\hat{H}a}$, no matter how large $\hat{H}a$ is.

In lattice field theory the Hamiltonian is defined through the time evolution operator, or “transfer matrix” [2]. Therefore, sect. 5.1 starts by reviewing and extending the construction of ref. [24] to the actions $S_{0}$ and $S=S_{0}+S_{B}+S_{E}$. A by-product of this analysis is the demonstration that there is no need to improve the temporal derivative in eq. (2.3). This feature is familiar from the static and nonrelativistic formulations. It is a special blessing here, because a temporal next-nearest-neighbor interaction would introduce unphysical states [6], and at large $m_{0}a$ the physical and unphysical levels cross. With the transfer matrix in hand, sects. 5.2 and 5.3 develop an expansion in $a$ for the Hamiltonian itself.

5.1 Construction of the transfer matrix

The transfer-matrix construction with two hopping parameters differs little from the usual case [24]. The transfer matrix acts as an integral operator in the space of gauge fields; in the $U_{0}=1$ axial gauge a wave functional $\Omega_{t}(U)$ at time $t$ evolves to

[TABLE]

at time $t+1$. The wave functional $\Omega_{t}(U)$ is also a vector in the fermion Hilbert space. For the standard gauge action the kernel may be written

[TABLE]

The factors arising from the fermion action $\hat{{\cal T}}_{F}^{(\dagger)}$ are operators in the fermion Hilbert space. The factors arising from the gauge field, ${\cal T}_{G}$ and ${\cal K}_{G}$, are given in ref. [24]; in the following, they do not play a crucial role, so we do not discuss them further.666Different from ref. [24] is the convention for the factors $(1\pm\gamma_{0})$ in the action (compare eq. (1.6) with eq. (2) of ref. [24]). With our convention it is natural for time-ordering to place later times to the left. Thus, the kernel ${\cal K}(V,U)$ transfers the field from $U$ at time $t$ to $V$ at time $t+1$.

The fermion operator for action $S_{0}$ can be written

[TABLE]

where (cf. ref. [24])

[TABLE]

[TABLE]

in a matrix notation in which $\Psi$ and $\bar{\Psi}$ are vectors and ${\cal B}_{U}$, \mbox{\boldmathD}_{U}, and ${\cal M}_{U}$ are matrices depending on gauge field $U$. The vectors and matrices of this notation are labeled by spin, color, and space. The covariant difference operator $D$ is as in eq. (2.5) and

[TABLE]

[TABLE]

The operators $\hat{\Psi}$ and $\hat{\bar{\Psi}}=\hat{\Psi}^{\dagger}\gamma_{0}$ obey canonical anti-commutation relations

[TABLE]

where $m$ and $n$ label spatial sites and $a$ and $b$ are multi-indices for spin and color. The fields corresponding to these operators are related to the original fields by

[TABLE]

This discrepancy in normalization between the integration variables in the functional integral and the canonical operators in Hilbert space demonstrates again that the normalization convention for the field $\psi(x)$, cf. eq. (2.1), is arbitrary. On the other hand, the propagator of $\Psi$ has unit residue at tree level, and a perturbative series $Z_{2\Psi}=1+g_{0}^{2}Z_{2\Psi}^{[1]}(m_{0}a)+\cdots$ beyond tree level.

The generalization of eqs. (5.2)–(5.9) to include the chromoelectric interactions suffers from a technical difficulty. Usually one uses the “four-leaf clovers” in eq. (2.12) as the lattice approximants to the chromomagnetic and chromoelectric fields. For the chromomagnetic interaction, this choice poses no problem, because $B$ involves link variables from one timeslice only. For the chromoelectric interaction, however, the time-space four-leaf clover involves link variables from three timeslices. In that case, the construction of the gauge-field transfer matrix is more complicated, and, if the improved gauge action is any indication, it may no longer be positive [25].

To avoid this complication one can define a chromoelectric field on only two timeslices. Consider

[TABLE]

where

[TABLE]

is defined on a two-leaf clover. The projection operators \mbox{\small\frac{1}{2}}(1\pm\gamma_{0}) in eq. (5.10) are chosen by analogy with the Wilson time derivative, cf. eq. (2.3), and as a result the standard transfer-matrix construction goes through with minor modifications. The two-leaf version $S_{E2}$ differs from the four-leaf version $S_{E}$ by an interaction of dimension six, so it should not alter the tree-level tuning of $c_{E}$.

Extending the transfer-matrix construction to $S_{0}+S_{B}+S_{E2}$ (eqs. (1.6), (1.7), and (5.10)), one finds the following changes. The chromomagnetic interaction modifies the matrices ${\cal B}$ and ${\cal M}$ to

[TABLE]

[TABLE]

Except for the new ${\cal B}$, eqs. (5.8) and (5.9) still hold. The chromoelectric interaction, eq. (5.10), modifies the fermion operator $\hat{{\cal T}}_{F}$ so that it depends on initial and final gauge fields $U$ and $V$:

[TABLE]

with $\hat{H}_{0}$ as in eq. (5.4) and ${\cal M}$ from eq. (5.13), but

[TABLE]

where the subscripts on $D$ and $E$ specify the spatial link fields, out of which they are constructed. The sign of the chromoelectric term in $\hat{H}_{I}^{\dagger}$ can be checked as follows: in our sign and $i$ conventions ${{t^{a}}^{\dagger}=-t^{a}}$ and Euclidean electric fields are anti-Hermitian operators in the gauge-field Hilbert space, {{\hat{\mbox{\boldmathE}\mkern 3.0mu}\mkern-3.0mu}{}^{a}}^{\dagger}=-{\hat{\mbox{\boldmathE}\mkern 3.0mu}\mkern-3.0mu}{}^{a}.

Comparing eqs. (5.13) and (5.15) with eqs. (5.7) and (5.5), respectively, one notices a pattern emerging. Interactions with block-diagonal Dirac matrices append to $e^{\cal M}$, whereas those with block–off-diagonal Dirac matrices modify $H_{I}$. This pattern depends only on the special Wilson time derivative and the technical assumption that all interactions live on only one or two timeslices. It proves that there is no need to alter the temporal derivative in eq. (2.3): higher-dimension “spatial” interactions are enough to achieve on-shell improvement, as asserted at the end of sect. 3.

5.2 Small $a$ expansion (general considerations)

From the transfer matrix one would like to derive the exact lattice Hamiltonian $\hat{H}=-\log\hat{{\cal K}}$. Of course, with this definition the Hamiltonian cannot be represented by a finite number of local operators. According to the Symanzik philosophy, however, one ought to expand it in powers of the lattice spacing. After obtaining the transfer matrix, (higher) time derivatives are no longer a concern, so the lattice-spacing expansion will hold if the quantities \mbox{\boldmathD}a, \mbox{\boldmathB}a^{2}, \mbox{\boldmathE}a^{2}, …, are small.

One can anticipate the expansion by enumerating the terms allowed by symmetry:

[TABLE]

where the subscript “cont” refers to an underlying continuum gauge field; below we usually suppress this subscript, for brevity. The coefficients $b_{i}$ depend on $m_{0}a$, and since eq. (5.16) is to be interpreted as an expansion in $a$ (rather than $1/m_{0}$), the $b_{i}$ for small $m_{0}a$ must be ${\rm O}((m_{0}a)^{p})$, with $p$ nonnegative. The coefficients for the action $S_{0}+S_{B}+S_{E}$, given in sect. 5.3, satisfy this requirement.

The general objective is to adjust the couplings so that eq. (5.16) takes the relativistic Dirac form, i.e. $b_{0}=b_{1}=1$ and $b_{2}=b_{B}=b_{E}=\cdots=0$. But, based on the considerations of sect. 3, there must be some leeway in the redundant directions. In the operator formalism adopted here, unitary changes of variables are possible, and these play the role of the isospectral transformation, eq. (3.2). Under a change of variables the Hamiltonian becomes

[TABLE]

where $\partial_{t}$ is a derivative with respect to imaginary time, and $\hat{{\cal U}}$ is the unitary operator implementing the change of variables in Hilbert space. Consider, for example, the following transformation:

[TABLE]

for which

[TABLE]

Such transformations are familiar from studies of the nonrelativistic limit of the Dirac equation, where they are called Foldy-Wouthuysen-Tani transformations [20]. Their characteristic feature is that the exponent is always a block–off-diagonal Dirac matrix.

The transformed Hamiltonian $\hat{H}^{\prime}$ has an expansion of the same form as in eq. (5.16), but with transformed coefficients:

[TABLE]

In light of the transformations, it is, therefore, enough to adjust $m_{0}a$, $\zeta$, $r_{s}$, $c_{B}$, and $c_{E}$, so that for some (hidden) value of $\xi_{1}$ the transformed Hamiltonian takes the Dirac form \hat{H}^{\prime}=\hat{\bar{\Psi}}(m_{q}+\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD})\hat{\Psi}. That means that one wants $b^{\prime}_{0}=b^{\prime}_{1}=1$ and $b^{\prime}_{2}=b^{\prime}_{B}=b^{\prime}_{E}=\cdots=0$.

The Foldy-Wouthuysen-Tani parameter $\xi_{1}$ drops out of on-shell quantities. It is preferable, therefore, to parameterize the redundant direction by one of the couplings. To this end, it is efficient to note that the following combinations of the $b_{i}$’s do not depend on $\xi_{1}$:

[TABLE]

A Hamiltonian unitarily equivalent to the Dirac Hamiltonian is then obtained whenever

[TABLE]

Eqs. (5.21) do not contain an invariant corresponding to $b_{E}$. This is analogous to the result, eq. (4.11), that the general mass dependence of $\zeta$ can only be determined hand-in-hand with $r_{s}$. In the present language, that connection arises as follows. Consider truncating eq. (5.16) at dimension four. Then only $b_{0}$ and $b_{1}$ remain. The Foldy-Wouthuysen-Tani transformation is {\rm O}(\mbox{\boldmathp}a), and superficially not worth considering. If one introduces it anyway, one sees immediately that $b_{1}$ changes, and in the transformation law, eq. (5.20), one power of $a$ has combined with the fermion mass to give $m_{q}a$. At $m_{q}a=0$ the condition $b_{1}=1$ is enough to determine $\zeta(0)$ unambiguously. But to obtain fully the mass dependence $\zeta(m_{0})$ one must consider simultaneously the interactions \bar{\Psi}\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}\Psi and \bar{\Psi}(\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD})^{2}\Psi.

A similar fate awaits the coefficient $b_{E}$ and its coupling $c_{E}$. Consider a two-parameter Foldy-Wouthuysen-Tani transformation

[TABLE]

The new parameter $\xi_{E}$ introduces changes that are superficially {\rm O}(\mbox{\boldmathp}^{2}a^{2}). The other coefficients are unaffected by $\xi_{E}$, but

[TABLE]

Again, one power of $a$ has combined with the fermion mass to give $m_{q}a$. Thus, the condition $b_{E}-\xi_{1}=0$ is enough to determine only $c_{E}(0)$. The full mass dependence of $c_{E}$ can only be revealed by considering simultaneously \bar{\Psi}\mbox{\boldmath\alpha}\!\cdot\!\mbox{\boldmathE}\Psi and the dimension-six interaction \bar{\Psi}\{\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD},\mbox{\boldmath\alpha}\!\cdot\!\mbox{\boldmathE}\}\Psi. This analysis is deferred to Appendix A.

The next subsection adjusts $m_{0}a$, $\zeta$, $r_{s}$, and $c_{B}$ to ensure eq. (5.22), and $c_{E}(0)$ to ensure $b^{\prime}_{E}(0)=0$. For heavy-light systems, the resulting lattice theory has cutoff artifacts of ${\rm O}(\Lambda_{\rm QCD}^{2}a^{2})$ and, only when $m_{Q}a\gg 1$, (yet smaller) artifacts of ${\rm O}(\Lambda_{\rm QCD}^{2}a/m_{Q})$ and of ${\rm O}(\Lambda_{\rm QCD}^{2}/m_{Q}^{2})$ as well. See sect. 6 for details. Moreover, for quarkonia, the lattice theory is similarly correct through ${\rm O}(v^{2})$.

5.3 Small $a$ expansion for $S_{0}+S_{B}+S_{E}$

Combining eqs. (5.2) and (5.14), and omitting factors that depend only on the gauge field, the fermion Hamiltonian of the lattice theory is

[TABLE]

where $H_{0}$ is specified by eqs. (5.4) and (5.13), and $H_{I}$ is specified by eq. (5.15). To derive an expression for the fermion Hamiltonian, one must coalesce the four exponents in eq. (5.25) into one. Owing to nontrivial commutators between $\hat{H}_{0}$, $\hat{H}_{I}$, and $\hat{H}_{I}^{\dagger}$, this is too difficult in general. But through order \mbox{\boldmathD}^{2}a^{2}, Appendix B achieves the desired result by a trick. There the field theory is mimicked by a toy model with the same algebraic structure but only two degrees of freedom. In the toy model one needs only to take the logarithm of a two-by-two transfer matrix, and expand the result in powers of $a$.

For small \mbox{\boldmathD}a the Hamiltonian becomes

[TABLE]

where

[TABLE]

The rest mass $M_{1}$ and the terms in parentheses come from expanding ${\cal M}$, and the functions $f_{i}$ are extracted from the toy model:

[TABLE]

Note that $f_{1}(0)=1$ and f_{2}(0)=\mbox{\small\frac{1}{2}}.

In the spirit of an underlying continuum gauge field one can identify $\triangle^{(3)}$ with \mbox{\boldmathD}^{2}_{\rm cont}, \mbox{\boldmathD}\times\mbox{\boldmathD} with \mbox{\boldmathB}_{\rm cont}, and [\partial_{t},\mbox{\boldmathD}] with \mbox{\boldmathE}_{\rm cont}. With these identifications one can cast eq. (5.26) into the form of eq. (5.16). Thus, the Hamiltonian of the action $S_{0}+S_{B}+S_{E}$ has coefficients

[TABLE]

and the invariants $B_{i}$ are

[TABLE]

The masses $M_{1}$ and $M_{2}$ are as before, and

[TABLE]

After imposing eq. (5.22), $M_{1}$, $M_{2}$, and $M_{B}$ all equal the physical mass.

The mass dependence of the couplings follows immediately from eqs. (5.30) and (5.22). The requirement $B_{1}=1$ implies

[TABLE]

precisely as in eq. (4.10). The requirement $B_{B}=0$ implies

[TABLE]

Finally, for small mass the chromoelectric coupling should be tuned to

[TABLE]

to enforce $b^{\prime}_{E}(0)=0$. With the axis-interchange invariant boundary condition $r_{s}(0)=1$, one thus recovers the action of ref. [6], with $\zeta(0)=r_{s}(0)=c_{E}(0)=c_{B}(0)=1$.

Our analysis has not yet specified the relevant couplings $g_{0}^{2}$ and $m_{0}a$. They, of course, are fixed not by theory but by experiment. In the Hamiltonian language, the bare mass $m_{0}a$ is adjusted so that $B_{0}=1$, i.e. $M_{1}=m_{q}$. Then the improvement conditions, eqs. (5.32)–(5.33), guarantee that $m_{q}=M_{2}=M_{B}$ also.

There is a special case of eqs. (5.32)–(5.34) that is of at least passing interest, namely the one for which the Foldy-Wouthuysen-Tani parameter $\xi_{1}=0$. This is obtained by choosing $r_{s}$ so that (untransformed) $b_{2}=0$:

[TABLE]

Then the condition $b_{1}=1$ requires

[TABLE]

and the condition $b_{E}=0$ requires

[TABLE]

and, as before, the condition $b_{B}=0$ requires $c_{B}=r_{s}$. After substituting eq. (5.35) into eq. (5.32) one re-obtains the right-hand side of eq. (5.36). Appendix A shows that—with $r_{s}$ and $\zeta$ from eqs. (5.35)–(5.36)—$c_{E}=1$ can be maintained for arbitrary $m_{0}a$.

6 Truncation Criteria Revisited

This section reexamines criteria for truncating a cutoff theory, with some emphasis on the errors left over after truncation. The analysis of the previous sections takes the scaling dimension of the interaction as a guide. For massless quarks that is certainly correct. But the most appropriate organization may vary when the same cutoff theory is applied to different physical systems. Thus, conclusions about the accuracy of a massive-fermion action must be refined, after deciding whether the action is to be applied to heavy-light systems or to quarkonia.

After the couplings have been adjusted to some practical accuracy, the Hamiltonian (possibly after a Foldy-Wouthuysen-Tani transformation) is

[TABLE]

the Coulomb potential appears if one transforms to a gauge without $A_{0}=0$. A lattice artifact $\delta\hat{H}_{\rm lat}$ remains, because one cannot exactly incorporate infinitely many terms into eq. (1.3).

One can estimate the errors induced by $\delta\hat{H}_{\rm lat}$ by treating it as a perturbation. There is an advantage to estimating cutoff effects from the Hamiltonian. In the action formalism, eq. (1.5), it may not be clear how the time discretization trickles down to physical quantities. But by proceeding through the transfer matrix these effects are treated exactly.

From the line of argument leading to eq. (5.16), one expects that $\delta\hat{H}_{\rm lat}$ consists of operators multiplied by mass-dependent coefficients

[TABLE]

where the power $s_{n}=\dim H_{n}-4$, and $L_{n}$ is the number of loops already under control. One can determine the effect of $\delta\hat{H}_{\rm lat}$ on a physical quantity from order-of-magnitude estimates for the operators $\hat{H}_{n}$ and general properties of the coefficients $b_{n}^{[l]}(m_{0}a)$. While the former depend on the physical process under study, the latter are process independent.

The dimension-five, tree-level coefficients have two important properties, which we believe are generic. First, at asymptotically large $m_{0}a$ the tree-level coefficients either approach a constant or fall as a power of $1/(m_{0}a)$. An analysis of higher-order Feynman diagrams (sect. 8) shows that tree-level patterns persist to all orders in perturbation theory. Indeed, the asymptotic behavior is presumably a consequence of the heavy-quark symmetries obeyed by all lattice actions under consideration. Second, the coefficients always contain the recurring ingredients ${1+m_{0}a}$, $m_{0}a(2+m_{0}a)$, and $\log(1+m_{0}a)$ in a way that makes implausible any combination that would blow up at an intermediate value of $m_{0}a$. Indeed, all evidence suggests that the functions $b(m_{0}a)$ are smaller than their low-order Taylor expansions, once m_{0}a\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{>}}1.

Let us now discuss the typical size of the operators in the Hamiltonian. Table 2 gives ballpark estimates for the

dimension-three, -four, and -five interactions for three systems: those in which all quarks are light, those with one heavy quark, and quarkonia. The row labeled $E_{0}$ in Table 2 gives the non-trivial dynamical scales, to which artifacts should be compared. In all-light and heavy-light systems, the estimates start from naive dimensional analysis, but heavy-quark bilinears with an off-diagonal Dirac matrices are $\Lambda_{\rm QCD}/m_{Q}$ times smaller still. In quarkonia, the estimates are those of ref. [12], with $v$ denoting the typical velocity of the heavy (anti-)quark in the bound state ($v\sim\alpha_{s}(m_{Q})$).

A conservative estimate of the artifact is then as follows: Choose a system, multiply by $a^{s_{n}}$, and compare to $E_{0}$. The coefficient $b(m_{0}a)$ is a number of order 1 (or less), for any value of $m_{0}a$, so its numerical value does not affect the (conservative) conclusion. If, after suitable adjustment of the couplings, one finds $b\sim m_{q}a$ for $m_{q}a\ll 1$, or $b\sim 1/(m_{Q}a)$ for $m_{Q}a\gg 1$, the artifact might be even smaller.

Consider the chromoelectric interaction as an example. For the sake of argument, suppose that the rest mass $M_{1}$ and the kinetic mass $M_{2}$, and hence $m_{0}a$ and $\zeta$, have been adjusted nonperturbatively. If $c_{E}$ is not adjusted correctly, then the (transformed) coefficient $b^{\prime}_{E}$ of the chromoelectric term in the Hamiltonian does not vanish. Then, relative to the corresponding $E_{0}$, there are artifacts of ${\rm O}(\Lambda_{\rm QCD}a)$ for all-light, ${\rm O}(\Lambda_{\rm QCD}^{2}a/m_{Q})$ for heavy-light, and ${\rm O}(m_{Q}av^{2})$ for heavy-heavy. If instead $c_{E}(m_{0}a)$ is adjusted to $c_{E}(0)$ in eq. (5.34), the artifacts in all-light systems fall to ${\rm O}(m_{q}\Lambda_{\rm QCD}a^{2})$. With heavy quarks the estimates depend on $m_{Q}a$. If $a$ is so tiny that $m_{Q}a\ll 1$, then $b^{\prime}_{E}\sim m_{Q}a$, and the chromoelectric artifact is reduced to ${\rm O}(\Lambda_{\rm QCD}^{2}a^{2})$ for heavy-light and to ${\rm O}(m_{Q}^{2}a^{2}v^{2})$ for heavy-heavy. But if m_{Q}a\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{>}}1, it turns out that $b^{\prime}_{E}$ either remains constant or falls as $1/(m_{Q}a)$, depending on the mass dependence of the redundant coupling $r_{s}$. The artifacts are then either ${\rm O}(\Lambda_{\rm QCD}^{2}a/m_{Q})$ for heavy-light and ${\rm O}(m_{Q}av^{2})$ for heavy-heavy, or $1/(m_{Q}a)$ times smaller.

The appearance of $1/(m_{Q}a)$ in coefficients, in addition to the $\Lambda_{\rm QCD}/m_{Q}$ in heavy-light dynamics, makes the error analysis of heavy-light systems somewhat delicate. Since the $1/(m_{Q}a)$ behavior arises only if $m_{Q}a\gg 1$, it leads only to errors that are smaller than the usual discretization errors, relative to $E_{0}$. On the other hand, occasionally one is interested in effects that are subleading in the heavy-quark expansion. For a given lattice action, such quantities may have a larger relative error. For example, even with $c_{E}(0)$ adjusted correctly, the fine structure of the heavy-light spectrum, which is ${\rm O}(\Lambda_{\rm QCD}^{2}/m_{Q})$, suffers a relative error of order $(\Lambda_{\rm QCD}^{3}a/m_{Q})/(\Lambda_{\rm QCD}^{2}/m_{Q})\sim\Lambda_{\rm QCD}a$.

Similar comments apply to quarkonia. Though the chromomagnetic and chromoelectric interactions are of order $m_{Q}v^{4}/m_{Q}v^{2}\sim v^{2}$ smaller than the spin-independent kinetic energy, they introduce relative errors on spin-dependent structure of order $m_{Q}v^{4}/m_{Q}v^{4}\sim 1$. A full ${\rm O}(v^{4})$ analysis requires a few dimension-six and -seven interactions, which we consider in Appendix A.

Once the dimension-five couplings $c_{B}$ and $c_{E}$ have been properly adjusted, lattice artifacts remain from dimension six and higher. Table 3 lists bilinear operators that can appear in the

Hamiltonian. The conservative estimate of the absolute errors caused by these operators is to multiply Table 3 by $a^{2}$. When $m_{Q}a\gg 1$, however, some of the contributions may be, as before, a factor of $1/(m_{Q}a)$ or $1/(m_{Q}a)^{2}$ smaller. But, again, this subtlety is only crucial when quantities subleading in the heavy-quark expansion are at issue.

The four-fermion interactions, listed in Table 4,

are also of dimension six. To generalize the analysis of sect. 3 to encompass these operators, one must simultaneously treat dimension-six gauge-field interactions [6]. The result is that $(\bar{\psi}\gamma_{0}t^{a}\psi)^{2}$ and $(\bar{\psi}\gamma_{i}t^{a}\psi)^{2}$ are redundant [6], even without axis-interchange symmetry. The other four-fermion interactions arise first at the one-loop level.

Let us summarize the main points of this section for heavy-light spectroscopy with action ${S_{0}+S_{B}+S_{E}}$. After the tree-level adjustments of sect. 5.3 have been applied, the largest remaining lattice artifacts are ${\rm O}(g_{0}^{2}\Lambda_{\rm QCD}a)$ from the one-loop maladjustment of $S_{B}+S_{E}$ and ${\rm O}(\Lambda_{\rm QCD}^{2}a^{2})$ from unadjusted dimension-six interactions. The mass dependence of the artifacts is solely in the coefficients $b(m_{q}a)$, which is a number of order unity at any mass.

7 Electroweak Perturbations

This section extends the formalism of the previous sections to the two- and four-quark operators of the electroweak Hamiltonian, which may be treated as a first-order perturbation to QCD. The construction of the renormalized (or continuum-limit) operator is analogous to the construction of the renormalized trajectory. Let ${\cal O}$ denote the continuum operator. Then

[TABLE]

where the sum runs over all lattice operators $O_{n}$ with the same quantum numbers as ${\cal O}$.777By convention, the zeroth lattice operator $O_{0}$ has the same dimension as the continuum-limit operator ${\cal O}$. The role of the other $C_{n}O_{n}$ is to remove terms suppressed (or enhanced!) by a power of $a$. The role of $Z_{\cal O}$ is to convert to a preferred renormalization convention. Like the couplings in the action, the coefficients $Z_{\cal O}$ and $C_{n}$ are functions of the relevant couplings, all fermion masses $\{m_{0}a\}$ and the gauge coupling $g_{0}^{2}$. Eq. (7.1) is general, but we again consider perturbative expansions in $g_{0}^{2}$,

[TABLE]

and focus on tree level. Previous work [7] applied to small masses, but we treat the mass dependence of $Z_{\cal O}^{[l]}(\{m_{0}a\})$ and $C_{n}^{[l]}(\{m_{0}a\})$ exactly. We also do not impose axis-interchange invariance in classifying the lattice operators $O_{n}$.

The coefficients $Z_{\cal O}$ and $C_{n}$ can be determined from low-momentum matrix elements of all $O_{n}$, analogously to sect. 4. In perturbation theory it is enough to compute matrix elements between quark, anti-quark, and gluon states, both with lattice and ${\overline{\rm MS}}$ regulators. It is essential to impose consistent normalization conditions. Appendix C derives external-state rules for lattice perturbation theory. There one finds that the contraction of $\psi(x)$ with a normalized fermion state corresponds to a factor u_{\rm lat}(\xi,\mbox{\boldmathp}){\cal N}(\mbox{\boldmathp}), where $u_{\rm lat}$ is a normalized spinor on the lattice mass shell. The factor

[TABLE]

where (for $S_{0}$) \mu(\mbox{\boldmathp})=1+m_{0}a+\mbox{\small\frac{1}{2}}r_{s}\zeta\hat{\mbox{\boldmathp}}^{2}a^{2}. A relativistic theory has instead u_{\rm rel}(\xi,\mbox{\boldmathp})\sqrt{m_{q}/E}, where $u_{\rm rel}$ and $E$ comply with the relativistic mass shell.

Consider the bilinear operator ${\cal J}_{\Gamma}^{fg}$ that creates flavor $f$ and annihilates flavor $g$ with spin coupling $\Gamma$. At tree level its matrix elements should be

[TABLE]

where $E_{f}^{(\prime)}$ is the energy of flavor $f$ with momentum \mbox{\boldmathp}^{(\prime)} and mass $m_{f}$. Note that the relativistic spinors $u_{\rm rel}$ and $v_{\rm rel}$ appear on the right-hand side.

With the right-hand side of eq (7.4) as a target, we now consider lattice operators $O_{n}$. The simplest lattice bilinear with the correct dimension and quantum numbers is

[TABLE]

which corresponds to $O_{0}$ in eq. (7.1). Recall that $\psi(x)$ is the field appearing in the mass form of the action, eq. (2.3). At tree level the matrix elements are

[TABLE]

where {\cal N}_{f}(\mbox{\boldmathp}) is the normalization factor of flavor $f$, cf. eq. (7.3). Note that the lattice spinors $u_{\rm lat}$ and $v_{\rm lat}$ appear on the right-hand side.

Setting \mbox{\boldmathp}=\mbox{\boldmathp}^{\prime}={\bf 0}, the matrix elements differ only because of the factors ${\cal N}({\bf 0})$. Thus $Z_{\Gamma}J_{\Gamma}$ has the same zero-momentum matrix elements as the target ${\cal J}_{\Gamma}$, in all four channels, if the (re)normalization factor

[TABLE]

This is a tree-level result, but the mass dependence shown here remains dominant to all orders, cf. sect. 8.

Further terms in the three-momentum expansion cannot be matched without considering higher-dimension terms in eq. (7.1). At tree level one sees the differences between eqs. (7.4) and (7.6) in the factors ${\cal N}\neq\sqrt{m/E}$ and spinors $u_{\rm lat}\neq u_{\rm rel}$. Eq. (7.1) can therefore be extended to higher dimension by introducing an improved field. To first order in \mbox{\boldmathp}a consider

[TABLE]

with flavor labels implied. Then

[TABLE]

is the target operator of interest, through first order in \mbox{\boldmathp}a, if $d_{1}$ is adjusted properly. Comparing the bracket in eq. (7.8) with those in eqs. (C.25) and (C.27), one finds

[TABLE]

identifying $m_{q}=M_{2}$.888The substitution of the kinetic mass $M_{2}$ for the rest mass $M_{1}$ is done so that the expression remains valid under a nonrelativistic interpretation explained in sect. 9.

For small mass one finds $d_{1}\propto m_{0}a$; the only ${\rm O}(a)$ improvement needed is the normalization factor $e^{M_{1}a/2}$. At large mass, however, the rotation of eq. (7.8) becomes important. Analogously to the Hamiltonian coefficients discussed in sect. 6, when $m_{0}a\gg 1$, one has $d_{1}\approx 1/(2m_{q})$. Consequently, the contribution of d_{1}\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD} is essential for computing the $1/m_{q}$ correction to the static limit of matrix elements of ${\cal J}_{\Gamma}$. Similarly, higher-dimension generalizations of eq. (7.8) are needed to obtain $1/m_{q}^{2}$ and corrections of higher order in $1/m_{q}$..

The improved field $\Psi_{\rm I}(x)$ in eq. (7.8) coincides, through {\rm O}(\mbox{\boldmathp}), with the one denoted by $\Psi$ in sect. 5. Combining eqs. (5.9) and (5.18), the Foldy-Wouthuysen-Tani transformed field is

[TABLE]

where $\xi_{1}$ parameterizes the solution of the tuning conditions. This expression is (numerically) cumbersome, but one may expand consistently the exponentials in $a$. This exercise identifies $\xi_{1}$ with $d_{1}$. Indeed, solving $b^{\prime}_{1}=1$ for $\xi_{1}$ yields the right-hand side of eq. (7.10), after replacing the rest mass $M_{1}$ with the kinetic mass $M_{2}$.88footnotemark: 8

The special role of $\Psi$ should not be too surprising, because it possesses two important properties. First, it satisfies canonical anti-commutation relations and is thus properly normalized. Second, its dynamics are given by the Dirac Hamiltonian \hat{H}=\hat{\bar{\Psi}}(m_{q}+\gamma_{0}A_{0}+\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD})\hat{\Psi}—at least at tree level and up to {\rm O}(\mbox{\boldmathp}^{2}). Therefore, any operator built out of the transformed field yields the desired matrix elements, also at tree level and up to {\rm O}(\mbox{\boldmathp}^{2}).

Let us conclude this section with some comments on two other Ansätze for the currents. A formal argument based on the Ward identity suggests that a conserved current999Both “Noether” and “gauge” currents are conserved; they differ by $\sigma_{\mu\nu}$ terms. is especially suited to the determination of form factors of a vector current or the decay amplitude of a vector meson. But although the Ward identity implies a certain universality in radiative corrections, it does not imply any special mass dependence at tree (or any other) level.

With standard Feynman rules and Appendix C, straightforward algebra yields the tree-level on-shell matrix elements. To {\rm O}(\mbox{\boldmathp}) the (conserved) gauge current $V^{\rm G}_{\mu}$ has matrix elements

[TABLE]

[TABLE]

through {\rm O}(\mbox{\boldmathp}a), where $M_{1}$, $M_{2}$, and $M_{B}$ are the tree-level masses,101010For the Noether current the terms proportional to $c_{B}$ (implicitly in $1/M_{B}$) and $c_{E}$ (in $1/M_{G}$) in eqs. (7.12) and (7.13) would not appear. but

[TABLE]

The Ward identity asserts that these tree-level masses all renormalize in a coherent way. But although the “forward-scattering” matrix elements in eq. (7.12) are correct (assuming $M_{2}=M_{B}=m_{q}$), the “annihilation” matrix elements in eq. (7.13) are not (unless $m_{q}a\ll 1$). We conclude, therefore, that $V^{\rm G}_{\mu}$ is not useful for determining the decay constant of a massive vector meson.

Ref. [7] suggests using a “(four-dimensional) rotated current”

[TABLE]

To ascertain if $J_{\Gamma,\rm rot}^{fg}$ matches the target continuum operator, one must evaluate matrix elements, as above. The timelike translations in $D_{0}$ greatly change the mass dependence. One finds that $Z_{\Gamma,\rm rot}J_{\Gamma,\rm rot}^{fg}$ has correctly normalized matrix elements only if

[TABLE]

where $Z_{\Gamma}$ is the normalization factor of the unrotated bilinear, eq. (7.7). Moreover, when $m_{q}a\not\ll 1$ the rotation of eq. (7.15) must be supplemented à la eq. (7.8), with the same $d_{1}$ as in eq. (7.10). Thus, mass-dependent improvement of eq. (7.15) is analogous to improvement of eq. (7.5), but the latter is simpler.

In summary, the mass dependence of electroweak operators is tractable, if one proceeds as follows. First, start with a simple operator $O_{0}$ and expand its on-shell matrix elements in external, spatial momenta small in lattice units. As usual, there is no need to expand in $m_{0}a$. Second, add additional terms $C_{n}O_{n}$ to correct the momentum dependence of the matrix elements. At least to tree level this step can be accomplished by field rotations, as in eqs. (7.8) and (A.17). Finally, normalize $O_{0}+\sum_{n}C_{n}O_{n}$ to obtain the fully renormalized operator ${\cal O}$ in the desired renormalization scheme. For example, through {\rm O}(\mbox{\boldmathp}a) the renormalized bilinear ${\cal J}_{\Gamma}^{fg}$ is given by eqs. (7.9), (7.8), and (7.10).

8 Beyond Tree Level

In the previous sections, the $m_{0}a$ dependence of the couplings in the action is derived at tree level. This section considers what happens beyond tree level.

In perturbation theory the expressions for the masses introduced previously become power series in $g_{0}^{2}$. For example,

[TABLE]

and

[TABLE]

where $M_{1}^{[0]}$ and $M_{2}^{[0]}$ are given by eqs. (4.8) and (4.9), respectively. After calculating the self-energy to $l$ loops, one can extract the coefficients $M_{i}^{[l]}$ as functions of $m_{0}a$. The requirement $M_{1}=M_{2}=m_{q}$ subsequently yields the perturbative power series for the couplings $r_{s}$ and $\zeta$. (Based on the arguments of sects. 3 and 5, the Wilson term’s coupling $r_{s}$ should be redundant to all orders in $g_{0}^{2}$.) In the same vein, the on-shell fermion-gluon vertex function to $l$ loops yields $c_{B}^{[l]}(m_{0}a)$ and $c_{E}^{[l]}(m_{0}a)$, and electroweak matrix elements to $l$ loops yield $Z_{\cal O}^{[l]}(m_{0}a)$ and $d_{i}^{[l]}(m_{0}a)$.

A complete derivation of one-loop corrections is beyond the scope of this paper. It is easy, however, to assess two qualitative features: the mass dependence of loop diagrams (sect. 8.1) and the expected size of corrections from tadpole diagrams (sect. 8.2).

8.1 Mass dependence of loop diagrams

This subsection shows that the mass dependence of loop diagrams is benign. Although we focus on the specific action $S=S_{0}+S_{B}+S_{E}$, the conclusions hold for any action with the Wilson time derivative and arbitrary spatial interactions. Actions with next-nearest-neighbor interactions in time are problematic starting at tree level [6], so they are not considered here.

Let us first consider vacuum polarization. At one loop it is easy to see that the lattice-regulated Feynman integrals for vacuum polarization are smooth functions of the fermion mass. Moreover, for large fermion mass the integrals vanish as $(m_{0}a)^{-2}$; rigorously so, because the momentum-dependent terms in the fermion propagator are bounded. The behavior is the same for a closed fermion loop with any number of gluons attached. Hence, internal heavy-fermion loops decouple precisely as expected.

The self-energy and vertex corrections are less trivial, because the external momenta are set on shell. Fig. 1 shows the one- and two-gluon vertices. For the action $S=S_{0}+S_{B}+S_{E}$ (with the four-leaf clover for $S_{E}$)

[TABLE]

[TABLE]

where $p$ and $k$ are the incoming fermion and gluon momenta, respectively.

The expression for $X^{ab}_{\mu\nu}$ is not needed, except to note that its mass dependence is qualitatively the same as $\Lambda_{\mu}^{[0]}$.

A formal way of going to the mass shell is to put $p_{0}=iE\gamma_{0}$, with $E$ from eq. (4.4). The $\gamma_{0}$ in the analytic continuation is not rigorous when applied under an integral, but the $m_{0}a$ dependence comes out right. The temporal vertex $\Lambda_{0}^{[0]}$ is proportional to $e^{M_{1}^{[0]}a}=1+m_{0}a$. The spatial vertex $\Lambda_{i}^{[0]}$ is proportional to $\zeta$, which, when it is tuned so that $M_{1}=M_{2}$, satisfies $r_{s}\zeta\approx(1+m_{0}a)/M_{1}a$ for large mass. Similar behavior holds for quark–multi-gluon vertices. Finally, the inverse propagator is also proportional to $1+m_{0}a$, for $p$ close to the mass shell.

Consider any process with an external fermion line. Loop diagrams can be built up from the tree diagram by adding more gluons. Each additional vertex on the external line requires an additional fermion propagator. The dominant mass dependence of the propagator-and-vertex combination is $(1+m_{0})/(1+m_{0})$ or $\zeta/(1+m_{0})$, and thus cancels always.

For example, all diagrams in the self-energy are proportional to $1+m_{0}a$. After summing the geometric series and integrating over $p_{0}$ one finds

[TABLE]

[TABLE]

where $M_{1}^{[1]}(m_{0}a)$ and $Z_{2}^{[1]}(m_{0}a)$ depend mildly on the mass, varying smoothly from the value obtained for massless fermions to the value in the static formulation.

The same happens to the fermion-gluon vertex. The gauge-coupling renormalization factor is defined through the fermion-gluon vertex via

[TABLE]

where \Lambda_{\mu}(\mbox{\boldmathp}^{\prime},\mbox{\boldmathp}) is the full vertex function, including leg contributions.111111At tree level one verifies $Z_{g}^{[0]}=1$ from $\Lambda_{0}^{[0]}$, and also from $\Lambda_{i}^{[0]}$ if $m_{q}=M_{2}$. In perturbation theory one usually organizes the calculation by treating the legs and the proper vertex separately. By gauge invariance

[TABLE]

where $Z_{3}$ ($Z_{2}$) and $Z_{1}$ ($Z_{\rm 1F}$) are the gluon (fermion) wavefunction and proper vertex renormalization factors. The strong mass dependence of $Z_{2}$ must, therefore, cancel against $Z_{\rm 1F}$. (The residual mass dependence of $Z_{\rm 1F}/Z_{2}$ should be the same as $Z_{1}/Z_{3}$ to satisfy the expectations of decoupling.) Indeed, at tree level the temporal vertex provides the asymptotic factor $1+m_{0}a$, and, by the general argument, the full proper vertex has the same (dominant) mass behavior to all orders. Hence,

[TABLE]

where $Z_{\rm 1F}^{[1]}(m_{0}a)$ again depends only mildly on the mass. With the spatial vertex $\Lambda^{[0]}_{i}$ the factor $\zeta$ compensates for the missing factor of ${1+m_{0}a}$ to ensure that $1/M_{2}$ appears, so the ${1+m_{0}a}$ counting is the same.

For electroweak currents and four-quark operators, the analysis of the mass dependence is similar. Again, loop diagrams have the same leading mass dependence as tree diagrams for the same process. For example, the bilinear $J_{\Gamma}^{fg}$, defined in eq. (7.5), has renormalization constant

[TABLE]

The mass dependence of the loop corrections $Z^{[l]}_{\Gamma}(m_{0f}a,m_{0g}a)$ smoothly connects massless and static results. Such behavior is borne out in sect. 10’s nonperturbative check of the local vector current, for which $\Gamma=\gamma_{0}$.

The considerations of this subsection argue that the large-mass limit of actions described by eq. (1.3) is well-behaved in perturbation theory. More generally, the physical masses and, hence, the couplings could depend on the gauge coupling in a nonperturbative way. But because the origin of the gauge-coupling dependence is the region of momentum space near the cutoff, it seems unlikely that nonperturbative contributions would overwhelm the perturbative contribution, at least once the cutoff is large enough. Should perturbation theory prove inadequate, however, a nonperturbative renormalization group could, in principle, substitute for perturbative calculations.121212For example, to tune $\zeta$ nonperturbatively, compute the energy of a meson and imposing $M_{q}=M_{2}$. Nevertheless, it seems implausible that nonperturbative effects are more worrisome at large mass than at small. Thus, the main conclusion, that the large-mass behavior of interacting fermions is benign, is probably valid nonperturbatively.

8.2 Mean field theory

To estimate the one-loop corrections, recall that the dominant contributions come (in Feynman gauge) from tadpole diagrams, which originate from higher-order terms in the expansion of the link matrix U_{\mu}=1+g_{0}A_{\mu}+\mbox{\small\frac{1}{2}}g_{0}^{2}A_{\mu}^{2}+\cdots. It is possible to make this observation more systematic [23]. Wherever the gauge field appears, substitute

[TABLE]

where $u_{0}$ is a gauge-invariant average of the link matrices. The substitution should be understood in the following sense: The second factor $[U_{\mu}/u_{0}]$ produces perturbative series with small coefficients. The first factor $u_{0}$, which has a nasty tadpole-dominated perturbative series, should be absorbed into the couplings $c_{n}$ and into renormalization factors $Z_{\cal O}$. When a numerical value for $u_{0}$ is needed, for example in a Monte Carlo calculation, it should be taken from the Monte Carlo itself.

With this prescription the hopping-parameter form of $S_{0}$ remains as in eq. (1.6), but with $U_{\mu}\to U_{\mu}/u_{0}$ and

[TABLE]

The mass form of $S_{0}$ is given by eqs. (2.3) but with difference operators defined with

[TABLE]

instead of $T_{\pm\mu}$, and mass

[TABLE]

instead of $m_{0}$. Finally, an overall factor of $u_{0}$ multiplies each term in the action.

The clover-leaf construction used to define the chromomagnetic and chromoelectric fields contains products of four $U$ matrices. If one replaces the gauge fields $B$ and $E$ with tadpole-improved clovers, the interactions $S_{B}$ and $S_{E}$ are given by eqs. (1.7) and (1.8), respectively, but with

[TABLE]

instead of $c_{B}$ and $c_{E}$, and $\tilde{\kappa}_{s}$ instead of $\kappa_{s}$. The fourth factor of $u_{0}$ corresponds to the overall factor mentioned above.

After these rearrangements one can immediately generalize the expressions in sects. 4 and 5 to the mean-field level. They remain the same as before, but with $m_{0}\to\widetilde{m}_{0}$, $c_{B}\to\tilde{c}_{B}$, and $c_{E}\to\tilde{c}_{E}$. Consequently, the couplings $\zeta$, $\tilde{c}_{B}$, and $\tilde{c}_{E}(0)$ should be adjusted to the right-hand sides of eqs. (5.32)–(5.34), but with $m_{0}\to\widetilde{m}_{0}$. The resulting conditions represent a set of mean-field-theory predictions at $g_{0}^{2}\neq 0$, given a nonperturbative input for $u_{0}$. One-loop calculations with $m_{0}\neq 0$, $c_{B}\neq 0$, and $c_{E}\neq 0$ will test and correct mean-field theory estimates.

At currently accessible lattice spacings ref. [23] has shown that, with this mean-field reorganization and a sensible choice of expansion parameter, the bare perturbative series converges quickly in many cases. Calculations [18] in one-loop perturbation theory of Feynman diagrams needed to determine the $c_{n}^{[1]}(m_{0}a)$ show a smooth transition from the massless to the static limits.131313For tadpole and scale-choice improvement [23] of the static limit and of nonrelativistic QCD, see refs. [26, 27]. One therefore expects the essential concepts of ref. [23] to apply to the $\tilde{c}_{n}$ and to the coefficients in eq. (7.1) too. Indeed, in the one case for which a nonperturbative check is unambiguous, the normalization of the vector current, there is excellent agreement with mean-field theory, cf. fig. 4 in sect. 10.

9 The Nonrelativistic Limit

It is illuminating to adapt the methods of sect. 5 to the nonrelativistic and static limits. Rather than adjusting the couplings to obtain the Dirac Hamiltonian, one could instead aim for the nonrelativistic Pauli Hamiltonian (and generalizations thereof). An advantage of this avenue is that it provides a useful physical picture even when the couplings are maladjusted, in particular when $\zeta=1$. Many Monte Carlo studies have used actions with $\zeta=1$, and it would be helpful to have a framework for interpreting their data in the heavy-quark regime. Indeed, the analysis of this section shows that $\zeta=1$ is acceptable for nonrelativistic fermions. Even the Wilson action (somewhat crudely) approximates the properties of nonrelativistic or heavy-quark effective theory, provided $m_{0}a$ is adjusted correctly. Similarly, the Sheikholeslami-Wohlert action is a better approximation.

The Hamiltonian of the action $S_{0}+S_{B}+S_{E}$ can be brought to the nonrelativistic Pauli form with the Foldy-Wouthuysen-Tani transformation. Imagine transforming the \mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD} and \mbox{\boldmath\alpha}\!\cdot\!\mbox{\boldmathE} terms away completely. Afterwards, the transformed Hamiltonian reads

[TABLE]

where $M_{1}$, $M_{2}$, and $M_{B}$ are as in sects. 4 and 5. The new mass $M_{E}$ reduces to $M_{2}$ with suitable mass dependence of $c_{E}$ (cf. Appendix A), or as $m_{0}a\to 0$. The specific expression is not needed here. The Pauli form of eq. (9.1) has no coupling between the upper (particle) and lower (anti-particle) components of $\Psi$, as in the explicitly nonrelativistic formulations [11, 12]. Here, however, eq. (9.1) is derived within the lattice theory, rather than being an Ansatz for an effective lattice theory.

Let us discuss the physics of each term in eq. (9.1). The first three are the rest mass, Coulomb potential, and kinetic energy141414Because of this physical interpretation the quantity $M_{2}$, defined in eq. (4.7), is called the kinetic mass. of the fermion. The \mbox{\boldmath\Sigma}\!\cdot\!\mbox{\boldmathB} term, as one recalls from atomic physics, produces the hyperfine splitting. The last term can be rewritten

[TABLE]

The two parentheses give the (non-Abelian) spin-orbit and Darwin interactions, respectively.

The Pauli Hamiltonian is quantitatively useful only if the fermion is nonrelativistic. Given nonrelativistic velocities, however, eq. (9.1) remains applicable even when the various masses are unequal. Fig. 2 is a sketch of the quarkonium spectrum, illustrating how the masses affect the spectrum.

The interesting gross feature of the spectrum is not the overall mass gap—close to $2M_{1}$—but the pattern of radial and orbital excitations, e.g. $m_{\rm 2S}-m_{\rm 1S}$ or $m_{\rm 1P}-m_{\rm 1S}$. These splittings are dictated by the kinetic mass $M_{2}$. Following the analysis of ref. [12] they are of order $M_{2}v^{2}$, where $v$ is the typical velocity of a heavy quark in quarkonium. ($v\sim 0.3$ for charmonium, and $v\sim 0.1$ for bottomonium.) Further application of the velocity counting in ref. [12] to eq. (9.1) shows that the hyperfine splittings are \mbox{\boldmath\Sigma}\!\cdot\!\mbox{\boldmathB}/M_{B}\sim M_{2}^{2}v^{4}/M_{B}, and the spin-orbit splittings are [\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD},\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathE}]/M_{E}^{2}\sim M_{2}^{3}v^{4}/M_{E}^{2}.

The preceding paragraph merely reviews the well-known argument that the rest mass of a nonrelativistic particle decouples from the interesting dynamics. In our formalism the reasoning suggests the following strategy: forget about $M_{1}$ and adjust the bare mass so that the kinetic mass $M_{2}$ takes the physical value. Meanwhile, choose the coupling $\zeta$ by convenience. The obvious example is to take $\zeta=1$, as in the Wilson and Sheikholeslami-Wohlert actions.

Since the Wilson and Sheikholeslami-Wohlert actions represent viable nonrelativistic field theories, it makes sense to compare them to the explicitly nonrelativistic theories. The (tree-level) masses for the Wilson action are plotted as a function of $m_{0}a$ in fig. 3.

Assuming $m_{0}a$ is chosen so that $M_{2}=m_{q}$, the other masses satisfy $M_{1}<m_{q}$, $M_{B}^{-1}<m_{q}^{-1}$, and $M_{E}^{-1}<m_{q}^{-1}$. The simplest form of nonrelativistic QCD [11] has Hamiltonian \hat{H}_{\rm NR}=\hat{\bar{\Psi}}\mbox{\boldmathD}^{2}\hat{\Psi}/(2m_{q}). Thus, in our notation, $\hat{H}_{\rm NR}$ has $M_{1}=0$ and $M_{B}^{-1}=M_{E}^{-1}=0$. Thus, the Hamiltonians of the Wilson and simplest nonrelativistic theories make the same errors qualitatively. For example, in both one expects the fine and hyperfine splittings to be too small. Similarly, for the Sheikholeslami-Wohlert action one finds $M_{B}=M_{2}$, and thus good hyperfine splittings, but $M_{E}^{-1}>m_{q}^{-1}$, so the splittings between $\chi_{J}$ states ought to be too large. To obtain the correct spin-orbit splittings, one needs the mass dependence of $c_{E}$, cf. Appendix A.

A parallel set of remarks applies to heavy-light systems. The Hamiltonian of the lattice theory satisfies the usual heavy-quark symmetries as $m_{q}\to\infty$, no matter what $M_{2}$, $M_{B}$, and $M_{E}$ are. On the other hand, the lattice theory possesses the right151515Many phenomenological applications require matrix elements of operators of the electroweak Hamiltonian. These operators must also be constructed to the appropriate order in $1/m_{Q}$, cf. sect. 7. In particular, to first order in $1/m_{Q}$ the coefficient $d_{1}$ must be chosen according to eq. (7.10). $1/m_{Q}$ corrections only if $M_{2}=M_{B}=m_{Q}$. A computation with the Wilson action and $m_{Q}=M_{2}$ obtains spin-averaged features correctly, but underestimates the chromomagnetic $1/m_{Q}$ corrections. Compared to corrections from the kinetic energy, the spin dependent effects are thought to be small [17], so again the most essential adjustment is $M_{2}=m_{Q}$. A better computation with the Sheikholeslami-Wohlert action and $m_{Q}=M_{2}=M_{B}$ obtains all1515footnotemark: 15 $1/m_{Q}$ features correctly.

Despite the similarity between previous nonrelativistic field theories [8, 9, 10, 11, 12] and the view adopted in this section, there are significant technical differences. The four-component approach explicitly includes the terms $m_{0}\bar{\psi}(x)\psi(x)$ and terms that couple upper and lower components, such as {\bar{\psi}(x)\mbox{\boldmath\gamma}\!\cdot\!\mbox{\boldmathD}\psi(x)} and {\bar{\psi}(x)\mbox{\boldmath\alpha}\!\cdot\!\mbox{\boldmathE}\psi(x)}. The program of Lepage, et al, [12] omits these interactions in practice, though perhaps not in principle. There are advantages to leaving out the Dirac block–off-diagonal interactions. Fermion propagators are the solution of a (one-sweep) initial-value problem, whereas they are otherwise the solution of a boundary-value problem, solvable only by iteration; with fewer interactions, perturbation theory is easier [27]. On the other hand, these interactions are necessary to take $a\to 0$ (by brute force) without the scourge of power-law divergences, or to reach into the semirelativistic regime.

10 Numerical Tests and Applications

With a few examples this section tests the results of the previous sections with Monte Carlo data. All data were generated with the axis-interchange symmetric Wilson or Sheikholeslami-Wohlert actions, so as $m_{q}a$ increases, we rely on the nonrelativistic interpretation of sect. 9. The tests verify the most important lessons. The bare mass should be adjusted until the kinetic mass $M_{2}$, defined in eq. (4.7), takes the desired value. In particular, in an extrapolation to vanishing lattice spacing, one ought to hold the kinetic mass161616This can be done nonperturbatively with a meson instead of a quark state. fixed. On the other hand, the dynamically irrelevant rest mass may deviate from $m_{q}$. For matrix elements, it is also important to use the improved field $\Psi_{\rm I}$ of eq. (7.8) or (A.17). The factor $e^{M_{1}a/2}$ is more important than the bracket in eq. (7.8), because it guarantees a smooth approach to the static limit.171717Neglecting the bracket introduces only {\rm O}(m_{0}\hbox to0.0pt{\hskip 0.50003ptp\hss}p\hskip 0.50003pta^{2}) lattice artifacts at $m_{0}a\ll 1$, but {\rm O}(\hbox to0.0pt{\hskip 0.50003ptp\hss}p\hskip 0.50003pt/m_{q}) at m_{0}a\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{>}}1.

First, consider the mass spectrum, in particular the hyperfine splitting in heavy-light systems. By heavy-quark spin symmetry the vector-pseudoscalar mass difference ${m_{V}-m_{P}}$ is expected to be proportional to $1/m_{Q}$. Obviously the leading term in the sum ${m_{V}+m_{P}}$ is proportional to $m_{Q}$, so the combination ${m_{V}^{2}-m_{P}^{2}}$ should be nearly independent of $m_{Q}$. Numerical work [28, 29] found, however, that ${m_{V}^{2}-m_{P}^{2}}$ decreases for increasing $m_{Q}$, with lattice spacing $a$ held fixed. These analyses take $m_{V}$ and $m_{P}$ from the rest mass. From sect. 9 the rest mass $M_{1}$ of the quark governs the rest mass of the mesons, while the chromomagnetic mass $M_{B}$ governs the hyperfine splitting. The computed lattice quantity, therefore, is proportional to $M_{1}/M_{B}$, which decreases for increasing quark mass, cf. fig. 3. Given $M_{1}$, $M_{B}$ is not as large with the Sheikholeslami-Wohlert action as with the Wilson action. Numerical data with $c_{B}=1$ show behavior [29] qualitatively similar to $c_{B}=0$.

To improve the determination of ${m_{V}-m_{P}}$ one should tune to the kinetic mass instead of the rest mass and use the mean-field or one-loop estimate of $c_{B}$. The chosen value of $c_{B}$ could be tested in quarkonia. Then in heavy-light systems one could verify two predictions of heavy-quark symmetry, as applied to the lattice theory: a falling $M_{1}/M_{B}$ behavior when using the mesons’ rest masses, and a flat $M_{2}/M_{B}$ behavior when using the mesons’ kinetic masses.

Next, consider the improvement and normalization of multi-quark operators from sect. 7. The normalization of the vector current can be checked nonperturbatively. The fermion number

[TABLE]

counts the number of $h$-flavored fermions in $|H\rangle$. If $|H\rangle$ has one and only one $h$ in it, the condition $N_{h}=1$ defines the factor $2\kappa_{t}Z_{V}$.181818In the notation of sect. 8, $Z_{V}=Z_{\gamma_{\mu}}$ and $Z_{A}=Z_{\gamma_{\mu}\gamma_{5}}$. Fig. 4 compares this nonperturbative definition

of $2\kappa Z_{V}$ with the mean-field-improved, tree-level perturbative approximation. The symbols are from Monte Carlo calculations [30] of eq. (10.1), with $|H\rangle$ a meson with a spectator anti-quark of different flavor, and the solid curve is the mean-field improved, tree-level approximation $2\kappa(1+\widetilde{m}_{0h})$. Fig. 4 exhibits several interesting features. The solid curve accurately tracks the dominant mass dependence from $\widetilde{m}_{0}=0$ to $\widetilde{m}_{0}=\infty$. From eq. (8.10) one expects a subdominant mass dependence from loop corrections $Z^{[l]}_{V}(m_{0h}a)$. Indeed, near $\widetilde{m}_{0}=0$ the massless one-loop correction [31, 32] accounts quantitatively for the discrepancy, and near $\widetilde{m}_{0}=\infty$ the discrepancy becomes smaller, in accord with a Ward identity, which requires $2\kappa Z_{V}=1$ at infinite mass [33]. Neglecting the dominant mass dependence, as in the dashed curve, is obviously completely wrong for \widetilde{m}_{0}\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{>}}1.

Finally, consider the decay constant of a heavy-light pseudoscalar meson, computed with the local axial current $J^{ub}_{\mu 5}(x)=\bar{\psi}^{u}(x)\gamma_{\mu}\gamma_{5}\psi^{b}(x)$. Fig. 5 shows Monte Carlo data [34, 35]

at $\beta=5.7$ (for which $a^{-1}\approx 1$ GeV) for

[TABLE]

where $v_{\mu}$ is the meson’s four-velocity. (The vacuum and one-meson states are normalized to unity.) We have deliberately chosen a largish lattice spacing to enhance lattice artifacts, and thus test our control over them. We analyze the data two different ways. The lower set of points takes the meson mass from the rest mass and neglects the factor1818footnotemark: 18 $Z_{A}^{ub}=e^{M_{1b}a/2}$ in eq. (7.8). As suggested by the curve, the neglectful analysis would produce a locus of points that approaches zero in the static limit. The upper set of points uses the normalization factor, and—just as important—it defines the meson mass through a mean-field approximation to the kinetic mass.191919Ref. [34] provides hopping parameters and rest masses only, but the mean-field approximation is adequate for illustrative purposes. Fig. 5 shows how crucial both refinements are, if the Wilson-action data are to approach the static limit smoothly.

An important application of a plot like fig. 5 is to compute the slope in $1/m_{q}$. In the heavy-quark effective theory the slope is seen to arise from three sources: the kinetic energy, the chromomagnetic interaction, and a correction to the infinite-mass current [17]. The lattice theory has direct analogs: the kinetic energy requires tuning $M_{2}=m_{q}$, the chromomagnetic contribution requires tuning $M_{B}=m_{q}$, and the local lattice current requires a correction, given most compactly by eqs. (7.8) and (7.10). All three ingredients are needed to obtain the correct slope [36].

11 Conclusions

This paper (and conference reports [37, 18, 36] anticipating it) provides the foundations of a theory of lattice fermions, valid at any mass—large or small. The action starts with a set of interactions that encompasses those of both light-fermion actions [13, 6] and heavy-fermion actions [10, 11, 12]. The couplings of such a general action are then tuned in successive approximations to the renormalized trajectory. In applying renormalization-group techniques to analyze and reduce cutoff effects we do not, however, expand in either $m_{q}a$ or $\Lambda_{\rm QCD}/m_{q}$.

Although there are several methods for tuning the action (for example, refs. [38, 39]), our analysis is based on Symanzik-like, on-shell improvement criteria. This entails the computation of on-shell correlation functions, or, equivalently, of the Hamiltonian. Enforcing continuum-limit behavior—for example the relativistic mass shell—yields conditions on couplings of higher-dimension interactions. In practice, we here compute on-shell quantities in (tadpole-improved) perturbation theory.

An examination of the lattice theory’s Hamiltonian, derived from the transfer matrix, is especially illuminating. It shows that it is unnecessary to improve Wilson’s discretization of the time derivative. Instead higher-dimension interactions can be built from fields on one (or, in some cases, two adjacent) timeslices. Thus, our class of improved lattice actions automatically has an easy-to-construct transfer matrix. The actions, consequently, all automatically satisfy heavy-quark symmetries in the limit of large mass.

The Hamiltonian is a useful tool for examining lattice artifacts. A term of dimension $s+4$ factorizes as

[TABLE]

with dependence on the theory’s relevant couplings, $m_{q}a$ and $g^{2}$, in the coefficient $b$. From this expression it is plain how the artifacts behave as the mass increases: the absolute error induced by eq. (11.1) is \langle\delta\hat{H}_{\rm lat}\rangle\sim(\mbox{\boldmathp}a)^{s}\mbox{\boldmathp}, where $p$ is the typical three-momentum of the process. At small mass it is standard that the associated coefficient is a benign numerical factor $b(0)$. At large mass the lattice action’s heavy-quark symmetry implies that it is a (generally different) numerical factor $b(\infty)$. Our explicit results at tree level and our analysis of higher orders show that $b(m_{q}a)$ is a smooth, gentle function connecting the two extremes. In a nutshell, therefore, the characteristic measures of cutoff artifacts are $\Lambda_{\rm QCD}a$ and \mbox{\boldmathp}a, but never $m_{q}a$.

In general our actions have two hopping parameters. Then it is possible to maintain equality between the rest mass (energy at zero momentum, eq. (4.6)) and the kinetic mass (inertial response, eq. (4.7)). In nonrelativistic systems, however, there is a noteworthy simplification. Embracing the philosophy of the static [8] or nonrelativistic [10, 11] theories, one can ignore the rest mass and, hence, forgo one of the hopping parameters. The obvious application is to set them equal, as in the Wilson [3] and Sheikholeslami-Wohlert [6] actions. Therefore, the correct interpretation of numerical data generated with these actions at $m_{q}a\sim 1$ is a nonrelativistic one. In particular, the hopping parameter must be adjusted so that the kinetic mass agrees with the physical mass.

If one knows a priori that a quark is nonrelativistic, it is computationally cheaper to use a two-component formalism [12]. But there are many instances in which one would like to trace the mass dependence from the static limit down to, say, the strange quark. One example is fig. 5, which, with reliable calculations, should indicate how and where the heavy-quark expansion deteriorates.

The results of this paper can be extended in several ways. The couplings have been computed at tree level with mean-field improvement. One-loop calculations are desirable, and better still would be a nonperturbative determination, perhaps in a mass-dependent generalization of ref. [40]. Once one is confident that ${\rm O}(a)$ artifacts are under control, one can extend the analysis to dimension-six interactions. Tree-level, ${\rm O}(a^{2})$ improvement should be manageable; beyond tree level the bothersome four-fermion interactions enter the fray.

One would like to use the actions presented here in Monte Carlo calculations of QCD. If one uses the ${\rm O}(v^{4})$-improved action to compute the spectrum of charmonium and bottomonium, then, without more tuning, one could calculate properties of $D$ and $B$ mesons, including the electroweak matrix elements needed to determine the unknown elements of the Cabibbo-Kobayashi-Maskawa matrix.202020For relevant reviews, see ref. [41, 21]. Note that for these matrix elements, as well as for the quark mass in the ${\overline{\rm MS}}$ scheme, a small lattice spacing is helpful to reduce perturbative corrections. Once the lattice spacing is small enough so that m_{Q}a\mathop{\raisebox{-1.72218pt}{\hbox to0.0pt{\sim\hss}}\raisebox{1.72218pt}{<}}1, our formulation is especially advantageous. The two-component nonrelativistic theory breaks down as $m_{Q}a$ gets smaller [11], yet the old small-mass theory would have leading lattice artifacts of order $\alpha_{s}m_{Q}a$ and $(m_{Q}a)^{2}$. Our improved action, on the other hand, remains viable for any mass, and its cutoff effects are small, of order $\alpha_{s}\Lambda_{\rm QCD}a$ and $(\Lambda_{\rm QCD}a)^{2}$. To obtain comparable accuracy through brute force in the old theory, one would have to reduce the lattice spacing by a factor of $m_{Q}/\Lambda_{\rm QCD}$—about five for the charm quark. Even for a perfect algorithm the savings in computer time is, therefore, a factor of $5^{4}$.

Acknowledgements

We would like to thank Peter Lepage for numerous discussions and Bart Mertens for collaboration on the one-loop calculations mentioned in sect. 8. During this work we have also had useful conversations with Claude Bernard, Estia Eichten, Jim Labrenz, Martin Lüscher, Maarten Golterman, and Jim Simone. AXK would like to thank the Fermilab Theory Group for hospitality. Fermilab is operated by Universities Research Association, Inc., under contract DE-AC02-76CH03000 with the U.S. Department of Energy.

Appendix A Quarkonium to ${\rm O}(v^{4})$

This appendix extends the analysis of the main text to incorporate interactions that contribute in quarkonium through ${\rm O}(v^{4})$. Naively, this would entail close scrutiny of all interactions through {\rm O}(\mbox{\boldmathp}^{4}), i.e. up to dimension seven. Some dimension-seven bilinear interactions are listed in Table 5,

with their magnitude in quarkonium estimated by the velocity-counting rules of ref. [12]. Together with Table 3, one sees that not all dimension-six and -seven interactions are necessary to ${\rm O}(v^{4})$; this Appendix considers only the entries that are.

One must consider the dimension-six interaction

[TABLE]

with coupling $c_{\rm so}$, and the dimension-six and -seven interactions

[TABLE]

with couplings $c_{4D}$ and $c_{4L}$. Further dimension-six and -seven interactions contribute in ${\rm O}(v^{6})$ or higher [12]. We discuss the adjustment of $c_{\rm so}$ in eq. (A.1) in sect. A.1, and the adjustment of $c_{4D}$ and $c_{4L}$ in eq. (A.2) in sect. A.2.

The discretization of the covariant difference $D$ and Laplacian $\triangle^{(3)}$ also must be improved, to remove $\bar{\Psi}\gamma_{i}D_{i}^{3}\Psi$ and $\bar{\Psi}D_{i}^{4}\Psi$, respectively. These interactions break rotational invariance, and we treat them in sect. A.2.

A.1 Chromoelectric interactions to ${\rm O}(v^{4})$

It is easiest to treat the “spin-orbit” interaction $S_{\rm so}$ in the Hamiltonian formalism of sect. 5. The coupling $c_{\rm so}$ can only appear in a dimension-six term implied by the ellipsis in eq. (5.16). The spin-orbit Hamiltonian is

[TABLE]

Under the two-parameter Foldy-Wouthuysen-Tani transformation of eq. (5.23)

[TABLE]

With $b_{\rm so}$ it is possible to give an invariant involving $b_{E}$:

[TABLE]

and one wants $B_{E}=1$. Just as the redundancy associated with the Foldy-Wouthuysen-Tani parameter $\xi_{1}$ intertwines the mass dependence of $\zeta$ and $r_{s}$, the redundancy associated with the other Foldy-Wouthuysen-Tani parameter $\xi_{E}$ intertwines the mass dependence of $c_{E}$ and $c_{\rm so}$ too.

Assuming a discretization of $S_{\rm so}$ that resides on two timeslices only, it is straightforward to generalize the transfer-matrix construction to the action $S_{0}+S_{B}+S_{E}+S_{\rm so}$. After expanding the transfer matrix in powers of $a$ one finds

[TABLE]

Combining eqs. (A.5), (5.29), and (A.6) and setting $B_{E}=1$ yields the mass-dependent condition

[TABLE]

Here the rest mass $M_{1}$ has been eliminated in favor of the kinetic mass $M_{2}$ as appropriate to the nonrelativistic interpretation of sect. 9.

The redundant direction associated to $\xi_{E}$ permits a free choice of $c_{\rm so}$. In our framework, which stresses a smooth matching to the massless limit, the most convenient choice is probably $c_{\rm so}=0$. But for purely nonrelativistic applications ref. [12] would choose $c_{E}=0$ and $\zeta c_{\rm so}\propto m_{q}^{-2}$. Other possibilities correspond to the special choice for $r_{s}$ in eq. (5.35). Then eq. (A.7) reduces to

[TABLE]

The further special case corresponding to $\xi_{E}=0$ is $c_{E}=1$ (independent of $m_{0}a$) and $c_{\rm so}=0$.

A.2 Kinetic energy to ${\rm O}(v^{4})$

The interaction $S_{4}$ produces corrections to the kinetic energy. It is easiest to analyze from the energy-momentum relation, as in sect. 4. Expanding eq. (4.4) to {\rm O}(\mbox{\boldmathp}^{4}) yields

[TABLE]

where

[TABLE]

and

[TABLE]

The relativistic mass shell satisfies $M_{4}=M_{2}=M_{1}$, and a nonrelativistic mass shell with leading relativistic correction satisfies $M_{4}=M_{2}$. In both cases rotational invariance requires $w_{4}=0$.

A straightforward way is enforce $w_{4}=0$ is to take an improved covariant difference

[TABLE]

and an improved covariant Laplacian

[TABLE]

The coefficients are chosen so that the Fourier transforms have no $p_{i}^{3}$ or $\sum_{i}p_{i}^{4}$ terms, respectively. Then (at tree-level) $w_{4}=0$ automatically.

Since $S_{4}$ contains no higher time derivatives, the transfer-matrix construction proceeds as usual. After deriving $M_{4}$ for $S_{0}+S_{4}$ ($M_{1}$ and $M_{2}$ are unchanged), one finds $M_{4}=M_{2}$ at tree level if the couplings $c_{4D}$ and $c_{4L}$ obey

[TABLE]

This result holds whether $\zeta$ is tuned so that $M_{1}=M_{2}$ or not. As with $\zeta$ and $c_{E}$, only the massless limit of $c_{4D}$ is unambiguous. For the full mass dependence the redundant $c_{4L}$ must be specified, for example $c_{4L}=0$.

It is instructive to look explicitly at the consequences of omitting $S_{4}$ from Monte Carlo calculations. With unimproved $D$ and $\triangle^{(3)}$ and no $S_{4}$

[TABLE]

and

[TABLE]

The rotational-invariance breaking artifact is, thus, ${\rm O}(p_{i}^{4}a^{2}/m_{q})$ for $m_{q}a$ large and small. The rotationally invariant, relativity-breaking artifact is {\rm O}(\mbox{\boldmathp}^{4}a^{2}/m_{q}) at small $m_{q}a$, and {\rm O}(\mbox{\boldmathp}^{4}a^{1}/m_{q}^{2}) at large $m_{q}a$.

A.3 Electroweak operators

For electroweak decays of quarkonia to ${\rm O}(v^{4})$, one needs a higher-dimensional generalization of eq. (7.8),

[TABLE]

The $d$’s are easiest to derive from the Foldy-Wouthuysen-Tani transformed field. Combining eqs. (5.9) and (5.23)

[TABLE]

where $\xi_{1}$ and $\xi_{E}$ parameterize the solution of the tuning conditions. Expanding the cumbersome exponentials to {\rm O}(a^{2}[\mbox{\boldmathp}^{2},\mbox{\boldmathB},\mbox{\boldmathE}]) and eliminating $\xi_{1}$ and $\xi_{E}$ in favor of the couplings $\zeta$, $r_{s}$, $c_{B}$, and $c_{E}$, one finds

[TABLE]

and $d_{1}$ as in eq. (7.10). In eqs. (A.19) the kinetic mass $M_{2}$ has been substituted for the rest mass $M_{1}$. Thus, these formulae remain valid under the nonrelativistic interpretation explained in sect. 9.

Appendix B Combining Exponents

This appendix presents a proof of eq. (5.26), i.e. that the functions $f_{1}$ and $f_{2}$ are given by the expressions in eq. (5.28). If $H_{0}$, $H_{I}$, and $H_{I}^{\dagger}$ were to commute, it would be trivial to combine the exponents. They do not, so the combined exponent depends on their commutators as well. The commutators are

[TABLE]

To obtain the last commutator we have written \mbox{\boldmathD}_{U}=(\mbox{\boldmathD}-\mbox{\small\frac{1}{2}}\mbox{\boldmathE})_{\rm cont} and \mbox{\boldmathD}_{V}=(\mbox{\boldmathD}+\mbox{\small\frac{1}{2}}\mbox{\boldmathE})_{\rm cont}, and we have neglected higher powers of \mbox{\boldmathE}_{\rm cont}. The operators $\Theta$ and $\triangle^{(3)}$ carry one and two powers of \mbox{\boldmathp}a, respectively; thus {\cal M}=M_{1}-\mbox{\small\frac{1}{2}}r_{s}\zeta e^{-M_{1}}\triangle^{(3)}+{\rm O}(\mbox{\boldmathp}^{4}a^{4}). Note that although $[\hat{H}_{I},\hat{H}_{I}^{\dagger}]$ is {\rm O}(\mbox{\boldmathp}^{2}a^{2}), further commutators such as $[\hat{H}_{0},[\hat{H}_{I},\hat{H}_{I}^{\dagger}]]$ are at least {\rm O}(\mbox{\boldmathp}^{4}a^{4}).

Instead of solving the field theory, it is enough to consider a toy model with two degrees of freedom, a fermion (annihilated by $\hat{a}$) and an anti-fermion (annihilated by $\hat{b}$). With discrete time the action is

[TABLE]

The transfer matrix has the same form as eq. (5.25) with

[TABLE]

where $A=(1+m)^{1/2}a$, and $B=(1+m)^{1/2}b$. With the identification of $\vartheta$ with $\zeta\Theta$ and $m$ with $e^{\cal M}-1$, the operators $\hat{H}_{0}$ and $\hat{H}_{I}^{(\dagger)}$ of the toy model and the field theory have the same algebraic structure.

This model has only four states, the vacuum, fermion, anti-fermion, and a fermion–anti-fermion state. The strategy is to work out the transfer matrix elements explicitly, and then take the logarithm. These steps are easier in the Grassman-number approach, where the matrix elements of the transfer matrix are the coefficients of monomials in (Grassman numbers) $A$, $A^{\dagger}$, $B$, and $B^{\dagger}$, when

[TABLE]

is expressed as a polynomial. Up to factors analogous to $\det(2\kappa_{t}B)$, which we can drop without loss, the transfer matrix in the neutral sector is

[TABLE]

Writing ${\cal T}=VDV^{\dagger}$, where $D$ is diagonal, the Hamiltonian is $H=-V\log(D)V^{\dagger}$. Expanding the result to ${\rm O}(\vartheta^{2})$ one finds

[TABLE]

where $e^{m_{1}}=1+m$. Expressed in terms of Fock-space operators

[TABLE]

Substituting $m_{1}$ and $\vartheta$ for ${\cal M}$ and $\zeta\Theta$ completes the derivation of eq. (5.26).

Appendix C Spinors, Creation and Annihilation Operators

This Appendix gives the construction of spinors and of creation and annihilation operators in $d=4$ space-time dimensions. These are needed to calculate amplitudes of on-shell fermions via Feynman diagrams.

Consider an arbitrary bilinear fermion action

[TABLE]

with an implied sum over $\mu$. We assume that $\tilde{K}_{\mu}$ and $\tilde{L}$ are translation invariant. With parity $\tilde{L}(x,y)$ is symmetric, and $\tilde{K}_{\mu}(x,y)$ anti-symmetric, under interchange of $x$ and $y$. The field $\psi(x)$ has the following equation of motion

[TABLE]

In such a free theory, one searches for solutions of the form \psi(x)=e^{ip_{0}t+i\hbox to0.0pt{\hskip 0.35002pt\scriptstyle p\hss}\scriptstyle p\hskip 0.35002pt\cdot\hbox to0.0pt{\hskip 0.35002pt\scriptstyle x\hss}\scriptstyle x\hskip 0.35002pt}u(p). The four-component spinor $u(p)$ must satisfy

[TABLE]

The Fourier transforms $K_{\mu}(p)$ and $L(p)$ are real functions of $p$ by parity and translation invariance. Multiplying by $-i\gamma_{\mu}K_{\mu}(p)+L(p)$ one sees that solutions exist only if

[TABLE]

This mass shell coincides with the one derived from the propagator, as in sect. 4.

The actions that we consider all have the Wilson time derivative. Writing p=(p_{0},\mbox{\boldmathp}), one then has

[TABLE]

K_{i}(\mbox{\boldmathp}) is an odd, and \mu(\mbox{\boldmathp}) an even, function of $p$. Thus, solutions exist only when

[TABLE]

where

[TABLE]

It is convenient to label the solutions of eq. (C.3) by the sign in eq. (C.6). For each sign there are two solutions, u(\pm 1,\mbox{\boldmathp}) and u(\pm 2,\mbox{\boldmathp}). Setting \mbox{\boldmathp}={\bf 0} the equation of motion simplifies to

[TABLE]

where $M_{1}=E({\bf 0})=\log[\mu({\bf 0})]$. Choosing $\gamma_{0}$ as in eq. (2.6) the four solutions at \mbox{\boldmathp}={\bf 0} are

[TABLE]

where the subscript is the Dirac index, and all other components are zero. Direct substitution verifies that

[TABLE]

solves eq. (C.3) for \mbox{\boldmathp}\neq{\bf 0}, if $\sin p_{0}=i\mathop{\rm sign}\xi\sinh E$ and L=\mu(\mbox{\boldmathp})-\cosh E(\mbox{\boldmathp}). The denominator yields the normalization convention in eq.(C.14).

The two solutions with “negative energy” ($\xi<0$) correspond to anti-particle states. As usual we introduce

[TABLE]

The spinors $v$ obey the equation of motion

[TABLE]

which is solved by

[TABLE]

now with $\sin p_{0}=+i\sinh E$. From now on we shall use $u$ and $v$ with $\xi\in\{1,2\}$ and $\sin p_{0}=+i\sinh E$ only.

The spinors obey the conventional orthonormality properties

[TABLE]

where $\bar{u}=u^{\dagger}\gamma_{0}$ and $\bar{v}=v^{\dagger}\gamma_{0}$. Moreover,

[TABLE]

In a relativistic theory $E/m$ would appear here.

The general solution to eq. (C.2) is a linear superposition

[TABLE]

with $\sin p_{0}=i\sinh E$. The normalization factor {\cal N}(\mbox{\boldmathp}) is fixed below, after invoking this expansion for Hilbert-space operators. The operator-valued expansion coefficients \hat{b}^{\dagger}(\xi,\mbox{\boldmathp}) and \hat{d}^{\dagger}(\xi,\mbox{\boldmathp}) create particle and anti-particle states respectively:

[TABLE]

where $|0\rangle$ is the Fock state annihilated by all \hat{b}(\xi,\mbox{\boldmathp}) and \hat{d}(\xi,\mbox{\boldmathp}). Assuming the vacuum is normalized to $\langle 0|0\rangle=1$, the fermion states are normalized to

[TABLE]

and similarly for the anti-fermion state |\bar{q}(\xi,\mbox{\boldmathp})\rangle, if and only if the anti-commutator

[TABLE]

and similarly for \{\hat{d}(\xi^{\prime},\mbox{\boldmathp}^{\prime}),\hat{d}^{\dagger}(\xi,\mbox{\boldmathp})\}.

The transfer-matrix construction provides the anti-commutation relation for \hat{\psi}(t,\mbox{\boldmathx}) and {\hat{\bar{\psi}\mkern 5.0mu}\mkern-5.0mu}{}(t,\mbox{\boldmathx}). Eq. (5.9) becomes

[TABLE]

After inverting the Fourier series and evaluating the anti-commutators one finds

[TABLE]

The convention $\phi=1$ is the most convenient.212121Another candidate is the relativistic convention, which is not at all natural in nonperturbative calculations. How can one normalize to \phi=\sqrt{m^{2}+\hbox to0.0pt{\hskip 0.50003ptp\hss}p\hskip 0.50003pt^{2}}/m without solving the theory? Even here, in a footnote to an appendix, one may not forget that the aim of these perturbative calculations is to understand, interpret, and improve the nonperturbative calculations. Then

[TABLE]

Note that ${\cal N}({\bf 0})=e^{-M_{1}/2}$.

With all this machinery we can now state the main result of this appendix. The Feynman rules for vertices in standard references (e.g. ref. [42] for the Wilson action) are derived from the functional integral, i.e. using $\psi(x)$. To obtain on-shell matrix elements, these rules must be supplemented by rules for contractions between $\psi(x)$ and conventionally normalized external states. They are

[TABLE]

multiplied by the sign appropriate to the anti-commutation implied by the $\cdots\,$. For all states the momentum flow is physical, i.e. with (against) the charge flow for particles (anti-particles).

For the specific action discussed in this paper, $S_{0}$, one finds (restoring $a$)

[TABLE]

The chromomagnetic and chromoelectric interactions do not modify these functions, but the kinetic corrections in Appendix A do.

It is useful to record the small $p$ expansion of the external line factor here:

[TABLE]

where the subscript “lat” abbreviates “lattice,” and $M_{X}$ is an “external line mass.” For $S_{0}$

[TABLE]

For a unified treatment of fermions and anti-fermions in initial and final states, it is handy to note that eq. (C.25) holds for positive and negative $\xi$. The analogous expression for {\cal N}(\mbox{\boldmathp})v_{\rm lat}(\xi,\mbox{\boldmathp}) then follows from eq. (C.11).

Unless $m_{0}a\ll 1$ the lattice external line factor {\cal N}(\mbox{\boldmathp})u_{\rm lat}(\xi,\mbox{\boldmathp}) deviates from the relativistic one. With our normalization conventions the relativistic analog of eq. (C.25) is

[TABLE]

where the subscript “rel” abbreviates “relativistic.” The rotations in sects. 7 and A.3 are needed to convert the bracket of eq. (C.25) into the bracket of eq. (C.27), assuming $m_{q}=M_{2}$.