# UTCCP-P-61 Comparative Study of full QCD Hadron Spectrum and Static Quark Potential with Improved Actions

## I Introduction

With the progress over the last few years of quenched simulations of QCD, it has become increasingly clear that the quenched hadron spectrum shows deviations from the experiment if examined at a precision better than 5–10%. For light hadrons the first indication was that the strange quark mass cannot be set consistently from pseudo scalar and vector meson channels in quenched QCD . For heavy quark systems calculations both with relativistic  and non-relativistic  quark actions have shown that the fine structure of quarkonium spectra can not be reproduced on quenched gluon configurations. Most recently an extensive calculation by the CP-PACS collaboration found a systematic departure of both the light meson and baryon spectra from experiment . These results raise the question as to whether the discrepancies can be accounted for by the inclusion of dynamical sea quarks. It is therefore timely to study more thoroughly the effects of full QCD in order to answer this question.

Full QCD simulations are, however, computationally much more expensive than those of quenched QCD. Simple scaling estimates coupled with past experience place a hundred-fold or more increase in the amount of computations for full QCD compared to that of quenched QCD with current algorithms. Since $`32^3\times 64`$ is a typical maximal lattice size for quenched QCD which can be simulated with high statistics on computers with a speed in the 10 GFLOPS range , reliable full QCD results are difficult to obtain on lattice sizes exceeding $`32^3\times 64`$ even with TFLOPS-class computers such as CP-PACS  and QCDSP . Recalling that a physical lattice size of $`L2.5`$–3.0 fm is needed to avoid finite-size effects , the smallest lattice spacing one can reasonably reach at present is therefore $`a^12`$ GeV. Hence lattice discretization errors have to be controlled through simulations carried out at inverse lattice spacings below this value, e.g. in the range $`a^112`$ GeV. It is, however, known that with the standard plaquette and Wilson quark actions discretization errors are already of order 10% even for $`a^12`$ GeV. These observations suggest the use of improved actions for simulations of full QCD.

Studies of improved actions have been widely pursued in the last few years. Detailed tests of improvement for the hadron spectrum, however, have been carried out mostly within quenched QCD  with only a few full QCD attempts . In particular, a systematic investigation of how gauge and quark action improvement, taken separately, affects light hadron observables has not been carried out in full QCD. Prior to embarking on a large scale simulation, we examine this question as the first subject of the full QCD program on the CP-PACS computer.

For a systematic comparison of action improvement we employ four possible types of action combinations, the standard plaquette or a renormalization-group improved action  for the gauge part, and the standard Wilson or the improvement of Sheikholeslami and Wohlert  for the quark part. Since effects of improvement are clearer to discern at coarser lattice spacings, we carry out simulations at an inverse lattice spacing of $`a^11`$ GeV with quark masses in the range corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.7`$–0.9. Results for the four action combinations are used for comparative tests of improvement on the light hadron spectrum and the static quark potential.

Another limiting factor for full QCD simulations is how close one can approach the chiral limit with present computing power. To investigate this question we take the action in which both gauge and quark parts are improved, and carry out simulations down to a quark mass corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$. In addition to exploring the chiral behavior of hadron masses, this simulation allows an examination of signs of string breaking in the static quark-antiquark potential.

In this article we present results of our study on the two questions discussed above, expounding on the preliminary accounts reported in Refs. . We begin with discussions on our choice of actions for our comparative studies in Sec. II. Details of the full QCD configuration generation procedure and measurements of hadron masses and potential are described in Sec. III. Results for the hadron masses are discussed in Sec. IV where, after a description of the chiral extrapolation or interpolation of our data, we examine the effects of action improvement for the scaling behavior of hadron mass ratios. In Sec. V we turn to discuss the static potential. The influence of action improvement on the restoration of rotational symmetry of the potential is examined, and the consistency of the lattice spacing determined from the vector meson mass and the string tension is discussed. In Sec. VI we report on our effort to approach the chiral limit, where our attempt to observe a flattening of the potential at large distances due to string breaking is also presented. We end with a brief conclusion in Sec. VII. Detailed numerical results on run performances, hadron masses and string tensions are collected at the end in Appendices A, B and C.

## II Choice of action

The discretization error of the standard plaquette gauge action is $`O(a^2)`$ while that of the Wilson quark action is $`O(a)`$. In principle one would only need to improve the quark action to the same order as the gauge action. On the other hand, violations of rotational invariance have been found to be strong for the plaquette gauge action at coarse lattice spacings . Hence improving the gauge action is still advantageous for coarse lattices. In this spirit we employ (besides the standard actions) improved actions both in the gauge and quark sectors in the forms specified below.

Let us denote the standard plaquette gauge action by P. Improving this action requires the addition of Wilson loops with a perimeter of six links or more. The number, the precise form and the coefficients of the added terms differ depending on the principle one follows for the improvement . In this study we test the action determined by an approximate block-spin renormalization group analysis of Wilson loops, denoted by R in the pursuant, which is given by 

$$S_g^𝐑=\frac{\beta }{6}\left(c_0W_{1\times 1}+c_1W_{1\times 2}\right),$$
(1)   
where the $`1\times 2`$ rectangular shaped Wilson loop $`W_{1\times 2}`$ has the coefficient $`c_1=0.331`$ and from the normalization condition defining the bare coupling $`\beta =6/g_0^2`$ follows $`c_0=18c_1=3.648`$.

The discretization error of the R action is still $`O(a^2)`$. The coefficients of $`O(a^2)`$ terms in physical quantities, however, are expected to be reduced from those of the plaquette action. Indeed, the quenched static quark potential calculated with this action was found to exhibit good rotational symmetry and scaling already at $`a^11`$ GeV , and so does the scaling of the ratio $`T_c/\sqrt{\sigma }`$ of the critical temperature of the pure gauge deconfining phase transition and the string tension $`\sigma `$ . The degree of improvement is similar to those observed for tadpole-improved and fixed point actions .

To improve the quark action we adopt the clover improvement proposed by Sheikholeslami and Wohlert , denoted by C in the following and defined by

$$D_{xy}^𝐂=D_{xy}^𝐖\delta _{xy}c_{\mathrm{SW}}K\underset{\mu <\nu }{}\sigma _{\mu \nu }F_{\mu \nu },$$
(2)   
where $`D_{xy}^𝐖`$ is the standard Wilson quark matrix given by

$$D_{xy}^𝐖=\delta _{xy}K\underset{\mu }{}\{(1\gamma _\mu )U_{x,\mu }\delta _{x+\widehat{\mu },y}+(1+\gamma _\mu )U_{x,\mu }^{}\delta _{x,y+\widehat{\mu }}\}$$
(3)   
and $`F_{\mu \nu }`$ is the lattice discretization of the field strength,

$$F_{\mu \nu }=\frac{1}{8i}(f_{\mu \nu }f_{\mu \nu }^{}),$$
(4)   
where $`f_{\mu \nu }`$ is the standard clover-shaped combination of gauge links.

The complete removal of $`O(a)`$ errors requires a non-perturbative tuning of the clover coefficient $`c_{\mathrm{SW}}`$. This has been carried out for the plaquette gauge action in both quenched  and two-flavor full QCD . A similar analysis for the R gauge action is yet to be made, however. In this study we compare three different choices:

* The tree level value $`c_{\mathrm{SW}}=1`$.
* The mean-field (MF) improved value  $`c_{\mathrm{SW}}=P^{3/4}`$ with $`P`$ the self-consistently determined plaquette average.
* A perturbative mean-field (pMF) improved value $`c_{\mathrm{SW}}=P^{3/4}`$ with the plaquette $`P`$ calculated in one-loop perturbation theory. For the R gauge action $`P=10.8412\beta ^1`$ .

For all three choices the leading discretization error in physical quantities is $`O(g_0^2a)`$. The magnitude of the coefficients of this term should be reduced in the cases of (b) and (c) as compared to (a). The one-loop value of $`c_{\mathrm{SW}}`$ has been recently reported to be $`c_{\mathrm{SW}}=1+0.678(18)/\beta `$ . This value is close to the pMF value $`c_{\mathrm{SW}}^{\mathrm{pMF}}=1+0.631/\beta +\mathrm{}`$. We also find that the one-loop value of $`P`$ reproduces the measured values from simulations within 10% for the R action. Hence the pMF value of the clover coefficient is similar to the MF value employed in (b). The advantage of the pMF choice is that it does not require a self-consistent tuning of $`c_{\mathrm{SW}}`$ for each choice of $`\beta `$ and $`K`$.

We carry out simulations employing either the plaquette (P) or rectangular action (R) for gluons, combining it with either the Wilson (W) or clover action (C) for quarks.

## III Simulations

### A choice of simulation parameters

We choose the coupling constant $`\beta `$ so that it gives an inverse lattice spacing of $`a^11`$ GeV. For each action combination we choose at least two values of $`\beta `$ to allow us to interpolate (or extrapolate) to a desired common lattice spacing.

Simulations are generally carried out at three values of the hopping parameter $`K`$ corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.7`$–0.9. The lattice size employed is $`12^3\times 32`$.

In Table I we give an overview of the calculations performed for the action comparison. Details of the simulation parameters at each run are collated in Appendix A. Our procedure for estimating the critical hopping parameter $`K_c`$, and the physical scale of lattice spacing either from the $`\rho `$ meson mass ($`a_\rho `$) or the string tension ($`a_\sigma `$) will be discussed in Sec. IV A and Sec. V C.

We take the RC<sub>pMF</sub> action at $`\beta =1.9`$ to explore how close one can take the calculation towards the chiral limit. For this study we employ two lattice sizes $`12^3\times 32`$ and $`16^3\times 32`$. In Table II we list the main features of these two runs whereas details can be found in Appendix A.

### B configuration generation and matrix inversion

Simulations are carried out for two flavors of dynamical quarks using the hybrid Monte Carlo (HMC) algorithm. The integration of molecular dynamics (MD) equations is made with the standard leapfrog scheme and with a step size $`\mathrm{\Delta }\tau `$ chosen to yield an acceptance ratio of 70–90% for trajectories of unit length. The actual values chosen for $`\mathrm{\Delta }\tau `$ in each case and the measured acceptance are given in Appendix A.

For the inversion of the fermion matrix we employed the minimal residue (MR) algorithm for our early simulations but switched later to BiCGStab. In both cases we use an even-odd preconditioning of the quark matrix $`D`$. $`D`$ can be decomposed into

$$D(K)=MK(D_{eo}+D_{oe}),$$
(5)   
where $`M`$ is only defined on single sites and the remaining connects neighboring sites. For the Wilson quark action $`M`$ is a unit matrix, whereas for the clover action it is non-trivial in color and Dirac space. The even-odd preconditioning consists of solving the equation $`AG_e=B_e^{}`$ where $`A=1K^2M_e^1D_{eo}M_o^1D_{oe}`$ and $`B_e^{}=M_e^1\left(B_e+KD_{eo}M_o^1B_o\right)`$ instead of the equation $`D(K)G=B`$. As an initial guess for the solution vector on even sites, the right-hand-side vector $`G_e=B_e^{}`$ is used. The preconditioning requires the inversion of the local matrix $`M`$, which is trivial for the Wilson quark action. For the clover quark action we precalculate $`M^1`$ and store it before the solver starts.

As a stopping condition for the matrix inversion during the fermionic force evaluation we generally use, on the $`12^3\times 32`$ lattice, the criterion

$$r_1=DGB^210^{10}$$
(6)   
which we found to be approximately equivalent to the condition

$$r_2=DGB/G10^8.$$
(7)   
The actual stopping conditions chosen for each run and the number of iterations needed to reach this condition are listed in Appendix A. For the evaluation of the Hamiltonian we choose stricter stopping criteria for $`r_1`$ between $`10^{14}`$ and $`10^{18}`$.

A necessary condition for the validity of the HMC algorithm is the reversibility of the MD evolution . The CP-PACS computer, on which the present work is made, employs 64 bit arithmetic for floating point operations. Flipping the sign of momenta after a unit trajectory, with the stopping condition (7) above, we checked that, (i) the gauge link and conjugate momenta variables return to the starting values within a relative error of less than $`10^7`$ on the average, and (ii) the relative error in the evaluation of the Hamiltonian is less than $`10^{10}`$ (absolute error better than $`10^4`$ for the $`16^3\times 32`$ lattice where the check was made) so that the effects in the accept/reject procedure are far below the level of statistical fluctuations.

At each simulated parameter we first run for 100–200 HMC trajectories of unit length for thermalization and then generate 500–1500 trajectories for measurements. Hadron propagators are measured on configurations separated by 5 trajectories. The static quark potential is measured on a subset of the configurations separated by either 5 or 10 trajectories. The detailed numbers are again given in Appendix A.

### C hadron mass measurement

We calculate quark propagators for the hopping parameter equal to that for the dynamical quarks used in the configuration generation. Two quark propagators are prepared for each configuration, one with the point source and the other with an exponentially smeared source with the smearing function $`\psi (r)=A\mathrm{exp}(Br)`$. For the latter we fix the gauge configuration to the Coulomb gauge. The choice of the smearing parameters $`A`$ and $`B`$ is guided by previous quenched results for the pion wave function , readjusted by hand so that hadron effective masses reach a plateau as soon as possible.

Hadron propagators are constructed by combining quark propagators for the point (P) or the smeared (S) sources in various ways, but always adopting the point sink. For example, PS represents a meson propagator calculated with the point source for quark and the smeared source for antiquark. In Fig. 1 we show a typical example of effective masses for a variety of source combinations.

In most cases the effective masses for the SS (SSS for baryons) propagators comes from below, show the best plateau behavior, and have the smallest statistical errors estimated with the jackknife procedure. We therefore determine hadron masses with a fit to SS (SSS) hadron propagators, supplementing results for other source combinations as a guide in choosing the fit range.

Hadron masses are extracted from propagators by employing a single hyperbolic cosine fit for mesons and a single exponential fit for baryons. We use uncorrelated fits and determine the error with the jackknife method. While our runs of at most 1500 HMC trajectories are not really long enough to carry out detailed autocorrelation analysis, examining the bin size dependence of the estimated error indicates a bin size of 5 configurations or 25 HMC trajectories to be a reasonable choice, which we adopt for all of our error analyses.

The hadron mass results for all our runs are collected in Appendix B.

### D potential measurement

We measure Wilson loops $`W(r,t)`$ both in the on- and off-axis directions in space. The spatial paths of $`W(r,t)`$ are formed by connecting one of the following spatial vectors repeatedly,

$$(1,0,0),(1,1,0),(2,1,0),(1,1,1),(2,1,1),(2,2,1).$$
(8)   
We measure $`W(r,t)`$ up to $`r6`$ and $`t8`$ on the $`12^3\times 32`$ lattice, while we enlarge the largest spatial size to $`r4\sqrt{3}`$ on the $`16^3\times 32`$ lattice in order to investigate the large distance behavior of the potential. The smearing procedure of Ref.  is applied to the link variables, up to 6 times on the $`12^3\times 32`$ lattice and up to 8 times on the $`16^3\times 32`$ lattice, respectively. The Wilson loop is measured at every smearing step in order to choose the optimal smearing number for each value of $`r`$.

We extract the potential $`V(r)`$ and the overlap function $`C(r)`$ by a fully correlated fit of the Wilson loop to the form

$$W(r,t)=C(r)\mathrm{exp}\left[V(r)t\right].$$
(9)   
The optimum smearing number at each $`r`$ is determined by the condition that the overlap $`C(r)`$ takes the largest value smaller than 1.

Typical results for the effective mass defined by

$$m_{\mathrm{eff}}=\mathrm{ln}\left[W(r,t)/W(r,t+1)\right],$$
(10)   
are shown in Fig. 2. We find that noise generally dominates over the signal for $`t>4`$. Thus we set the upper limit of the fitting range to $`t_{\mathrm{max}}=4`$. Since choosing the lower limit $`t_{\mathrm{min}}=1`$ leads to an increase of $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}`$ by 3–10 times compared to the choice $`t_{\mathrm{min}}=2`$ for most values of $`r`$ and simulation parameters, we fix the fitting range to be $`t=2`$–4.

The statistical error of $`V(r)`$ is estimated by the jackknife method. We find that a bin size of 30 HMC trajectories is generally sufficient to ensure stability of errors against bin size. We therefore adopt this bin size for all of our error estimates with potential data.

## IV Hadron spectrum

### A chiral fits

A basic parameter characterizing the chiral behavior of hadron masses is the critical hopping parameter $`K_c`$ at which the pseudo scalar meson mass $`m_{\mathrm{PS}}a`$ vanishes. Results for $`(m_{\mathrm{PS}}a)^2`$ exhibit deviations from a linear function in $`1/K`$, and hence we extract $`K_c`$ by assuming

$$(m_{\mathrm{PS}}a)^2=B_{\mathrm{PS}}\left(\frac{1}{K}\frac{1}{K_c}\right)+C_{\mathrm{PS}}\left(\frac{1}{K}\frac{1}{K_c}\right)^2.$$
(11)   
The fitted values of the critical hopping parameter are listed in Table I and II.

Another important parameter is the vector meson mass $`m_\mathrm{V}a`$ in the chiral limit $`m_{\mathrm{PS}}a=0`$, which allows us to set the physical lattice spacing. We determine this quantity by a chiral fit of the vector meson mass in terms of the pseudo scalar meson mass, both of which are measured quantities. Our results for this relation show curvature (see Fig. 8 in Sec. VI A for an example), and hence for the fitting function we employ

$$m_\mathrm{V}a=A_\mathrm{V}+B_\mathrm{V}(m_{\mathrm{PS}}a)^2+C_\mathrm{V}(m_{\mathrm{PS}}a)^3,$$
(12)   
where the cubic term is inspired by chiral perturbation theory.

A practical problem with this fit is that we have only three data points for most of our runs. We estimate systematic uncertainties in the extrapolation by repeating the fit without the cubic term to the two points of data for lighter quark masses. Results for the vector meson mass in the chiral limit, translated into the lattice spacing through $`a_\rho =A_\mathrm{V}/768\mathrm{M}\mathrm{e}\mathrm{V}`$, are listed in Table I and II.

Results for the nucleon and $`\mathrm{\Delta }`$ also show curvature in terms of $`m_{\mathrm{PS}}a`$. We therefore fit them employing a cubic polynomial without the linear term (12) as for the vector meson mass.

### B scaling of mass ratios

We show in Fig. 3 a compilation of our hadron mass results for the four action combinations in terms of the mass ratios $`m_\mathrm{N}/m_\mathrm{V}`$ and $`m_\mathrm{\Delta }/m_\mathrm{V}`$ as a function of $`(m_{\mathrm{PS}}/m_\mathrm{V})^2`$. In order to avoid overcluttering of points, we include results for only two values of $`\beta `$ per action combination. Furthermore, for the PC action combination the results with $`c_{\mathrm{SW}}=`$ MF are displayed whereas for the RC action results for $`c_{\mathrm{SW}}=`$ pMF are shown.

We observe two features in this figure. In the first instance, for each action combination the baryon to vector meson mass ratio decreases as the coupling decreases. This is a well-known trend of scaling violation for Wilson-type quark actions. Secondly, the magnitude of scaling violation, measured by the distance from the phenomenological curve (solid line in Fig. 3)  has an order where $`\mathrm{𝐏𝐖}>\mathrm{𝐑𝐖}>\mathrm{𝐏𝐂}>\mathrm{𝐑𝐂}`$. In particular the results for the PC and RC cases show a significant improvement over those for the PW and RW cases in that they lie close to the phenomenological curve even though the lattice spacing is as large as $`a_\rho ^11`$–1.3 GeV (see Tables I and II).

A point of caution, however, is that the lattice spacings for the data sets displayed in Fig. 3 do not exactly coincide. In order to disentangle effects associated with action improvement from those of a finer lattice spacing for each action, we need to plot results at the same lattice spacing.

One way to make such a comparison is to take a cross section of Fig. 3 at a fixed value of $`m_{\mathrm{PS}}/m_\mathrm{V}`$ and plot the resulting value of $`m_{\mathrm{N},\mathrm{\Delta }}/m_\mathrm{V}`$ as a function of $`m_\mathrm{V}a`$ at that value of $`m_{\mathrm{PS}}/m_\mathrm{V}`$. This requires an interpolation of hadron mass results, for which we employ the cubic chiral fits described in Sec. IV A and the jackknife method for error estimation.

In Fig. 4 we show results of this analysis for $`m_\mathrm{N}/m_\mathrm{V}`$ and $`m_\mathrm{\Delta }/m_\mathrm{V}`$ at $`m_{\mathrm{PS}}/m_\mathrm{V}=0.8`$. It is interesting to observe that the PW and RW results lie almost on a single curve, while the PC and RC results, respectively using the MF and pMF value of $`c_{\mathrm{SW}}`$, fall on a different, much flatter curve. This clearly shows that the improvement of the gauge action has little effect on decreasing the scaling violation in the baryon masses. The improvement is due to the use of the clover quark action for the PC and RC cases. An apparently better behavior of RW results in Fig. 3 compared to those for the PW case is merely an effect of the finer lattice spacing of the former.

We have commented in Sec. II that the values of $`c_{\mathrm{SW}}`$ for the MF and pMF cases are similar. This would explain why results for the PC action with the MF value of $`c_{\mathrm{SW}}`$ and those for the RC action with the pMF value of $`c_{\mathrm{SW}}`$ lie almost on a single curve. For both MF and pMF choices, the magnitude of $`c_{\mathrm{SW}}`$ is significantly larger than the tree-level value $`c_{\mathrm{SW}}=1`$. As is shown in Fig. 4 with open symbols, the degree of improvement with the tree-level $`c_{\mathrm{SW}}`$ is substantially less than those for the MF and pMF choices.

## V Static quark potential

### A restoration of rotational symmetry

In Fig. 5, we plot our potential data for the four action combinations at a quark mass corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.8`$ and $`a^11`$ GeV. We find a sizable violation of rotational symmetry in the PW case at this coarse lattice spacing. Looking at the potential for the PC case, we cannot observe any noticeable restoration of the symmetry. In contrast, a remarkable restoration of rotational symmetry is apparent in the RW and RC cases.

In order to quantify the violation of rotational symmetry and its improvement depending on the action choice, we consider the difference between the on-axis and off-axis potential at a distance $`r=3`$ defined by

$$\mathrm{\Delta }V=\frac{V\left(r=(3,0,0)\right)V\left(r=(2,2,1)\right)}{V\left(r=(3,0,0)\right)+V\left(r=(2,2,1)\right)}.$$
(13)   
We find that the value of $`\mathrm{\Delta }V`$ monotonously decreases as the sea quark mass decreases for most cases. We ascribe this trend to the fact that one effect of dynamical sea quarks is to renormalize the coupling toward a smaller value, and hence reduces violation of rotational symmetry.

In order to make a comparison at the same quark mass, we estimate $`\mathrm{\Delta }V`$ at $`m_{\mathrm{PS}}/m_\mathrm{V}=0.8`$ by an interpolation as a linear function of $`(m_{\mathrm{PS}}a)^2`$. In Fig. 6 we plot results for $`\mathrm{\Delta }V`$ obtained in this way against the value of $`m_\mathrm{V}a`$ at $`m_{\mathrm{PS}}/m_\mathrm{V}=0.8`$. This figure confirms the qualitative impression from Fig. 5. Rotational symmetry is badly violated for the PW and PC cases, which is significantly improved by changing the gauge action as demonstrated by the small values of $`\mathrm{\Delta }V`$ for the RW and RC results. In contrast the effect of quark action improvement on the restoration of rotational symmetry appears to be small. This may not be surprising since dynamical quarks affect the static potential only indirectly through vacuum polarization effects.

### B string tension

The static potential in full QCD is expected to flatten at large distances due to string breaking. None of our potential data, which typically extends up to the distance of $`r1`$ fm, show signs of such a behavior, but rather increase linearly. As we discuss in more detail in Sec. VI this is probably due to a poor overlap of the Wilson loop operator with the state of a broken string. This suggests that we can extract the string tension from the present data for the potential $`V(r)`$ by assuming the form

$$V(r)=V_0\frac{\alpha }{r}+\sigma r.$$
(14)   
In practice we find that the Coulomb coefficient $`\alpha `$ is difficult to determine from the fit, even if we introduce the tree-level correction term corresponding to the one lattice gluon exchange diagram. This may be due to the fact that our potential data taken at coarse lattice spacings do not have enough points at short distance to constrain the Coulomb term. As an alternative we test a two-parameter fitting with a fixed Coulomb term coefficient $`\alpha _{\mathrm{fixed}}=0.1`$, 0.125, …, 0.475, and 0.5, using the fitting range $`r_{\mathrm{min}}`$$`r_{\mathrm{max}}`$ with $`r_{\mathrm{min}}=1`$, $`\sqrt{2}`$, $`\sqrt{3}`$ and $`r_{\mathrm{max}}=5`$–6. We find that the value of $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}`$ takes its minimum value around $`\alpha _{\mathrm{fixed}}=0.3`$$`0.4`$ for most fitting ranges and simulation parameters.

Based on this result, we extract the string tension by fitting the potential at large distances, where a linear behavior dominates, to the form (14) with a fixed Coulomb coefficient $`\alpha _{\mathrm{fixed}}=0.35`$. The shift of the fitted $`\sigma `$ over the range $`\alpha =0.3`$$`0.4`$ is taken into estimates of the systematic error.

The result for the string tension $`\sigma `$ with this two-parameter fit is quite stable against variations of $`r_{\mathrm{max}}`$. It does depend more on $`r_{\mathrm{min}}`$, however. This leads us to repeat the two-parameter fit with $`\alpha _{\mathrm{fixed}}=0.35`$ over the interval of $`r_{\mathrm{min}}`$ listed in Appendix C, and determine the central value of $`\sigma `$ by the weighted average of the results over the ranges. The variance over the ranges are included into the systematic error of $`\sigma `$. We collate the final results for the string tension $`\sigma `$ in Appendix C.

### C consistency in lattice spacings

The scaling violation in the ratio $`m_\rho /\sqrt{\sigma }`$ leads to an inconsistency in the lattice spacings determined from the $`\rho `$ meson mass $`a_\rho `$ and the string tension $`a_\sigma `$ in the chiral limit. Thus, examination of this consistency provides another test of effectiveness of improved actions. For the physical value we use $`m_\rho =768`$ MeV and $`\sqrt{\sigma }=440`$ MeV. We should note that the latter value is uncertain by about 5–10% since the string tension is not a directly measurable quantity by experiment.

The chiral extrapolation of the vector meson mass was already discussed in Sec. IV A. We follow a similar procedure for the chiral extrapolation of the string tension. Namely we fit results to a form

$$\sigma a^2=A_\sigma +B_\sigma (m_{\mathrm{PS}}a)^2+C_\sigma (m_{\mathrm{PS}}a)^3.$$
(15)   
In most cases we find a quadratic ansatz ($`C_\sigma =0`$) to be sufficient, which we then adopt for all data sets. Results for the string tension in the chiral limit, converted to the physical scale of lattice spacing $`a_\sigma `$, are listed in Table I and II.

In Fig. 7 we plot $`m_\mathrm{V}a/768\mathrm{M}\mathrm{e}\mathrm{V}`$ and $`\sqrt{\sigma }a/440\mathrm{M}\mathrm{e}\mathrm{V}`$ as a function of $`(m_{\mathrm{PS}}a)^2`$ for the four action combinations with a similar lattice spacing $`a_\rho ^11`$$`1.3`$ GeV determined from the vector meson mass. A distinctive difference between the results for the Wilson and the clover quark action is clear; while results for $`m_\mathrm{V}`$ and $`\sqrt{\sigma }`$ cross each other at heavy quark masses where $`m_{\mathrm{PS}}/m_\mathrm{V}0.75`$$`0.8`$ for the PW and RW cases, leading to a mismatch of $`a_\rho `$ and $`a_\sigma `$ in the chiral limit, the two sets of physical scales converge well toward the chiral limit for the PC and RC cases.

We expect the large discrepancy for the Wilson quark action to disappear closer to the continuum limit. This is supported by the results obtained at $`\beta =5.5`$ with $`a^12`$ GeV in Ref.. Our results show that the clover term helps to improve the consistency between $`a_\rho `$ and $`a_\sigma `$ already at $`a^11`$ GeV.

## VI Approaching the chiral limit

The analyses presented so far show that the RC action has the best scaling behavior for hadron masses and static quark potential among the four action combinations we have examined. We then take this action and attempt to lower the quark mass as much as possible.

Two runs are made at $`\beta =1.9`$, one on a $`12^3\times 32`$ lattice down to $`m_{\mathrm{PS}}/m_\mathrm{V}0.5`$, and the other on a $`16^3\times 32`$ lattice down to $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$. We discuss results from these runs below.

### A hadrons with small quark masses

In Fig. 8 we plot the results of hadron masses as functions of $`(m_{\mathrm{PS}}a)^2`$. The existence of a curvature is observed, necessitating a cubic ansatz for extrapolation to the chiral limit. The lattice spacing determined from $`m_\rho =768`$ MeV equals $`a_\rho =0.20(2)`$ fm using mass results from the larger lattice. Hence the spatial size equals 2.4 fm ($`12^3\times 32`$) and 3.2 fm ($`16^3\times 32`$) for the two lattice sizes employed.

Finite-size effects are an important issue for precision determinations of the hadron mass spectrum. Our results in Fig. 8 do not show clear signs of such effects down to the second lightest mass, which corresponds to $`m_{\mathrm{PS}}/m_\mathrm{V}0.5`$. We feel, however, that it is premature to draw conclusions with the present low statistics of approximately 1000 trajectories.

The results for mass ratios are plotted in Fig. 9. While errors are large, and may even be underestimated because of the shortness of the runs, we find it encouraging that the ratios exhibit a trend of following the phenomenological curve toward the experimental points as the quark mass decreases. If we use the chiral extrapolation described above for the results on the $`16^3\times 32`$ lattice, we obtain $`m_\mathrm{N}/m_\mathrm{V}=1.342(25)`$ and $`m_\mathrm{\Delta }/m_\mathrm{V}=1.700(33)`$ at the physical ratio $`m_{\mathrm{PS}}/m_\mathrm{V}=0.1757`$, which are less than 10% off the experimentally observed ratios of 1.223 and 1.603, respectively, despite the coarse lattice spacing of $`a0.2`$ fm.

### B static potential at large distances

We have mentioned in Sec. V that our results for the static potential do not show signs of flattening, indicative of string breaking up to the distance of $`r1`$ fm. Similar results have been reported by other groups . A possible reason for these results is that potential data do not extend to large enough distances where string breaking becomes energetically favorable. Another related possibility is that the dynamical quark masses, which in most cases correspond to $`m_{\mathrm{PS}}/m_\mathrm{V}=0.7`$$`0.9`$, are too heavy. With our runs on the $`16^3\times 32`$ lattice we can examine these points up to the distance of $`r2`$ fm and for quark masses down to $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$.

In Fig. 10 we plot our potential data obtained on the $`16^3\times 32`$ lattice at the lightest sea quark mass corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$. We find that the potential increases linearly up to $`r2`$ fm, without any clear signal of flattening. The situation is similar for our data at heavier sea quark masses.

An interesting and crucial question here is whether the Wilson loop operator has sufficient overlap with the ground state at large $`r`$ so that the potential in that state is reliably measured there . In Fig. 11 we compare results for the overlap function $`C(r)`$ for the full QCD run at $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$ with that obtained in a quenched run with the R gauge action on a $`9^3\times 18`$ lattice at $`\beta =2.1508`$ ($`a^11`$ GeV) . For the quenched run the overlap $`C(r)`$ of the smeared Wilson loop operator with the ground string state is effectively 100 % at all distances. For full QCD, on the other hand, $`C(r)`$ significantly decreases as $`r`$ increases. Such a behavior of $`C(r)`$ is observed in all of our data including those taken with action choices other than RC. These results may be taken as a tantalizing hint that the Wilson loop operator develops mixings with states other than a single string, possibly a pair of static-light mesons in full QCD. We leave further investigations of this interesting question for future studies.

### C computer time

An important practical information in full QCD is the computer time needed for the approach to the chiral limit. In Table III we assemble the relevant numbers for our runs on the $`16^3\times 32`$ lattice. These runs have been performed on a partition of 256 nodes, which is 1/8 of the CP-PACS computer. For a partition of this size, our full QCD program, written in FORTRAN with the matrix multiplication in the quark solver hand-optimized in the assembly language, sustains about 37% of the peak speed of 75 GFLOPS. Adding the CPU time per trajectory of Table III, we find that accumulating 5000 trajectories for each of the 6 hopping parameters for this lattice size would take about 160 days with the full use of the CP-PACS computer. Carrying out such a simulation is certainly feasible. For larger lattice sizes such as $`24^3\times 48`$, however, we would have to stop at $`m_{\mathrm{PS}}/m_\mathrm{V}0.5`$ since the run at $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$ alone increases the computer time by a factor two. Let us add that the CPU time for a unit of HMC trajectory increases roughly proportional to $`(1/K1/K_c)^{1.6}`$ for the 4 smallest quark masses.

## VII Conclusions

In this paper we have presented a detailed investigation of the effect of improving the gauge and the quark action in full QCD. We have found that the consequence of improving either of the actions is different depending on the observable examined.

For the light hadron spectrum the clover quark action with a mean-field improved coefficient drastically improves the scaling of hadron mass ratios. Improving the gauge action, on the other hand, has almost no influence in this aspect. The SW-clover action also has the good property that the physical scale determined from the vector meson mass and the string tension in the chiral limit of the sea quark are consistent already at scales $`a^11`$ GeV, which is not the case with the Wilson quark action.

We have also confirmed that the use of improved gauge actions leads to a significant decrease of the breaking of rotational symmetry of the static quark potential.

Finally, we have made an exploratory simulation toward the chiral limit employing a renormalization group improved gauge and clover improved quark actions.

The results obtained in the present study suggest that a significant step toward a systematic full QCD simulation can be made with the present computing power using improved gauge and quark actions at relatively coarse lattice spacings of $`a^11`$–2 GeV.

###### Acknowledgements.

This work was supported in part by the Grants-in-Aid of the Ministry of Education (Nos. 08NP0101, 08640349, 08640350, 08640404, 08740189, 08740221, 09304029, 10640246 and 10640248). Two of us (GB, RB) were supported by the Japan Society for the Promotion of Science.    

## A Run Parameters

In this appendix we assemble information about our runs. An overview of the runs has been given in Table I. For the inversion of the quark matrix either the MR algorithm (M) or the BiCGStab algorithm (B) is used with the stopping condition $`r_1`$ stop defined through Eq.(6). During the HMC update $`D^{}D`$ has to be inverted. We do this in two steps, first inverting $`D^{}`$ and then $`D`$. In the tables we quote the number of iterations $`N_{\mathrm{inv}}`$ needed for the first inversion $`D^{}`$. Finally we also quote the statistics, giving the number of configurations for spectrum and potential measurements separately. Configurations for the hadron spectrum are separated by 5 HMC trajectories, whereas for the potential the separation is either 5 or 10 trajectories. Unless stated otherwise the lattice size is $`12^3\times 32`$.

## B Hadron Masses

In this appendix we assemble the results of our hadron mass measurements. We quote numbers for pseudo scalar and vector mesons, nucleons and $`\mathrm{\Delta }`$ baryons together with mass ratios against vector mesons. Additionally we quote numbers for the bare quark mass based on the axial Ward identity defined by

$$m_qa=m_{\mathrm{PS}}a\underset{t\mathrm{}}{lim}\frac{_\stackrel{}{x}A_4(\stackrel{}{x},t)P}{_\stackrel{}{x}P(\stackrel{}{x},t)P},$$
(B1)   
where $`A_4`$ is the local axial current and $`P`$ is the pseudo scalar density. Masses are extracted with an uncorrelated fit to the propagator and the errors are determined with the jackknife method with bin size 5.

## C String Tension