Light quark masses from unquenched lattice QCD
UTCCS-P-31, UTHEP-542, KEK-CP-193, HUPD-0702, RBRC-666
Light quark masses from unquenched lattice QCD
T. Ishikawa,1, 2 S. Aoki,3, 2 M. Fukugita,4 S. Hashimoto,5, 6 K-I. Ishikawa,7
N. Ishizuka,1, 3 Y. Iwasaki,3 K. Kanaya,3 T. Kaneko,5, 6 Y. Kuramashi,1, 3 M. Okawa,7
Y. Taniguchi,1, 3 N. Tsutsui,5 A. Ukawa,1, 3 N. Yamada,5, 6 and T. Yoshié1, 3
(CP-PACS and JLQCD Collaborations)
1 Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan
2 RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973, USA
3 Graduate School of Pure and Applied Sciences, University of Tsukuba, Tsukuba 305-8571, Japan
4 Institute for Cosmic Ray Research, University of Tokyo, Kashiwa 277-8582, Japan
5 High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan
6 School of High Energy Accelerator Science, The Graduate
University for Advanced Studies (Sokendai), Tsukuba 305-0801, Japan
7 Department of Physics, Hiroshima University, Higashi-Hiroshima 739-8526, Japan
(Dated: October 31, 2018)
We calculate the light meson spectrum and the light quark masses by lattice QCD simulation,
treating all light quarks dynamically and employing the Iwasaki gluon action and the nonperturba-
tively O(a)-improved Wilson quark action. The calculations are made at the squared lattice spacings
at an equal distance a2 ≃ 0.005, 0.01 and 0.015 fm2, and the continuum limit is taken assuming an
O(a2) discretization error. The light meson spectrum is consistent with experiment. The up, down
and strange quark masses in the MS scheme at 2 GeV are m = (mu +md)/2 = 3.55
+0.65
−0.28 MeV and
ms = 90.1
+17.2
−6.1 MeV where the error includes statistical and all systematic errors added in quadra-
ture. These values contain the previous estimates obtained with the dynamical u and d quarks
within the error.
The masses of light quarks are fundamental parame-
ters of QCD. They cannot be measured experimentally
since quarks are confined in hadrons. Lattice QCD en-
ables calculations of hadron masses as functions of quark
masses, and hence allows a determination of the quark
masses from the experimental hadron masses. This ap-
proach has been successfully applied, first in quenched
QCD [1] and then in Nf = 2 QCD where degenerate
up (u) and down (d) quarks are treated dynamically [2].
These studies have revealed that the light quark mass val-
ues are significantly reduced by dynamical u and d quark
effects. In this article, we present our attempt to deter-
mine the quark masses in Nf = 2 + 1 QCD where the
heavier strange (s) quark is also treated dynamically. We
wish to examine to what extent the dynamical s quark
affects the light quark masses. We determine the quark
masses in the continuum limit and estimate all possible
systematic errors. We also calculate the prerequisite light
meson spectrum. A similar attempt has been made by
the MILC Collaboration [3].
We adopt the Iwasaki RG gauge action [4] and the
clover quark action with the improvement coefficient cSW
determined nonperturbatively for the RG action [5]. The
choice of the gauge action is made to avoid a first-order
phase transition (lattice artifact) observed for the plaque-
tte gauge action [6]. We employed the Wilson quark for-
malism because we prefer an unambiguous quark-flavor
interpretation over the computational ease of the stag-
gered formalism adopted by the MILC collaboration [7].
Configurations are generated at three values of the cou-
pling β ≡ 6/g2 = 2.05, 1.90 and 1.83 corresponding to the
squared lattice spacing a2 ≃ 0.005, 0.01 and 0.015 fm2,
with the physical volume fixed to about (2.0fm)3. At
each β, we perform simulations for 10 quark mass com-
binations using a combined algorithm [8] of the Hybrid
Monte Carlo (HMC) for the degenerate u and d quarks
and the polynomial Hybrid Monte Calro (PHMC) for the
s quark. Table I summarizes the simulation parameters.
The meson and quark masses at the simulation
points are determined from single exponential correlated
χ2 fits to the correlators 〈P (t)P (0)〉, 〈V (t)V (0)〉 and
〈A4(t)P (0)〉, where P , V and Aµ denote pseudoscalar,
vector and nonperturbatively O(a)-improved [9] axial-
vector current operators, respectively. We use an ex-
ponentially smeared source and a point sink, and mea-
surements are made at every 10 HMC trajectories in the
Coulomb gauge. For the pseudoscalar sector, 〈P (t)P (0)〉
and 〈A4(t)P (0)〉 are fitted simultaneously ignoring cor-
relations among them. Errors are estimated by the jack-
knife method with a bin size of 100 HMC trajectories;
errors do not increase for larger bin sizes.
Chiral fits are made to the light-light (LL), light-
strange (LS) and strange-strange (SS) meson masses
simultaneously ignoring their correlations, using a
quadratic polynomial function of the sea quark masses
(mu,md,ms) and valence quark masses (mval1,mval2) in
mesons;
f(Ms,Mv) (1)
= A+BStrMs +BV trMv +DSV trMstrMv
+CS1trM
s + CS2(trMs)
2 + CV 1trM
v + CV 2(trMv)
where f = m2PS or vector meson mass mV , MS =
http://arxiv.org/abs/0704.1937v2
TABLE I: Simulation parameters; L3×T is the lattice size, (κud, κs) is the hopping parameter combination, 1/δτ is the number
of molecular dynamics steps in one trajectory, Npoly is the PHMC polynomial order, and traj. is analyzed trajectory length.
Pseudoscalar vector mass ratios
are also listed for light-light (LL) and strange-strange (SS) mass combinations.
β = 1.83, L3 × T = 163 × 32, cSW = 1.761
κud κs δτ Npoly traj.
(SS) κud κs δτ Npoly traj.
0.13655 0.13710 1/80 80 7000 0.7772(13) 0.7522(15) 0.13655 0.13760 1/90 110 7000 0.7769(14) 0.7235(19)
0.13710 1/85 80 7000 0.7524(21) 0.7524(21) 0.13710 1/100 110 8600 0.7448(14) 0.7128(16)
0.13760 1/100 100 7000 0.7076(18) 0.7414(17) 0.13760 1/110 120 8000 0.7033(18) 0.7033(18)
0.13800 1/120 110 8000 0.6629(22) 0.7365(16) 0.13800 1/120 130 8100 0.6525(23) 0.6941(20)
0.13825 1/140 120 8000 0.6213(24) 0.7343(15) 0.13825 1/150 150 8100 0.6083(32) 0.6884(21)
β = 1.90, L3 × T = 203 × 40, cSW = 1.715
κud κs δτ Npoly traj.
(SS) κud κs δτ Npoly traj.
0.13580 0.13580 1/125 110 5000 0.7673(15) 0.7673(15) 0.13580 0.13640 1/125 140 5200 0.7667(16) 0.7211(21)
0.13610 1/125 110 6000 0.7435(18) 0.7647(17) 0.13610 1/125 140 8000 0.7444(15) 0.7182(17)
0.13640 1/140 110 7600 0.7204(19) 0.7687(15) 0.13640 1/140 140 9000 0.7145(16) 0.7145(16)
0.13680 1/160 110 8000 0.6701(27) 0.7673(17) 0.13680 1/160 140 9200 0.6630(21) 0.7127(17)
0.13700 1/180 110 7900 0.6390(22) 0.7691(15) 0.13700 1/180 140 7900 0.6243(28) 0.7102(20)
β = 2.05, L3 × T = 283 × 56, cSW = 1.628
κud κs δτ Npoly traj.
(SS) κud κs δτ Npoly traj.
0.13470 0.13510 1/175 200 6000 0.7757(26) 0.7273(29) 0.13470 0.13540 1/175 250 6000 0.7790(23) 0.6821(32)
0.13510 1/195 200 6000 0.7316(24) 0.7316(24) 0.13510 1/195 250 6000 0.7341(29) 0.6820(39)
0.13540 1/225 200 6000 0.6874(30) 0.7395(23) 0.13540 1/225 250 6000 0.6899(34) 0.6899(34)
0.13550 1/235 200 6500 0.6611(34) 0.7361(25) 0.13550 1/235 250 6500 0.6679(45) 0.6899(43)
0.13560 1/250 200 6500 0.6337(38) 0.7377(28) 0.13560 1/250 250 6500 0.6361(47) 0.6852(46)
0 0.02 0.04 0.06 0.08
Ks=0.13580
Ks=0.13640
β=1.90, L
xT=20
x40, polynomial
0 0.02 0.04 0.06 0.08
Ks=0.13580
Ks=0.13640
β=1.90, L
xT=20
x40, polynomial
FIG. 1: Chiral fits of meson masses with mAWIq at β = 1.90.
diag(mu,md,ms), MV = diag(mval1,mval2), and “tr”
means the trace of matrices. In the fits, we use
the axial-vector Ward identity quark mass mq =
limt→∞〈∂4A4(t)P (0)〉/(2〈P (t)P (0)〉) and set A = BS =
CS1 = CS2 = 0 for m
PS . These fits reproduce mea-
sured data well, as illustrated in Fig.1, with reasonable
χ2/d.o.f. of at most 1.36.
The physical quark mass point and the lattice spac-
ing are determined from the experimental values of π0,
ρ0 and K (K-input) or π0, ρ0 and φ (φ-input) meson
masses. Taking the ρ0 mass as input may cause a large
systematic error, because the ρ → ππ decay mode is
not open for our mass range (the lightest pion mass in
this simulation ∼ 620 MeV) and hence chiral extrapo-
TABLE II: Lattice spacings in fm units.
β K-input φ-input [π,K, φ]-input
1.83 0.1174(23) 0.1184(26) 0.1095(25)
1.90 0.0970(26) 0.0971(25) 0.0936(33)
2.05 0.0701(29) 0.0702(28) 0.0684(41)
lation of mV for lighter quarks may be quite different
from our fits. In order to estimate this uncertainty, we
also check another combination [π0,K, φ]. We assume
the ideal mixing for the vector isosinglets. Since our sim-
ulation is made with degenerate u and d quarks, we con-
sider the isospin averages m
= {(m2
)/2}1/2
and m
= (mK∗± + mK∗0)/2 and predict the aver-
age light quark mass m = (mu + md)/2. The electro-
magnetic (EM) effects, not included in our simulations,
are removed from the mK± above using Dashen’s the-
orem [10] (m2
− m2K0)EM = (m
− m2π0)EXP. The
isospin breaking effects and the EM effects for other
mesons we consider are expected to be small and thus are
not considered. The experimental values we use are taken
from the PDG booklet [11]; mπ0 = 0.1350GeV, mπ± =
0.1396GeV, mK0 = 0.4976GeV, mK± = 0.4937GeV,
mρ0 = 0.7755GeV, mK∗0 = 0.8960GeV, mK∗± =
0.8917GeV and mφ = 1.0195GeV. Lattice spacings (Ta-
ble II) for theK- and φ- inputs are consistent, while those
for the [π,K, φ]-input are slightly smaller by at most 7%.
An agreement of the meson spectrum with experiment
is a necessary condition for a reliable estimate of the
quark masses. To confirm this, we extrapolate the me-
TABLE III: Meson masses in the continuum limit (in MeV
units), compared to experiment. The EM effect is subtracted
using Dashen’s theorem.
K-input φ-input [π,K, φ]-input EXP.
K̂ - 491(19) - 495.0
ρ0 - - 761(32) 775.5
K̂∗ 900.5(9.9) 898.0(1.4) 891(16) 893.9
φ 1025(19) - - 1019.5
1.000
1.050
0.860
0.880
0.900
0.880
0.890
0.900
0 0.01 0.02 0.03 0.04 0.05
0.450
0.500
0.550
Nf=0     (IW+   PT clover)
Nf=2     (IW+   PT clover)
Nf=2+1 (IW+NPT clover)
experiment
φ (K-input)
K* (K-input)
K* (φ-input)
K (φ-input)
FIG. 2: Continuum extrapolation of meson masses for Nf =
2 + 1 QCD (circles), compared to experiment (stars) and re-
sults in Nf = 2 (squares) and Nf = 0 (triangles) QCD [2].
son masses linearly in a2, because our action is O(a)
improved and data are well fitted, as shown in Fig.2,
with small χ2/d.o.f ≤ 1.4. The masses in the continuum
limit, summarized in Table III, are consistent with exper-
iment with at most 2.9σ deviation. The K̂∗ mass turns
out to be slightly heavier than experiment, though the
supplemental [π,K, φ]-input gives consistent results with
experiment with large statistical error. Possible origin
of the deviation is due to uncertainty of chiral fits. In
fact, an alternative fit based on chiral perturbation the-
ory (χPT) we discuss later yields m
= 894(12) MeV
(K-input). In Fig.2 we overlay the previous results of
meson masses [2] in the Nf = 2 and quenched (Nf = 0)
QCD with tadpole improved one-loop cSW . The dynam-
ical u and d quarks significantly reduce the O(10%) de-
viation of the quenched spectrum from experiment. We
find no further dynamical s quark effect beyond statisti-
cal errors.
The quark masses are evaluated for the MS scheme
at the scale µ = 2GeV using the tadpole improved one-
loop matching [12] at µ = a−1 with an improved cou-
pling determined from plaquette and rectangular loop
and four-loop renormalization group equation. In the
continuum extrapolation of the quark masses, we assume
0 0.01 0.02 0.03 0.04 0.05
Nf=0     (K, IW+   PT clover)
Nf=2     (K, IW+   PT clover)
Nf=2+1 (K, IW+NPT clover)
0 0.01 0.02 0.03 0.04 0.05
Nf=0    (K, IW+  PT clover)
Nf=0    (φ, IW+   PT clover)
Nf=2    (K, IW+  PT clover)
Nf=2    (φ, IW+   PT clover)
Nf=2+1(K, IW+NPT clover)
Nf=2+1(φ, IW+NPT clover)
FIG. 3: Continuum extrapolations of the up, down and
strange quark masses. For comparison, results for Nf = 0
and Nf = 2 QCD [2] are overlaid.
the O(g4ma) contributions are small and neglect it. As
Fig. 3 shows, the quark masses are well described by
a linear function in a2, and the values determined for
either the K- or the φ-inputs, while different at finite
lattice spacings, extrapolate to a common value in the
continuum limit. Therefore the continuum limit is esti-
mated from a combined linear fit with the K- and the
φ-inputs. We obtain mMS(µ = 2GeV) = 3.55(19) MeV
and mMSs (µ = 2GeV) = 90.1(4.3) MeV with a suffi-
ciently small χ2/d.o.f. < 0.42. Note that the supplemen-
tal [π,K, φ]-input gives larger statistical error and hence
is not used to estimate central values.
We now turn to estimates of possible systematic errors.
Finite size effect (FSE) — The meson masses at
the infinite volume are estimated at β = 1.90 using
data on a V ∼ (2.0fm)3 lattice and those from our
exploratory study on a V ∼ (1.6fm)3 lattice [13], and
assuming a strong volume dependence of (mhad,V −
mhad,V=∞)/mhad,V=∞ ∝ 1/V [14]. The chiral fits to the
infinite volume values lead to less than a 4% change for
the meson masses at the physical point. For the quark
masses, however, we find a larger shift of 12.2% from a
V ∼ (2.0fm)3 lattice to V = ∞ for m with φ-input and
8.1% for ms with K-input (differences are smaller for the
other cases). Assuming that FSE is independent of lat-
tice spacing, we take the differences as estimates of FSE
for the quark masses in the continuum limit.
Chiral extrapolation — In addition to the polynomial
chiral fits, we fit the meson masses using χPT formu-
laes modified for the Wilson quark action (WχPT) [15].
Namely, we fit mπ, mK̂ , mρ and mK̂∗ using the NLO
Nf = 2+1 QCD WχPT formulae for the O(a) improved
theory [16]. Since the formula in Ref. [16] is not appli-
cable for the φ meson, we estimate the effect only for
K-input. In the fits we obtain m to be 3.1% smaller and
ms to be 1.2% larger than those of the polynomial fit. We
note that our WχPT fits to data do not exhibit a clear
chiral logarithm, probably because u and d quark masses
in our simulation are not sufficiently small. Further pos-
sible systematic error from a long chiral extrapolation for
ρ0, mentioned above, is estimated by the supplemental
[π,K, φ]-input, which gives 3.0% larger for m and 3.4%
larger value for ms than the central one. For an estimate
of systematic errors from chiral fits, differences of the two
alternative fits from the central value are added linearly.
Renormalization factor — Uncertainty of the one-
loop calculation of the renormalization factor is esti-
mated by shifting the matching scale from µ = 1/a to
µ = π/a and also using an alternative tadpole improved
coupling [2].
Continuum extrapolation — Possible O(a3) effects
are investigated by performing the continuum extrapola-
tion adding an O(a3) term to the fit function.
Electromagnetic (EM) effects — Systematic error
due to uncertainty of the EM effects is estimated fol-
lowing extensive arguments [7, 17, 18] to Dashen’s
theorem [10]. Namely, we estimate the effects by a
further mass shift of our input m
using a relation
)EM = (1 + ∆E)(m
)EXP assum-
ing the EM effects for other mesons are negligeble. We
vary the ∆E in range [−1,+1] as our estimate of the EM
effects, and we find a quite small change in ms and no
change in m.
Isospin breaking effects — Isospin breaking effects
are estimated by chiral fits with Eq. (1) for mu 6= md
and taking mπ0 , mρ0 , mK± and mK0 as inputs. We
find that mu/md = 0.577(25), and that m and ms
have no change from the K̂ input result. We note that
mu/md strongly depends on an estimate of the EM ef-
fects; mu/md =0.663–0.498 for ∆E = [−1,+1], though
m and ms almost do not.
Finally we obtain
mMS(µ = 2GeV)
= 3.55(19)(+43
−0 )(
−11)(
−17)(
−0 )(
−0), (2)
mMSs (µ = 2GeV)
= 90.1(4.3)(+7.3
−0 )(
−0 )(
−4.3)(
+12.8
−0 )(
−0.2)(
−0), (3)
in MeV units, where the errors are statistical, systematic
due to FSE, chiral extrapolation, renormalization factor,
continuum extrapolation, EM effect and isospin break-
ing effect, respectively. Adding the errors in quadrature
yields the values quoted in the abstract. These values
agree well with the latest report from the MILC Collab-
oration [3] m = 3.3 ± 0.3 MeV and ms = 90 ± 6 MeV
where we added the quoted errors in quadrature. They
also include the Nf = 2 values [2] within the error.
Scaling violation in the quark masses is unexpectedly
large, while that for the meson masses are reasonably
bounded at a percent level at a ≈ 0.1 fm. To gain a
better control over systematic uncertainties, a significant
reduction in the simulated light quark masses on a
correspondingly larger lattice is needed. An attempt is
underway to meet these challenges [19].
This work is supported by the Epoch Making Simula-
tion Projects of Earth Simulator Center, the Large Scale
Simulation Program No.132 (FY2005) of High Energy
Accelerator Research Organization (KEK), the Large
Scale Simulation Projects of Academic Computing and
Communications Center of University of Tsukuba, Inter
University Services of Super Computers of Information
Technology Center of University of Tokyo, Super Sinet
Projects of National Institute of Informatics, and also
by the Grant-in-Aid of the Ministry of Education (Nos.
13135204, 13135216, 15540251, 16540228, 16470147,
17340066, 17540259, 18104005, 18540250, 18740130).
[1] S. Aoki et al. [CP-PACS Collaboration], Phys. Rev. Lett.
84, 238 (2000); Phys. Rev. D 67, 034503 (2003).
[2] A. Ali Khan et al. [CP-PACS Collaboration], Phys.
Rev. D 65, 054505 (2002) [Erratum-ibid. D 67, 059901
(2003)]; Phys. Rev. Lett. 85, 4674 (2000) [Erratum-ibid.
90, 029902 (2003)]; S. Aoki et al. [JLQCD Collabora-
tion], Phys. Rev. D 68, 054502 (2003).
[3] C. Bernard et al. [MILC Collaboration], PoS LAT2006,
163 (2006).
[4] Y. Iwasaki, Nucl. Phys. B 258, 141 (1985); Univ. of
Tsukuba Report No. UTHEP-118, 1983.
[5] S. Aoki et al. [CP-PACS/JLQCD Collaboration], Phys.
Rev. D 73, 034501 (2006).
[6] S. Aoki et al. [JLQCD Collaboration], Phys. Rev. D 72,
054510 (2005).
[7] C. Aubin et al. [HPQCD Collaboration, MILC Collab-
oration and UKQCD Collaboration], Phys. Rev. D 70,
031504(R) (2004); Phys. Rev. D 70, 114501 (2004).
[8] S. Aoki et al. [JLQCD Collaboration], Phys. Rev. D 65,
094507 (2002).
[9] T. Kaneko et al. [CP-PACS/JLQCD and ALPHA Col-
laborations], arXiv:hep-lat/0703006.
[10] R. F. Dashen, Phys. Rev. 183, 1245 (1969).
[11] W. M. Yao et al. [Particle Data Group], J. Phys. G 33,
1 (2006).
[12] S. Aoki et al., Phys. Rev. D 58, 074505 (1998).
[13] T. Kaneko et al. [CP-PACS/JLQCD Collaboration],
Nucl. Phys. Proc. Suppl. 129, 188 (2004).
[14] M. Fukugita et al., Phys. Lett. B 294, 380 (1992).
[15] S. R. Sharpe and R. L. . Singleton, Phys. Rev. D 58,
074501 (1998); S. Aoki, Phys. Rev. D 68, 054508 (2003).
[16] S. Aoki et al., Phys. Rev. D 73, 014511 (2006); S. Aoki et
al., Phys. Rev. D 73, 094501 (2006); S. Takeda, Doctral
thesis at Univ. of Tsukuba (2005).
[17] J. Bijnens and J. Prades, Nucl. Phys. B 490, 239 (1997);
J. Bijnens and P. Gosdzinsky, Phys. Lett. B 388, 203
(1996).
[18] D. R. Nelson, arXiv:hep-lat/0212009.
[19] Y. Kuramashi et al. [PACS-CS Collaboration], PoS
LAT2006, 029 (2006)); A. Ukawa et al. [PACS-CS Col-
laboration], PoS LAT2006, 039 (2006).
http://arxiv.org/abs/hep-lat/0703006
http://arxiv.org/abs/hep-lat/0212009
