%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: <EM> Phys.Rev.  </EM> <B> D78:   </B> (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}
