text:
# Light hadron spectroscopy with two flavors of 𝑂(𝑎)-improved dynamical quarks

## I Introduction

Lattice QCD calculation of the light hadron mass spectrum has witnessed significant progress in recent years Spectrum.review.lat00 ; Spectrum.review.lat01 ; KS-Spectrum.review.lat01 . In the quenched approximation in which the quark vacuum polarization effects are ignored, the CP-PACS Collaboration performed a precise calculation of hadron masses, in which the estimated accuracy reached the level of a few percent in the continuum limit Spectrum.Nf0.CP-PACS . They found that the quenched spectrum shows a significant and systematic deviation from experiment; the $`K^{}`$$`K`$ hyperfine splitting is smaller by about 10 % than experiment. The decuplet baryon mass splittings are also small, and the octet baryon masses are themselves smaller than experiment.

Since this work, the focus of efforts has shifted toward full QCD simulations including vacuum polarization effects of dynamical quarks. A number of simulations now exist, pursued by the SESAM-T$`\chi `$Spectrum.Nf2.SESAM ; Spectrum.Nf2.TchiL ; Spectrum.Nf2.SESAM-TchiL , UKQCD Spectrum.Nf2.UKQCD.csw176 ; Spectrum.Nf2.UKQCD , CP-PACS mq.Nf2.CP-PACS ; Spectrum.Nf2.CP-PACS  and QCDSF-UKQCD Spectrum.Nf2.QCDSF  Collaborations for two flavors using the Wilson-type quark action, and by the MILC Collaboration Spectrum.Nf2.MILC ; Spectrum.Nf3.MILC  for two and three flavors using the Kogut-Susskind(KS) quark action. In particular the CP-PACS calculation mq.Nf2.CP-PACS ; Spectrum.Nf2.CP-PACS  made a first attempt toward execution of the chiral and continuum extrapolations within the single set of simulations, as pioneered by the GF11 Collaboration Spectrum.Nf0.GF11  in their quenched spectrum study. The chief finding of this work was that the $`K^{}`$$`K`$ hyperfine splitting agrees much better with experiment in two-flavor full QCD than in quenched QCD, and that light quark masses decrease by about 25% by the inclusion of dynamical $`u`$ and $`d`$ quarks.

A subtle point with the CP-PACS results is that the dynamical sea quark effects become manifest only after the continuum extrapolation. Further studies are required to consolidate effects of dynamical sea quarks. The CP-PACS simulation used a renormalization group (RG) improved gauge action RGaction , but the $`O(a)`$-improved Wilson quark action SWaction  with only tadpole-improved tadpole\_improvement  clover coefficient. This leaves scaling violation of $`O(g^2a)`$. The use of non-perturbatively determined clover coefficient removing all of $`O(a)`$ errors should be much better to control the continuum extrapolation NPimprovement . Studies along this direction were previously carried out by the UKQCD and QCDSF Collaborations using the plaquette gauge action. However, sea quark effects are not clear in their results of hadron masses, albeit encouraging evidence is seen in the static quark potential Spectrum.Nf2.UKQCD ; Spectrum.Nf2.QCDSF .

In the present work, we explore sea quark effects in hadron and quark masses in two-flavor QCD using the plaquette gauge action and the non-perturbatively $`O(a)`$-improved Wilson quark action. Our simulations are performed at a single lattice spacing $`a^12`$ GeV at $`\beta =5.2`$ using a $`20^3\times 48`$ lattice. We also carry out calculations in quenched QCD with the same action and similar simulation parameters to those in full QCD in order to make a direct comparison between full and quenched QCD. Preliminary results of these calculations have been reported in Refs. Spectrum.Nf2.JLQCD.lat00 ; Spectrum.Nf2.JLQCD.lat01 ; Spectrum.Nf2.JLQCD.lat02.Kaneko ; Spectrum.Nf2.JLQCD.lat02.Hashimoto .

We pay particular attention to two points which are important for an unambiguous identification of sea quark effects. One is the finite spatial volume effect whose magnitude is believed to be more pronounced in full QCD simulations than in quenched QCD FSE.Z3 . An increase of hadron masses due to this effect could mimic sea quark effects. There are only a few studies of finite-size effects in full QCD for the Wilson-type quark action Spectrum.Nf2.TchiL ; Spectrum.Nf2.UKQCD.csw176 . This leads us to perform a systematic investigation of finite-size effects employing $`12^3\times 48`$, $`16^3\times 48`$ and $`20^3\times 48`$ lattices.

Another point is the chiral extrapolation. With currently available computer power and simulation algorithms, the sea quark mass which can be explored with the Wilson-type quark action is limited to values corresponding to $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.6`$. The long extrapolation to the physical $`u`$ and $`d`$ quark masses may involve sizable systematic errors, potentially blurring, or artificially enhancing, sea quark effects. This can be avoided if one examines sea quark effects at the quark masses actually simulated. In this study, therefore, sea quark effects are examined in detail not only at the physical quark mass but also at our simulation points.

We have also made efforts to accumulate high statistics of 12000 trajectories each at five values of sea quark masses. Our implementation of the symmetric preconditioning of the lattice clover-Dirac operator Even-Odd2 ; PHMC.JLQCD  speeded up the configuration generation by a factor two by allowing a doubly larger step size over the even-odd preconditioning.

This paper is organized as follows. We describe details of configuration generation in full and quenched QCD in Sec. II. Method of measurement of hadron masses and the static quark potential is explained in Sec. III. Finite-size effects on hadron masses are discussed in Sec. IV. Section V is devoted to detailed description of the chiral extrapolation of our hadron mass data. We examine sea quark effects in light hadron masses in Sec. VI. Results of the decay constants and quark masses are presented in Secs. VII and VIII. Our conclusion is given in Sec. IX.

## II Simulation Method

### II.1 Simulation parameters and algorithm

We carry out numerical simulations of lattice QCD with two flavors of degenerate dynamical quarks which are identified with the up and down quarks. We use the standard plaquette action for gauge fields defined by

$`S_g`$ $`=`$ $`{\displaystyle \frac{\beta }{6}}{\displaystyle \underset{x,\mu \nu }{}}\text{Tr }U_{x,\mu \nu },`$ (1)
where $`U_{x,\mu \nu }`$ is the product of gauge link variables $`U_{x,\mu }`$ around the plaquette given by

$`U_{x,\mu \nu }`$ $`=`$ $`U_{x,\mu }U_{x+\widehat{\mu },\nu }U_{x+\widehat{\nu },\mu }^{}U_{x,\nu }^{}.`$ (2)
The $`O(a)`$-improved Wilson action SWaction  defined by

$`S_q`$ $`=`$ $`{\displaystyle \underset{x,y}{}}\overline{q}_xD_{xy}q_y`$ (3)
$`D_{xy}`$ $`=`$ $`\delta _{xy}K{\displaystyle \underset{\mu }{}}\left\{\left(1\gamma _\mu \right)U_{x,\mu }\delta _{x+\widehat{\mu },y}+\left(1+\gamma _\mu \right)U_{x,\mu }^{}\delta _{x,y+\widehat{\mu }}\right\}{\displaystyle \frac{1}{2}}Kc_{\mathrm{SW}}\sigma _{\mu \nu }F_{x,\mu \nu }\delta _{xy}`$ (4)
is used for the quark part. The field strength tensor on the lattice is defined by

$`F_{x,\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{8i}}\left\{\left(U_{x,\mu \nu }+U_{x,\mu \nu }+U_{x,\mu \nu }+U_{x,\mu \nu }\right)\left(\text{h.c.}\right)\right\},`$ (5)
where $`\left(\text{h.c.}\right)`$ denotes the hermitian conjugate of the preceding bracket, and $`\sigma _{\mu \nu }=(i/2)[\gamma _\mu ,\gamma _\nu ]`$. The clover coefficient $`c_{\mathrm{SW}}`$ is set to the non-perturbative value determined by the ALPHA Collaboration NPcsw.Nf2.ALPHA .

Our simulations are performed at a single value of $`\beta =5.2`$. The lattice spacing fixed from $`m_\rho `$ at the physical sea quark mass is found to be 0.0887(11) fm. Our value of $`\beta `$ is slightly off the range $`\beta =12.0`$ – 5.4 where the ALPHA Collaboration carried out a non-perturbative determination of $`c_{\mathrm{SW}}`$. We set $`c_{\mathrm{SW}}=2.02`$ by extrapolating their parametrization formula of the non-perturbative $`c_{\mathrm{SW}}`$ as a function of the bare coupling. We performed an independent non-perturbative determination of $`c_{\mathrm{SW}}`$ at $`\beta =5.2`$ and confirmed that our preliminary result $`c_{\mathrm{SW}}=1.98(7)`$ is consistent with 2.02 within the error NPcsw.Nf2.JLQCD .

We employ three lattice sizes that differ in spatial volumes, $`N_s^3\times N_t=12^3\times 48`$, $`16^3\times 48`$ and $`20^3\times 48`$. The hadron spectrum and quark masses at the physical point are calculated using the data on the largest lattice. The data on the two smaller lattices are used to investigate finite-size effects on hadron masses.

On each lattice size, we adopt five values of the sea quark mass corresponding to the hopping parameter $`K_{\mathrm{sea}}=0.1340`$, 0.1343, 0.1346, 0.1350 and 0.1355. This choice covers $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}=0.6`$ – 0.8, and enables us to extrapolate our data to the physical sea quark mass. These simulation parameters are summarized in Table 1.

We note that the UKQCD Collaboration also performed a set of simulations using the same lattice action at $`a^12`$ GeV Spectrum.Nf2.UKQCD . There are, however, some differences in the choice of $`\beta `$ and $`K_{\mathrm{sea}}`$: The UKQCD simulations shift $`\beta `$ with the sea quark mass keeping the Sommer scale $`r_0/a`$ r0  fixed, while our simulations are performed at fixed $`\beta `$. Another difference is the range of the sea quark mass covered in the two simulations. We explore light sea quark masses down to $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.6`$, whilst UKQCD’s lightest point is around $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.7`$. Although the UKQCD Collaboration made another simulation at a smaller sea quark mass $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.6`$ at a spatial extent of $`N_s=16`$ ($`N_sa1.6`$ fm), finite-size effects seem to be significant there (see discussion in Sec. IV).

Gauge configurations are generated using the Hybrid Monte Carlo (HMC) algorithm HMC.Duane ; HMC.Gottlieb . We use simulation programs with three variants of HMC for the $`O(a)`$-improved Wilson action:

* HMC with the even/odd preconditioning Even-Odd.D  only for the inversion of the quark matrix $`D_{xy}`$. This algorithm is used in the simulations on the $`16^3\times 48`$ lattice.
* HMC with the asymmetric preconditioning for the lattice action (A-HMC) Even-Odd1 ; Even-Odd2 ; PHMC.JLQCD . Whole simulations on the $`12^3\times 48`$ lattice are performed with this algorithm.
* HMC with the symmetric preconditioning for the action (S-HMC) Even-Odd2 ; PHMC.JLQCD , which shows the best performance among the three algorithms.

Our main simulation on the $`20^3\times 48`$ lattices is initially started with the A-HMC algorithm, but is later switched to the S-HMC to speed up the calculations. The trajectory length in each HMC step is fixed to the unit length. We use the conventional leap-frog integration scheme for the molecular dynamics equation. The step size $`\mathrm{\Delta }\tau `$ is chosen to achieve an acceptance of 60 – 80%.

The even/odd preconditioned BiCGStab algorithm BiCGStab  is used for the quark matrix inversion to solve the equation $`D_{xy}G_y=B_x`$. We take the stopping condition of the form $`DGB<\mathrm{\Delta }`$ in the HMC program. A modified form $`DGB/B<\mathrm{\Delta }`$ is used in the A-HMC and S-HMC programs. The value of $`\mathrm{\Delta }`$ in the evaluation of the fermionic force is determined so that the reversibility over unit length is satisfied to a relative level better than $`10^{13}`$ for the Hamiltonian. We use a stricter stopping condition in the calculation of the Hamiltonian in the Metropolis accept/reject test. Table 1 shows our choice of $`\mathrm{\Delta }`$ together with the average number of the BiCGStab iteration in the quark matrix inversion for the force calculation, $`N_{\mathrm{inv}}`$.

We accumulate 12000 HMC trajectories at each sea quark mass on the $`20^3\times 48`$ lattice. The statistics on smaller lattices are 3000 trajectories. Measurements of light hadron masses and the static quark potential are carried out at every 10 HMC trajectories. Details of the measurement method will be described in the next section.

All simulations are performed on the parallel computer HITACHI SR8000 model F1 installed at KEK. This machine consists of 100 nodes and has a peak speed of 1.2 TFLOPS and 448 GB of main memory in total. The CPU time needed per unit HMC trajectory on the full machine is listed in Table 1. The total time for configuration generation on each lattice size is 8.6 days on $`12^3\times 48`$, 58 days on $`16^3\times 48`$, and 130 days on $`20^3\times 48`$ lattices. Additional 100 days are spent for the measurement of the hadron masses and the static potential.

### II.2 Simulation in quenched QCD

While many calculations of the hadron spectrum have been performed in quenched QCD, comparisons between our full QCD results and quenched results from other simulations may be subject to systematic uncertainties due to the difference in the simulation details. We therefore carry out a set of quenched calculations of the hadron spectrum using the same lattice actions and simulation parameters as those for full QCD runs.

Our simulations are performed at $`\beta =6.0`$, where the lattice spacing fixed from $`m_\rho `$ equals 0.1074(14) fm. We take $`c_{\mathrm{SW}}=1.769`$ which is the value determined non-perturbatively by the ALPHA Collaboration NPcsw.Nf0.ALPHA . Three lattice sizes $`12^3\times 48`$, $`16^3\times 48`$ and $`20^3\times 48`$ are employed in order to investigate finite-size effects.

Gauge configurations are generated with a combination of the heat-bath and the over-relaxation algorithms. We call four heat-bath sweeps with a succeeding over-relaxation step an iteration. We accumulate statistics of 60000 iterations on each lattice size. Hadron masses and the static potential are calculated at every 200 iterations.

## III Measurement

### III.1 Hadron Masses

In measurements in full QCD, we use six values of the valence quark mass corresponding to the hopping parameter $`K_{\mathrm{val},i}(i=1,\mathrm{}6)=0.1340`$, 0.1343, 0.1346, 0.1350, 0.1355 and 0.1358, which cover the range of $`m_{\mathrm{PS},\mathrm{val}}/m_{\mathrm{V},\mathrm{val}}0.5`$ – 0.8. At each sea quark mass, therefore, there is one value of $`K_{\mathrm{val},i}`$, which equals $`K_{\mathrm{sea}}`$ and is identified as the light quark mass. Other five values of $`K_{\mathrm{val},i}`$ correspond to the mass of strange quarks treated in the quenched approximation. In the following, we use the abbreviation “diagonal data” to represent hadron correlators or masses with a quark mass combination in which all valence quark masses are equal to the sea quark mass.

We employ meson operators defined by

$`M(x)=\overline{q}_x^{(f)}\mathrm{\Gamma }q_x^{(g)},\mathrm{\Gamma }=I,\gamma _5,\gamma _\mu ,\gamma _5\gamma _\mu ,`$ (6)
where $`f`$ and $`g`$ are flavor indices and $`x`$ is the coordinates on the lattice. Meson correlators $`M(x)M(0)^{}`$ are calculated for the following eleven combinations of valence quark masses

$`\begin{array}{c}(K_{\mathrm{val},i},K_{\mathrm{val},i})(i=1,\mathrm{}6),\hfill \\ (K_{\mathrm{sea}},K_{\mathrm{val},i})(i=1,\mathrm{}6,K_{\mathrm{sea}}K_{\mathrm{val},i}).\hfill \end{array}`$ (9)
The former is identified with a degenerate light or strange meson and the latter with a non-degenerate light-strange meson. This choice of the valence quark masses enables us to calculate the full spectrum of strange and non-strange mesons.

For baryons, we use the same operators as those employed in Ref. Spectrum.Nf2.CP-PACS . Namely, the octet baryon operator is defined as

$`O^{fgh}(x)=ϵ^{abc}\left(q_x^{(f)aT}C\gamma _5q_x^{(g)b}\right)q_x^{(h)c},`$ (10)
where $`a,b,c`$ are color indices and $`C=\gamma _4\gamma _2`$ is the charge conjugation matrix. We measure baryon correlators with two types of flavor structure ($`\mathrm{\Sigma }`$\- and $`\mathrm{\Lambda }`$-like baryons),

$`\mathrm{\Sigma }`$ $`:`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(O^{[fh]g}+O^{[gh]f}\right),`$ (11)
$`\mathrm{\Lambda }`$ $`:`$ $`{\displaystyle \frac{1}{\sqrt{6}}}\left(O^{[fh]g}O^{[gh]f}2O^{[fg]h}\right),`$ (12)
where $`O^{[fg]h}=O^{fgh}O^{gfh}`$. Decuplet baryon correlators are calculated using an operator defined by

$`D^{fgh}(x)=ϵ^{abc}\left(q_x^{(f)aT}C\gamma _\mu q_x^{(g)b}\right)q_x^{(h)c}`$ (13)
with symmetrized flavor structure

$`\begin{array}{c}O^{fff},\frac{1}{\sqrt{3}}\left(D^{ffg}+D^{fgf}+D^{gff}\right),\hfill \\ \frac{1}{\sqrt{6}}\left(D^{fgh}+D^{hfg}+D^{ghf}+D^{fhg}+D^{gfh}+D^{hgf}\right).\hfill \end{array}`$ (16)
We take quark mass combinations of $`(K_{\mathrm{val},i},K_{\mathrm{val},i},K_{\mathrm{val},i})(i=1,\mathrm{}6)`$, $`(K_{\mathrm{sea}},K_{\mathrm{val},i},K_{\mathrm{val},i})`$ and $`(K_{\mathrm{sea}},K_{\mathrm{sea}},K_{\mathrm{val},i})(i=1,\mathrm{}6,K_{\mathrm{sea}}K_{\mathrm{val},i})`$ for the baryon correlators.

In order to construct the smeared hadron operators, we measure the wave function of the pseudoscalar (PS) meson

$`\varphi (𝐫)`$ $`=`$ $`{\displaystyle \frac{_𝐱\overline{q}(𝐱,t)\gamma _5q(𝐱+𝐫,t)P(\mathrm{𝟎},0)^{}}{_𝐱P(𝐱,t)P(\mathrm{𝟎},0)^{}}},`$ (17)
where $`P`$ is the PS meson operator, Eq. (6), with $`\mathrm{\Gamma }=\gamma _5`$ and $`t`$ fixed to 12. The measurement is performed at each sea quark mass and lattice size using a subset of gauge configurations (30 configurations every 100 trajectories). We parameterize $`\varphi (𝐫)`$ using a polynomial approximation $`\varphi (𝐫)=1+_{n=1,8}c_n|𝐫|^n`$ and use it as the smearing function. We employ three types of the meson operator : i) local operator, ii) smeared operator $`M(x)=_𝐫\varphi (𝐫)\overline{q}(𝐱,t)\mathrm{\Gamma }q(𝐱+𝐫,t)`$, iii) doubly smeared operator $`M(x)=_{𝐫,𝐫^{}}\varphi (𝐫)\varphi (𝐫^{})\overline{q}(𝐱+𝐫,t)\mathrm{\Gamma }q(𝐱+𝐫^{},t).`$ Additionally, we use “triply smeared operator”

$`O^{fgh}(x)`$ $`=`$ $`{\displaystyle \underset{𝐫_1,𝐫_2,𝐫_3}{}}\varphi (𝐫_1)\varphi (𝐫_2)\varphi (𝐫_3)ϵ^{abc}\left(q^{aT}(𝐱+𝐫_1,t)C\mathrm{\Gamma }q^b(𝐱+𝐫_2,t)\right)q^c(𝐱+𝐫_3,t)`$ (18)
for baryons. Hadron correlators are measured with a) point source and sink operators, b) smeared source and point sink, and c) smeared source and sink operators. We fix configurations to the Coulomb gauge, since b) and c) are not gauge invariant.

We observe that, when valence quarks are lighter than sea quarks, the hadron correlator takes an exceptionally large value on a small number of configurations. This might be caused by a fluctuation of the lowest eigenvalue of the Dirac operator of the $`O(a)`$-improved Wilson action. If the PS meson correlator on the $`i`$-th gauge configuration takes a value larger than 20 times the statistical average, which is evaluated without that configuration, at a certain time slice

$`P(x)P(0)^{}_i>{\displaystyle \frac{20}{N_{\mathrm{conf}}1}}{\displaystyle \underset{k=1,ki}{\overset{N_{\mathrm{conf}}}{}}}P(x)P(0)^{}_k`$ (19)
where $`N_{\mathrm{conf}}`$ is the total number of configurations, we consider it as an exceptional configuration and remove it from the following analysis. The number of the removed configurations is given in Table 1.

In order to reduce the statistical fluctuation of hadron correlators on the $`20^3\times 48`$ lattice, we repeat the measurement for two choices of the location of the hadron source, $`t_{\mathrm{src}}=1`$ and $`N_t/2+1(=25)`$, and take the average over the two sources:

$`{\displaystyle \frac{1}{2}}\left(M(t_{\mathrm{src}}+t)M(t_{\mathrm{src}})^{}_{t_{\mathrm{src}}=1}+M(t_{\mathrm{src}}+t)M(t_{\mathrm{src}})^{}_{t_{\mathrm{src}}=25}\right).`$ (20)
We find that this procedure reduces the statistical error of hadron correlators by typically 20 %, which suggests that the statistics is increased effectively by a factor of 1.5. For further reduction of the statistical fluctuation, we take the average over three polarization states for vector mesons, two spin states for octet baryons and four spin states for decuplet baryons.

Figures 1 – 4 show examples of effective mass plots. We find that the best plateau of the effective mass is obtained from hadron correlators with the point sink and the doubly smeared source for mesons and the triply smeared one for baryons. Therefore, hadron masses are extracted from these types of correlators.

We carry out $`\chi ^2`$ fits to hadron correlators by taking account of correlations among different time slices. A single hyperbolic cosine form is assumed for mesons, and a single exponential form for baryons. The lower cut of the fit range $`t_{\mathrm{min}}`$ is determined by inspecting stability of the fitted mass. The upper cut ($`t_{\mathrm{max}}`$) dependence of the fit results is small and, therefore, we fix $`t_{\mathrm{max}}`$ to $`N_t/2`$ for all hadrons. Our choice of fit ranges and resulting hadron masses are summarized in Tables LABEL:tab:meas:had:meson\_12x48 – LABEL:tab:meas:had:baryon\_20x48 in Appendix A. Statistical errors of hadron masses are estimated with the jack-knife procedure. We adopt the bin size of 100 trajectories by inspecting the bin size dependence of the jack-knife error as discussed in Sec.III.3.

In Fig. 5, we test double exponential fits to extract hadron masses at $`K_{\mathrm{sea}}=K_{\mathrm{val}}=0.1355`$. While these fits are unstable and lead to a large error for the mass of the first excited state, the result for the ground state mass is consistent with that from the single exponential fit. The situation is similar at other sea and valence quark masses. This suggests that the hadron masses in Tables LABEL:tab:meas:had:meson\_12x48 – LABEL:tab:meas:had:baryon\_20x48 and the light hadron spectrum calculated from these results have small contamination from excited states.

Hadron correlators in quenched QCD are calculated in an analogous manner. We use six values of $`K_{\mathrm{val}}`$, 0.13260, 0.13290, 0.13331, 0.13384, 0.13432 and 0.13465, corresponding to $`m_{\mathrm{PS},\mathrm{val}}/m_{\mathrm{V},\mathrm{val}}0.50`$ – 0.80 and the hadron operators and smearing procedure same as those in the full QCD study. A difference is that we can take more combinations of valence quark masses than in full QCD, since any value of the six valence quark masses can be identified with either light or strange quark mass. We take all combinations $`(K_{\mathrm{val},i},K_{\mathrm{val},j})(i,j=1,..,6)`$ for mesons, and somewhat restricted choices $`(K_{\mathrm{val},i},K_{\mathrm{val},i},K_{\mathrm{val},j})(i,j=1,..,6)`$ for baryons. Statistical errors are estimated with the jack-knife procedure with bin size of 200 iterations. The exceptional configurations are discarded with the same criterion as defined in Eq. (19). Results of hadron masses are collected in Tables LABEL:tab:meas:had:meson\_q12x48 – LABEL:tab:meas:had:baryon\_q20x48 in Appendix A.

### III.2 Static quark potential

We calculate the static quark potential in order to determine the Sommer scale r0  which we use in our analysis of hadron masses. For this purpose, the temporal Wilson loops $`W(r,t)`$ up to $`t=16`$ and $`r=(\sqrt{3}N_s/2)`$ are measured both in full and quenched QCD simulations. We apply the smearing procedure of Ref. Potential.Nf0.Bali  up to twelve steps and the measurements are carried out every four steps.

The static quark potential $`V(r)`$ is determined from the correlated fit of the form

$`W(r,t)`$ $`=`$ $`C(r)\mathrm{exp}\left[V(r)t\right].`$ (21)
We take the fit range $`[t_{\mathrm{min}},t_{\mathrm{max}}]=[3,7]`$ in all simulations in full and quenched QCD by inspecting the $`t`$ dependence of the effective potential

$`V_{\mathrm{eff}}(r,t)=\mathrm{ln}\left[W(r,t)/W(r,t+1)\right].`$ (22)
Examples of $`V_{\mathrm{eff}}`$ are plotted in Fig. 6. For each $`r`$, the number of smearing steps is fixed to its optimum value at which the overlap to the ground state $`C(r)`$ takes the largest value.

As shown in Fig.7, we do not observe any clear indication of the string breaking. Therefore $`V(r)`$ is fitted to the form

$`V(r)`$ $`=`$ $`V_0{\displaystyle \frac{\alpha }{r}}+\sigma r+\delta V(r),`$ (23)
where $`\delta V(r)`$ is the lattice correction to the Coulomb term calculated perturbatively from one lattice-gluon exchange diagram deltaV

$`\delta V(r)`$ $`=`$ $`g_c\left(G(𝐫){\displaystyle \frac{1}{r}}\right),`$ (24)
$`G(𝐫)`$ $`=`$ $`4\pi {\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{cos}\left[\mathrm{𝐤𝐫}\right]}{4_{i=1}^3\mathrm{sin}^2\left[k_i/2\right]}}.`$ (25)
The Sommer scale $`r_0`$ defined through r0

$`r_0^2{\displaystyle \frac{dV(r)}{dr}}|_{r=r_0}=1.65`$ (26)
is then determined from the parametrization of the corrected potential $`V(r)\delta V(r)`$:

$`r_0`$ $`=`$ $`\sqrt{{\displaystyle \frac{1.65\alpha }{\sigma }}}.`$ (27)
The lower cut of the fit range in Eq. (23), $`r_{\mathrm{min}}`$, is determined by inspecting the $`r_{\mathrm{min}}`$ dependence of $`r_0`$. We observe that $`r_0`$ is relatively stable for $`r_{\mathrm{min}}[\sqrt{2},2\sqrt{2}]`$ as shown in Fig.8. With $`r_{\mathrm{min}}<\sqrt{2}`$, $`\chi ^2/\mathrm{dof}`$ takes an unacceptably large value due to the violation of rotational symmetry, while $`\alpha `$ becomes ill-determined with $`r_{\mathrm{min}}>2\sqrt{2}`$. We therefore take $`r_{\mathrm{min}}=\sqrt{5}`$. While the $`r_{\mathrm{max}}`$ dependence of $`r_0`$ is rather mild, the covariance matrix becomes ill-determined with $`r_{\mathrm{max}}`$ greater than $`9\sqrt{2}`$ on $`20^3\times 48`$, $`7\sqrt{2}`$ on $`16^3\times 48`$ and $`6\sqrt{2}`$ on $`12^3\times 48`$. We therefore fix $`r_{\mathrm{max}}`$ to these values.

We repeat the fits, Eqs. (21) and (23), with other choices of the range: $`t_{\mathrm{min}}=4`$ or $`r_{\mathrm{min}}[\sqrt{2},2\sqrt{2}]`$. The largest deviations in the fit parameters and $`r_0`$ are included into their systematic errors. Other systematic errors due to the choice of $`t_{\mathrm{max}}`$, the optimum number of the smearing steps and $`r_{\mathrm{max}}`$ are small and ignored. Fit parameters in Eq. (23) and $`r_0`$ are summarized in Table 2 for full QCD, and in Table 3 for quenched QCD.

### III.3 Autocorrelation

The autocorrelation in our full QCD data is studied by calculating the cumulative autocorrelation time

$`\tau _𝒪^{\mathrm{cum}}(\mathrm{\Delta }t_{\mathrm{max}})={\displaystyle \frac{1}{2}}+{\displaystyle \underset{\mathrm{\Delta }t=1}{\overset{\mathrm{\Delta }t_{\mathrm{max}}}{}}}\rho _𝒪(\mathrm{\Delta }t),`$ (28)
where $`\rho _𝒪(t)`$ is the autocorrelation function

$`\rho _𝒪(\mathrm{\Delta }t)={\displaystyle \frac{\mathrm{\Gamma }_𝒪(\mathrm{\Delta }t)}{\mathrm{\Gamma }_𝒪(0)}},`$ (29)
$`\mathrm{\Gamma }_𝒪(\mathrm{\Delta }t)=\left(𝒪(t)𝒪\right)\left(𝒪(t+\mathrm{\Delta }t)𝒪\right)`$ (30)
and we take $`\mathrm{\Delta }t_{\mathrm{max}}=200`$.

In Table 4, we summarize $`\tau _𝒪^{\mathrm{cum}}`$ in the A-HMC and S-HMC simulations on the $`20^3\times 48`$ lattices for three quantities: i) the plaquette which is measured at every trajectory, ii) the PS meson propagator at $`t=12`$, iii) the temporal Wilson loop with $`(r,t)=(5,4)`$. The results do not show any systematic differences in $`\tau _𝒪^{\mathrm{cum}}`$ between the A-HMC and S-HMC runs. The plaquette shows the largest autocorrelation with $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}=10`$ – 30, which is similar to those found in the UKQCD simulation Spectrum.Nf2.UKQCD  using the same lattice action and similar simulation parameters. We obtain smaller values of $`\tau _𝒪^{\mathrm{cum}}`$ for the other two quantities. This is contrary to a naive expectation that these long-distance observables have a longer autocorrelation than the local quantity like the plaquette. This suggests that the size of noise arising from short correlation modes is larger than that of the longest mode in these observables and our statistics are not sufficient to extract $`\tau _𝒪^{\mathrm{cum}}`$ of the longest but weak mode.

The statistical error including the effect of autocorrelation is given by $`\sqrt{2\tau _𝒪^{\mathrm{cum}}}`$ times the naive error. Therefore, the above observation tells us that the bin size in the jack-knife procedure of 60 HMC trajectories or larger is a safe choice to take account of the autocorrelation in our data.

The bin size dependence of the jack-knife error of hadron masses and the static potential is plotted in Figs. 9 and 10. We use errors obtained from uncorrelated fits because, with large bin sizes, the number of bins would not be sufficiently large to determine the covariance matrix reliably. For both hadron masses and the static potential, the jack-knife error reaches its plateau at bin size of 50–100 trajectories, which is roughly consistent with the above estimate from $`\tau _𝒪^{\mathrm{cum}}`$. The situation is similar on smaller volumes $`16^3\times 48`$ and $`12^3\times 48`$. We therefore take the bin size of 100 trajectories in the error analysis in full QCD.

We also investigate the bin size dependence of the jack-knife error in quenched QCD. As shown in Figs. 11 and 12, the bin size of 200 iterations is reasonable.

Another point of interest is the sea quark mass dependence of the autocorrelation. A natural expectation is that smaller sea quark mass leads to a larger correlation length and hence a longer autocorrelation. This expectation is supported by the CP-PACS observation in Ref. Spectrum.Nf2.CP-PACS , where they used the RG-improved gauge and clover quark actions. However, our result of $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ in the S-HMC simulations, which is determined more precisely than that for the A-HMC due to the higher statistics, shows the contrary sea quark mass dependence: $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ decreases as the sea quark mass decreases. This is consistent with the UKQCD’s observation in Ref. Spectrum.Nf2.UKQCD . We also note that $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ in our simulations is much larger than in the CP-PACS’s runs particularly at the heaviest sea quark masses $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.8`$ ($`K_{\mathrm{sea}}=0.1340`$).

In our determination of non-perturbative $`c_{\mathrm{SW}}`$ at $`\beta =5.2`$ NPcsw.Nf2.JLQCD , we find that the expectation value of the plaquette varies rapidly around $`K_{\mathrm{sea}}0.132`$, where the plaquette shows the strongest autocorrelation in the investigated region $`K[0.100,0.136]`$. Since such a behavior, somewhat similar to a phase transition, is not observed at higher $`\beta `$, we consider the unexpected behavior of the plaquette to be an artifact due to finite lattice spacing. This artifact is probably absent or well suppressed with the CP-PACS’s choice of the improved actions. At sufficiently small lattice spacings, we then expect that $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ shows the natural sea quark mass dependence, namely larger $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ for lighter sea quark masses. We also expect that, even at $`\beta =5.2`$, $`\tau _{\mathrm{plaq}}^{\mathrm{cum}}`$ will increase if the sea quark mass becomes sufficiently small compared to that corresponding to $`K_{\mathrm{sea}}=0.132`$.

## IV Finite-size Effects

Finite-size effects (FSE) are one of the major sources of systematic errors in lattice calculations. Since our largest volume size $`(\text{1.8 fm})^3`$ is still not so large, it is important to check FSE in our data. We discuss how much FSE is present in our data on the largest lattice using data on three spatial volumes $`12^3\times 48`$, $`16^3\times 48`$ and $`20^3\times 48`$.

In Figs. 13 – 15, we plot diagonal data of hadron masses as a function of the spatial volume inverse. For $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.7`$, including the quenched case, hadron masses obtained on the $`16^3\times 48`$ and $`20^3\times 48`$ lattices are consistent with each other within two standard deviations. On the other hand, hadron masses decrease monotonously up to $`V=20^3`$ at the lightest sea quark mass corresponding to $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.6`$.

The magnitude of FSE also depends on the valence quark mass. Figure 16 shows the valence quark mass dependence of the relative mass shift between the two larger lattices for PS mesons and octet baryons. We observe that, except at the lightest sea quark mass, the mass shift is at most a few percent level in the whole range of the simulated valence quark mass. The situation is similar for vector mesons and decuplet baryons. Therefore, we conclude that the size of FSE on our largest lattice is small over our range of valence quark masses down to the second lightest sea quark mass.

The mass shift is non-negligible at the lightest sea quark mass. While the magnitude is of the order of a few percent for the heaviest valence quarks, it clearly increases as the valence quark mass decreases.

We consider that the observed FSE is caused by valence quarks wrapping around the lattice in spatial directions (namely squeezing of hadrons into the small box) rather than wrapping of virtual pions. As shown in Fig. 13, the magnitude of FSE caused by the effects of virtual pions (long dashed line) FSE.expL  given by $`m_{\mathrm{PS}}(L)m_{\mathrm{PS}}(L=\mathrm{})\mathrm{exp}\left[m_{\mathrm{PS}}(L=\mathrm{})L\right]`$ with $`L=N_sa`$ is too small compared to observed effects.

A qualitative understanding of the observed FSE is as follows. The wrapping of valence quarks is suppressed by the center $`Z(3)`$ symmetry in quenched QCD FSE.Z3 . In full QCD, $`Z(3)`$ symmetry is broken by the wrapping of sea quarks in the spatial directions, whose magnitude increases toward lighter sea quark. A possible reason why FSE is significant only at our lightest sea quark mass would be that the $`Z(3)`$ breaking turns on rather quickly around the lightest sea quarks.

The enhancement of FSE toward the lighter valence quarks leads to a decrease of the slope $`dm_{\mathrm{had}}/dm_{\mathrm{PS}}^2`$ in Fig. 17 and, hence, underestimation of hadron mass splittings, such as the $`K^{}`$-$`K`$ hyperfine splitting. The mass splittings are expected to be increased by sea quark effects, since these are underestimated in quenched QCD as well established in Ref. Spectrum.Nf0.CP-PACS . Therefore, FSE makes sea quark effects less clear. It is crucial to check how large FSE is in our hadron mass data at the lightest sea quark mass on the largest lattice.

Figures 13 and 14 show that the volume dependence of our data is well described by a power law

$`m_{\mathrm{had}}(L)=m_{\mathrm{had}}(L=\mathrm{})+c/L^3`$ (31)
as found in Ref. FSE.Vinv  using the KS fermion. The relative size of FSE on the largest lattice

$`\mathrm{\Delta }m={\displaystyle \frac{m(L=20a)m_{\mathrm{had}}(L=\mathrm{})}{m(L=20a)}}`$ (32)
is estimated from this ansatz and is plotted in Fig. 18. We find that, for PS and vector mesons, $`\mathrm{\Delta }m`$ is about 5 % for diagonal data and is reduced to a few percent at $`K_{\mathrm{val}}=0.1350`$, which roughly corresponds to the strange quark mass $`m_s`$. It is expected, however, that the volume dependence (31) turns into a milder form $`\mathrm{exp}[m_{\mathrm{PS}}L]`$ for sufficiently large volumes. The actual size of FSE should be smaller than the above estimation, say, a few percent. Since this is smaller than the typical size of quenching errors, which is 5 – 10%, we consider that the examination of sea quark effects is feasible in the meson sector, particularly in strange meson masses.

Finite-size effects are more pronounced for baryon masses as observed in Fig. 18. For diagonal data, $`\mathrm{\Delta }m`$ is roughly comparable with typical quenching errors in the baryon spectrum of the order of 5 – 10 %. Sea quark effects in the light baryon masses, such as $`m_N`$ and $`m_\mathrm{\Delta }`$, may become unclear by the contamination of FSE. We note, however, that $`\mathrm{\Delta }m`$ decreases for heavier valence quark masses. The examination of sea quark effects becomes more feasible for strange baryon masses like $`m_\mathrm{\Xi }`$ and $`m_\mathrm{\Omega }`$.

Figure 19 shows $`\mathrm{\Delta }m`$ for the diagonal data of the octet baryon mass as a function of $`N_s`$. We find that the size of $`N_s30`$, which corresponds to $`L2.7`$ fm, is required to suppress $`\mathrm{\Delta }m`$ to a few percent level. The required size becomes slightly smaller, $`L2.4`$ fm, for the valence quark mass around $`m_s`$. These sizes are larger than our largest spatial lattice size $`L1.8`$ fm. Further simulations on such large lattices will be needed to obtain a definite conclusion on sea quark effects in the baryon spectrum.

Tables 2 and 3 show that FSE in $`r_0`$ are much smaller than in hadron masses. While the central value on $`20^3`$ is systematically higher than that on $`16^3`$ in both full and quenched QCD, the difference is about 1 % and not significant with the accuracy of our data. The size of FSE is small, namely a few percent level, even on $`12^3`$. Therefore FSE in $`r_0`$ on the largest lattice can be safely neglected both in full and quenched QCD.

## V Chiral Extrapolation

The hadron spectrum in full QCD is calculated using hadron masses measured on the $`20^3\times 48`$ lattice. This requires a parametrization of the mass data as a function of sea and valence quark masses in order to extrapolate (up-down) or interpolate (strange) quark masses to their physical values. We make this parametrization by combined fits to masses of a given hadron at all sea quark masses. We test the following two methods for the combined fit:

1. The effective lattice spacing, determined from $`r_0`$ for instance, may vary as a function of the sea quark mass. In order to separate this effect from the physical quark mass dependence, we carry out the chiral extrapolation using dimensionless quantities such as $`(r_0(K_{\mathrm{sea}})m_{\mathrm{had}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2}))`$, where $`m_{\mathrm{had}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})`$ represents the measured hadron mass composed of valence quark masses corresponding to $`K_{\mathrm{val},1}`$ and $`K_{\mathrm{val},2}`$ on the gauge configurations generated at $`K_{\mathrm{sea}}`$. We refer to this way as “method-A”.
2. It is also possible to fit hadron masses in lattice units, as was done by the SESAM Spectrum.Nf2.SESAM  and the CP-PACS Spectrum.Nf2.CP-PACS  Collaborations. We call this “method-B”.

A detailed description of the two methods will be given in Secs. V.1 and V.2. They should yield a consistent hadron spectrum, since fit forms in method-B can be reproduced from those in method-A by expanding $`r_0`$ as a function of sea quark mass. This consistency is examined in Sec. V.2.

### V.1 chiral extrapolation using $`r_0`$ (method-A)

Chiral perturbation theory (ChPT) ChPT  provides a guide to obtain a controlled chiral limit of hadron mass data. For the quark mass dependence of diagonal data of the PS meson mass, ChPT predicts the presence of logarithmic singularities. At the one-loop level, the ChPT prediction reads

$`{\displaystyle \frac{m_{\mathrm{PS}}^2}{2B_0m_q}}`$ $`=`$ $`1+{\displaystyle \frac{1}{N_f}}y\mathrm{ln}\left[y\right]+Ay`$ (33)
where $`y=2B_0m_q/(4\pi f)^2`$ and $`A`$ is a linear combination of the low energy constants $`\alpha _i`$ of ChPT Lagrangian : $`A=\left(2\alpha _8\alpha _5\right)+N_f\left(2\alpha _6\alpha _4\right)`$. The mass ratio on the left hand-side, $`m_{\mathrm{PS}}^2/m_\mathrm{q}`$, is plotted as a function of $`m_qy`$ in Fig. 20. As we already reported Spectrum.Nf2.JLQCD.lat01 ; Spectrum.Nf2.JLQCD.lat02.Hashimoto , our data show no hint of the curvature predicted by the chiral logarithm. The fit of Eq. (33), assuming $`f`$ to be a free parameter, gives $`f6`$ GeV, which is much larger than its experimental value 93 MeV. On the other hand, the fit gives an unacceptably large $`\chi ^2/\mathrm{dof}=O(100)`$, if we fix $`f`$ to the experimental value.

A similar test using formulae from partially quenched ChPT (PQChPT) PQChPT  also shows that the coefficient of the chiral logarithm term obtained from our data is much smaller than the prediction from PQChPT Spectrum.Nf2.JLQCD.lat01 ; Spectrum.Nf2.JLQCD.lat02.Hashimoto ; ChPT\_test.Nf2.JLQCD . A possible reason for the absence of the chiral logarithm is that the sea quark mass in our simulations is still too large and higher order corrections of ChPT should be included to describe the data.

In this study, therefore, we use simple polynomial fitting forms in terms of the quark mass for the chiral extrapolation. The systematic error due to the chiral extrapolation is estimated by testing several different polynomial forms. However, the inconsistency between our data and ChPT suggests that the extrapolation may have larger uncertainty than this estimation. This point will be examined in detail in a separate paper ChPT\_test.Nf2.JLQCD .

Since the sea quark mass in our simulations is not so small as discussed above, it is important to check the convergence property of the polynomial expansion of hadron masses in our range of the sea quark mass. We carry out both quadratic and cubic chiral fits to diagonal data of PS and vector meson masses:

$`(r_0(K_{\mathrm{sea}})m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}}))^2`$ $`=`$ $`B_{\mathrm{diag}}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}+C_{\mathrm{diag}}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}^2+D_{\mathrm{diag}}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}^3,`$ (34)
$`r_0(K_{\mathrm{sea}})m_\mathrm{V}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}})`$ $`=`$ $`A_{\mathrm{diag}}^\mathrm{V}+B_{\mathrm{diag}}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}+C_{\mathrm{diag}}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}^2+D_{\mathrm{diag}}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}^3,`$ (35)
where $`\mu _{\mathrm{q},\mathrm{diag}}`$ and $`\mu _{\mathrm{PS},\mathrm{sea}}`$ are the quark mass defined through the vector Ward identity (VWI) and the PS meson mass normalized by $`r_0`$:

$`\mu _{\mathrm{q},\mathrm{sea}}`$ $`=`$ $`r_0(K_{\mathrm{sea}})m_{\mathrm{q},\mathrm{sea}},`$ (36)
$`m_{\mathrm{q},\mathrm{sea}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{K_{\mathrm{sea}}}}{\displaystyle \frac{1}{K_c}}\right),`$ (37)
$`\mu _{\mathrm{PS},\mathrm{sea}}`$ $`=`$ $`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS},\mathrm{sea}}\right)^2,`$ (38)
$`m_{\mathrm{PS},\mathrm{sea}}`$ $`=`$ $`m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}}).`$ (39)
Fit parameters and $`\chi ^2/\mathrm{dof}`$ are collected in Tables 5 and 6. The coefficient of the cubic term is small and consistent with zero for both Eqs. (34) and (35). Consequently, the quadratic and cubic fits show a good consistency with each other in the whole range of the quark mass and toward the chiral limit, as seen in Fig. 21. These observations suggest that the polynomial expansion up to the quadratic order is sufficient to describe the quark mass dependence of our data of meson masses in the method-A.

We carry out a combined fit to PS meson masses as a function of sea and valence quark masses using the quadratic form

$`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})\right)^2`$ (40)
$`=`$ $`B_\mathrm{s}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}+B_\mathrm{v}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{val}}+C_\mathrm{s}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}^2+C_\mathrm{v}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{val}}^2+C_{\mathrm{sv}}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{sea}}\mu _{\mathrm{q},\mathrm{val}},`$
where $`\mu _{\mathrm{q},\mathrm{sea}}`$ is defined in Eq. (36) and

$`\mu _{\mathrm{q},\mathrm{val}}`$ $`=`$ $`r_0(K_{\mathrm{sea}})m_{\mathrm{q},\mathrm{val}},`$ (41)
$`m_{\mathrm{q},\mathrm{val}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m_{\mathrm{q},\mathrm{val},1}+m_{\mathrm{q},\mathrm{val},2}\right),`$ (42)
$`m_{\mathrm{q},\mathrm{val},i}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{K_{\mathrm{val},i}}}{\displaystyle \frac{1}{K_c}}\right)(i=1,2).`$ (43)
The presence of the monomial term in $`m_{\mathrm{q},\mathrm{sea}}`$ means that the PS meson mass does not vanish in the chiral limit $`m_{\mathrm{q},\mathrm{val}}=0`$ for non-zero values of $`m_{\mathrm{q},\mathrm{sea}}`$. This is because the value of $`K_{\mathrm{val}}`$ where the PS meson mass vanishes depends on the sea quark mass due to explicit violation of chiral symmetry with the Wilson-type quark action.

We employ uncorrelated fits in the combined chiral extrapolations although the data with the same sea quark mass are expected to be correlated. Therefore, the obtained $`\chi ^2/\mathrm{dof}`$ can be considered only as a guide to judge the quality of the fit. Figure 22 shows that this fit form describes our data well. Parameters of the fit are summarized in Table 7. We note that $`K_c`$ determined from the diagonal fit Eq. (34) and the combined fit Eq. (40) are consistent with each other, as they should be.

The most general quadratic fit ansatz for the PS meson masses should include an additional cross term

$`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})\right)^2`$ $`=`$ $`\text{“r.h.s of Eq. (}\text{40}\text{)”}+C_{\mathrm{vv}}^{\mathrm{PS}}\mu _{\mathrm{q},\mathrm{val},1}\mu _{\mathrm{q},\mathrm{val},2},`$ (44)
where $`\mu _{\mathrm{q},\mathrm{val},i}=r_0(K_{\mathrm{sea}})m_{\mathrm{q},\mathrm{val},i}`$. However, the coefficient is small as shown in fit parameters in Table 7, and hence does not change the hadron spectrum. We use this fit to estimate the systematic error due to the choice of the fitting function.

For the vector meson, we find that the following form describes our data well

$`r_0(K_{\mathrm{sea}})m_\mathrm{V}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})`$ $`=`$ $`A^\mathrm{V}+B_\mathrm{s}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}+B_\mathrm{v}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{val}}+C_{\mathrm{sv}}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val}},`$ (45)
where $`\mu _{\mathrm{PS},\mathrm{sea}}`$ is defined in Eq. (38) and

$`\mu _{\mathrm{PS},\mathrm{val}}`$ $`=`$ $`(r_0(K_{\mathrm{sea}})m_{\mathrm{PS},\mathrm{val}})^2,`$ (46)
$`m_{\mathrm{PS},\mathrm{val}}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m_{\mathrm{PS},\mathrm{val},1}^2+m_{\mathrm{PS},\mathrm{val},2}^2\right),`$ (47)
$`m_{\mathrm{PS},\mathrm{val},i}`$ $`=`$ $`m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},i},K_{\mathrm{val},i}).`$ (48)
For a more general fit of form

$`r_0(K_{\mathrm{sea}})m_\mathrm{V}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})`$ $`=`$ $`\text{“r.h.s of Eq. (}\text{45}\text{)”}+C_\mathrm{s}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{sea}}^2+C_\mathrm{v}^\mathrm{V}\mu _{\mathrm{PS},\mathrm{val}}^2,`$ (49)
the additional parameters $`C_\mathrm{s}^\mathrm{V}`$ and $`C_\mathrm{v}^\mathrm{V}`$ are not well-determined as seen in Table 8. We use the former fit, which is shown in Fig. 23, to calculate the hadron spectrum, and the latter to estimate systematic error of the chiral extrapolation.

We also carry out a partially quenched fit to vector meson masses at each sea quark mass. We use a linear form, which is obtained from Eq. (45) by dropping all terms describing the sea quark mass dependence

$`m_\mathrm{V}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})=A_{\mathrm{PQ}}^\mathrm{V}+B_{\mathrm{PQ}}^\mathrm{V}m_{\mathrm{PS},\mathrm{val}}^2.`$ (50)
Parameters given in Table 9 are used to calculate the $`J`$ parameter at each sea quark mass.

The chiral extrapolation of octet baryon masses is carried out using a quadratic form based on the leading order prediction of ChPT ChPT.octet , which was also used in Ref. Spectrum.Nf2.CP-PACS . We carry out the simultaneous fit to the $`\mathrm{\Sigma }`$\- and $`\mathrm{\Lambda }`$-like octet baryon masses using the functions

$`r_0(K_{\mathrm{sea}})m_{\mathrm{oct},\mathrm{\Sigma }}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$ (52)
$`=`$ $`A^\mathrm{O}+B_\mathrm{s}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{sea}}+(F_\mathrm{v}^\mathrm{O}D_\mathrm{v}^\mathrm{O})\mu _{\mathrm{PS},\mathrm{val},1}+2F_\mathrm{v}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{val},2}+C_\mathrm{s}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{sea}}^2+(C_\mathrm{v}^\mathrm{O}+C_\mathrm{v}^{\mathrm{O},\mathrm{\Sigma }})\mu _{\mathrm{PS},\mathrm{val},1}^2`$
$`+(C_\mathrm{v}^\mathrm{O}C_\mathrm{v}^{\mathrm{O},\mathrm{\Sigma }})\mu _{\mathrm{PS},\mathrm{val},2}^2+(C_{\mathrm{sv}}^\mathrm{O}+C_{\mathrm{sv}}^{\mathrm{O},\mathrm{\Sigma }})\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val},1}+(C_{\mathrm{sv}}^\mathrm{O}C_{\mathrm{sv}}^{\mathrm{O},\mathrm{\Sigma }})\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val},2},`$
$`r_0(K_{\mathrm{sea}})m_{\mathrm{oct},\mathrm{\Lambda }}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$
$`=`$ $`A^\mathrm{O}+B_\mathrm{s}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{sea}}+\left(F_\mathrm{v}^\mathrm{O}+{\displaystyle \frac{D_\mathrm{v}^\mathrm{O}}{3}}\right)\mu _{\mathrm{PS},\mathrm{val},1}+2\left(F_\mathrm{v}^\mathrm{O}{\displaystyle \frac{2}{3}}D_\mathrm{v}^\mathrm{O}\right)\mu _{\mathrm{PS},\mathrm{val},2}`$
$`+C_\mathrm{s}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{sea}}^2+(C_\mathrm{v}^\mathrm{O}+C_\mathrm{v}^{\mathrm{O},\mathrm{\Lambda }})\mu _{\mathrm{PS},\mathrm{val},1}^2+(C_\mathrm{v}^\mathrm{O}C_\mathrm{v}^{\mathrm{O},\mathrm{\Lambda }})\mu _{\mathrm{PS},\mathrm{val},2}^2+(C_{\mathrm{sv}}^\mathrm{O}+C_{\mathrm{sv}}^{\mathrm{O},\mathrm{\Lambda }})\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val},1}`$
$`+(C_{\mathrm{sv}}^\mathrm{O}C_{\mathrm{sv}}^{\mathrm{O},\mathrm{\Lambda }})\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val},2},`$
where $`\mu _{\mathrm{PS},\mathrm{val},i}=(r_0(K_{\mathrm{sea}})m_{\mathrm{PS},\mathrm{val},i})^2`$.

The decuplet baryon masses are well described by the following form

$`r_0(K_{\mathrm{sea}})m_{\mathrm{dec}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})=A^\mathrm{D}+B_\mathrm{s}^\mathrm{D}\mu _{\mathrm{PS},\mathrm{sea}}+B_\mathrm{v}^\mathrm{D}\mu _{\mathrm{PS},\mathrm{val}}+C_{\mathrm{sv}}^\mathrm{D}\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val}},`$ (53)
where $`\mu _{\mathrm{PS},\mathrm{val}}`$ stands for the average of three valence quark masses

$`\mu _{\mathrm{PS},\mathrm{val}}`$ $`=`$ $`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS},\mathrm{val}}\right)^2,`$ (54)
$`m_{\mathrm{PS},\mathrm{val}}^2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left(m_{\mathrm{PS},\mathrm{val},1}^2+2m_{\mathrm{PS},\mathrm{val},2}^2\right).`$ (55)
Figures 24 and 25 show the fit for octet and decuplet baryon masses. Parameters are summarized in Tables 10 and 11. We also test the following forms to estimate systematic error of the baryon spectrum due to the choice of the fitting form

$`r_0(K_{\mathrm{sea}})m_{\mathrm{oct}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$
$`=`$ $`\text{“r.h.s of Eqs. (}\text{52}\text{) and (}\text{52}\text{)”}+C_{\mathrm{vv}}^\mathrm{O}\mu _{\mathrm{PS},\mathrm{val},1}\mu _{\mathrm{PS},\mathrm{val},2},`$
$`r_0(K_{\mathrm{sea}})m_{\mathrm{dec}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$
$`=`$ $`\text{“r.h.s of Eq. (}\text{53}\text{)”}+C_\mathrm{s}^\mathrm{D}\mu _{\mathrm{PS},\mathrm{sea}}^2+C_\mathrm{v}^\mathrm{D}\mu _{\mathrm{PS},\mathrm{val}}^2.`$ (57)
We carry out the chiral extrapolation of $`r_0(K_{\mathrm{sea}})`$ in order to determine $`r_0`$ at the physical sea quark mass, which is required to calculate the hadron spectrum in the method-A. We use a linear form

$`{\displaystyle \frac{1}{r_0}}=A_{r_0}+{\displaystyle \frac{B_{r_0}}{K_{\mathrm{sea}}}}.`$ (58)
As seen in Fig. 26 and Table 12, this fit describes our data well and gives a reasonable value of $`\chi ^2/\mathrm{dof}1.5`$.

### V.2 chiral extrapolation in lattice units (method-B)

In order to study the convergence properties of polynomial fit forms in the method-B, we carry out quadratic and cubic diagonal fits to PS and vector meson masses

$`m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}})^2`$ $`=`$ $`B_{\mathrm{diag}}^{\mathrm{PS}}m_{\mathrm{q},\mathrm{sea}}+C_{\mathrm{diag}}^{\mathrm{PS}}m_{\mathrm{q},\mathrm{sea}}^2+D_{\mathrm{diag}}^{\mathrm{PS}}m_{\mathrm{q},\mathrm{sea}}^3,`$ (59)
$`m_\mathrm{V}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}})`$ $`=`$ $`A_{\mathrm{diag}}^\mathrm{V}+B_{\mathrm{diag}}^\mathrm{V}m_{\mathrm{PS},\mathrm{sea}}^2+C_{\mathrm{diag}}^\mathrm{V}m_{\mathrm{PS},\mathrm{sea}}^4+D_{\mathrm{diag}}^\mathrm{V}m_{\mathrm{PS},\mathrm{sea}}^6.`$ (60)
Fit parameters are summarized in Tables 13 and 14.

Fit curves of the quadratic and cubic fits to vector meson masses are shown in Fig. 27. While quadratic and cubic fits describe our data reasonably well at quark masses used in the simulation, they develop a deviation toward the chiral limit and for heavy quarks.

Figure 28 compares the relative magnitude of the linear and quadratic terms in the quadratic diagonal fit to PS meson mass in the method-A and B. As the quark mass increases, the magnitude of the quadratic contribution in the method-B increases more rapidly than in the method-A; it is no longer a small correction at the simulated quark masses ($`m_{\mathrm{q},\mathrm{diag}}0.015`$ – 0.055). A similar situation is observed in the chiral fit to vector meson masses as shown in Fig. 29.

We come to conclude that the chiral expansions of meson masses in lattice units, Eqs. (59) and (60), have poor convergence properties compared to those in unit of $`r_0`$, Eqs. (34) and (35) in the method-A, and the cubic term should not be ignored in the method-B. We directly confirm this point in Fig. 30, where the fit results for vector mesons from the method-A are converted to lattice units using Eq. (58) and compared with the fits of method-B. The cubic fit in the method-B shows a good consistency with the quadratic fit in the method-A, while the quadratic fit in method-B does not.

The combined chiral fit including cubic terms is not very stable because it contains a number of free parameters. In this study, therefore, we do not use the method-B to extract the physical hadron spectrum.

Before we turn to details of the determination of the hadron spectrum, let us make additional comments on the failure of the method-B with our data. Figure 31 shows the chiral fit of $`r_0`$ as a function of $`m_{\mathrm{PS},\mathrm{sea}}^2`$

$`r_0(K_{\mathrm{sea}})=A_{r_0}^{}+B_{r_0}^{}m_{\mathrm{PS},\mathrm{sea}}^2+C_{r_0}^{}m_{\mathrm{PS},\mathrm{sea}}^4+D_{r_0}^{}m_{\mathrm{PS},\mathrm{sea}}^6.`$ (61)
Fit parameters in Table 15 show that a large contribution of higher order terms is present also in this fit. By substituting this parametrization of $`r_0`$ to Eqs. (34) and (35) (diagonal fits in the method-A), the large higher order corrections appear in Eqs. (59) and (60) (diagonal fits in the method-B). Conversely, why the method-A works well is that large contributions of higher order terms in hadron masses and $`r_0`$ cancel with each other at least partially.

We note that the method-B works well in the CP-PACS’s study Spectrum.Nf2.CP-PACS , where they took similar simulation parameters but with different lattice actions, namely the RG-improved gauge action and the tadpole improved clover quark action. We compare the CP-PACS data of $`r_0m_\mathrm{V}`$ at $`\beta =2.2`$ and ours in Fig. 32. A good consistency in the whole range of the quark mass suggests that two groups’ data are in the scaling region. However, we find that the CP-PACS data of $`r_0`$ in lattice units show much milder dependence on the sea quark mass than ours. This is the reason why the method-B works well in the CP-PACS study, but does not with our data. It is, of course, not surprising that different lattice actions lead to different sea quark mass dependences of hadron masses and $`r_0`$ in the lattice units. However, as discussed above, the much stronger dependence with our choice of the lattice action is practically problematic, if one carries out the chiral extrapolation in lattice units.

### V.3 calculation of hadron spectrum

The hadron spectrum at the physical quark mass is determined as follows. The pion and $`\rho `$ meson masses normalized by $`r_0`$ are determined by tuning their ratio $`(r_0m_\pi )/(r_0m_\rho )`$ to its experimental value, i.e., by solving the equation,

$`{\displaystyle \frac{(r_0m_\pi )}{A^\mathrm{V}+\left(B_\mathrm{s}^\mathrm{V}+B_\mathrm{v}^\mathrm{V}\right)(r_0m_\pi )^2+C_{\mathrm{sv}}^\mathrm{V}(r_0m_\pi )^4}}={\displaystyle \frac{m_{\pi ,\mathrm{exp}}}{m_{\rho ,\mathrm{exp}}}},`$ (62)
where we denote the experimental value of hadron mass by $`m_{\mathrm{had},\mathrm{exp}}`$, and $`r_0`$ represents the Sommer scale at the physical sea quark mass. The hopping parameter corresponding to the physical light quark mass, $`K_{ud}`$, is fixed by solving

$`\left\{\text{ r.h.s of Eq. (}\text{40}\text{) with }K_{\mathrm{sea}}=K_{\mathrm{val},1}=K_{\mathrm{val},2}=K_{ud}\right\}`$ $`=`$ $`(r_0m_\pi )^2.`$ (63)
Then we determine $`r_0`$, which is required to convert $`r_0m_\pi `$ and $`r_0m_\rho `$ to $`m_\pi `$ and $`m_\rho `$, from $`K_{ud}`$ and Eq. (58).

We test two meson mass inputs to fix the strange quark mass:

* In the first method, we use the kaon mass as input. The hopping parameter corresponding to the strange quark mass, $`K_s`$, is determined by solving

$`{\displaystyle \frac{\sqrt{\text{ r.h.s of Eq. (}\text{40}\text{)}}}{r_0m_\rho }}={\displaystyle \frac{m_{K,\mathrm{exp}}}{m_{\rho ,\mathrm{exp}}}},`$ (64)
where we set $`K_{\mathrm{sea}}=K_{\mathrm{val},1}=K_{ud}`$ and $`K_{\mathrm{val},2}=K_s`$ in the r.h.s of Eq. (40). Then the mass of the “$`\eta _s`$” meson, that is an unphysical $`\overline{s}s`$ PS meson, is determined from Eq. (40), and used to calculate strange vector meson and strange baryon masses. We refer to this meson mass input as $`K`$-input.
* In the second method, we use the $`\varphi `$ meson mass as input assuming that it is a pure $`\overline{s}s`$ vector meson. The $`\eta _s`$ meson mass is fixed from

$`{\displaystyle \frac{\text{ r.h.s of Eq. (}\text{45}\text{)}}{r_0m_\rho }}={\displaystyle \frac{m_{\varphi ,\mathrm{exp}}}{m_{\rho ,\mathrm{exp}}}},`$ (65)
where we set $`m_{\mathrm{PS},\mathrm{sea}}=m_\pi `$ and $`m_{\mathrm{PS},\mathrm{val},1}=m_{\mathrm{PS},\mathrm{val},2}=m_{\eta _s}`$. We determine $`K_s`$ from $`m_{\eta _s}`$ and Eq. (40). This input is called $`\varphi `$-input.

The full spectrum of non-strange and strange hadrons is determined by substituting $`K_{ud}`$, $`m_\pi `$, $`K_s`$, $`m_{\eta _s}`$ and $`r_0`$ to Eqs. (40), (45), (52), (52) and (53). We use the lattice spacing $`a`$ determined from $`m_\rho `$ to convert the hadron masses in lattice units to those in physical units. We note that this estimate of the scale $`a`$ is subject to a systematic uncertainty due to the use of the polynomial fitting forms for the chiral extrapolation. However, if we use $`r_0`$ as the input to set the scale, we obtain a consistent result for $`a`$ within errors. The results of $`K_{ud}`$, $`K_s`$ and $`a^1`$ are collected in Table 16.

We repeat the above analysis using each of the alternative fit forms Eqs. (44), (49) (V.1) and (57). The largest deviation in the hadron spectrum among these analyses is taken as the systematic error due to the choice of the chiral fit forms.

For the chiral extrapolation of $`r_0`$, we find that an alternative form

$`r_0(K_{\mathrm{sea}})=A_{r_0}^{}+B_{r_0}^{}\mu _{\mathrm{PS},\mathrm{sea}}^2`$ (66)
also describes our data well. However, the hadron spectrum calculated using this fit is completely consistent with those using Eq. (58). We therefore ignore the systematic error due to the choice of the fit form Eq. (58).

The systematic error of the measured value of $`r_0(K_{\mathrm{sea}})`$ leads to an additional uncertainty in the result of the hadron spectrum. We perform the calculation of the spectrum with $`r_0(K_{\mathrm{sea}})`$ shifted by its systematic error at one value of $`K_{\mathrm{sea}}`$. This calculation is repeated for all $`K_{\mathrm{sea}}`$ and the largest deviation in the spectrum is included into the systematic error.

### V.4 chiral extrapolation in quenched QCD

The chiral extrapolations in quenched QCD are performed using fit forms which are obtained from those used in the full QCD analysis by dropping all terms describing the sea quark mass dependence. Namely, fitting forms for meson masses are

$`m_{\mathrm{PS}}(K_{\mathrm{val},1},K_{\mathrm{val},2})^2`$ $`=`$ $`B_\mathrm{q}^{\mathrm{PS}}m_{\mathrm{q},\mathrm{val}}+C_\mathrm{q}^{\mathrm{PS}}m_{\mathrm{q},\mathrm{val}}^2,`$ (67)
$`m_\mathrm{V}(K_{\mathrm{val},1},K_{\mathrm{val},2})`$ $`=`$ $`A_\mathrm{q}^\mathrm{V}+B_\mathrm{q}^\mathrm{V}m_{\mathrm{PS},\mathrm{val}}^2.`$ (68)
The following fitting forms are used for baryon masses

$`m_{\mathrm{oct},\mathrm{\Sigma }}(K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$ $`=`$ $`A_\mathrm{q}^\mathrm{O}+(F_\mathrm{q}^\mathrm{O}D_\mathrm{q}^\mathrm{O})m_{\mathrm{PS},\mathrm{val},1}^2+2F_\mathrm{q}^\mathrm{O}m_{\mathrm{PS},\mathrm{val},2}^2+(C_\mathrm{q}^\mathrm{O}+C_\mathrm{q}^{\mathrm{O},\mathrm{\Sigma }})m_{\mathrm{PS},\mathrm{val},1}^4`$ (69)
$`+(C_\mathrm{q}^\mathrm{O}C_\mathrm{q}^{\mathrm{O},\mathrm{\Sigma }})m_{\mathrm{PS},\mathrm{val},2}^4,`$
$`m_{\mathrm{oct},\mathrm{\Lambda }}(K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$ $`=`$ $`A_\mathrm{q}^\mathrm{O}+\left(F_\mathrm{q}^\mathrm{O}+{\displaystyle \frac{D_\mathrm{q}^\mathrm{O}}{3}}\right)m_{\mathrm{PS},\mathrm{val},1}^2+2\left(F_\mathrm{q}^\mathrm{O}{\displaystyle \frac{2}{3}}D_\mathrm{q}^\mathrm{O}\right)m_{\mathrm{PS},\mathrm{val},2}^2`$ (70)
$`+(C_\mathrm{q}^\mathrm{O}+C_\mathrm{q}^{\mathrm{O},\mathrm{\Lambda }})m_{\mathrm{PS},\mathrm{val},1}^4+(C_\mathrm{q}^\mathrm{O}C_\mathrm{q}^{\mathrm{O},\mathrm{\Lambda }})m_{\mathrm{PS},\mathrm{val},2}^4,`$
$`m_{\mathrm{dec}}(K_{\mathrm{val},1},K_{\mathrm{val},2},K_{\mathrm{val},2})`$ $`=`$ $`A_\mathrm{q}^\mathrm{D}+B_\mathrm{q}^\mathrm{D}m_{\mathrm{PS},\mathrm{val}}^2.`$ (71)
These forms fit to our data very well as shown in Fig. 33. Fit parameters are summarized in Table 17. The hadron spectrum is calculated in an analogous way to that for full QCD. Resulting values of $`K_{ud}`$, $`K_s`$ and $`a^1`$ are summarized in Table 18.

## VI Sea Quark Effects on Hadron Spectrum

### VI.1 Sea quark effects at simulated quark mass

Figure 34 compares the valence quark mass dependence of the vector meson mass at each sea quark mass in full QCD and in quenched QCD. We observe that quenched data have a significantly smaller slope than experimental points. This leads to the underestimation of the $`K^{}`$$`K`$ hyperfine splitting in quenched QCD.

The slopes in full QCD data are clearly larger than in quenched QCD, and increase for decreasing sea quark mass. This is reflected in a negative value of $`C_{\mathrm{sv}}^\mathrm{V}`$ for Eq. (45) in Table 8. This sea quark effect leads to a better agreement of the meson spectrum in full QCD with experiment than in quenched QCD.

The $`J`$ parameter J  defined by

$`J=m_\mathrm{V}{\displaystyle \frac{dm_\mathrm{V}}{dm_{\mathrm{PS}}^2}}|_{m_{\mathrm{PS}}/m_\mathrm{V}=m_K/m_K^{}}`$ (72)
is useful to quantify the sea quark effect. Numerical results of $`J`$ calculated from the partially quenched chiral fit Eq. (50) are given in Table 19. In Fig. 35, we plot $`J`$ in full QCD (filled circles) as a function of the sea quark mass together with the quenched result (open circle in the right panel). We observe that $`J`$ in full QCD is close to the quenched value at heavy sea quark masses corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.75`$ and increases as the sea quark masses decreases.

In the same figure, we also plot $`J`$ reproduced from the combined chiral fit Eq. (45) (dashed lines). The result is consistent with $`J`$ from the partially quenched fit, as it should be, and shows a similar sea quark mass dependence. We observe that $`J`$ extrapolated to the physical sea quark mass is closer to the phenomenological value J

$`J=m_\mathrm{K}^{}{\displaystyle \frac{m_K^{}m_\rho }{m_K^2m_\pi ^2}}=0.48(2)`$ (73)
than in quenched QCD.

Figure 35 also shows $`J`$ calculated from our results of the meson spectrum (see Table 20) using the above alternative definition Eq. (73) (filled square). The result is in good agreement with other determinations, showing the magnitude of the sea quark effect in $`J`$ to be stable against the definition of $`J`$.

In Sec. IV, we pointed out that FSE decreases the slope $`dm_\mathrm{V}/dm_{\mathrm{PS}}^2`$. This is confirmed numerically in $`dm_\mathrm{V}/dm_{\mathrm{PS}}^2`$ determined from the partially quenched chiral fit Eq. (50) at the lightest sea quark mass: $`dm_\mathrm{V}/dm_{\mathrm{PS}}^2`$ is 0.906(14) on $`20^3`$, 0.814(88) on $`16^3`$ and 0.68(30) on $`12^3`$. The slope would be larger if we increase the spatial size beyond $`20^3`$. Therefore, the observed sea quark mass dependence of the slope and $`J`$ is a genuine effect of dynamical quarks and not an artifact of FSE.

A similar effect of sea quarks can be found in decuplet baryon masses as shown in Fig. 36. However, a significant deviation still exists in the slope between full QCD data and the experimental spectrum. We consider that a larger slope in full QCD is still partly masked by FSE on the $`20^3`$ volume; volume as large as $`30^3`$ would be needed to reduce FSE to a few % level as discussed in Sec. IV.

We emphasize that the evidence of sea quark effects observed in this subsection does not suffer from possibly large systematic errors due to the chiral extrapolation and the choice of inputs to fix the scale and quark masses: the increase of the slopes, $`dm_\mathrm{V}/dm_{\mathrm{PS}}^2`$ and $`dm_{\mathrm{dec}}/m_{\mathrm{PS}}^2`$, is observed without any chiral extrapolation and inputs. The sea quark mass dependence of $`J`$ is obtained by a short extrapolation or an interpolation to a relatively heavy valence quark mass corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}=m_K/m_K^{}`$.

### VI.2 hadron spectrum

The meson spectrum in full and quenched QCD is summarized in Table 20. Since our fitting functions to vector meson masses, Eqs. (45) and (68), are linear in terms of the valence quark mass, $`m_K^{}`$ with $`\varphi `$-input equals $`\left(m_{\rho ,\mathrm{exp}}+m_{\varphi ,\mathrm{exp}}\right)/2`$ in both full and quenched QCD. The deviation of this value from the experimental mass $`m_{K^{},\mathrm{exp}}`$ is only 0.2 %.

A clear difference between full and quenched QCD is observed in other meson masses as shown in Fig. 37. While the quenched meson spectrum shows a significant deviation from experiment, sea quark effects reduce the deviation by about 40%. This closer agreement of the meson spectrum with experiment is a consequence of the sea quark effects observed in the previous subsection.

In Fig. 38, $`m_K^{}`$ with the $`K`$-input in full QCD is compared with the CP-PACS results obtained with the RG-improved gauge and clover quark actions Spectrum.Nf2.CP-PACS . We observe that our $`m_K^{}`$ is consistent with the CP-PACS’ result at a similar lattice spacing, and is at the lower edge of their estimate in the continuum limit. In the same figure, we also make a comparison in quenched QCD with the CP-PACS results obtained with the plaquette gauge and the Wilson quark actions in Ref. Spectrum.Nf0.CP-PACS . Two groups’ results show a good agreement with each other. These observations suggest that the scaling violation is small in our data both in full and quenched QCD and, hence, the closer agreement of the meson spectrum in full QCD with experiment is a genuine effect of sea quarks.

The baryon masses in full and quenched QCD are listed in Table 21. These masses are compared with experiment in Fig. 39. For heavier baryons, such as $`\mathrm{\Sigma }`$, $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$, full QCD results show a closer agreement with experiment than in quenched QCD. The sea quark effect is, however, less clear for lighter baryons. This is partly due to FSE in full QCD data which is more pronounced for lighter valence quarks.

In Fig. 40, $`m_N`$ and $`m_\mathrm{\Xi }`$ with the $`K`$-input are compared with the CP-PACS results Spectrum.Nf0.CP-PACS ; Spectrum.Nf2.CP-PACS . While the full QCD results of two groups show a reasonable agreement with each other, the CP-PACS results in the continuum limit in quenched QCD are systematically smaller than ours. This suggests that our quenched data has non-negligible scaling violation, which is another source making sea quark effects less clear. Therefore, further investigations of FSE in full QCD and scaling violation in quenched QCD are required to obtain a clear conclusion on sea quark effects in the baryon spectrum.

We now turn to theoretical predictions which can be derived from our data. The first is the mass of the $`\eta _s`$ meson, for which our full QCD data predict $`m_{\eta _s}=0.6948(3)(+8/1)(+2)`$ GeV with the $`K`$-input, and $`0.7381(46)(+57)(+40/46)`$ GeV with the $`\varphi `$-input, where the first error is statistical, and the second and third ones are due to the choice of the fitting form and a systematic uncertainty of $`r_0`$. These results are to be compared with those in quenched QCD, $`0.6988(9)`$ GeV ($`K`$-input) and $`0.7719(58)`$ GeV ($`\varphi `$-input). While the values themselves do not differ by going from quenched to full QCD, the difference between the two inputs is reduced by about 40% in full QCD. This reflects the closer agreement of the meson spectrum in full QCD.

Another interesting prediction is the physical value of $`r_0`$. Our full QCD simulation gives $`r_0=0.497(6)(9)(+11/12)`$ fm, where the meaning of three errors are the same to those of $`m_{\eta _s}`$. We note that this is close to the phenomenological estimate in the original paper r0 , $`r_0=0.49`$ fm. The quenched simulation gives 0.5702(75)(50) fm, where the first error is statistical and the second comes from the systematic uncertainty of the measurement. About a $`14`$ % difference between full and quenched QCD arises from the following two sea quark effects. One is the difference of the physical value of $`r_0`$ itself due to the change of the shape of the static quark potential. The other source is the reduction of the quenching error in $`m_\rho `$ in full QCD, which is used to fix the lattice scale.

## VII Decay constants

The PS meson decay constants are calculated using the fourth component of the improved axial vector current, which is defined by

$`A_4^{\mathrm{imp}}`$ $`=`$ $`A_4+c_A\mathrm{\Delta }_4P`$ (74)
with the symmetric lattice derivative $`\mathrm{\Delta }_4`$. We extract the amplitude of the $`A_4^{\mathrm{imp}}(t)P(0)^{}`$ correlator, $`C_A^{\mathrm{LS}}`$, by the correlated fit of the form

$`A_4^{\mathrm{imp}}(t)P(0)^{}`$ $`=`$ $`C_A^{\mathrm{LS}}\left\{\mathrm{exp}\left[m_{\mathrm{PS}}t\right]\mathrm{exp}\left[m_{\mathrm{PS}}(N_tt)\right]\right\}`$ (75)
with $`m_{\mathrm{PS}}`$ fixed to the results given in Appendix A. We use the local operator for $`A_4^{\mathrm{imp}}`$, while the double smearing is applied to $`P`$. The amplitude of the $`P(t)P(0)^{}`$ operator with the doubly smeared source and the sink operators, $`C_P^{\mathrm{SS}}`$, is extracted assuming a single hyperbolic cosine form.

The renormalized decay constant is calculated by

$`f_{\mathrm{PS}}`$ $`=`$ $`2KZ_A\left(1+b_Am_{\mathrm{q},\mathrm{val}}\right)C_A^{\mathrm{LS}}\sqrt{{\displaystyle \frac{2}{m_{\mathrm{PS}}C_P^{SS}}}}.`$ (76)
Since non-perturbatively determined values for $`Z_A`$, $`b_A`$ and $`c_A`$ are not available for two-flavor QCD, we adopt one-loop perturbative values in Refs. 1loop\_cSWcA.ALPHA ; 1loop\_Z.Roma ; 1loop\_b.ALPHA ; 1loop\_Zb.Tsukuba ; 1loop\_Zbc.Tsukuba ; 1loop\_bc.Aoki  with the tadpole improvement. We calculate $`\alpha _P(3.40/a)`$ from the plaquette average $`U_P`$ according to alpha\_P.JLQCD ; alpha\_P

$`\mathrm{ln}\left[U_P\right]`$ $`=`$ $`{\displaystyle \frac{4\pi }{3}}\alpha _P(3.40/a)\left\{1\left(1.19050.2266N_f\right)\alpha _P\right\}.`$ (77)
Then, $`\alpha _P(3.40/a)`$ is evolved to the optimum scale ($`q_{Z_A}^{}=1.803/a`$ for $`Z_A`$, $`q_{b_A}^{}=2.289/a`$ for $`b_A`$ and $`q_{c_A}^{}=2.653/a`$ for $`c_A`$ q\_star ) using the universal two-loop beta function and is used as the expansion parameter of tadpole improved perturbation theory.

The consistent chiral extrapolation of the decay constant should include the chiral logarithmic term as predicted by ChPT ChPT . However, our data do not show the characteristic curvature of the chiral logarithm as discussed in Sec.V.1 (and also in Refs. Spectrum.Nf2.JLQCD.lat01 ; Spectrum.Nf2.JLQCD.lat02.Hashimoto ). We therefore use the following polynomial form for the chiral extrapolation leaving the problem of the chiral logarithm and associated uncertainty for future publication ChPT\_test.Nf2.JLQCD

$`r_0(K_{\mathrm{sea}})f_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})=A^\mathrm{f}+B_\mathrm{s}^\mathrm{f}\mu _{\mathrm{PS},\mathrm{sea}}+B_\mathrm{v}^\mathrm{f}\mu _{\mathrm{PS},\mathrm{val}}+C_{\mathrm{sv}}^\mathrm{f}\mu _{\mathrm{PS},\mathrm{sea}}\mu _{\mathrm{PS},\mathrm{val}}.`$ (78)
The fit is plotted in Fig. 41 with parameters summarized in Table 22. Pion and kaon decay constants, $`f_\pi `$ and $`f_K`$, are calculated by tuning $`\mu _{\mathrm{PS},\mathrm{sea}}`$ to $`(r_0m_\pi )^2`$ and substituting $`(r_0m_\pi )^2`$ or $`(r_0m_{\eta _s})^2`$ for $`\mu _{\mathrm{PS},\mathrm{val},i}`$ $`(i=1,2)`$ in $`\mu _{\mathrm{PS},\mathrm{val}}`$. The systematic errors due to the choice of the fitting function and the uncertainty of $`r_0`$ are estimated in a way similar to those described in Sec. V.3. In the estimation of the former error, we use

$`r_0(K_{\mathrm{sea}})f_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})`$ $`=`$ $`\text{“r.h.s of Eq. (}\text{78}\text{)”}+C_\mathrm{s}^\mathrm{f}\mu _{\mathrm{PS},\mathrm{sea}}^2+C_\mathrm{v}^\mathrm{f}\mu _{\mathrm{PS},\mathrm{val}}^2`$ (79)
and Eqs. (44) and (49) as the alternative fitting functions for the chiral extrapolation.

In quenched QCD, we use the chiral extrapolation form

$`f_{\mathrm{PS}}(K_{\mathrm{val},1},K_{\mathrm{val},2})`$ $`=`$ $`A_\mathrm{q}^\mathrm{f}+B_\mathrm{q}^\mathrm{f}m_{\mathrm{PS},\mathrm{val}}^2.`$ (80)
and obtain parameters summarized in Table 23. For $`Z_A`$, $`b_A`$ and $`c_A`$, we test the one-loop perturbative value and the non-perturbative one in Ref. NP\_Zbc.LANL .

Our results of the decay constants are summarized in Table 24. A comparison between full and quenched QCD results obtained by the one-loop matching is made in Fig. 42. We observe that $`f_\pi `$ is consistent with the experimental value within two standard deviations in both full and quenched QCD. While $`f_K`$ in quenched QCD is significantly smaller than the experimental value, the deviation is reduced by sea quark effects and the full QCD result becomes consistent with experiment.

The results obtained with one-loop renormalization factors are subject to higher order corrections. However, as shown in Table 24 for the quenched results, the difference between the perturbative and non-perturbative matchings is not large. This is because the $`O(a)`$ correction to the improved current in Eq. (74) is not large, and the non-perturbative values for $`Z_A`$ and $`b_A`$ are close to those in tadpole improved perturbation theory. We may therefore expect that the uncertainty due to the perturbative matching is small also in the full QCD results.

It is expected that various systematic uncertainties, including the scaling violation, would partially cancel in the ratio $`f_K/f_\pi `$. This expectation is supported by the good agreement of the quenched results between the perturbative and non-perturbative matching, which suggests that higher order corrections to the renormalization factors almost cancel in the ratio. Therefore the ratio is useful to discuss sea quark effects. Comparison of this quantity shows that the full QCD values are significantly closer to the experimental value $`1.22`$ by about two standard deviations than in quenched QCD.

## VIII Quark masses

We calculate the up-down and strange quark masses through the axial vector Ward identity (AWI). The bare quark mass at simulation points is obtained by

$`m_q^{\mathrm{AWI}}`$ $`=`$ $`{\displaystyle \frac{m_{\mathrm{PS}}C_A^{\mathrm{LS}}}{2C_P^{\mathrm{LS}}}},`$ (81)
where $`C_A^{\mathrm{LS}}`$ and $`C_P^{\mathrm{LS}}`$ are the amplitudes of $`A_4^{\mathrm{imp}}(t)P(0)^{}`$ and $`P(t)P(0)^{}`$ with the doubly smeared source and the local sink operators.

We then carry out the chiral fit of the PS meson mass as a function of the AWI bare quark mass. The fitting function is obtained from Eq. (40) with the replacement of the VWI masses with the AWI ones. We also drop all monomial terms in the sea quark mass, since the PS meson mass vanishes in the chiral limit $`m_{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}=0`$ even for non-zero sea quark masses PQChPT . The adopted form is

$`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})\right)^2`$ (82)
$`=`$ $`B_\mathrm{v}^{\mathrm{PS},\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}+C_\mathrm{v}^{\mathrm{PS},\mathrm{AWI}}(\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}})^2+C_{\mathrm{sv}}^{\mathrm{PS},\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{sea}}^{\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}},`$
where

$`\mu _{\mathrm{q},\mathrm{sea}}^{\mathrm{AWI}}`$ $`=`$ $`r_0(K_{\mathrm{sea}})m_\mathrm{q}^{\mathrm{AWI}}(K_{\mathrm{sea}};K_{\mathrm{sea}},K_{\mathrm{sea}}),`$ (83)
$`\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}`$ $`=`$ $`r_0(K_{\mathrm{sea}})m_{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}},`$ (84)
$`m_{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(m_{\mathrm{q},\mathrm{val},1}^{\mathrm{AWI}}+m_{\mathrm{q},\mathrm{val},2}^{\mathrm{AWI}}\right),`$ (85)
$`m_{\mathrm{q},\mathrm{val},i}^{\mathrm{AWI}}`$ $`=`$ $`m_\mathrm{q}^{\mathrm{AWI}}(K_{\mathrm{sea}};K_{\mathrm{val},i},K_{\mathrm{val},i}).`$ (86)
Our data and fit are shown in Fig. 43.

We adopt the fit (82) because it is consistent with that in terms of the VWI quark mass (Eq. (40)). However, a function with fewer terms

$`\left(r_0(K_{\mathrm{sea}})m_{\mathrm{PS}}(K_{\mathrm{sea}};K_{\mathrm{val},1},K_{\mathrm{val},2})\right)^2=B_\mathrm{v}^{\mathrm{PS},\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}+C_{\mathrm{sv}}^{\mathrm{PS},\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{sea}}^{\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}`$ (87)
also gives an acceptable $`\chi ^2/\mathrm{dof}`$. We use this as an alternative fit in our estimation of the systematic error due to the choice of the fitting function (see below). Parameters of these two fits are summarized in Table 25.

The bare AWI masses of the up-down and strange quarks are fixed in a way analogous to that described in Sec. V.3 by using Eq. (82) instead of Eq. (40). The matching to the $`\overline{\mathrm{MS}}`$ scheme is made at the scale $`\mu =2/a`$ using the one-loop renormalization constant 1loop\_Z.Roma ; 1loop\_bc.Aoki ; 1loop\_Zbc.Tsukuba ; 1loop\_b.ALPHA ; 1loop\_Zb.Tsukuba  with the tadpole improvement. We use $`\alpha _P(q_{Z_A}^{})`$, $`\alpha _P(q_{b_A}^{})`$ as the expansion parameter in the one-loop expression of $`Z_A`$ and $`b_A`$, while we set $`q^{}=2/a`$ for other coefficients for which $`q^{}`$ is not known. The $`\overline{\mathrm{MS}}`$ quark mass is evolved to $`\mu =2`$ GeV using the four-loop beta function 4loop.running.1 ; 4loop.running.2 .

The VWI quark mass may differ from the AWI one because of explicit violation of chiral symmetry at finite lattice spacings. The difference between the AWI and VWI masses, therefore, gives insight into the size of scaling violation in our results. This leads us to repeat the calculation of quark masses using the VWI definition. The bare quark mass is calculated from $`K_c`$, $`K_{ud}`$ and $`K_s`$ in Tables 7 and 16. The $`\overline{\mathrm{MS}}`$ mass is obtained by the one-loop matching at $`\mu =2/a`$ and the four-loop running to $`\mu =2`$ GeV. The resulting AWI and VWI masses are summarized in Table 26.

For a calculation of quark masses in quenched QCD, the chiral extrapolation is carried out using a quadratic fit form obtained from Eq. (82) by dropping the third term which represents the sea quark mass dependence

$`m_{\mathrm{PS}}(K_{\mathrm{val},1},K_{\mathrm{val},2})^2=B_\mathrm{q}^{\mathrm{PS},\mathrm{AWI}}\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}}+C_\mathrm{q}^{\mathrm{PS},\mathrm{AWI}}(\mu _{\mathrm{q},\mathrm{val}}^{\mathrm{AWI}})^2.`$ (88)
Obtained parameters are listed in Table 27. We use either one-loop or non-perturbative value in Refs. NPcsw.Nf0.ALPHA ; NP\_ZA.ALPHA ; NP\_ZP.ALPHA ; NP\_bA-bP.ALPHA  for the renormalization factors. Numerical results are summarized in Table 28.

In Figure 44, our full QCD results are compared with estimates by the CP-PACS Spectrum.Nf2.CP-PACS ; mq.Nf2.CP-PACS  and the QCDSF-UKQCD Collaborations mq.Nf2.QCDSF . We observe good agreement for the AWI masses among the three groups even at the finite lattice spacing of $`a^12`$ GeV. These results are consistent also with the CP-PACS’s result extrapolated to the continuum limit. This suggests that various uncertainties, such as scaling violation and higher order corrections to renormalization factors, are likely to be canceled in the ratio defining the AWI mass (Eq. (81)).

On the other hand, such a cancellation is not expected in the VWI mass. Indeed there is a sizable difference between our AWI and VWI results. We also observe that our and the CP-PACS results of VWI $`m_{ud}`$ show a large discrepancy of about 18 % (6 standard deviations). These observations suggest that the scaling violation in our results of the VWI masses is not small.

In both full and quenched QCD, therefore, we quote the AWI masses as the central value. We adopt $`K`$-input for $`m_s`$. The difference between $`m_s`$ with $`K`$\- and $`\varphi `$-inputs is treated as the systematic error due to the choice of the meson mass input to fix $`m_s`$.

Additional systematic errors due to the choice of the chiral fit form and the uncertainty of the measured value of $`r_0`$ are included in our final results in full QCD. These errors are estimated in a similar way to that described in Sec. V.3 by using Eqs. (44), (49) and (87) as alternative fit forms.

Adding all errors in quadrature, we obtain

$`m_{ud}^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`3.223({}_{0.069}{}^{+0.046})`$ (89)
$`m_s^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`84.5({}_{1.7}{}^{+12.0})`$ (90)
$`m_s/m_{ud}`$ $`=`$ $`26.13({}_{0.02}{}^{+3.65})`$ (91)
in two-flavor QCD and

$`m_{ud}^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`4.020\left(0.077\right)`$ (92)
$`m_s^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`104.1({}_{1.6}{}^{+24.1})`$ (93)
$`m_s/m_{ud}`$ $`=`$ $`25.90({}_{0.13}{}^{+5.98}).`$ (94)
in quenched QCD using the one-loop matching. The scaling violation is expected to be small and hence is ignored here. This point, however, should be checked in future studies.

There is an additional uncertainty arising from the use of the perturbative value for $`c_A`$, $`Z_{A,P}`$, and $`b_{A,P}`$. Comparing the quenched result of Eqs. (92) – (94) with those obtained by the non-perturbative matching given by

$`m_{ud}^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`3.522\left(0.66\right)\text{MeV}`$ (95)
$`m_s^{\overline{\mathrm{MS}}}(\text{2 GeV})`$ $`=`$ $`91.9({}_{1.4}{}^{+21.3})\text{MeV}`$ (96)
$`m_s/m_{ud}`$ $`=`$ $`26.08({}_{0.13}{}^{+6.02}),`$ (97)
we observe a systematic error of about 13 %. In quenched QCD the non-perturbative estimate of $`Z_A/Z_P=1.19`$ is very close to that in one-loop perturbation theory 1.22, since higher order corrections in $`Z_A`$ and $`Z_P`$ partially cancel with each other. The non-perturbative value $`b_Ab_P=0.171`$ deviates significantly from that at one-loop $`0.011`$. The $`O(am_q)`$ term, however, is a small correction in our data. Most of the 13 % difference originates from the large deviation between the non-perturbative value $`c_A=0.083`$ in Ref. NPcsw.Nf0.ALPHA  and its one-loop value $`0.013`$. Therefore, a non-perturbative determination of $`c_A`$ in full QCD is an important task toward a more precise calculation of the quark masses in future studies.

In Fig. 45, we compare the quark masses in full and quenched QCD. The chief observation is that sea quark effects reduce the light and strange quark masses by about 20 %. The magnitude of the sea quark effect is roughly consistent with the CP-PACS observation in Refs. Spectrum.Nf2.CP-PACS ; mq.Nf2.CP-PACS .

The quark mass ratio $`m_s/m_{ud}`$ in full QCD is consistent with the quenched value, because the sea quark effects in $`m_{ud}`$ and $`m_s`$ almost cancel with each other in the ratio. We note that $`m_s/m_{ud}26`$ is in a good agreement with the estimate of one-loop chiral perturbation theory 24.4(1.5) mq.ChPT .

Another important observation is that the deviation in $`m_s`$ between $`K`$\- and $`\varphi `$-inputs is reduced by effects of sea quarks: the deviation is 21 % ($`24`$ MeV) in quenched QCD and 13 % ($`12`$ MeV) in full QCD. This reflects the closer agreement of the meson spectrum in full QCD with experiment. The remaining deviation may be attributed to quenching of strange quarks and scaling violation.

## IX Conclusions

In this paper we have presented a high statistics study of the hadron spectrum and quark masses in two-flavor QCD using the plaquette gauge action and the fully $`O(a)`$-improved Wilson quark action. We find firm evidence of sea quark effects at the simulated quark masses: the slopes $`dm_\mathrm{V}/dm_{\mathrm{PS}}^2`$ and $`dm_{\mathrm{dec}}/dm_{\mathrm{PS}}^2`$ are larger than in quenched QCD and the $`J`$ parameter increases for lighter sea quarks. These findings do not suffer from systematic errors arising from the chiral extrapolation with respect to the sea quark mass, which is a major uncertainty particularly in recent studies with the Wilson-type quark action. Note that the use of a volume $`La1.8`$ fm at smaller sea quark masses $`m_{\mathrm{PS},\mathrm{sea}}/m_{\mathrm{V},\mathrm{sea}}0.6`$ – 0.7 is an important factor to control finite-size errors and in reaching our observations.

The sea quark effect observed at the simulated quark masses means that the strange meson and baryon masses in full QCD show a better agreement with experiment than in quenched QCD. A similar reduction of quenching errors is also observed in the ratio $`f_K/f_\pi `$. We also find that the sea quark effects lead to about 20 % reduction of quark masses.

For baryons finite-size effects are large for the volume we used, which render sea quark effects unclear for lighter baryons. Further investigations on larger spatial volumes, of the order of 3 fm at the lightest sea quark mass, are needed to observe sea quark effects in the light baryons.

The present work is carried out at a single lattice spacing of $`a^1=2.221(28)`$ GeV. The $`O(a)`$-improved Wilson quark action we employed is designed to have reduced scaling violation, and experiences in quenched QCD scaling.Nf0.UKQCD  support this expectation. Nonetheless scaling study of both the hadron spectrum and quark masses with this action is needed to establish the sea quark effects in the spectral quantities on a quantitative basis.

Another important subject in future is simulations at much lighter sea quark mass, in particular below the $`\rho \pi \pi `$ threshold. Such simulations will lead to a better control of the chiral extrapolation. This would also give insight into the chiral logarithmic singularity in the PS meson mass and decay constant, which we have not observed in our data.

###### Acknowledgements.

This work is supported by the Supercomputer Project No.79 (FY2002) of High Energy Accelerator Research Organization (KEK), and also in part by the Grant-in-Aid of the Ministry of Education (Nos. 11640294, 12640253, 12740133, 13135204, 13640259, 13640260, 14046202, 14740173). N.Y. is supported by the JSPS Research Fellowship.

## Appendix A Hadron masses

Measured hadron masses are summarized in Tables LABEL:tab:meas:had:meson\_12x48 – LABEL:tab:meas:had:baryon\_20x48 for full QCD, and Tables LABEL:tab:meas:had:meson\_q12x48 – LABEL:tab:meas:had:baryon\_q20x48 for quenched QCD. Our choice of the fitting range and resulting value of $`\chi ^2/\mathrm{dof}`$ are also shown in these tables.