%Title: High-Precision Charm-Quark Mass from Current-Current Correlators in Lattice and Continuum QCD
%Author: I. Allison, E. Dalgic, C.T.H. Davies, E. Follana, R.R. Horgan, K. Hornbostel, G.P. Lepage, J. Shigemitsu, H. Trottier, R.M. Woloshyn, K.G. Chetyrkin, J.H. Kuhn, M. Steinhauser, C. Sturm
%Published: * Phys.Rev. * ** D78: ** (2008.) 054513.
%arxiv/0805.2999
\documentclass[aps,prd,superscriptaddress,amsmath,groupedaddress,
twocolumn,showpacs]{revtex4}
\usepackage{amssymb,graphicx}
\usepackage{dcolumn}
\newcommand{\msb}{{\rm\overline{MS}}}
\newcommand{\eq}[1]{Eq.~(\ref{#1})}
\newcommand{\order}{{\cal O}}
\newcommand{\xv}{{\bf x}}
\newcommand{\psib}{\overline{\psi}}
\newcommand{\g}{\gamma}
\newcommand{\lmc}{l_{ m_c}}
\begin{document}
\title{High-Precision Charm-Quark Mass
from Current-Current Correlators\\ in Lattice and Continuum
QCD}
\author{I.\ Allison}
\affiliation{TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T 2A3,
Canada}
\author{E.\ Dalgic}
\affiliation{Physics Department, Simon Fraser University, Vancouver,
British Columbia, Canada}
\author{C.\ T.\ H.\ Davies}
\affiliation{Department of Physics and Astronomy, University of
Glasgow, Glasgow G12 8QQ, UK}
\author{E.\ Follana}
\affiliation{Physics Department, The Ohio State University, Columbus,
Ohio 43210, USA}
\author{R.\ R.\ Horgan}
\affiliation{
Department of Applied Mathematics and Theoretical Physics, Cambridge
University, Wilberforce Road, Cambridge CB3 0WA, UK}
\author{K.\ Hornbostel}
\affiliation{Southern Methodist University, Dallas, Texas 75275, USA}
\author{G.\ P.\ Lepage}
\email{g.p.lepage@cornell.edu}
\affiliation{Laboratory for Elementary-Particle Physics, Cornell
University, Ithaca, NY 14853, USA}
\author{J.\ Shigemitsu}
\affiliation{Physics Department, The Ohio State University, Columbus,
Ohio 43210, USA}
\author{H.\ Trottier}
\affiliation{Physics Department, Simon Fraser University, Vancouver,
British Columbia, Canada}
\author{R.\ M.\ Woloshyn}
\affiliation{TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T 2A3,
Canada}
\collaboration{HPQCD Collaboration}
\noaffiliation
\author{K.\ G.\ Chetyrkin}
\author{J.\ H.\ K\"uhn}
\author{M.\ Steinhauser}
\affiliation{Institut f\"ur Theoretische Teilchenphysik,
Universit\"at Karlsruhe, D-76128 Karlsruhe, Germany}
\author{C.\ Sturm}
\affiliation{Physics Department,
Brookhaven National Laboratory,
Upton, New York 11973, U.S.A.
}
\date{April 24, 2008}
\pacs{11.15.Ha,12.38.Aw,12.38.Gc}
\preprint{ BNL-HET-08/11}
\preprint{ TTP08-18}
\preprint{ SFB/CPP-08-28}
\begin{abstract}
We use lattice QCD simulations, with MILC configurations and HISQ
$c$-quark propagators,
to make very precise determinations of moments
of charm-quark pseudoscalar, vector and axial-vector correlators.
These moments are combined with new four-loop results from continuum
perturbation
theory to obtain several new determinations of the $\msb$ mass of the
charm quark. We find
$m_c(3\,\mathrm{GeV})=0.984\,(16)$\,GeV, or, equivalently,
$m_c(m_c)=1.266\,(14)$\,GeV.
This agrees well with results from
continuum analyses of the vector correlator using experimental data for
$e^+e^-$~annihilation (instead of using lattice QCD simulations). These
lattice and
continuum results are the most accurate determinations to date of this
mass. We
also obtain a new result for the QCD coupling:
$\alpha_\msb^{(n_f=4)}(3\,\mathrm{GeV}) =
0.230\,(18)$, or, equivalently, $\alpha_\msb^{(n_f=5)}(M_Z) =
0.113\,(4)$.
\end{abstract}
\maketitle
\section{Introduction}
Precise values for the charm quark's mass~$m_c$ are important for
high-precision tests of the Standard Model. Some of the most accurate
determinations currently come from zero-momentum moments of
current-current correlators built from the $c$~quark's electromagnetic
current (see, for example, \cite{Kuhn:2001dm,Kuhn:2007vp}). Low moments
are perturbative and have long been known through three-loop
order~\cite{Chetyrkin:1995ii,Chetyrkin:1996cf,Chetyrkin:1997mb}. New
techniques have recently extended these results to much higher
moments~\cite{Boughezal:2006uu,Maier:2007yn} and, in some cases, to
four-loop order~\cite{Chetyrkin:2006xg,Boughezal:2006px}. These moments
can be estimated nonperturbatively, using dispersion relations, from
experimental data for the electron-positron annihilation cross section,
$\sigma(e^+ e^-\to\gamma^*\to X)$. The $c$~quark's mass is extracted by
comparing the perturbative and experimental determinations.
In this paper we show how to compute such moments directly using
accurately tuned, highly realistic numerical simulations of QCD in the
lattice approximation~\cite{Bochkarev:1995ai}. As we will show, the
correlator moments obtained nonperturbatively from such simulations can
be used in place of data from $e^+ e^-$~annihilation to obtain new,
few-percent accurate determinations of the $c$~quark's mass. With
lattice QCD, it is also possible to replace the electromagnetic current
in the correlator by the pseudoscalar operator $m_c\psib\gamma_5\psi$,
thereby providing a completely new set of mass determinations and an
important cross check on the entire methodology. In fact our most
accurate results come from the pseudoscalar correlators.
In Section~\ref{corr-from-latt} we describe how to compute pseudoscalar
correlators and their moments using lattice QCD. We discuss techniques
for reducing lattice artifacts in Section~\ref{sec:red-mom}, and
present new determinations of~$m_c(\mu)$ from our lattice ``data'' in
Section~\ref{sec:mc}. Here we also present a new determination of the
QCD coupling. In Section~\ref{sec:other} we extend our analysis to
include vector and axial-vector correlators. We summarize our main
results in Section~\ref{concl}. In the Appendix we review the continuum
perturbation theory needed for this analysis, including new four-loop
results for the pseudoscalar case.
\section{Lattice QCD and Pseudoscalar
Correlators}\label{corr-from-latt}
Few-percent accurate QCD simulations have only become possible in the
last few years (see, for example,~\cite{Davies:2003ik,Mason:2005zx}),
and accurate simulations of relativistic $c$~quarks only in the past
year\,---\,with the new Highly Improved Staggered Quark (HISQ)
discretization of the quark
action~\cite{Follana:2006rc,Follana:2007uv}, which we use here. A
lattice QCD simulation proceeds in two steps. First the QCD
parameters\,---\,the bare coupling constant and bare quark masses in
the Lagrangian\,---\,must be tuned. Then the tuned simulation is used
to compute vacuum matrix elements of various quantum operators from
which physics is extracted. An obvious approach to the tuning is to
choose a lattice spacing~$a$, and then tune each of the QCD parameters
so that the simulation reproduces the experimental value for a
corresponding physical quantity that is well measured. It is more
efficient, however, to first choose a value for the bare coupling and
then adjust the lattice spacing and bare masses to give physical
results.
In the simulations used here, we set the lattice spacing to reproduce
the correct $\Upsilon^\prime-\Upsilon$~meson mass difference in the
simulations, while we tuned the $u/d$, $s$, $c$ and~$b$ masses to give
correct values for $m_\pi^2$, $2m^2_K-m_\pi^2$, $m_{\eta_c}$, and
$m_\Upsilon$, respectively. (For efficiency we set $m_u=m_d$; this
leads to negligible errors in the analysis presented here.) The
important parameters for the particular simulations used in this paper
are listed in Table~\ref{tab:param}; further details can be found
in~\cite{Follana:2007uv,Davies:2003ik}. Once these parameters are set,
there are no further physics parameters, and the simulation will
accurately reproduce QCD physics for momenta much smaller than the
ultraviolet~(UV) cutoff~($\Lambda\sim\pi/a$). We have tested these
simulations extensively (see, for example,
\cite{Davies:2003ik,Mason:2005zx}) and, in particular, we have done
very precise tests for the charm-quark physics most relevant to this
work. These demonstrate, for example, that our simulations reproduce
the low-lying spectrum, including spin structure, of both charmonium
and heavy-light mesons ($D$ and $D_s$) to within our simulation
uncertainties (a few percent or
less)~\cite{Follana:2006rc,Follana:2007uv}.
\begin{table}
\caption{Parameters for the QCD simulations used in this paper. The
inverse lattice spacing $a^{-1}$ is in units of~GeV; $L$ and $T$
are the
spatial and temporal size of the lattices used. The $u$~and $d$~masses
are set
equal to $m_{u/d}$. The configurations used here were generated by
the MILC
collaboration~\cite{MILC}. $N_\mathrm{cfg}$ is the number of
configurations
used in each case; we used multiple time origins and, in some cases,
random wall sources to
render statistical errors negligible.}
\label{tab:param}
\begin{center}
\begin{ruledtabular}\begin{tabular}{ccccccc}
$a^{-1}$ & $am_{0u/d}$ & $am_{0s}$ & $am_{0c}$ & $L/a$ & $T/a$ &
$N_\mathrm{cfg}$ \\ \hline
1.31 & 0.013 & 0.066 & 0.850 & 16 & 48 & 631 \\
1.60 & 0.014 & 0.055 & 0.660 & 20 & 64 & 595 \\
2.26 & 0.007 & 0.037 & 0.430 & 28 & 96 & 566 \\
\end{tabular}\end{ruledtabular}
\end{center}
\end{table}
Given a tuned simulation, it is straightforward to calculate
correlators of the sort used to determine~$m_c$. The simplest of these
is for the $c$~quark's pseudoscalar density,
$j_5\equiv{\psib}_{c}\gamma_5\psi_c$:
\begin{equation}
G(t) \equiv a^6\,\sum_\xv (am_{0c})^2 \langle0| j_{5}(\xv,t)
j_{5}(0,0)
|0\rangle
\end{equation}
where $m_{0c}$ is the $c$~quark's bare mass (in the lattice
Lagrangian). Here time~$t$ is euclidean, and the sum over spatial
position~$\xv$ sets the total three momentum to zero. Note that
$G(t)=G(T-t)=G(T+t)$ where $T$ is the temporal length of the lattice.
We include two factors of $am_{0c}$ in the definition of $G(t)$ so that
$G(t)$ becomes independent of the UV cutoff as~$a\to 0$. Consequently
the lattice and continuum $G(t)$s become equal in this limit. Moments
$G_n$ are trivially computed:
\begin{equation} \label{Gn-def}
G_n \equiv \sum_t (t/a)^n G(t),
\end{equation}
where, on our periodic lattice~\cite{footnote:0},
\begin{equation}
t/a \in \{0,1,2\,\ldots\,T/2a-1,0,-T/2a+1\,\ldots\,-2,-1\}.
\end{equation}
The cutoff independence of $G(t)$ implies that
\begin{equation}\label{gn-def}
G_n = \frac{g_n(\alpha_\msb(\mu),\mu/m_c)}{(am_c(\mu))^{n-4}}
\end{equation}
for $n\ge4$, where $m_c(\mu)$ is the $\msb$ mass at scale~$\mu$ and
$g_n$ is dimensionless. The $c$~mass can be determined from moments
with $n\ge6$ given $G_n$ from lattice simulations and $g_n$ from
perturbation theory (see Appendix). This assumes that perturbation
theory is applicable, which should be the case for small enough~$n$.
Note that here and elsewhere in this paper we omit annihilation
contributions from $c\bar{c}\to\mathrm{gluons}\to c\bar{c}$. This is
allowed provided the same contributions are
omitted from perturbation theory, which we do. Annihilation
contributions to the nonperturbative
part of our analysis would be negligible in any case (for example, they
shift the $\eta_c$~mass by approximately~2.4\,MeV, which is less
than~0.1\%~\cite{Follana:2006rc}).
\section{Reduced Moments}\label{sec:red-mom}
The two biggest challenges for lattice QCD in producing these moments
lie in controlling: 1) $\order((am_c)^n)$ errors caused by the lattice
approximation; and 2) tuning errors in the QCD parameters, and
especially in the lattice spacing and the $c$~quark's bare mass. Each
of these potential sources of error is reduced by replacing $G_n$ by a
reduced moment:
\begin{equation}\label{Rn-def}
R_n \equiv \left\{
\begin{aligned}
& G_4/G^{(0)}_4 & & \text{for $n=4$,} \\
&\frac{a m_{\eta_c}}{2a m_{0c}}
\left(G_n/G^{(0)}_n\right)^{1/(n-4)}
& & \text{for $n\ge6$,}
\end{aligned}\right.
\end{equation}
where $G^{(0)}_n$ is the $n^\mathrm{th}$~moment of the correlator to
lowest order in lattice QCD perturbation theory, $a m_{\eta_c}$ is the
$\eta_c$~meson mass (in lattice units) from the simulation, and $a
m_{0c}$ is the bare lattice $c$~mass (in lattice units) used in the
simulation. The reduced moments can again be written in terms of
continuum quantities:
\begin{equation}\label{Rn-eqn}
R_n \equiv \left\{
\begin{aligned}
& r_4(\alpha_\msb,\mu/m_c) & & \text{for $n=4$,} \\
& \frac{r_n(\alpha_\msb,\mu/m_c)}{2 m_c(\mu)/m_{\eta_c}}
& & \text{for $n\ge6$,}
\end{aligned}\right.
\end{equation}
where $r_n$ is obtained from $g_n$ (\eq{gn-def}) and its value,
$g_n^{(0)}$, in lowest-order continuum perturbation theory:
\begin{equation}\label{rn-def}
r_n = \begin{cases} g_4/g_4^{(0)} & \text{for $n=4$,}\\
\left(g_n/g_n^{(0)}\right)^{1/(n-4)} & \text{for
$n\ge6$.}
\end{cases}
\end{equation}
The $c$~mass is obtained from \eq{Rn-eqn} with $n\ge6$ using the
nonperturbative lattice QCD (LQCD) value for $R_n$, the perturbative
QCD (PQCD) estimate for $r_n$, and the experimental value for
$m_{\eta_c}$,~2.980\,GeV:
\begin{equation}\label{extract-mc}
m_c(\mu) = \frac{m_{\eta_c}^\mathrm{exp}}{2}\,
\frac{r_n^\mathrm{PQCD}}{R_n^\mathrm{LQCD}}.
\end{equation}
The $G^{(0)}_n$~factor in $R_n$ removes the lattice spacing~$a$ from
the reduced moment, and cancels all tree-level $\order(a^n)$~errors,
thereby reducing finite-$a$ errors overall by about a factor of
three~\cite{footnote:1}. The $am_{\eta_c}$~factor reduces $R_n$'s
sensitivity to slight mistunings of the $c$ quark's bare mass, again by
a factor of three. The sensitivity to~$m_{0c}$ is reduced because $m_c$
now enters in the ratio $m_c/m_{\eta_c}$, so small errors in the
simulation parameter $m_{0c}$ are mostly cancelled by corresponding
shifts in the simulation value for~$m_{\eta_c}$~\cite{footnote:2}.
\begin{figure}
\includegraphics[scale=1.0]{Rn}
\caption{Reduced moments $R_n$ from lattice simulations with
different lattice spacings~$a$. The dashed lines show the functions
used to fit the lattice results. These extrapolation functions were
used to obtain the $a=0$~results shown in the plot.}
\label{fig:rn}
\end{figure}
\begin{table}
\caption{Simulation results for $R_n(a)$ for different lattice
spacings~$a$
and moments~$n$. The inverse lattice spacing is in~GeV.
Extrapolations to $a=0$
are given for each~$n$, together with the
corresponding value (in GeV) for $m_c(\mu)$ when $\mu=3$\,GeV.}
\label{tab:rn}
\begin{center}
\begin{ruledtabular}\begin{tabular}{cccc|cc}
$a^{-1}$ & 1.31 & 1.60 & 2.26 & $R_n(a=0)$ & $m_c(\mu)$ \\\hline
$R_{4}$ & 1.163(1) & 1.187(1) & 1.223(1) & 1.248(25) & -- \\
$R_{6}$ & 1.424(3) & 1.481(3) & 1.510(3) & 1.522(17) & 0.990(16) \\
$R_{8}$ & 1.351(2) & 1.378(2) & 1.373(2) & 1.367(11) & 0.980(16) \\
$R_{10}$ & 1.309(2) & 1.319(2) & 1.305(2) & 1.298(8) & 0.979(19) \\
$R_{12}$ & 1.274(2) & 1.276(2) & 1.262(2) & 1.256(6) & 0.978(23) \\
$R_{14}$ & 1.245(2) & 1.244(2) & 1.231(2) & 1.226(5) & 0.978(28) \\
$R_{16}$ & 1.221(2) & 1.219(2) & 1.207(2) & 1.202(4) & 0.979(33) \\
$R_{18}$ & 1.201(2) & 1.199(2) & 1.187(2) & 1.183(4) & 0.980(39) \\
\end{tabular}\end{ruledtabular}
\end{center}
\end{table}
Our simulation results for $R_n(a)$ are listed for different
moments~$n$ and lattice spacings~$a$ in Table~\ref{tab:rn}. The
uncertainties quoted for each of these values come from two sources.
First, for sake of efficiency, we used $u/d$~sea-quark masses that were
about 4--5~times too large in our simulations
($m_{u/d}/m_s\approx0.2$). Simulations with sea masses that were half
and twice this value shift our results by less than our $0.2$\%
simulation errors~\cite{footnote:2.1}. This error in the moments can
also be analyzed using perturbation theory. It implies a fractional
error of order $\alpha_s^2(\mu) (m_{u/d}/m_c)^2$, which is much smaller
than~0.1\%\,---\,much too small to justify the effort of extrapolating
in the sea masses.
A related concern is the contribution from $c$-quark vacuum
polarization, which we omit from the simulation. The importance of
$c$-quark loops in radiative corrections to the moments can be
estimated using perturbation theory. We find that the $c$-loops add
0.7\% to~$R_4$ but less than 0.1\% to $R_6$ and still less for higher
moments~\cite{footnote:2.2}. We have corrected our $R_4$ values to
account for $c$ loops, and we include an uncertainty of 0.1\% in our
other reduced moments to allow for any residual uncertainty due to
sea-quark masses or the absence of $c$~quarks in the sea.
The second major source of uncertainty, significant for $n\ge6$ only,
comes from residual errors due to mistuning of~$m_{0c}$ in the
simulation. This results in relative errors of order $0.2\%$, and an
overall tuning error of order~$0.5\%$~\cite{footnote:3}. We include the
former in the errors reported in the Table~\ref{tab:rn} for~$R_n(a)$,
but include the latter only in our final results (since it affects all
$R_n(a)$s by the same amount). Other sources of uncertainty in
the~$R_n$, like statistics, are negligible compared to these errors.
The reduced moments from the simulation are, as expected, only weakly
dependent (few percent) upon the lattice spacing for the lattice
spacings we use. To correct for this variation, we extrapolate each
moment to $a=0$ by fitting the simulation values to a function of the
form
\begin{align}
R_n(a) = R_n(0)\left(1\right. &+ c_0(ap_\mathrm{av})^2\alpha_s
+ c_1(am_c)^2(ap_\mathrm{av})^2\alpha_s\nonumber\\ \nonumber
&+ c_2 (am_c)^4\alpha_s + c_3 (am_c)^6\alpha_s \\
& \left. + c_4 (am_c)^8\alpha_s + \cdots\right)
\end{align}
where $p_\mathrm{av}\approx7m_c/n$ is the typical three momentum
carried by the $c$~quark in the lowest-order perturbation theory
diagram for $R_n$~\cite{footnote:4}.
The extrapolated results depend only weakly on the exact functional
form used here provided reasonable Bayesian priors are used for the
coefficients~$c_i$. We use the same Gaussian prior, centered at zero
with width~$\sigma_c$, for all of the $c_i$ and all values of $n$; we
tune $\sigma_c$ to roughly twice the optimal value indicated by the
empirical Bayes procedure described in~\cite{Lepage:2001ym}. (This
tuning gives~$\sigma_c=1/6$.) This is quite conservative, and leads to
uncertainties in our final extrapolated results, $R_n(a=0)$ (see
Table~\ref{tab:rn}), that are of order the amount of the extrapolation
from the finest lattice spacing to~$a=0$. Our extrapolations for the
$n=6,8,10$~moments are shown in Fig.~\ref{fig:rn}.
\section{Extracting $m_c(\mu)$}\label{sec:mc}
To convert the extrapolated reduced moments into $c$~masses, we require
perturbative expansions for the~$r_n$ in \eq{extract-mc}. These are
easily computed from the expansions for
$g_n$~\cite{Chetyrkin:1995ii,Chetyrkin:1996cf,Chetyrkin:1997mb,Boughezal:2006uu,Maier:2007yn,Chetyrkin:2006xg,Boughezal:2006px} using \eq{rn-def}; details can be found in the~Appendix. The perturbative expansions have the form
\begin{equation} \label{rn-exp}
r_n = 1 + r_{n,1}\alpha_\msb(\mu) + r_{n,2}\alpha_\msb^2
+ r_{n,3}\alpha_\msb^3 + \ldots
\end{equation}
where we set the renormalization scale~$\mu$
to~3\,GeV~\cite{footnote:4.5}.
The full third-order coefficients for the $n=4,6$ moments were computed
for this analysis and are presented in the Appendix. The third-order
coefficients for moments with $n\ge8$ are only partially complete: our
analysis includes all $\mu$-dependent terms (that is, $\log^n(\mu/m_c)$
terms), but the constant parts have not yet been computed. Consequently
we take the truncation uncertainty in $r_n$ to be of
order~\cite{footnote:5}
\begin{equation}\label{rn-err}
\sigma_{r_n} = \begin{cases}
r_n^\mathrm{max}\alpha_\msb^4(\mu) & \text{for $n=4,6$,}
\\[2ex]
r_n^\mathrm{max}\alpha_\msb^3(\mu) & \text{for $n\ge8$,}
\end{cases}
\end{equation}
where
\begin{equation}
r_n^\mathrm{max} = \mathrm{max}\left(|r_{n,1}|,|r_{n,2}|,|r_{n,3}|
\right).
\end{equation}
In this paper we use $\alpha_\msb(3\,\mathrm{GeV})=0.252\,(10)$ for the
($n_f=4$) coupling constant. We derived this from the current Particle
Data Group average for the $n_f=5$~coupling at $\mu=M_Z$, which
is~0.1176\,(20)~\cite{pdg}.
Note that coefficients in the $r_n$~expansion, \eq{rn-exp}, depend
upon~$m_c(\mu)$ through scale-dependent
logarithms,~$\log^n(\mu/m_c(\mu))$. Consequently \eq{extract-mc} is an
implicit equation for~$m_c(\mu)$, since the mass appears on both sides
of the equation. The $m_c(\mu)$-dependence on the right-hand side,
however, is suppressed by~$\alpha_\msb(\mu)$, and therefore the
equation is easily solved numerically.
Our final results for the $c$-quark's mass, $m_c(\mu)$ at
$\mu=3\,\mathrm{GeV}$ in the $\msb$~scheme, are listed in
Table~\ref{tab:rn}, and plotted in the upper-left panel of
Figure~\ref{fig:mc}. As is clear from the figure, all moments agree on
the mass although the higher moments may be less trustworthy
(see~\cite{Kuhn:2007vp}). We obtain our final result from the first
three moments ($n=6,8,10$), which average to give:
\begin{equation} \label{eq:ps-mc}
m_c(3\,\mathrm{GeV}) = 0.984\,(16)\,\mathrm{GeV}
\end{equation}
Evolving down to scale $\mu=m_c(\mu)$ using fourth-order
evolution~\cite{van
Ritbergen:1997va,Czakon:2004bu,Chetyrkin:1997dh,Vermaseren:1997fq},
this is equivalent to~\cite{mc-footnote}
\begin{equation}
m_c(m_c) = 1.266\,(14)\,\mathrm{GeV}.
\end{equation}
\begin{figure}
\includegraphics[scale=1.0]{pp-mc}\includegraphics[scale=1.0]{aa-mc}\\
\includegraphics[scale=1.0]{vv-mu-mc}\includegraphics[scale=1.0]{vv-1-mc}
\caption{$m_c(\mu=3\,\mathrm{GeV})$ from different moments of
correlators built
from four different lattice operators. The gray band
is our final result for the mass, 0.984\,(16)\,GeV, which comes
from the first three
moments of the pseudoscalar correlator (upper-left panel).}
\label{fig:mc}
\end{figure}
The leading sources of uncertainty in $m_c(\mu)$ are listed in
Table~\ref{tab:err} for each of the four lowest moments. Most of the
systematic changes with $n$ in these uncertainties are due to the fact
that the moments become increasingly infrared as $n$~increases.
Consequently, we are not surprised to see uncertainties related to the
convergence of perturbation theory grow as $n$~grows, while
uncertainties due to the finite lattice spacing~$a$ decrease. Similarly
nonperturbative effects in the continuum part of our analysis are
expected to be negligible for the lower moments ($n\le10$ at
least~\cite{Kuhn:2007vp}), but could become important for larger~$n$ as
confinement turns on.
To test for nonperturbative effects, we examined the leading
gluon-condensate contribution to our
analysis~\cite{Novikov:1977dq,Broadhurst:1994qj,Kuhn:2007vp}. We did
not correct the central
values of our masses for condensate contributions because the
condensate's size is not well known; but we did add a contribution to
the uncertainty for each mass that covers the current range of possible
condensate contributions~\cite{cond-footnote}. This has negligible
effect on the moments with~$n\le12$ but becomes quite large as
$n$~increases.
There are also uncertainties due to the finite spatial volume of our
lattices; our lattices were approximately 2.5\,fm~across. While our
simulations showed no measurable volume dependence~\cite{footnote:2.1},
lattice perturbation theory shows finite-volume sensitivity for the
higher (more infrared) moments. This is negligible for lower moments
but grows with~$n$. The finite-volume sensitivity is mostly an artifact
of perturbation theory; confinement significantly reduces finite-volume
effects. Consequently we assign a finite-volume error to our
perturbative factors that is equal to the entire finite-volume
correction in perturbation theory.
\begin{table}
\caption{Sources of uncertainty in the determinations
of $m_c(\mu=3\,\mathrm{GeV})$ from different reduced moments
$R_n$ of the pseudoscalar correlator. The uncertainties
listed are percentages of the final result~0.984\,(16)\,GeV.}
\label{tab:err}
\begin{center}
\begin{ruledtabular}\begin{tabular}{lddddd}
& \multicolumn{1}{c}{$R_6$}
& \multicolumn{1}{c}{$R_8$}
& \multicolumn{1}{c}{$R_{10}$}
& \multicolumn{1}{c}{$R_{12}$} \\ \hline
$a^2$ extrapolation & 1.3\% & 0.9\% & 0.7\% &
0.5\%\\
pert'n theory & 0.4 & 0.9 & 1.3 & 1.6\\
$\alpha_\msb$ uncertainty & 0.3 & 0.6 & 1.0 & 1.3\\
gluon condensate & 0.3 & 0.0 & 0.3 & 0.7\\
statistical errors & 0.0 & 0.0 & 0.0 & 0.0\\
relative scale errors & 0.4 & 0.4 & 0.3 & 0.3\\
overall scale errors & 0.6 & 0.6 & 0.7 & 0.7\\
sea quarks & 0.2 & 0.2 & 0.1 & 0.1\\
finite volume & 0.1 & 0.1 & 0.3 & 0.4\\
\hline
Total & 1.6\% & 1.6\% & 2.0\% & 2.4\%
\end{tabular}\end{ruledtabular}
\end{center}
\end{table}
One check on the reliability of our analysis comes from the $n=4$
moment, which is dimensionless. We compared our simulation result for
this moment with perturbation theory in order to extract a new value
for the QCD coupling. We obtained
\begin{equation}
\alpha_\msb^{(n_f=4)}(3\,\mathrm{GeV}) = 0.230\,(18)
\end{equation}
which is equivalent to $\alpha_\msb^{(n_f=5)}(M_Z) = 0.113(4)$. This
agrees well with the Particle Data Group's world average of
0.1176\,(20) for the result at $\mu=M_Z$~\cite{pdg} but is far from
being the most accurate determination from either the lattice or the
continuum. The bulk of the uncertainty in our new result for the
coupling constant comes from uncertainties in the $a^2$~extrapolation.
Our coupling is three to four times more sensitive to such effects than
are our mass determinations since the coupling is determined from
radiative corrections (which are suppressed by a power of~$\alpha_s$).
\section{Other Correlators}\label{sec:other}
The close agreement between different moments is important evidence
that we understand our systematic errors since these enter quite
differently in different moments. To further check this we repeated our
analysis for three different correlators, which we formed by replacing
the pseudoscalar operator $m_{0c}j_5$ with each of the following
$c$-quark currents on the lattice:
\begin{align}
j_{\mu}^{(1)} &\equiv \psib_c(x+a\hat\mu)\gamma_\mu\psi_c(x),\\
j_{\mu}^{(\mu)} &\equiv \psib_c(x)\gamma_\mu\psi_c(x), \\
j_{5\mu}^{(5\mu)} &\equiv \psib_c(x)\gamma_5\gamma_\mu\psi_c(x).
\end{align}
The first two currents are different lattice discretizations of the
vector current and were evaluated for space-like $\mu$s; and the first
of these was evaluated in Coulomb gauge. The third current is a lattice
discretization of the axial vector current and was evaluated for
time-like~$\mu$. The superscript on each $j$ labels the ``taste''
carried by that operator, using the notation presented in the
Appendices of~\cite{Follana:2006rc}. Taste is a spurious quantum
number, analogous to flavor, that is an artifact of staggered-quark
lattice discretizations like the HISQ formalism. Taste should not
affect physical results and therefore operators carrying different
taste here should give identical results in the $a\to0$~limit. By
studying these different currents, we not only test for conventional
systematic errors, but also verify that HISQ-specific taste effects are
negligible~\cite{footnote:6}.
\begin{table}
\caption{Simulation results for the reduced moments $R_n^{(j)}$,
extrapolated to~$a=0$,
from correlators of local axial-vector and vector lattice currents,
and a point-split lattice
vector current. Corresponding values for $m_c(\mu=3\,\mathrm{GeV})$
(in GeV) are also given.}
\label{tab:rnj}
\begin{ruledtabular}
\begin{tabular}{c|cc|cc|cc} \multicolumn{1}{c}{}
& \multicolumn{2}{c}{$j_{5\mu}^{(5\mu)}$} &
\multicolumn{2}{c}{$j_{\mu}^{(\mu)}$} &
\multicolumn{2}{c}{$j_{\mu}^{(1)}$} \\
\multicolumn{1}{c}{n}
& $R_n^{(j)}$ & \multicolumn{1}{c}{$m_c(\mu)$}
& $R_n^{(j)}$ & \multicolumn{1}{c}{$m_c(\mu)$}
& $R_n^{(j)}$ & $m_c(\mu)$ \\ \hline
6 & 1.243(13) & 0.95(3) & 1.256(15) & 0.98(3) & 1.268(15) & 0.96(3) \\
8 & 1.167(11) & 0.98(4) & 1.156(13) & 1.02(4) & 1.174(13) & 1.00(4) \\
10 & 1.131(10) & 0.98(4) & 1.127(12) & 1.02(5) & 1.137(12) & 1.01(5) \\
12 & 1.107(10) & 0.98(5) & 1.113(11) & 1.01(7) & 1.116(11) & 1.01(7) \\
14 & 1.086(10) & 0.99(7) & 1.096(11) & 1.01(8) & 1.097(11) & 1.01(8) \\
16 & 1.068(9) & 0.99(9) & 1.079(11) & 1.01(11) & 1.079(11) & 1.01(11)
\\
18 & & & 1.063(11) & 1.01(14) & 1.063(11) & 1.01(14)
\\
\end{tabular}
\end{ruledtabular}
\end{table}
A complication in our lattice analysis of these vector (or
axial-vector) correlators is that none of the currents is conserved (or
partially conserved) on the lattice. Consequently, each lattice current
is related to its corresponding continuum operator by a renormalization
constant:
\begin{align}
j_\mathrm{cont} &= Z^{(j)}\,j \,+\, \order(a^2) \\
Z^{(j)} &\equiv Z^{(j)}(\alpha_\msb(\pi/a),am_{0c}) \nonumber
\end{align}
where $j$ is one of the lattice currents $j_\mu^{(1)}$,
$j_\mu^{(\mu)}$, or $j_{5\mu}^{(5\mu)}$, and $j_\mathrm{cont}$ is the
continuum current $j_\mu=\psib\gamma_\mu\psi$ for the first two $j$s
and $j_{5\mu}=\psib\gamma_5\gamma_\mu\psi$ for the last. Consequently
moments of the correlators of these lattice currents have the form
\begin{equation} \label{Gnv-def}
G_n^{(j)} = \frac{1}{{Z^{(j)}}^2}\,
\frac{g_n^{(j_\mathrm{cont})}(\alpha_\msb(\mu),\mu/m_c)}{(am_c(\mu))^{n-2}}.
\end{equation}
where ${Z^{(j)}}^2 G_n^{(j)}$ is the continuum result for $n\ge4$. To
cancel the renormalization factor we redefine the reduced moments for
these correlators to be
\begin{align}
R_n^{(j)} &\equiv \frac{am^{(j)}}{2am_{0c}}\left(
\frac{G_n^{(j)}}{G_{n-2}^{(j)}} \,
\frac{G_{n-2}^{(j0)}}{G_n^{(j0)}}
\right)^{1/2} \\ \label{rnj-def}
&\equiv
\frac{r_n^{(j_\mathrm{cont})}(\alpha_\msb,\mu/m_c)}{2m_c(\mu)/m^{(j)}}
\end{align}
where $n\ge6$, and $m^{(j)}$ is the $\psi$~mass for the vector currents
(which couple to the $\psi$) and the $\eta_c$~mass for the axial-vector
current. Again we divide each moment $G_n^{(j)}$ by its value
$G_n^{(j0)}$ in leading-order lattice perturbation theory in order to
minimize finite-lattice-spacing errors. And again the perturbative
expansion for $r_n^{(j_\mathrm{cont})}$ can be obtained from continuum
perturbation theory expansions for the $g_n^{(j_\mathrm{cont})}$ (see
the Appendix).
Our simulation results for~$R_n^{(j)}$, extrapolated to lattice
spacing~$a=0$, are given in Table~\ref{tab:rnj} for different
moments~$n$ and each of the three currents. Perturbative coefficients
for the vector-current~$r_n^{(j_\mathrm{cont})}$s are discussed in the
Appendix; the coefficients for the temporal axial-vector current can be
derived from the pseudoscalar coefficients (also in the Appendix) using
Ward identities which imply~\cite{footnote:7}:
\begin{equation}
r_n^{(j_{5\mu})} = \left(r_{n+2}^{n-2}/r_{n}^{n-4}\right)^{1/2}.
\end{equation}
By combining perturbative with nonperturbative results, we obtain the
values for $m_c(\mu)$, with $\mu=3$\,GeV, that are listed in
Table~\ref{tab:rnj} and plotted in the top-right and bottom panels of
Figure~\ref{fig:mc}. Results from all moments agree with each other and
with the pseudoscalar result (the gray band in the plots), although
here the errors are about twice as large for the smaller moments.
\begin{figure}
\includegraphics[scale=1.0]{high-pp-mc}\\
\includegraphics[scale=1.0]{high-vv-mu-mc}
\caption{$m_c(\mu=3\,\mathrm{GeV})$ from large-$n$ moments of
the pseudoscalar and the (local) vector correlators. The gray band
is our final result for the mass, 0.984\,(16)\,GeV. The
perturbative part of the
analysis was evaluated at $\mu=m_c(\mu)$ using formulas
from~\cite{Maier:2007yn}, and
the results evolved to $\mu=3$\,GeV using fourth-order evolution.
Uncertainties due
to the gluon condensate are not included (see text).}
\label{fig:mc-largen}
\end{figure}
\section{Conclusions}\label{concl}
Our results are by far the most accurate determination of the $c$-quark
mass from lattice QCD~\cite{old-lattice-mc}. Such precision is possible
because the matching between lattice parameters and continuum
parameters here relies upon continuum perturbation theory, which is
much simpler than lattice QCD perturbation theory. Consequently
perturbation theory can be pushed to much higher orders.
The agreement between masses from different moments, and from different
correlators\,---\,27 determinations in all\,---\,is an important check
on systematic errors of all sorts since these enter in very different
ways in each calculation. Note that the different reduced moments in
our analysis vary in value by as much as 43\% (from 1.06 to 1.52), and
yet they all agree on the value of~$m_c$ to within a few percent once
we account for differences in the perturbative parts.
One surprising feature of our results is that even the higher moments
give correct values for the quark mass, albeit with larger errors.
Nonperturbative effects grow with~$n$ but our results show no
systematic deviation until very large~$n$, as is evident in
Fig.~\ref{fig:mc-largen} which shows results for~$20\le n\le62$. We
have not included potential errors due to the gluon condensate in this
figure. The error bars would have been much larger had we done so. For
example, they would have been about 5~times larger at $n=40$ in the
pseudoscalar plot (16\%~rather than~3\%). This might suggest that the
condensate is smaller than we allowed for\,---\,say $\langle \alpha_s
G^2/\pi\rangle\le0.003\,\mathrm{GeV}^4$\,---\,but we have not analyzed
this carefully enough to make a strong statement. The error bars shown
in the plots start to grow rapidly just where it becomes clear that
perturbation theory is failing (because of large coefficients).
Our lattice result for the mass,
$m_c(3\,\mathrm{GeV})=0.984\,(16)$\,GeV, agrees well with the continuum
determination, which uses $e^+e^-$~annihilation data and gives
0.986\,(13)\,GeV~\cite{Kuhn:2007vp}. This provides further strong
evidence that the different systematic errors in each calculation are
understood. The lattice analysis will be improved as data becomes
available for smaller lattice spacings. Also a very accurate $c$-quark
mass will allow us to make similarly accurate determinations of the
$s$-quark mass~\cite{c-s-ratio}. This is because the ratio $m_c/m_s$
can be determined very accurately in lattice simulations where the
$s$~and $c$~quarks are analyzed using the same formalism, as here.
Finally four-loop perturbation theory for further moments would improve
both lattice and continuum determinations.
\section*{Acknowledgements}
We thank Craig McNeile for useful discussions. We are grateful to the
MILC Collaboration for the use of their gluon configurations. The
lattice QCD computing was done on UKQCD's QCDOCX cluster and at the
Ohio Supercomputer Center. The work was supported by grants from: the
Deutsche Forschungsgemeinschaft, SFB-Transregio 9; the Department of
Energy (DE-FG02-91-ER40690, DE-AC02-98CH10886(BNL)); the Leverhulme
Trust; the Natural Sciences and Engineering Research Council; the
National Science Foundation; and the Science and Technology Facilities
Council.
\newpage
\begin{widetext}
\begin{center}
\begin{table}[ht]
\begin{center}
{\begin{tabular}{r|rrrrrrrrrr}
\hline
\hline
$k$ &
$\bar{C}_k^{(0)}$ & $\bar{C}_k^{(10)}$ & $\bar{C}_k^{(11)}$ &
$\bar{C}_k^{(20)}$ & $\bar{C}_k^{(21)}$ & $\bar{C}_k^{(22)}$ &
$\bar{C}_k^{(30)}$ & $\bar{C}_k^{(31)}$ & $\bar{C}_k^{(32)}$ &
$\bar{C}_k^{(33)}$
\\
\hline
1&$ 1.3333$&$ 3.1111$&$ 0.0000$&$ 0.1154$&$ -6.4815$&$
0.0000$&$ -1.2224$&$ 2.5008$&$ 13.5031$&$ 0.0000$\\
2&$ 0.5333$&$ 2.0642$&$ 1.0667$&$ 7.2362$&$ 1.5909$&$
-0.0444$&$ 7.0659$&$ -7.5852$&$ 0.5505$&$ 0.0321$\\
3&$ 0.3048$&$ 1.2117$&$ 1.2190$&$ 5.9992$&$ 4.3373$&$
1.1683$&$$---$$&$ 7.3626$&$ 4.2523$&$ -0.0649$\\
4&$ 0.2032$&$ 0.7128$&$ 1.2190$&$ 4.2670$&$ 4.8064$&$
2.3873$&$$---$$&$ 14.7645$&$ 11.0345$&$ 1.4589$\\
5&$ 0.1478$&$ 0.4013$&$ 1.1821$&$ 2.9149$&$ 4.3282$&$
3.4971$&$$---$$&$ 16.0798$&$ 16.6772$&$ 4.4685$\\
6&$ 0.1137$&$ 0.1944$&$ 1.1366$&$ 1.9656$&$ 3.4173$&$
4.4992$&$$---$$&$ 14.1098$&$ 19.9049$&$ 8.7485$\\
7&$ 0.0909$&$ 0.0500$&$ 1.0912$&$ 1.3353$&$ 2.2995$&$
5.4104$&$$---$$&$ 10.7755$&$ 20.3500$&$ 14.1272$\\
8&$ 0.0749$&$ -0.0545$&$ 1.0484$&$ 0.9453$&$ 1.0837$&$
6.2466$&$$---$$&$ 7.2863$&$ 17.9597$&$ 20.4750$\\
\hline
1&$ 1.0667$&$ 2.5547$&$ 2.1333$&$ 2.4967$&$ 3.3130$&$
-0.0889$&$ -5.6404$&$ 4.0669$&$ 0.9590$&$ 0.0642$\\
2&$ 0.4571$&$ 1.1096$&$ 1.8286$&$ 2.7770$&$ 5.1489$&$
1.7524$&$$---$$&$ 6.7216$&$ 6.4916$&$ -0.0974$\\
3&$ 0.2709$&$ 0.5194$&$ 1.6254$&$ 1.6388$&$ 4.7207$&$
3.1831$&$$---$$&$ 7.5736$&$ 13.1654$&$ 1.9452$\\
4&$ 0.1847$&$ 0.2031$&$ 1.4776$&$ 0.7956$&$ 3.6440$&$
4.3713$&$$---$$&$ 4.9487$&$ 17.4612$&$ 5.5856$\\
5&$ 0.1364$&$ 0.0106$&$ 1.3640$&$ 0.2781$&$ 2.3385$&$
5.3990$&$$---$$&$ 0.9026$&$ 18.7458$&$ 10.4981$\\
6&$ 0.1061$&$ -0.1158$&$ 1.2730$&$ 0.0070$&$ 0.9553$&$
6.3121$&$$---$$&$ -3.1990$&$ 16.9759$&$ 16.4817$\\
7&$ 0.0856$&$ -0.2033$&$ 1.1982$&$ -0.0860$&$ -0.4423$&$
7.1390$&$$---$$&$ -6.5399$&$ 12.2613$&$ 23.4000$\\
8&$ 0.0709$&$ -0.2660$&$ 1.1351$&$ -0.0496$&$ -1.8261$&$
7.8984$&$$---$$&$ -8.6310$&$ 4.7480$&$ 31.1546$\\
\hline
\hline
\end{tabular}
}
\end{center}
\caption{\label{tab:1}
Moments of the pseudoscalar
(upper 8 lines) and the vector (lower 8 lines) correlators, where
currently unknown coefficients are denoted by a dash and set to zero in
our analysis.
The numbers correspond to QCD with one massive $c$-quark and three
massless $(u,d,s)$ quarks.
}
\end{table}
\end{center}
\end{widetext}
\section*{Appendix: Continuum perturbation
theory\label{app:PT}}\label{cont-pert-th}
The correlators of two pseudoscalar ($\imath\,
\bar{\psi_c}\g_5\psi_c$),
vector $\bar{\psi_c}\g_{\mu}\psi_c$, and
axial-vector $\bar{\psi_c}\g_{\mu}\g_5\psi_c$ currents are defined by
\begin{equation}
\label{eq:pseudoscalar}
q^2\*\Pi^{p}(q^2)=i\int\!dx e^{iqx}\langle0|Tj_5(x)j_5(0)|0\rangle,
\end{equation}
\begin{eqnarray}
\label{eq:axialvector}
(q_{\mu}q_{\nu}-q^2g_{\mu\nu})\*\Pi^{\delta}(q^2)
+q_{\mu}q_{\nu}\Pi^{\delta}_{L}(q^2)=\Pi^{\delta}_{\mu\nu}(q)
\nonumber\\
=i\int\!dx\,
e^{iqx}\langle0|Tj_{\mu}^{\delta}(x)j_{\nu}^{\delta}(0)|0\rangle
{}.
\end{eqnarray}
where $\delta=v$ and~$a$ for the vector and axial-vector cases,
respectively.
The polarization functions can be expanded in the variable
$z=q^2/(2\*m_c(\mu))^2$:
\begin{equation}
\label{eq:expand}
\bar{\Pi}^{\delta}(q^2)=\frac{3}{16\pi^2}\sum_{k=-1}^{\infty}\bar{C}^{\delta}_k\*z^k
{}.
\end{equation}
The polarization function and the $c$-quark mass $m_c(\mu)$ are
renormalized
in the $\overline{\mbox{MS}}$-scheme.
The longitudinal part of the axial-vector current (which is of interest
in the present context) and the pseudoscalar correlator are related by
the axial Ward-identity\cite{Broadhurst:1981jk}:
\begin{equation}
\label{eq:AWI}
q^\mu\*q^\nu\*\Pi^a_{\mu\nu}(q)=(2\*m)^2\*q^2\*\Pi^p(q^2)+\mbox{contact
term.}
\end{equation}
Comparing different orders in $z$,
\begin{equation}
\label{eq:AWIExp}
\bar{C}^p_{k+1}=\bar{C}^{a}_{L,k},\quad\mbox{for}\quad k\ge-1,
\end{equation}
allows us to extract moments of the pseudoscalar correlator from those
of the longitudinal part of the axial-vector correlator, and \emph{vice
versa}.
The contact term contributes only to $k=-2$.
The coefficients of the perturbative expansion depend logarithmically
on
$m_c$ and can be written in the form
\begin{eqnarray}
\bar{C}_k &=&
\bar{C}_k^{(0)}
+ \frac{\alpha_s(\mu)}{\pi}
\left( \bar{C}_k^{(10)} + \bar{C}_k^{(11)}\lmc \right) \nonumber \\
&+& \left(\frac{\alpha_s(\mu)}{\pi}\right)^2
\left( \bar{C}_k^{(20)} + \bar{C}_k^{(21)}\lmc
+ \bar{C}_k^{(22)}\lmc^2 \right)
\nonumber\\
&+& \left(\frac{\alpha_s(\mu)}{\pi}\right)^3
\left( \bar{C}_k^{(30)} + \bar{C}_k^{(31)}\lmc \right. \nonumber \\
&&\quad\quad\quad\quad\quad\left. + \bar{C}_k^{(32)}\lmc^2 +
\bar{C}_k^{(33)}\lmc^3 \right)
+ \ldots
\,.
\label{eq::cn}
\end{eqnarray}
where $l_{m_c}\equiv \log(m_c^2(\mu)/ \mu^2)$. The coefficients
$\bar C_k^{(ij)}$ up through $k=8$ are listed in Table~\ref{tab:1} for
both pseudoscalar
and vector correlators.
The four-loop coefficients $\bar C_1^{(30)}$ and $\bar C_2^{(30)}$
for the pseudoscalar correlator are new~\cite{Sturm:2008aa}; and
$\bar C_1^{(30)}$ for the vector correlator comes
from~\cite{Chetyrkin:2006xg,Boughezal:2006px}.
The order $\alpha_s^2$ terms up through $k=8$ are given in
\cite{Chetyrkin:1995ii,Chetyrkin:1996cf,Chetyrkin:1997mb},
while results for higher $k$-values are given in
\cite{Boughezal:2006uu,Maier:2007yn}. See also~\cite{Czakon:2007qi} for
$n_f$-dependent four-loop terms, and~\cite{Maier:2007yn} for the
pseudoscalar case. Throughout this paper we take the number of light
(massless) active quark flavors to be $n_l=3$, and
the number of heavy (massive) quarks is set to $n_h=1$. For numerical
work, we use~$\mu=3\,$GeV.
The expansion coefficients $\bar{C}_k$ are related to the
coefficients $g_n$ of \eq{gn-def} by:
\begin{equation}
\label{eq:CngnRelation}
\frac{g_{2k+2}}{g_{2k+2}^{(0)}} =
\frac{\bar{C}_{k}}{\bar{C}_{k}^{(0)}}
{},
\end{equation}
and the coefficients~$r_{n,i}$ in \eq{rn-exp}
are obtained through the series expansion of \eq{rn-def}.
For the vector correlator the coefficients $r^{(j_{\mu})}_{n,i}$
are defined through the series expansions of
\begin{equation}
r^{(j_{\mu})}_{2k+2} =
\left(
\frac{\bar{C}_k^{v} }{\bar{C}_k^{v,(0)}}
\,
\frac{\bar{C}_{k-1}^{v,(0)}}{\bar{C}_{k-1}^{v} }
\right)^{1/2}
{}.
\end{equation}
\bibliographystyle{plain}
\begin{thebibliography}{16}
\bibitem{Kuhn:2001dm}
J.~H.~Kuhn and M.~Steinhauser,
Nucl.\ Phys.\ ¬†B {\bf 619}, 588 (2001)
[Erratum-ibid.\ ¬†B {\bf 640}, 415 (2002)]
[arXiv:hep-ph/0109084].
\bibitem{Kuhn:2007vp}
J.~H.~Kuhn, M.~Steinhauser and C.~Sturm,
Nucl.\ Phys.\ B {\bf 778}, 192 (2007)
[arXiv:hep-ph/0702103].
\bibitem{Chetyrkin:1995ii}
K.~G.~Chetyrkin, J.~H.~Kuhn and M.~Steinhauser,
Phys.\ Lett.\ ¬†B {\bf 371}, 93 (1996)
[arXiv:hep-ph/9511430].
\bibitem{Chetyrkin:1996cf}
K.~G.~Chetyrkin, J.~H.~Kuhn and M.~Steinhauser,
Nucl.\ Phys.\ ¬†B {\bf 482}, 213 (1996)
[arXiv:hep-ph/9606230].
\bibitem{Chetyrkin:1997mb}
K.~G.~Chetyrkin, J.~H.~Kuhn and M.~Steinhauser,
Nucl.\ Phys.\ ¬†B {\bf 505}, 40 (1997)
[arXiv:hep-ph/9705254].
\bibitem{Boughezal:2006uu}
R.~Boughezal, M.~Czakon and T.~Schutzmeier,
Nucl.\ Phys.\ Proc.\ Suppl.\ ¬†{\bf 160}, 160 (2006)
[arXiv:hep-ph/0607141].
\bibitem{Maier:2007yn}
A.~Maier, P.~Maierhofer and P.~Marquard,
Nucl.\ Phys.\ ¬†B {\bf 797}, 218 (2008)
[arXiv:0711.2636 [hep-ph]].
\bibitem{Chetyrkin:2006xg}
K.~G.~Chetyrkin, J.~H.~Kuhn and C.~Sturm,
Eur.\ Phys.\ J.\ C {\bf 48}, 107 (2006)
[arXiv:hep-ph/0604234].
\bibitem{Boughezal:2006px}
R.~Boughezal, M.~Czakon and T.~Schutzmeier,
Phys.\ Rev.\ D {\bf 74}, 074006 (2006)
[arXiv:hep-ph/0605023].
\bibitem{Bochkarev:1995ai}
For an earlier exploratory analysis, using less advanced lattice and
continuum QCD results, see
A.~Bochkarev and P.~de Forcrand,
Nucl.\ Phys.\ B {\bf 477}, 489 (1996)
[arXiv:hep-lat/9505025].
\bibitem{Davies:2003ik}
C.~T.~H.~Davies {\it et al.} [HPQCD Collaboration],
Phys.\ Rev.\ Lett.\ {\bf 92}, 022001 (2004)
[arXiv:hep-lat/0304004].
\bibitem{Mason:2005zx}
Q.~Mason {\it et al.} [HPQCD Collaboration],
Phys.\ Rev.\ Lett.\ {\bf 95}, 052002 (2005)
[arXiv:hep-lat/0503005].
\bibitem{Follana:2006rc}
E.~Follana {\it et al.} [HPQCD Collaboration],
Phys.\ Rev.\ D {\bf 75}, 054502 (2007)
[arXiv:hep-lat/0610092].
\bibitem{Follana:2007uv}
E.~Follana, C.~T.~H.~Davies, G.~P.~Lepage and J.~Shigemitsu [HPQCD
Collaboration],
Phys.\ Rev.\ Lett.\ {\bf 100}, 062002 (2008)
[arXiv:0706.1726 [hep-lat]].
\bibitem{MILC} C. Aubin et al, Phys. Rev. D70:094505, 2004.
\bibitem{footnote:0} The value for $t/a$ at the middle of the lattice
was
set to~0 in order to preserve the symmetry between~$t$ and~$-t$. This
choice
causes odd moments to vanish, for example. The value chosen for this
point doesn't
matter in our analysis because contributions to our moments from this
point are
negligibly small (because of the exponential falloff in the correlators
at large~$t$).
\bibitem{footnote:1} The suppression is by a factor of $\alpha_s(1/a)$,
which explains the factor of three.
\bibitem{footnote:2} The ratio $m_c/m_{\eta_c}$ is independent of $m_c$
up to binding
corrections of order $(v_c/c)^2\approx 1/3$, which is why the
suppression
is again by roughly a factor of three. See C.~T.~H.~Davies {\it et
al.},
Phys.\ Rev.\ Lett.\ {\bf 73}, 2654 (1994)
[arXiv:hep-lat/9404012];
A.~Gray, I.~Allison, C.~T.~H.~Davies, E.~Dalgic, G.~P.~Lepage,
J.~Shigemitsu and M.~Wingate,
Phys.\ Rev.\ D {\bf 72}, 094507 (2005)
[arXiv:hep-lat/0507013].
\bibitem{footnote:2.1} We did additional simulations to test sea-quark
mass
and spatial volume dependence using $a^{-1}=1.6$\,GeV lattices with:
1) $am_{0u/d}=0.007$, $L/a=24$; and 2) $am_{0u/d}=0.028$, $L/a=20$.
\bibitem{footnote:2.2} We computed the contribution from $c$-quark
loops by comparing results from
\cite{Maier:2007yn} for $n_h=1$ (with $c$ quarks) with those for
$n_h=0$ (without $c$ quarks).
The reduced moments~$R_n$ are corrected, to include $c$~quarks, by
multiplying them
by $r_n(n_h=1)/r_n(n_h=0)$.
\bibitem{footnote:3} More specifically these errors come from a
0.5\%~uncertainty in our
ability to measure ratios of lattice spacings, and an overall,
correlated 1.5\%~uncertainty
in the measured values of all lattice spacings
(see~\cite{Follana:2007uv}). These errors are all reduced by a
factor~$1/3$
because of our use of reduced moments, as discussed in the text.
\bibitem{footnote:4} We estimate this by inserting a factor of
$1/(1+(ap)^2)$ into the Feynman integral for lowest order and
calculating the change, $\Delta G_n$, in the moment. The average
value of $(ap)^2$ is then approximated by $|\Delta G_n/G_n|$.
\bibitem{Lepage:2001ym}
G.~P.~Lepage, B.~Clark, C.~T.~H.~Davies, K.~Hornbostel,
P.~B.~Mackenzie, C.~Morningstar and H.~Trottier,
Nucl.\ Phys.\ Proc.\ Suppl.\ {\bf 106}, 12 (2002)
[arXiv:hep-lat/0110175].
\bibitem{footnote:4.5} The obvious choice, $\mu=m_c(m_c)$, gives
smaller perturbative
coefficients for the moments of interest, but $\alpha_\msb(\mu)$ is
substantially larger.
Consequently $\mu=3\,$GeV gives smaller errors than $\mu=m_c(m_c)$
despite having somewhat larger
expansion coefficients. Differences in the final results are small,
however.
\bibitem{footnote:5} This gives a perturbative error that is
approximately the same size as the error estimates obtained
in~\cite{Kuhn:2007vp} by varying the renormalization scale
between~$\mu=2$~and~4\,GeV.
The value for $R_6$, for example, varies from $+0.5$\% to $-0.4$\% of
the central
value as $\mu$ is varied over this range; the error estimate
from~\eq{rn-err} is~$\pm 0.4$\%.
\bibitem{pdg} W.-M. Yao {\it et al.} (Particle Data Group),
J.\ Phys.\ {\bf G 33}, 1 (2006). We use
the Particle Data Group's value of $m_b(m_b) = 4.25$\,GeV for the
$b$-quark mass when reducing from
$n_f=5$ to $n_f=4$.
\bibitem{van Ritbergen:1997va}
T.~van Ritbergen, J.~A.~M.~Vermaseren and S.~A.~Larin,
Phys.\ Lett.\ B {\bf 400}, 379 (1997)
[arXiv:hep-ph/9701390].
\bibitem{Czakon:2004bu}
M.~Czakon,
Nucl.\ Phys.\ B {\bf 710}, 485 (2005)
[arXiv:hep-ph/0411261].
\bibitem{Chetyrkin:1997dh}
K.~G.~Chetyrkin,
Phys.\ Lett.\ B {\bf 404}, 161 (1997)
[arXiv:hep-ph/9703278].
\bibitem{Vermaseren:1997fq}
J.~A.~M.~Vermaseren, S.~A.~Larin and T.~van Ritbergen,
Phys.\ Lett.\ B {\bf 405}, 327 (1997)
[arXiv:hep-ph/9703284].
\bibitem{mc-footnote} Note that $m_c(m_c)=1.266\,(14)$\,GeV implies
that
$m_c(1.266\,\mathrm{GeV})=1.266\,(20)$\,GeV. The two have the same
central value
but different uncertainties. This is because the scale in $m_c(m_c)$ is
not fixed but rather varies as $m_c(m_c)$ itself is varied across the
interval
$1.266\pm0.014$\,GeV.
\bibitem{Novikov:1977dq}
V.~A.~Novikov, L.~B.~Okun, M.~A.~Shifman, A.~I.~Vainshtein,
M.~B.~Voloshin and V.~I.~Zakharov,
Phys.\ Rept.\ ¬†{\bf 41}, 1 (1978).
\bibitem{Broadhurst:1994qj}
D.~J.~Broadhurst, P.~A.~Baikov, V.~A.~Ilyin, J.~Fleischer,
O.~V.~Tarasov and V.~A.~Smirnov,
Phys.\ Lett.\ B {\bf 329}, 103 (1994)
[arXiv:hep-ph/9403274].
\bibitem{cond-footnote}
We include the gluon condensate by multiplying $r_n$ by a factor of
the form $(1+d_n\langle \alpha_s G^2/\pi\rangle/(2m_c)^4))$
where, here, $m_c=m_c(m_c)$ and $d_n$ is computed through leading
order in $\alpha_\msb(m_c)$. We set
$\langle \alpha_s G^2/\pi\rangle = 0\pm0.012\,\mathrm{GeV}^4$,
which covers the range of most current values. For a recent analysis of
condensate values see B.~L.~Ioffe,
Prog.\ Part.\ Nucl.\ Phys.\ {\bf 56}, 232 (2006)
[arXiv:hep-ph/0502148].
\bibitem{footnote:6} A related problem is that correlators built from
these currents
are contaminated by contributions from correlators built from
opposite-parity
operators (see Appendix~G of~\cite{Follana:2006rc}). These lattice
artifacts
oscillate in sign, like $(-1)^{t/a}$, and so are
suppressed by $(am_c/\pi)^{n-2}$ in the $n^\mathrm{th}$~moment since
the
sum over~$t$ is dominated by~$E\approx0$. This makes them
negligible except possibly for $n=4$; the slightly lower masses we
obtain
from ratios of $n=6$~and 4~moments may be due to this effect. Such
contamination is absent in the case of the pseudoscalar correlators.
\bibitem{footnote:7} The continuum Ward identity relates the axial
vector current
to the pseudoscalar density~\cite{Broadhurst:1981jk}, implying that
$g_n^{(j_{5\mu})} \propto g_{n+2}$
where the $g_n$s are defined in equations~(\ref{Gnv-def})
and~(\ref{gn-def}).
\bibitem{old-lattice-mc} For older analyses that include sea quarks
see: A.~Dougall, C.~M.~Maynard and C.~McNeile,
JHEP {\bf 0601}, 171 (2006)
[arXiv:hep-lat/0508033];
and M.~Nobes and H.~Trottier,
PoS {\bf LAT2005}, 209 (2006)
[arXiv:hep-lat/0509128].
The precision of analyses that do not include sea quarks is hard to
quantify because the value of the $c$-mass depends upon how the lattice
spacing and the bare $c$-mass are tuned, and it is not possible to
obtain consistency from different choices. For example, tuning the bare
mass to give the correct $D_s$~mass gives very different results from
tuning to get the correct $D_s^*-D_s$ mass difference, or the correct
$\eta_c$ or $\psi$~mass. We avoid these problems here by including the
full effect of sea quarks. For determinations of the $c$-mass without
sea quarks with further discussion of these issues see:
J.~Rolf and S.~Sint [ALPHA Collaboration],
JHEP {\bf 0212}, 007 (2002) [arXiv:hep-ph/0209255]; and
G.~M.~de Divitiis, M.~Guagnelli, R.~Petronzio, N.~Tantalo and
F.~Palombi,
Nucl.\ Phys.\ B {\bf 675}, 309 (2003)
[arXiv:hep-lat/0305018].
\bibitem{c-s-ratio} In preparation.
\bibitem{Broadhurst:1981jk}
D.~J.~Broadhurst,
Phys.\ Lett.\ B {\bf 101}, 423 (1981).
\bibitem{Sturm:2008aa}
Details of the calculation, as well as analytical results, will
be presented in C.~Sturm, in preparation.
\bibitem{Czakon:2007qi}
M.~Czakon and T.~Schutzmeier,
arXiv:0712.2762 [hep-ph].
\end{thebibliography}
\end{document}