text:
# Chiral Logs in Quenched QCD

## I Introduction

One of the important goals of lattice QCD is to understand, from first principles, low-energy hadron phenomenology as a consequence of chiral symmetry. However, there have been problems with regard to chiral symmetry since the advent of the lattice formulation. The lack of chiral symmetry in traditional approaches (such as Wilson fermions and staggered fermions) leads to a plethora of conceptual and technical difficulties. In particular, due to critical slowing down and the existence of exceptional configurations, there is a practical limit on how low in quark (or pion) mass one can carry out a Monte Carlo simulation. The existing calculations with these fermions are typically limited to pion mass greater than $`300\mathrm{MeV}`$. As such, the extrapolation from the current lattice data with pion mass above $`300\mathrm{MeV}`$, sometimes well above, to the physical pion mass of $`137\mathrm{MeV}`$ is both a theoretical and a practical challenge. The primary concern is whether the quark mass region (of the lattice data) is inside the range where $`\chi `$PT applies. If not, one does not have a reliable analytic form in quark mass dependence to allow for an unbiased extrapolation. Furthermore, most of the quenched $`\chi `$PT calculations are carried out at one-loop order; this limits its range of applicability further, to even smaller quark masses. How small does the quark mass have to be in order to see the chiral behavior predicted by $`\chi `$PT? Is the strange quark in the radius of convergence of $`\chi `$PT? Is there a way to model the chiral extrapolation from well above a $`300\mathrm{MeV}`$ pion mass without introducing systematic errors? These were some of the issues addressed in the panel discussion during the Lattice 2002 conference ber02 .

Fortunately, recent progress in lattice fermions has led to formulations (such as domain-wall fermions, overlap fermions, and the fixed-point actions) with lattice chiral symmetry at finite cut-off. As a consequence, many chiral-symmetry relations neu01 ; nn93 ; nn95 ; neu98 ; lus98  and the quark propagator liu02  preserve the same structure as in the continuum. While a lot of progress has been made in checking the chiral symmetries in all these three formulations, we find that overlap fermions neu98  offer several distinct advantages. For example, it has been demonstrated that one can simulate near the physical quark-mass region with very gentle critical slowing down dll00 . Secondly, the overlap-fermion inverter can accommodate multiple masses in the calculation of quark propagators ehn99b  and thus is particularly well-suited for studying the details of the logarithmically-varying mass dependence of physical observables. Thirdly, we find that the $`O(a^2)`$ and $`O(m^2a^2)`$ errors are small liu02 ; dll00 ; ddh02 . (There are no $`O(a)`$ errors due to chiral symmetry knn97 .) For these reasons, in this paper we use overlap fermions in a quenched lattice calculation to study the chiral logs in the pion mass, the pseudoscalar decay constant $`f_P`$, the pion decay constant $`f_\pi `$, and the nucleon mass. We determine how small a quark mass or pion mass needs to be for one-loop $`\chi `$PT to be valid. In the present work, we only consider hadrons with degenerate quark masses. Quantities involving the strange quark such as the kaon mass and decay constant will be studied later.

## II Numerical details

Our calculation is done on a $`16^3\times 28`$ lattice with 80 quenched gauge configurations generated with the Iwasaki gauge action iwa85  at $`\beta =2.264`$ and the quark propagators are calculated with overlap fermions. The lattice spacing is $`0.200(3)\mathrm{fm}`$, as determined from the pion decay constant $`f_\pi (m_\pi )`$; the box size is $`3.2\mathrm{fm}`$.

The massive overlap Dirac operator neu98a  is defined so that at tree-level there is no mass or wavefunction renormalization ddh02 ,

$$D(m)=\rho +\frac{m}{2}+(\rho \frac{m}{2})\gamma _5ϵ(H).$$
(1)
Here $`ϵ(H)=H/\sqrt{H^2}`$ is the matrix sign function and $`H_W`$ is taken to be the hermitian Wilson-Dirac operator, i.e. $`H=H_W=\gamma _5D_\mathrm{W}`$. Here $`D_\mathrm{W}`$ is the usual Wilson fermion operator, except with a negative mass parameter $`\rho =1/2\kappa 4`$ in which $`\kappa _c<\kappa <0.25`$. We take $`\kappa =0.19`$ in our calculation which corresponds to $`\rho =1.368`$.

Throughout the paper, we shall use lattice units for dimensionful quantities, except that the lattice spacing $`a`$ will be explicit in figures.

### II.1 Zolotarev Approximation

The numerical aspects of the calculation have been given elsewhere dll00 ; ddh02 . The new ingredient is the adoption of the Zolotarev optimal rational approximation vfl02 ; chh02  of the matrix sign function.

To approximate the function $`\frac{1}{\sqrt{x}}`$, it is suggested vfl02  that one should consider the minimization of

$$1\sqrt{x}f(x)_{\mathrm{}}^{[\lambda _{min}^2,\lambda _{max}^2]},$$
(2)
with $`f(x)`$ being approximated by a rational polynomial $`\stackrel{~}{f}(x)R^{m1,m}=p(m1)/q(m)`$ where $`p(m1)`$ and $`q(m)`$ are polynomials with degree $`m1`$ and $`m`$, respectively. It was proven by Zolotarev zol1877  that $`x\stackrel{~}{f}(x^2)R^{2m1,2m}`$ is the $`||.||_{\mathrm{}}`$-optimal rational approximation for the sign function on $`[|\lambda _{max}|,|\lambda _{min}|][|\lambda _{min}|,|\lambda _{max}|]`$. By way of Zolotarev’s theorem zol1877 ; idk00 ; pp87 , there is an analytic solution

$$\stackrel{~}{f}(x)=A\frac{_{l=1}^{m1}(x+c_{2l})}{_{l=1}^m(x+c_{2l1})},$$
(3)
where the coefficients can be expressed in terms of Jacobian elliptic functions

$$c_l=\frac{sn^2(lK/2m;\kappa )}{1sn^2(lK/2m;\kappa )},l=1,2,\mathrm{},2m1,$$
(4)
where $`\sqrt{1\kappa ^2}=|\frac{\lambda _{max}}{\lambda _{min}}|`$ and $`K`$ is the complete elliptic integral. According to Ref. vfl02 , $`A`$ is determined by the condition

$$\underset{C[1,(\frac{\lambda _{\text{max}}}{\lambda _{\text{min}}})^2]}{\mathrm{max}}[1\sqrt{x}f(x)]=\underset{C[1,(\frac{\lambda _{\text{max}}}{\lambda _{\text{min}}})^2]}{\mathrm{min}}[1\sqrt{x}f(x)].$$
(5)
However, there is a worry when the sign function $`ϵ(H)`$ is not bounded by unity because of slight numerical imprecision. From the Cauchy-Schwartz inequality $`\gamma _5\stackrel{~}{ϵ}(H)<\gamma _5\stackrel{~}{ϵ}(H)`$, the overlap operator $`D(0)=1+\gamma _5ϵ(H)`$ with $`ϵ(H)`$ approximated by $`\stackrel{~}{ϵ}(H)`$ could have superfluous unphysical zero modes neu00b . In view of this, we fine-tune the parameter $`A`$ by reducing it a little from the value determined from Eq. (5) to minimize the positive residual. This doubles the worst residuals.

For the Wilson action kernel $`H_W=\gamma _5D_W`$ inside the sign function, the largest eigenvalue is around 2.43 and is fairly stable over the configuration space. Furthermore, we find that the eigenstates of $`H_W`$ to be projected out become dense at $`0.050.06`$. Thus, it is sufficient to specify $`\sqrt{1\kappa ^2}=|\frac{\lambda _{max}}{\lambda _{min}}|=80`$. With $`|\lambda _{max}|=2.5`$, this gives $`|\lambda _{min}|=0.0315`$ which is smaller than the largest eigenstates ($`0.050.06`$) we project out from the kernel $`H_W`$. Since the coefficients of the rational polynomial approximation are expressed in terms of the Jacobian elliptic functions, one can use as high a degree of the polynomials as the memory can hold in the multi-mass algorithm. In practice, we use 14th degree which is a sufficiently good approximation for our calculation. We plot the residual $`ϵ_r(x)=\frac{x}{\sqrt{x^2}}Z_{14}(x)`$ for $`x>0`$ in Fig. 1 as a function of $`x`$. Here $`Z_{14}(x)`$ is the 14th degree Zolotarev approximation $`x\stackrel{~}{f}(x^2)`$ from Eq. (2). As we can see, in the selected window \[0.031, 2.5\], the approximation is better than $`3.3\times 10^{10}`$. This is a good deal better than the earlier attempts which approximate $`g(x)`$ to $`\frac{1}{\sqrt{x}}`$ with the Remez algorithm in which case the resultant approximation to the sign function with $`xg(x^2)`$ is not optimal in the range $`[|\lambda _{max}|,|\lambda _{min}|][|\lambda _{min}|,|\lambda _{max}|]`$. We plot the residual $`ϵ_r(x)`$ for the 14th degree ehn99b  approximation $`R(14)`$ and 12th degree dll00  approximation $`R(12)`$ with the Remez method in comparison with the 14th degree Zolotarev approximation $`Z(14)`$. We see that the Remez method gives residuals in the order of $`10^5`$ in the relevant range \[0.031, 2.5\]. We note that according to Chebyshev’s theorem, the errors of the optimal rational approximation $`\stackrel{~}{f}(x)R^{m,n}=p(m)/q(n)`$ have $`m+n+2`$ alternating signs in the range of the approximation. Since the Zolotarev approximation to the sign function is optimal, it should satisfy this requirement. In other words, there should be 29 alternating maxima and minima (including the endpoints of the range) for the approximation with $`m=13,n=14`$. As seen in Fig. 1, this feature is preserved with the restriction allowing only positive residuals, as described above.

After partial fraction expansion, the approximated sign function in rational polynomial $`x\stackrel{~}{f}(x^2)R^{2m1,2m}`$ has the form

$$\frac{x}{\sqrt{x^2}}x\underset{i=1}{\overset{N=14}{}}\frac{c_i}{x^2+q_i}.$$
(6)
The coefficients $`c_i`$ and $`q_i`$ for our 14th degree Zolotarev approximation are listed in Table 1.

In order to improve the approximation of the matrix sign function as well as the convergence in the conjugate gradient inversion in the inner do loop ehn99b ; dll00 ; ddh02 , it is desirable to project out the lowest several eigenmodes of $`H_W`$. We use the Iwasaki improved gauge action which requires that a smaller number of eigenmodes be projected, compared to the Wilson and Lüscher-Weisz gauge actions eh02 . We calculated 140 small eigenvalues $`\lambda `$ and projected out the lowest $`110`$$`140`$ eigenmodes $`\mathrm{\Psi }_\lambda `$, based on the criterion that the residual of the eigenvalue, defined as $`H_W\mathrm{\Psi }_\lambda \lambda \mathrm{\Psi }_\lambda /\lambda \sqrt{\mathrm{dim}.(H_W)}`$, be less then $`10^7`$. The next smallest eigenvalue is $`0.05`$$`0.06`$, which is well within the window where the approximation for the sign function is better than $`3.3\times 10^{10}`$. The residuals of both the inner and the outer do loops for the inversion of the massive overlap Dirac operator in Eq. (1) are at the level of $`10^7`$.

## III Pion mass and zero modes

We look at the pion mass as calculated from the pseudoscalar correlator $`_\stackrel{}{x}P(\stackrel{}{x},t)P(0)`$ where $`P=\overline{\psi }i\gamma _5(1D/2)\psi `$. The pion masses $`m_\pi `$ obtained from this correlator are listed in Table 2 and plotted in Fig. 2 as a function of the quark mass $`m`$. Our smallest pion mass is $`182(8)\mathrm{MeV}`$.

Consider the contributions of the zero modes to the meson correlator of the interpolation field $`M=\overline{\psi }\mathrm{\Gamma }(1D/2)\psi `$,

$`{\displaystyle d^3xM(x)M^{}(0)}|_{zeromodes}=`$ (7)
$`{\displaystyle }d^3x[{\displaystyle \underset{i,j=zeromodes}{}}{\displaystyle \frac{\mathrm{𝚝𝚛}(\psi _j^{}(x)\gamma _5\mathrm{\Gamma }\psi _i(x))\mathrm{𝚝𝚛}(\psi _i^{}(0)\overline{\mathrm{\Gamma }}\gamma _5\psi _j(0))}{m^2}}`$
$`+2{\displaystyle \underset{i=0,\lambda >0}{}}{\displaystyle \frac{\mathrm{𝚝𝚛}(\psi _\lambda ^{}(x)\gamma _5\mathrm{\Gamma }\psi _i(x))\mathrm{𝚝𝚛}(\psi _i^{}(0)\overline{\mathrm{\Gamma }}\gamma _5\psi _\lambda (0))}{\lambda ^2+m^2}}],`$
where $`\mathrm{\Gamma }`$ is the spinorial matrix and $`\overline{\mathrm{\Gamma }}=\pm \mathrm{\Gamma }^{}`$. For the pseudoscalar case with $`\mathrm{\Gamma }=\gamma _5`$, the zero modes contribute to both the direct and cross terms in the correlator. They contribute also to the cross terms in the $`_\stackrel{}{x}A_4(\stackrel{}{x},t)P(0)`$ and $`_\stackrel{}{x}A_4(\stackrel{}{x},t)A_4(0)`$ correlators, where $`A_4=\overline{\psi }i\gamma _4\gamma _5(1D/2)\psi `$, but not to the direct terms. In view of the worry that the pseudoscalar correlator may be contaminated by the zero modes bcc00 ; ddh02  particularly at small volumes, we look at $`_\stackrel{}{x}A_4(\stackrel{}{x},t)A_4(0)`$ to check if the masses calculated are the same as those from the $`PP`$ correlator. We plot the ratios of the pion masses from the $`PP`$ correlator and the $`A_4A_4`$ correlator in Fig. 3. We see that the ratios are basically unity within errors all the way to the smallest quark mass. This suggests that at this volume, i.e. $`V=(3.2\mathrm{fm})^4`$, there is no detectable contamination due to the zero modes. However, this is not so at smaller volumes. We will show results on the $`12^3\times 28`$ lattice and discuss the volume dependence in Sec. V.

It has been shown that the correlator $`PPSS`$ can get rid of the zero mode contribution bcc00  for both the direct and cross terms. This is a good procedure for medium quark masses. For small quark masses, there is a complication due to the quenched $`\eta ^{}`$-$`\pi `$ ghost state contribution in the $`SS`$ bde02  which has a mass close to two times the pion mass. We have fitted the $`PPSS`$ correlator with the inclusion of the ghost state which has a negative weight. We show in Fig. 4 the ratio of pion masses determined from the $`PP`$ and the $`PPSS`$ correlators. We note that the errors from the $`PPSS`$ are larger for $`m_\pi 0.4438(36)`$ due to the presence of the ghost state which is close to the ground state pion mass. It is not clear to what extend the additional complication due to the ghost state is disentangled from that of the zero modes. In any case, we find that the pion masses obtained are consistent with those from the $`PP`$ and $`A_4A_4`$ correlators. Their central values differ by 2.5% at most. We will come back to this issue in the next section when we fit the quenched chiral log in $`m_\pi `$.

## IV Chiral log in $`m_\pi `$

It is predicted that in the quenched approximation, there are quenched chiral logs sha92 ; bg92  arising from the hairpin diagrams; the flavor-singlet pseudoscalar correlator does not yield a full $`\eta ^{}`$ mass due to the absence of dynamical fermion loops. Instead, the would-be $`\eta ^{}`$ propagator gives a double pole of the Goldstone boson. The predicted form for the one-loop formula in $`\chi `$PT is the following sha97 ; hsw00

$$m_\pi ^2=Am\{1\delta [\mathrm{ln}(Am/\mathrm{\Lambda }_\chi ^2)+1]\}+A_\alpha m^2[1+2\mathrm{ln}(Am/\mathrm{\Lambda }_\chi ^2)]+Bm^2,$$
(8)
where the coefficients $`A,A_\alpha `$ and $`B`$ are given in terms of the parameters $`\mathrm{\Sigma },f,\alpha _\mathrm{\Phi }`$ and $`2\alpha _8\alpha _5`$ in the quenched effective chiral lagrangian, i.e.

$`A`$ $`=`$ $`2\mathrm{\Sigma }/f^2,`$
$`A_\alpha `$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{\Phi }A^2}{3(4\pi f)^2}},`$
$`B`$ $`=`$ $`{\displaystyle \frac{(2\alpha _8\alpha _5)A^2}{(4\pi f)^2}},`$ (9)
and $`\mathrm{\Lambda }_\chi `$ is an unphysical cutoff scale which is usually taken in $`\chi `$PT to be in the range of $`0.5`$$`1.0\mathrm{GeV}`$. The parameter $`\alpha _\mathrm{\Phi }`$ is the singlet coupling in the quenched theory. The quenched chiral log parameter $`\delta `$ is

$$\delta =\frac{m_0^2}{16\pi ^2N_ff^2},$$
(10)
where $`m_0870\mathrm{MeV}`$ from the Veneziano model of the $`\eta ^{}`$ mass. From this, one estimates that $`\delta =0.183`$.

There are several calculations to extract the chiral log parameter $`\delta `$ abb00 ; bde00 ; qcdsf00 ; milc01 ; ko00 ; bcc00 ; aok01 ; gat02 ; ch02 . Save for Refs. gat02 ; ch02 , most of the calculations obtained values for $`\delta `$ in the range $`0.060.1`$ which are smaller than the expectation. We have summarized these results in Table 3. The present simulation adopts the overlap fermion action which has exact chiral symmetry on the lattice, and a highly accurate numerical approach (the matrix sign function is approximated to better than $`10^9`$). This should be a valuable tool to probe the relevant range of applicability of the one-loop $`\chi `$PT with pion mass as low as $`180\mathrm{MeV}`$ on a fairly large volume ($`3.2\mathrm{fm}`$). We plot in Fig. 5 the $`PP`$ results for $`m_\pi ^2/m`$ as a function of the bare quark mass $`m`$. We see clearly that below $`m=0.1295`$ ($`m_\pi 480\mathrm{MeV}`$), there is a dramatic divergent behavior. The ratio of the highest value at $`m=0.014`$ ($`m_\pi 180\mathrm{MeV}`$) to that at the trough at $`m=0.1295`$ ($`m_\pi 480\mathrm{MeV}`$) is $`1.33`$. Thus, the percentage excess, $``$, of the highest value at the minimum pion mass $`m_{\pi ,min}`$ relative to that at the trough at $`m_{\pi ,trou}`$ is 33%. Most of the calculations studying the quenched chiral log stopped above $`m_\pi 300\mathrm{MeV}`$; thus, they seem to be seeing the onset of this behavior with corresponding excess $`6\%11\%`$. In Table 3, we summarize the attempts made by various groups to calculate this quenched chiral log from $`m_\pi ^2/m`$ with different actions and lattices. Except for the present work, all the results are fitted to Eq. (8) without the $`A_\alpha `$ term.

We see from Table 3 that the results are not consistent with each other. There is a large spread in $`\delta `$ and $`m_{\pi ,trou}`$ which invokes many issues to ponder: To what extent does the fermion action matter? Is $`\chi `$PT valid in the range of pion masses calculated? If so, is the one-loop formula sufficient and how does one judge it? The parameters $`A`$, $`\delta `$, $`\mathrm{\Lambda }_\chi `$, and $`B`$ are not linearly independent in Eq. (8); how does one fit them? Practically all the calculations so far ignored $`A_\alpha `$; does it change the results much if included? How large are the finite volume effects for the range of pion masses in the calculation? (As will be demonstrated in Section V, the pion masses calculated from the $`PP`$ correlator in small volumes are contaminated by the zero modes even up to fairly large quark masses, i.e. the strange quark region. This could affect the position of $`m_{\pi ,trou}`$ and the result for $`\delta `$.) How does one control the finite volume and zero mode effects?

Before addressing these questions, we first point out that there is a fundamental problem in the way the lattice data are usually fitted to the form predicted by the $`\chi `$PT. The common practice is to fit the lattice data to the $`\chi `$PT result from a high pion mass down to cover the range of calculated masses and then extrapolate to the physical pion mass. This can force a fit but cannot answer the question as to whether the one-loop chiral formula is applicable for the range of masses considered. It is possible that the range of masses will require higher orders in $`\chi `$PT, or even that the mass range is outside the radius of convergence of $`\chi `$PT. Then one does not know the precise form to fit. Extrapolating down from the pion mass as high as $`1\mathrm{GeV}`$ is simply diametrically opposed to the philosophy of $`\chi `$PT which is based on the expansion of small quark masses and momenta.

So, how does one proceed? To be consistent with the $`\chi `$PT approach, we start from the smallest quark mass and ask where the one-loop formula applies note1 . Since we have come down to the pion mass as low as $`180\mathrm{MeV}`$, we will assume that it is in the range where one-loop $`\chi `$PT is valid. We fit the one-loop formula with the five smallest pion masses and then incrementally add higher pion masses one by one. If the formula is valid for the range of masses, then the parameters will remain constant, i.e. they are truly the low energy constants that we seek to evaluate. On the other hand, if the parameters begin to change as one includes more higher masses, we take it as an indication that the one-loop formula starts to lose its applicability. This could be due to the fact that higher orders are needed or that the mass is simply outside the range of validity of $`\chi `$PT.

We mentioned earlier that in the fitting of the one-loop formula, the parameters $`A,\delta `$, $`\mathrm{\Lambda }_\chi `$, and $`B`$ are not linearly independent. We can rewrite Eq. (8) as

$$m_\pi ^2=C_1m+C_{1L}m\mathrm{ln}(m)+C_2m^2+C_{2L}m^2\mathrm{ln}(m)$$
(11)
where

$`C_1`$ $`=`$ $`A[1\delta (\mathrm{ln}A2\mathrm{ln}\mathrm{\Lambda }_\chi +1)],`$
$`C_{1L}`$ $`=`$ $`A\delta ,`$
$`C_2`$ $`=`$ $`A_\alpha (1+2\mathrm{ln}A4\mathrm{ln}\mathrm{\Lambda }_\chi )+B,`$
$`C_{2L}`$ $`=`$ $`2A_\alpha ,`$ (12)
are the independent parameters. Therefore, we have one redundant variable which we take to be the cutoff scale $`\mathrm{\Lambda }_\chi `$. We adopt a Bayesian-based constrained curve-fitting method lcd01 ; mor01  with an adaptive procedure to obtain priors sequentially according to the relative importance of the parameters. The details of the constrained curve-fitting algorithm for the hadron masses will be given elsewhere ddh03 . In each of the quark mass ranges, the priors for $`C_1`$, $`C_{1L}`$, $`C_2`$, and $`C_{2L}`$ are obtained from the corresponding $`A_4A_4`$ correlators. After the $`C`$ parameters are determined from the fit, we solve Eq. (IV) to obtain $`A`$, $`\delta `$, $`A_\alpha `$, and $`B`$ for a given $`\mathrm{\Lambda }_\chi `$. We list their results evaluated at $`\mathrm{\Lambda }_\chi =0.8`$ together with $`C_1`$ and $`C_{1L}`$ in Table 4. We also plot $`A`$, $`A_\alpha `$, $`B`$, and $`\delta `$ as a function of the maximum pion mass $`m_{\pi ,max}`$ of the fitting range in Figs. 6, and 7. We should mention that due to the limited number of gauge configurations, we are not able to carry out correlated fits with a reasonable $`\chi ^2`$ for more than 17 quark masses. Therefore, for longer ranges, we choose a maximum of 17 representative masses to cover the range of the fit. For $`m_{max}`$ less than $`0.08983`$ ($`m_\pi =0.4050(38)`$ or $`400(7)\mathrm{MeV}`$), we find that the data are dominated by $`C_{1L}`$ and $`C_1`$. The contributions from $`C_2`$ and $`C_{2L}`$ cancel each other and the combined contribution is negligible compared to those of $`C_{1L}`$ and $`C_1`$. This is illustrated in Fig. 8 for a typical fit. As a result, it is practically impossible to obtain a reliable fit for $`C_2`$ and $`C_{2L}`$. We decided to fit these low mass data using a prior for $`C_{2L}`$ which is limited to less than 2 during the fitting procedure for the $`A_4A_4`$ correlators. We tried using the priors without such a limit and also tried dropping $`C_2`$ and $`C_{2L}`$ all together in these low mass ranges and found that the fitted $`C_1`$ and $`C_{1L}`$ are not changed within errors. The central values for $`\delta `$ are changed only by $`12`$%. In view of this inability to fit $`C_2`$ and $`C_{2L}`$, we do not quote the fitted $`A_\alpha `$ and $`B`$ below $`m_{max}=0.08983`$ ($`m_\pi =0.4050(38)`$) where the errors are larger than their central values. Hopefully, with better statistics, one will be able to extract the parameters $`A_\alpha `$ and $`B`$ in the range $`m_\pi <400(7)\mathrm{MeV}`$.

We see in Fig. 6 that values for $`A`$ are fairly stable throughout the range of the fit. $`A_\alpha `$ and $`B`$ are also stable in the range where $`m_{max}0.08983`$ and $`A_\alpha `$ is consistent with zero in this range. On the contrary, $`\delta `$ in Fig. 7 is constant up to $`m_{max}0.044`$ which corresponds to $`m_{\pi ,max}300\mathrm{MeV}`$. Beyond that, it has a sharp drop. We take the fact that $`\delta `$ is sensitive to the range of masses as an indication that the one-loop chiral formula in Eq. (8) is valid only up to $`m_\pi 300\mathrm{MeV}`$. Below this value, the parameters in the formula do not exhibit mass dependence and thus they are the low energy constants of the effective chiral theory at a certain chiral scale (e.g. $`\mathrm{\Lambda }_\chi =0.8`$).

Next, we simulate what has been commonly done in lattice calculations. We fix the maximum pion mass at $`m_\pi =0.8740(30)`$ ($`862(13)\mathrm{MeV}`$) and incrementally include lower pion masses to fit Eq. (11). We find that no matter how many low pion masses are included, the results are not compatible with those starting from $`m_{\pi ,min}=0.1844(71)`$ ($`182(8)\mathrm{MeV}`$) and ending at $`m_{\pi ,max}300\mathrm{MeV}`$. We illustrate this by plotting, in Fig. 7, the fitted $`\delta `$ (open circles) as a function of the minimum pion mass $`m_{\pi ,min}`$ for the range of the fit. We see that when pion masses down to $`300400\mathrm{MeV}`$ are included, the fitted $`\delta `$ is in the range of $`0.15`$ which is consistent with those found in Table 3 with a similar mass range and can explain why most of the calculations end up with a small $`\delta `$. The only exception is Ref. ch02  which obtains a large $`\delta 0.2`$ with $`m_{\pi ,min}440\mathrm{MeV}`$. We speculate that this may be due to the small volume used ($`L=1.5\mathrm{fm}`$). As will be explored further in Section V, we find that the zero modes make a sizable contribution to pion masses as high as $`800\mathrm{MeV}`$ on a lattice with $`L=2.4\mathrm{fm}`$, and likely higher at smaller volumes. It would be sensible to check if this is the case by comparing pion masses obtained from the $`PP`$ correlator to those from $`A_4A_4`$ as will be done in Section V.

In view of the fact that the rise of $`m_\pi ^2/m`$ left of the trough at $`m_\pi 480`$ MeV is mostly due to the chiral log while the rise right of the trough is mostly due to the $`m^2`$ term, we first fit the points around the trough (i.e. from $`m=0.07583`$ to 0.22633) with the form in Eq. (11) but without the $`A_\alpha `$ term and then extend the range by including smaller quark masses. The results of $`A,\delta `$ and $`B`$ turn out to be very close to the ones starting from the maximum pion mass at $`m_\pi =0.8740(30)`$ ($`862(13)\mathrm{MeV}`$). We show the results of $`\delta `$ (labeled $`\delta _{trou}`$) as a function of the minimum pion mass $`m_{\pi ,min}`$ in Fig. 7 and find them almost the same as $`\delta _{min}`$ at small quark masses and similarly smaller than those where the pion masses are restricted to lower than $`300\mathrm{MeV}`$.

For comparison, we show in Fig. 9 the fit of $`m_\pi ^2/m`$ for the forward fitting ending at $`m_\pi =0.269(5)`$ ($`m=0.03617`$) (solid line) and the backward fitting from $`m_\pi =0.641(2)`$ ($`m=0.22633`$) to the smallest quark mass (dash line).

As one can see from Eq. (IV), the parameters $`A`$, $`\delta `$, and $`B`$ depend logarithmically on $`\mathrm{\Lambda }_\chi `$. In Fig. 10 we show how $`A`$ and $`\delta `$ vary as a function of $`\mathrm{\Lambda }_\chi `$. The $`C_1`$ and $`C_{1L}`$ are taken from the fit for $`m_{max}=0.03617`$ ($`m_\pi =0.2693(48)`$ or $`266(6)\mathrm{MeV}`$) in Table 4. Incidentally, for a given $`C_1`$ and $`C_{1L}`$, there is a maximal value

$$\mathrm{\Lambda }_{\chi ,\mathrm{max}}=\sqrt{C_{1L}}e^{C_1/2C_{1L}}$$
(13)
(for which $`\delta =1`$ and $`A=C_1`$) beyond which no solution exists. But this is well beyond the range of applicability of Eq. (8). We should point out that the fact $`A`$ runs with $`\mathrm{ln}\mathrm{\Lambda }_\chi `$ is a peculiar feature of quenched QCD. This is not the case in full QCD where $`A`$ is basically the quark condensate where there is no quenched chiral $`\delta `$ sha92 . From Eq. (IV), we see that the combination $`A\delta `$ is independent of the chiral scale. This makes $`\delta `$ also run with $`\mathrm{ln}\mathrm{\Lambda }_\chi `$. This is not physical, but rather is due to the fact that Eq. (8) contains only the leading term in the expansion of $`\delta `$ sharpe03 . To see this, consider the example where the leading logarithm is re-summed for the “cactus diagrams” sha92  (one is not concerned with $`A_\alpha `$ and $`B`$ here) which leads to

$$m_\pi ^2=m^{\frac{1}{1+\delta }}(A\mathrm{\Lambda }_\chi ^{2\delta })^{\frac{1}{1+\delta }}.$$
(14)
We see that the $`m`$ and $`\mathrm{\Lambda }_\chi `$ dependence of $`\delta `$ are separable which implies that $`\delta `$ does not run with $`\mathrm{\Lambda }_\chi `$, and the unphysical $`\mathrm{\Lambda }_\chi `$ can be absorbed in $`A`$ with

$$A\mathrm{\Lambda }_\chi ^{2\delta }$$
(15)
which gives

$$\frac{d\mathrm{ln}A}{d\mathrm{ln}\mathrm{\Lambda }_\chi ^2}=\delta .$$
(16)
This is the same as obtained from Eq. (8) when one requires that $`m_\pi ^2`$ be independent of $`\mathrm{\Lambda }_\chi `$ and assumes that $`\delta `$ does not depend on $`\mathrm{\Lambda }_\chi `$.

We shall report our final results based on the fit of pion mass from $`0.184(7)`$ ($`m=0.014`$) to $`0.269(5)`$ ($`m=0.03617`$) with 9 quark masses. This choice is based on the observation that $`\delta `$ begins to decrease beyond $`m_{max}=0.03617`$ ($`m_{\pi max}=0.269(5)`$). Although this happens to coincide with the minimum in $`\chi ^2/dof`$ in Table 4, we should stress that a minimal $`\chi ^2/dof`$ is not and should not be the criterion for such a choice. As we pointed out earlier, the primary difficulty in chiral extrapolation is the uncertainty about the valid range of one-loop formula in $`\chi `$PT. A $`\chi ^2/dof`$ less than unity does not imply that the fitting formula is correct. One way to judge is the criterion that we adopt, namely the one-loop formula is not valid when the fitted parameters do not remain constant beyond certain $`m_{\pi max}`$.

As we said earlier, we are not able to obtain a reliable fit for $`A_\alpha `$, and $`B`$ at these low masses. As for the finite volume systematic errors, we note that the finite volume correction for the lowest $`m_\pi ^2`$ is estimated to be 5.4% and the second lowest $`m_\pi ^2`$ is estimated to be 3.6%. (This will be explained in more detail in Section V.) We thus estimate the systematic errors due to the finite volume to be at the few percent level which are smaller than our current statistical errors. We also dropped the last quark mass and fitted the 8 quark masses from $`m=0.01633`$ ($`m_\pi =0.195(7)`$) to $`m=0.03617`$ ($`m_\pi =0.269(5)`$) and found that $`A`$ and $`\delta `$ change only about 1%, which are much smaller than the statistical errors.

We quote our final results for the low-energy parameters of the quenched one-loop chiral perturbation theory at $`\mathrm{\Lambda }_\chi =0.8\mathrm{GeV}`$ close to the $`\rho `$ mass in Table 5. In view of the unphysical dependence of $`\delta `$ on $`\mathrm{\Lambda }_\chi `$ at the the order specified in Eq. (8), we shall estimate the systematic error to be the average deviation from the value at $`0.8\mathrm{GeV}`$ of the values at $`\mathrm{\Lambda }_\chi =1\mathrm{GeV}`$ and $`0.6\mathrm{GeV}`$. We see that the central value of $`\delta `$ is somewhat larger than the phenomenological value of $`0.183`$ in Eq. (10), but they are consistent within errors.

### IV.1 Fitting of Re-summed Cactus Diagrams

Since at small quark masses, the higher order results with $`(\delta \mathrm{ln}(m))^n`$ corrections are important, it is pertinent to fit the re-summed cactus diagram form sha92  for the chiral log which leads to a scale independent $`\delta `$. Thus, we shall fit with the power form

$$m_\pi ^2=Am^{\frac{1}{1+\delta }}+Bm^2.$$
(17)
Similar to fitting the one-loop formula in Eq. (11), we start from the smallest quark mass and fit upward first. Since at very small quark masses the $`B`$ contribution is very small, we put a very weak constraint at $`1.8\pm 1`$ which covers the range of $`B`$ as fitted in one-loop formula in Table 4. The resultant $`\delta ,A`$ and $`B`$ are presented in Figs. 11, 12, and 13. We see that, in contrast to the sharp transition at $`m_{\pi ,max}300\mathrm{MeV}`$ in the one-loop formula fit, the transition for $`\delta `$ is smoother. The behavior for $`A`$ and $`B`$ are similarly changed. This is to be expected, since the power form takes into account the re-summed cactus diagrams to all orders. As such, it should extend the range of applicability of the quenched $`\chi `$PT to beyond $`m_\pi 300\mathrm{MeV}`$ for the one-loop case, provided that it dominates the higher loop corrections. From the fitted results, it appears that the power form gives a stable fit to the trough of $`m_\pi ^2/m`$ and is perhaps valid to $`m_\pi 500600\mathrm{MeV}`$. Above $`m_\pi 640\mathrm{MeV}`$, all the parameters ($`\delta ,A`$ and $`B`$) are different from those below $`m_{\pi ,max}300400\mathrm{MeV}`$ by two sigmas. We shall quote the result of the scale independent $`\delta `$ at $`m=0.08083(m_\pi =400(7)\mathrm{MeV})`$ where $`\delta =0.20(3)`$. We have fitted also the pion mass from the $`PPSS`$ correlator and find that the changes of $`\delta `$ and $`A`$ are small, at a few percent level. In the range between $`m_\pi =0.2501(54)`$ and 0.4840(32), the change is at the 1% level. For example, $`\delta =0.20(4)`$ at $`m=0.08083(m_\pi =400(7)\mathrm{MeV})`$.

We also fitted the points around the trough ($`m=0.07583`$ to 0.22633) and then incrementally added smaller quark masses in the fit, as we did for the one-loop formula. The results are plotted in Figs. 11, 12, and 13 as a function of $`m_{\pi ,min}`$ and labeled $`\delta _{min},A_{min}`$ and $`B_{min}`$. We see again that they do not reproduce those parameters which are obtained from fitting the range from $`m=0.014`$ to 0.08983 ($`m_\pi =0.4050(38)`$).

As we are concerned about the finite volume effect (this will be addressed in Sec. V) for small quark masses, we drop the three smallest quark masses and then fit the pion mass upward. We plot the results in Figs. 11, 12, and 13 as a function of $`m_{\pi ,max}`$ and labeled $`\delta _3,A_3`$ and $`B_3`$. We see that in most cases they are only a few percent different from those ($`\delta _{max},A_{max}`$ and $`B_{max}`$) without dropping the last three quark masses. We have also performed the same test for the one-loop log formula fit and find that the parameters are also not changed much. We plot in Fig. 14 the results of $`\delta _3`$ from this fit together with $`\delta _{max}`$ as obatined from the one-loop log fit without dropping the last 3 points. We see that the characteristic transtion around $`m_\pi 300`$ MeV is still visible.

## V Finite volume effects

Intertwined in the discussion of the validity of the one-loop $`\chi `$PT for a certain mass range are the finite volume effects. There are two questions to consider: First of all, how does one know how large the zero mode contributions are, and how much do they affect the pion mass in the range of the fitting? Secondly, how large is the finite volume effect for the pion mass where the finite box is not larger than 4 times the Compton wavelength of the pion? We considered the zero mode contribution in Section III where we examined the ratio of pion masses from the $`PP`$ and $`A_4A_4`$ correlators for the $`16^3\times 28`$ lattice and found that the ratio is basically unity within errors for the range of quark masses that we are concerned with. We now look at the same ratio in a smaller lattice $`12^3\times 28`$ ($`La=2.4\mathrm{fm}`$ in this case). They are plotted in Fig. 15 together with those from the $`16^3\times 28`$ lattice. We see that the ratio for the $`12^3\times 28`$ lattice deviates from unity by two sigmas for pion mass in the range of 0.5 to 0.8, or $`500`$ to $`800\mathrm{MeV}`$. The central value of the ratio goes down to 98% at $`m_\pi 300\mathrm{MeV}`$ before it turns back to around unity for small pion mass due to the finite volume effect. This is a clear indication that the zero mode contribution to the $`PP`$ propagator is visible even for pion mass as heavy as $`800\mathrm{MeV}`$. Most of the calculations compiled in Table 3 have volumes similar to our $`12^3\times 28`$ lattice, therefore one needs to be concerned about the zero mode contamination which is a finite volume artifact. Since the zero mode contribution varies as $`\frac{1}{\sqrt{V}}`$ bcc00 ; ddh02 , the effect of the zero mode contribution will be larger for calculations with still smaller lattice sizes. It would be useful to compare the pion masses from $`PP`$ and $`A_4A_4`$ correlators to see how different they are in these cases.

To check the finite volume effect, we consider the pion mass from the $`A_4A_4`$ correlator which is less affacted by the zero mode contamination. We plot the ratio of pion masses between the $`12^3\times 28`$ and $`16^3\times 28`$ lattices in Fig. 16 as a function of $`m_\pi `$. We see that the ratio is basically unity for pion mass greater than $`300\mathrm{MeV}`$. Below this, the pion mass from the smaller volume, i.e. the $`12^3\times 28`$ lattice, becomes higher. The worst deviation, at the smallest pion mass of $`182(8)\mathrm{MeV}`$, is $`3.7\%`$. For comparison, we also plot in Fig. 16 the corresponding ratio of $`m_\pi `$ from the $`PP`$ correlator. Contrary to the $`A_4A_4`$ correlator, we see that the ratio has a dip between the $`m_\pi `$ range of 250 and $`800\mathrm{MeV}`$, again reflecting the fact that the pion mass calculated from the $`PP`$ correlator on the $`12^3\times 28`$ lattice has observable zero mode contributions.

To assess the finite volume correction, we first consider the leading finite volume correction calculated by replacing the infinite volume meson propagator with the finite volume counterpart. The correction bg92  for $`m_\pi ^2`$ is

$$\mathrm{\Delta }(m_\pi ^{1loop})^2=\frac{m_\pi ^2}{4\pi ^2f^2}(\mu ^2\alpha m_\pi ^2)\sqrt{\frac{2\pi }{m_\pi L}}e^{m_\pi L},$$
(18)
where $`\mu ^2=m_\eta ^{}^2m_\eta ^2/2m_{\pi ^0}^2/2=(871\mathrm{MeV})^2`$. With $`\alpha m_\pi ^2\mu ^2`$ bg92 , we can estimate the percentage difference of $`m_\pi ^2`$ between the $`12^3\times 28`$ and the $`16^3\times 28`$ lattices for the smallest $`m_\pi =182(8)\mathrm{MeV}`$ to be

$$\frac{\mathrm{\Delta }(m_\pi ^{1loop})^2(12)\mathrm{\Delta }(m_\pi ^{1loop})^2(16)}{m_\pi ^2}=0.13\%.$$
(19)
This is much smaller than the Monte Carlo calculation of 7.4% derived from the $`A_4A_4`$ correlators in Fig. 16. Next, we use the empirical $`L^3`$ dependence fuk92  to estimate the finite volume correction with

$$m_\pi (L)=m_\pi (L=\mathrm{})+\frac{a}{L^3}.$$
(20)
For our lowest $`m_\pi (L=16)=182(8)\mathrm{MeV}`$ where the corresponding $`m_\pi (L=12)`$ is $`3.7\%`$ higher, we obtain the finite volume correction for the $`m_\pi (L=16)=182(8)\mathrm{MeV}`$ to be $`2.7\%`$ and the correction for $`m_\pi ^2(L=16)=(182(8)\mathrm{MeV})^2`$ to be $`5.4\%`$. Similarly, for the second lowest pion mass at $`193(7)\mathrm{MeV}`$, the correction is $`3.6\%`$ for $`m_\pi ^2(L=16)`$. These are smaller than the statistical errors for the fitted low-energy parameters in Section IV.

We are not in a position to fully address the continuum limit extrapolation, as we only have results at one lattice spacing. However, we know there is no $`O(a)`$ error due to chiral symmetry knn97  and the $`O(a^2)`$ and $`O(m^2a^2)`$ errors are small liu02 ; dll00 ; ddh02 . We shall use the results from the small volume study in Ref. dll00  to gain some feeling for the magnitude of the $`O(a^2)`$ errors. We have calculated pion and rho masses on the $`6^3\times 12,8^3\times 16`$, and $`10^3\times 20`$ lattices with the Wilson gauge action at $`\beta =5.7,5.85`$, and 6.0 respectively. The figures of $`m_\pi a/\sqrt{\sigma }a`$ and $`m_\rho a/\sqrt{\sigma }a`$ ($`\sigma `$ is the string tension) vs. $`\sigma a^2`$ in Ref. dll00  show that both of them are quite flat at $`m_\pi /m_\rho =0.4,0.5`$, and 0.6. Since the lattice spacing at $`\beta =5.7`$, as measured with the Sommer $`r_0`$, is $`0.171`$ fm which is very close to our present lattice whose lattice spacing as determined from $`r_0`$ is 0.175 fm, we shall use the $`\beta =5.7`$ data as the estimate for the present lattice. After extrapolating the pion mass from $`\beta =5.7,5.85`$, and 6.0 in Ref. dll00  to the continuum limit, we find the central values of $`m_\pi `$ at $`\beta =5.7`$ are changed by 1 – 3 % for the cases with $`m_\pi /m_\rho =0.4,0.5`$, and 0.6. We thus tentatively suggest that our present pion masses are subject to the similar level of $`O(a^2)`$ errors due to the continuum limit extrapolation. This, of course, needs to be checked with calculations at smaller lattice spacings than what we have now.

## VI Chiral behavior of $`f_\pi `$ and $`f_P`$

We show in Fig. 17 our data on renormalized $`f_\pi `$ as obtained ddh02  from

$$f_\pi =\underset{t\mathrm{}}{lim}\frac{2m\sqrt{G_{PP}(\stackrel{}{p}=0,t)m_\pi }e^{m_\pi t/2}}{(m_\pi )^2}.$$
(21)
We see that if one fits a linear curve in quark mass from $`0.1085`$ ($`m_\pi =0.4438(36)`$ or $`438(7)\mathrm{MeV}`$) to $`0.4`$ ($`m_\pi =0.8740(30)`$ or $`862(13)\mathrm{MeV}`$) as is usually done in the literature, the fit gives $`f_\pi =0.0935(6)+0.110(2)m`$ with $`\chi ^2/NDF=1.0`$. In this case, the data below $`m`$ = 0.1085 shows a hump with a one-sigma deviation from the linear fit at the low mass end. Even though $`f_\pi `$ in full QCD with dynamical fermions has a leading order chiral log lp73 ; gl85  in the form $`m_\pi ^2\mathrm{log}(m_\pi ^2)`$ with a known coefficient, it is absent in the quenched theory sha92 . We interpret the bump as a consequence of the quenched chiral log in the pion mass. Thus, we plot the the $`f_\pi `$ as a function of $`m_\pi ^2`$ in Fig. 18 and fit with the following form

$$f_\pi =f_\pi (0)+C_1m_\pi ^2+C_2m_\pi ^4.$$
(22)
As in the fit of $`m_\pi ^2`$, we fix the minimum pion mass at 0.1844(71) and fit the data to the above formula with variable ranges. The results of $`f_\pi (0)`$ are plotted in Fig. 19 as a function of the maximum pion mass $`m_{\pi ,max}^2`$. We plot the same for $`C_1`$ and $`C_2/10`$ in Fig. 20. Even though there appear to be some changes in the region of $`m_{\pi ,max}`$ between $`0.3748(38)`$ and $`0.4840(32)`$ for $`f_\pi (0),C_1`$ and $`C_2`$, the errors are such that one cannot be as certain as in the case of the chiral log $`\delta `$. Thus, we shall quote the results for $`m_{\pi ,max}=0.3748(38)`$ or $`370(7)\mathrm{MeV}`$. They are given in Table 6. The fit for $`f_\pi `$ in the range up to $`m_{\pi ,max}=0.3748(38)`$ with the parameters in Table 6 is plotted in Fig. 18 as the solid lines for both the full figure and the insert. We also fit $`f_\pi `$ with only the quadratic pion mass term for the large pion mass region from $`m_\pi =0.4438(36)`$ or $`438(7)\mathrm{MeV}`$ to 0.8740(30) or $`862(13)\mathrm{MeV}`$) and obtain $`f_\pi =0.0952(6)+0.0551(13)m_\pi ^2`$. This fit is plotted as the dash line in Fig. 18.

According to the quenched chiral perturbation theory bde00 , $`C_1=\frac{4L_5}{f_\pi (0)}`$. From our data reported in Table 6, we obtain $`L_5=3.6(1.5)\times 10^3`$. This is comparable to $`L_5=2.5(5)\times 10^3`$ as obtained from the study of $`f_\pi `$ with the clover fermion action on Wilson gauge configurations at $`\beta =5.7`$. But it is substantially larger than the phenomenological value of $`1.4(5)\times 10^3`$ pic95 .

Using the experimental value $`f_\pi (m_\pi )=92.4\mathrm{MeV}`$ and $`m_\pi =137\mathrm{MeV}`$ to set the scale, we determine the lattice spacing $`a`$ to be $`0.200(3)\mathrm{fm}`$. We see that $`f_\pi (m_\pi )`$ is higher than $`f_\pi (0)`$ by $`5`$%.

The unrenormalized pseudoscalar decay constant of the pion is defined as

$$f_P^U=0|\overline{\psi }i\gamma _5(1D/2)\psi |\pi (\stackrel{}{p}=0)$$
(23)
From the chiral Ward identity, we have

$$Z_A_\mu A_\mu =2Z_mmZ_PP,$$
(24)
where $`A_\mu =\overline{\psi }i\gamma _\mu \gamma _5(1D/2)\psi `$ and $`P=\overline{\psi }i\gamma _5(1D/2)\psi `$. Since $`Z_m=Z_S^1`$ and $`Z_S=Z_P`$ due to the fact the scalar density $`S=\overline{\psi }(1D/2)\psi `$ and the pseudoscalar density $`P`$ are in the same chiral multiplet, $`Z_m`$ and $`Z_P`$ cancel in Eq. (24). Therefore, from the relation of the on-shell matrix elements

$$0|Z_A_4A_4|\pi (\stackrel{}{p}=0)=2m0|P|\pi (\stackrel{}{p}=0),$$
(25)
one obtains that

$$\frac{f_P^U}{f_\pi }=\frac{m_\pi ^2}{2m},$$
(26)
where $`f_\pi `$ is renormalized while $`f_P^U`$ and $`m`$ are unrenormalized. The relation also holds when both $`f_P`$ and $`m`$ are renormalized since $`Z_m=Z_P^1`$. It is clear from Eq. (26) that the chiral behavior of $`f_P^U`$ is the combination of those of $`m_\pi ^2/2m`$ and $`f_\pi `$.

## VII Quenched $`m_N`$

We first plot the nucleon mass as a function of $`m_\pi ^2`$ in Fig. 21. We see that it is fairly linear in $`m_\pi ^2`$ all the way down to $`m_\pi 0.4438(36)`$ ($`438(7)\mathrm{MeV}`$). Below that point, there appears to be a deviation from the linear behavior. Following the case of $`f_\pi `$, we shall first fit the $`N`$ mass only linear in $`m_\pi ^2`$ in the range of $`m_\pi ^2=0.1960.763`$ and obtain $`m_N=0.998(21)+0.932(74)m_\pi ^2`$. This fit is represented by a dashed line in Fig. 21. We see that the data (in the small $`m_\pi `$ region) in this case are systematically lower than the linear fit.

According to the quenched $`\chi `$PT for the baryons ls96 , the one-loop result predicts the presence of the non-analytic terms $`m_\pi `$ and $`m_\pi ^2\mathrm{log}(m_\pi )`$, both of which are proportional to the quenched $`\delta `$. Thus, we fit the nucleon data from the lowest pion mass, i.e. $`m_{\pi ,min}=182(8)\mathrm{MeV}`$ upward as done in previous sections with the form

$$m_N=m_N(0)+C_{1/2}m_\pi +C_1m_\pi ^2+C_{1L}m_\pi ^2\mathrm{log}(\frac{m_\pi }{\mathrm{\Lambda }_\chi })+C_{3/2}m_\pi ^3,$$
(27)
We should mention again that the fit does not depend on the value of $`\mathrm{\Lambda }_\chi `$, since there is an analytic term in $`m_\pi ^2`$. We take $`\mathrm{\Lambda }_\chi =1.16`$ for this case.

The fitted results of $`m_N(0)`$ are plotted in Fig. 22 as a function of the maximum pion mass $`m_{\pi ,max}`$ of the fitted range. We see that it is rather stable in the range fitted. Similarly, the coefficients $`C_{1/2}`$, $`C_{1L}`$, $`C_1`$ and $`C_{3/2}`$ are plotted in Fig. 23. We see that both $`C_{1/2}`$ and $`C_{1L}`$ have rather smooth behavior in $`m_{\pi ,max}`$ from $`0.3`$ to $`0.55`$ and approach zero around $`m_{\pi ,max}=0.6`$$`0.7`$. On the other hand, $`C_1`$ has more of a variation in the range $`m_{\pi ,max}=0.30.55`$. When $`m_{\pi ,max}`$ reaches the region $`0.6`$$`0.7`$, the error bars do not touch those at $`m_{\pi ,max}0.3`$. We note that $`C_{3/2}`$ is consistent with zero in the whole range of the fit. We have tried to include the next analytic term $`m_\pi ^4`$ leinweber03 . It turns out that $`m_N(0),C_{1/2}`$, and $`C_1`$ are not changed outside their errors in this case. However, $`C_{3/2}`$ and $`C_2`$ with the $`m_\pi ^4`$ term become ill-determined. They are consistent with zero but with large errors. Thus, the data seem to suggest that there is no evidence for a large $`m_\pi ^3`$ term in the quenched nucleon mass. We quote the fit evaluated at $`m_{\pi ,max}=0.3200(43)`$ ($`316(6)\mathrm{MeV}`$) in Table 7. This fit is shown as the solid curve in Fig. 24 for $`m_\pi ^2`$ up to $`0.4840(32)`$ ($`m_\pi =478(8)\mathrm{MeV}`$). Also shown in the figure as the dashed curve is the linear fit from the higher pion mass range of $`m_\pi =0.4438(36)`$ ($`438(7)\mathrm{MeV}`$) to $`0.8740(30)`$ ($`862(13)\mathrm{MeV}`$).

We should mention that even though the error is large, $`C_{1/2}`$ is consistent with that predicted by the one-loop $`\chi `$PT ls96 , for which $`C_{1/2}=\frac{3\pi }{2}(D3F)^2\delta `$. Using $`\delta =0.24(3)`$ from our fit in Sec. IV, it gives $`C_{1/2}=0.37(5)`$. This is consistent with the value in Table 7 which is $`0.358(234)`$. $`C_1=2.44(67)`$ from the low mass fit is much larger than that of the high mass fit where $`C_1=0.932(74)`$.

## VIII Conclusion

The chiral extrapolation has always been a major challenge in lattice QCD calculations, primarily due to the fact that one does not know if the formula used to fit is correct in the range of the available data. The idea of using $`\chi `$PT to guide the fit is useful and appropriate except it is complicated by the fact that one does not know the range in where the one-loop formula is valid. In this work, we attempt to answer this question by fitting the data from our smallest pion mass (close to the physical value) upward and seeing where the low-energy parameters change. This is commensurate with the $`\chi `$PT approach which expands around the small pion mass and momentum. The chiral log $`\delta `$ from fitting the pion mass to the quenched one-loop chiral formula clearly demonstrates a sharp transition around $`300400\mathrm{MeV}`$ in the pion mass. We take this as an evidence to reveal the upper limit for the applicability of the one-loop formula in quenched $`\chi `$PT. We have also fitted the pion mass with the cactus re-summed diagrams in quenched $`\chi `$PT and the results suggest that its range of applicability extends beyond the one-loop formula. We should point out that fitting with the cactus re-summed diagrams is preferred over the one-loop formula since it gives a scale-independent $`\delta `$. However, not all the physical quantities have cactus re-summed results in the $`\chi `$PT and in many cases, one would rely on one-loop formula to do the chiral fitting. Thus, it is essential to find out where the one-loop formula is valid. To be sure that this is not skewed by the finite volume effects, we compared our data on the $`16^3\times 28`$ lattice against those on the $`12^3\times 28`$ lattice and found that the finite volume errors are much smaller than the statistical errors. Although less certain due to larger errors, the analytic term in the $`N`$ mass and the non-leading chiral log in $`f_\pi `$ show variations in the pion mass range of $`300400\mathrm{MeV}`$, and the $`C_1`$ coefficients change by a factor $`23`$ compared with those from fits involving exclusively pion masses beyond this range. Thus we tentatively conclude that the one-loop quenched $`\chi `$PT is valid for pion mass less than $`300\mathrm{MeV}`$. It is claimed that using a different regularization scheme such as the finite range regularization lty03  will extend the range of validity for the one-loop quenched $`\chi `$PT from that of the dimensional regularization, but it remains to be checked by lattice Monte Carlo data at low pion mass.

One of the valuable lessons learned in this study is that it seems futile to extrapolate to the chiral region from pion mass beyond $`300\mathrm{MeV}`$. We tried to fit the pion mass with the one-loop formula by extrapolating from $`m_\pi =862(13)\mathrm{MeV}`$ and found that, no matter how small pion masses one includes, $`\delta `$ is inevitably different (by more than two sigmas) from that obtained with the data limited to $`m_\pi <250\mathrm{MeV}`$. One could try to model the functional form at higher masses, but it is not likely to succeed, because one does not know a priori how and when the modeled formula eases into that of the one-loop chiral $`\chi `$PT. Direct simulation at pion masses less than $`300\mathrm{MeV}`$ is nearly impossible (or at the very least an extremely difficult task) for lattice calculations employing fermion actions which suffer from exceptional configurations and/or critical slowing down. An overlap-fermion calculation is numerically intensive, but with approximately fixed cost – its critical slowing down is rather mild dll00 . The fact that it has chiral symmetry alleviates the worry that something may go wrong at small quark masses. Furthermore, the chiral behavior study is greatly enhanced by the incorporation of the multi-mass algorithm which facilitates the constrained curve-fitting method to discern the log structure. One does not know how different the chiral behavior in full QCD will be from that of the quenched approximation, but it is hard to imagine that there are different scales for the sea and valence quarks. If so, one will need to work harder to lower the dynamical quark masses such that the corresponding pion mass is below $`300\mathrm{MeV}`$.

## IX Acknowledgment

This work is partially supported by U.S. DOE Grants DE-FG05-84ER40154 and DE-FG02-95ER40907. We thank A. Alexandru, C. Bernard, M. Chanowitz, N. Christ, M. Golterman, A. Kennedy, D. Leinweber, S. Sharpe, A. Soni, S. Tamhankar, H. Thacker, and H. Wittig for useful discussions and suggestions. Thanks are also due to Y. Aoki, K. Holland, and C. Gattringer for providing us unpublished numbers.


format: text

text:
# Hubble Space Telescope observations of globular cluster systems along the Hubble sequence of spiral galaxies

## 1 Introduction

Recent observations with the Hubble Space Telescope (HST) and large-field ground-based CCD cameras have caused rapid advances in our knowledge of the formation and evolution of globular cluster (GC) systems of galaxies. However, ellipticals and S0 galaxies have received by far the largest amount of attention, mainly because early-type galaxies usually have much richer GC systems and also suffer less from internal extinction problems than do spirals. Consequently, our knowledge of the GC systems in spiral galaxies is still limited to a handful of galaxies (e.g., Harris 1991; Ashman & Zepf 1998; Kissler-Patig et al. 1999). However, it should be recognized that this scarcity of data for spiral galaxy GC systems constitutes an important limitation to the use of GC systems as probes of the formation and evolution of both early-type and late-type galaxies.

For instance, the ability to test predictions of scenarios for the formation of galaxies and their GCs depends on our knowledge of the ‘typical’ properties of GC systems of spirals. A good example is the ‘merger model’ (Schweizer 1987; Ashman & Zepf 1992), in which elliptical galaxies are formed by mergers of spiral galaxies. In this picture, the GC systems of elliptical galaxies are composite systems. One population of GCs is associated with the progenitor spirals (i.e., metal-poor GCs with a halo-like spatial distribution), while a second population of clusters forms in the merger event (i.e., metal-rich GCs with a bulge-like spatial distribution). Ashman and Zepf (1992; 1998) described several testable predictions arising from this scenario, most of which are however based on comparisons with properties of ‘typical’ spiral galaxy GC systems. A particular problem in this respect is that the characteristic specific frequency of GCs (i.e., the total number of GCs per unit galaxy luminosity) around spirals is poorly known, especially for late-type spirals (e.g., the compilation of Ashman & Zepf 1998 includes specific frequency measurements for only 5 spirals with Hubble type Sb or later). The absence of firm constraints on the number of GCs contributed by the progenitor spirals leads to uncertainties in many of the key predictions for elliptical galaxy GC systems that follow from the merger model. Furthermore, metal-rich ($`1.0\begin{array}{c}<\\ \end{array}\text{[Fe/H]}\begin{array}{c}<\\ \end{array}0.2`$) GCs are known to be associated with bulges of spiral galaxies like the Milky Way and M31 rather than with their discs or halos (Minniti 1995; Barmby, Holland & Huchra 2001; Forbes, Brodie & Larsen 2001), so that mergers of spirals containing significant bulges are likely to contribute metal-rich GCs to the merger remnant which are not formed during the merger. It is therefore very important to establish —using a significant sample of spiral galaxies— whether or not the properties of the Milky Way (MW) GC system are typical of spiral galaxy GC systems as a class, and to study the relationship between bulge luminosity and the number of metal-rich GCs in spiral galaxies.

Another important and timely area where knowledge of the properties of GC systems of spirals can yield significant progress is that of the formation and evolution of bulges of spiral galaxies. On the one hand, structural and dynamical properties of spiral bulges have long shown strong similarities with those of (low-luminosity) elliptical galaxies. Bulges follow the Fundamental Plane of elliptical galaxies (e.g., Bender, Burstein & Faber 1992), while their internal dynamics are consistent with oblate, isotropic models, just like low-luminosity ellipticals (e.g., Davies et al. 1983). Extinction-corrected optical and near-IR colours of bulges have been shown to be very similar to those of elliptical galaxies in the Coma cluster (Peletier et al. 1999), and spectroscopic metallicities and \[Mg/Fe\] ratios in bulges are similar to those of ellipticals at a given bulge luminosity (Goudfrooij, Gorgas & Jablonka 1999). As ellipticals in rich clusters most likely formed at redshift $`z>3`$ (e.g., Stanford, Eisenhardt & Dickinson 1997), these similarities argue for an early formation of bulges, in line with the original ‘monolithic collapse’ model by Eggen, Lynden-Bell & Sandage (1962). An alternative (and quite different) model is that bulges form from disc material through redistribution of angular momentum (Pfenniger & Norman 1990). In this scenario, large amounts of gas are driven into the central region of the galaxy by a stellar bar and trigger intense star formation. If enough mass is accreted, the bar itself will dissolve and the resulting galaxy will reveal a bigger bulge than before bar formation; galaxies would thus evolve from late to earlier types along the Hubble sequence (see also Pfenniger, Combes & Martinet 1994). Note that this scenario does not involve GCs, in that only disc stars would contribute to the secular building of bulges. Hence, there would be no reason for a relation between the bulge-to-total luminosity ratio and the number of GCs per galaxy in this scenario. This provides a testable prediction, and is one of the key tests performed in the present study.

This paper is built up as follows. Section 2 describes the sample selection and the HST observations, while the data analysis is described in section 3. The various results are discussed in section 4. Finally, section 5 summarizes the main conclusions of this study.

## 2 Approach

### 2.1 Galaxy sample selection

The sample consists of spiral galaxies in an edge-on configuration. We selected such galaxies in order to enable detection of GCs on both sides of the spiral discs and to minimize (spatially) the impact of dust absorption. In addition, the edge-on configuration allows one to assess radial number density distributions of the GC systems around these galaxies. Northern galaxies were selected from the Uppsala General Catalogue (UGC; Nilson 1973), while southern galaxies were selected from the Surface Photometry Catalogue of ESO-Uppsala Galaxies (ESO-LV; Lauberts & Valentijn 1989). Our two first-cut selection criteria were:

* the inclination of the galaxy $`i80^{}`$, and
* The Galactic latitude $`b`$ 30.

Inclinations were determined following Guthrie (1992), assuming an intrinsic flattening $`q_0(b/a)_0=0.11`$. Inclinations were then derived using Hubble’s (1926) formula

$$\mathrm{cos}^2i=(q^2q_0^2)/(1q_0^2)\text{,}$$

where $`q=b/a`$ is the observed (catalogued) axis ratio.

In order to get the highest observing efficiency out of the HST observations (see below), we further considered two main (and mutually counteracting) factors. On the one hand, a galaxy should be near enough to reach $`M_V=6.0`$, which is $``$ 1.5 magnitudes beyond the peak in the Galactic GC luminosity function (hereafter GCLF; Harris 1996), in a few HST orbits using WFPC2. On the other hand, a galaxy should be distant enough for a significant fraction of its GC system to fit in the field of view of one or two WFPC2 exposures.

From the remaining list of candidate spirals, we selected 2-3 galaxies of each main Hubble type (Sa, Sb, and Sc) to study the variation of the properties of GC systems as a function of Hubble type (i.e., of bulge-to-total luminosity ratio). Global properties of the final sample galaxies are given in Table 1.

### 2.2 Observations

Observations were made with the Wide Field and Planetary Camera 2 (WFPC2) aboard HST as part of General Observer program 6685. The data consist of multiple images through the F555W and F814W filters. Our program was supplemented by archival images of a few sample galaxies (NGC 4565 and NGC 4594), taken from other HST programs. The images of NGC 4565 were not reanalyzed; all relevant data on GC candidates in that galaxy was taken from Kissler-Patig et al. (1999). The HST observations are listed in Table 2 together with the exposure times for each galaxy. For a few galaxies, a subset of the images were spatially offset by 0$`.^{\prime \prime }`$5 from the others (corresponding to an approximately integer pixel shift in both PC and WF CCDs of the WFPC2) to enable a good correction for hot pixels. The locations of the WFPC2 fields are shown in Fig. 1, superposed onto grey-scale images from the Digital Sky Survey.

### 2.3 Data reduction and photometry

We used the iraf111iraf is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation, U.S.A. package for data reduction. After standard pipeline processing and alignment of the images, we combined the images using the IRAF/STSDAS task crrej, thereby effectively removing both cosmic rays and hot pixels. We also trimmed each image to exclude the obscured regions near the pyramid edges of WFPC2, yielding 751 $`\times `$ 751 usable pixels per CCD. Prior to performing source photometry, the strongly varying galaxy background was fitted and subtracted. The main reason for this is to minimize errors in the photometry due to any particular choice of object aperture and sky annulus. For the CCD chips covering a significant part of the galaxies’ bulges or dusty discs, we applied a median filter with a 30 $`\times `$ 30 pixel kernel to approximate the galaxy background. For the other CCDs, the smooth background gradient was fitted by a bi-cubic spline fit.

The photometry of star cluster candidates was carried out using the daophot-ii  package within iraf. The objects were selected by applying the daofind task to an image prepared by dividing the data image in question by the square root of the appropriate model image (created as described in the previous paragraph). This procedure ensures uniform shot noise characteristics over the whole image. We adopted fairly tight shape constraints \[$``$0.6 $`<`$ roundness $`<`$ 0.6; 0.2 $`<`$ sharpness $`<`$ 0.9\] in order to exclude extended background galaxies and faint objects distorted by noise or any residual bad pixels. The detection threshold was set at 4 sigma above the residual background. Although daofind returned with apparent point-like detections located within – or on the edge of – dust features in the galaxy discs, we decided to exclude those from further analysis due to the difficulty in judging the location of those sources relative to the dust features along the line of sight.

PSF photometry was performed using the ALLSTAR routine. For each chip, a PSF was constructed using, when possible, 10-15 bright, isolated GCs in the field. Aperture corrections were calculated using growth curve analysis on a subset of the PSF objects on each chip. One advantage of PSF photometry over the commonly-used two-pixel aperture photometry is that it avoids the systematic photometric offset due to underestimating aperture corrections for bright, resolved GCs.

One of the main sources of uncertainty in WFPC2 photometry (in general) is due to charge-transfer efficiency (CTE) problems of the CCDs for which correction recipes are available (Whitmore, Heyer & Casertano 1999; Dolphin 2000). We applied the Whitmore et al. formulae to our photometry. However, in the present case of GC photometry on a relatively high background (due to the diffuse galaxy light), the CTE correction typically affected the magnitudes by less than 0.02 mag and the colours by less than 0.01 mag, negligible for the purposes of this paper.

The measured WFPC2 STMAG magnitudes F547M, F555W, and F814W were converted into Johnson-Cousins $`V`$ and $`I`$ magnitudes using colour terms as derived by Holtzman et al.  and Goudfrooij et al. . Specifically, the following conversions were used:

$`V\text{F547M}`$ $`=`$ $`(0.009\pm 0.002)`$
$`+(0.005\pm 0.002)(\text{F547M}\text{F814W})`$
$`(0.008\pm 0.002)(\text{F547M}\text{F814W})^2`$
$`V\text{F555W}`$ $`=`$ $`(0.643\pm 0.001)`$
$`+(0.004\pm 0.002)(\text{F555W}\text{F814W})`$
$`(0.015\pm 0.002)(\text{F555W}\text{F814W})^2`$
$`I\text{F814W}`$ $`=`$ $`(1.266\pm 0.001)`$
$`+(0.018\pm 0.003)(\text{F555W}\text{F814W})`$
$`+(0.016\pm 0.002)(\text{F555W}\text{F814W})^2`$
Finally, the magnitudes of the point-like sources in the sample galaxies were corrected for Galactic foreground extinction (cf. Table 1). The $`A_V`$ values were converted to $`A_I`$ according to the Galactic extinction law of Rieke & Lebofsky .

### 2.4 Cluster candidate selection

To select candidate GCs, we adopted a fairly generous selection criterion for the colour cut: $`0.3VI2.0`$. The lower limit was set to 0.2 mag bluer than the bluest GC in the Milky Way, while the upper limit was set to prevent the exclusion of somewhat reddened objects (see the bottom right panel of Fig. 2). For NGC 4594, the GCs are clearly discernable on the colour-magnitude diagram (CMD; see below), so we were able to set a more stringent colour cut: $`0.5VI1.5`$. We also applied a brightness cut, eliminating GC candidates brighter than G1, the most luminous M31 GC ($`M_V=10.55`$, Rich et al. 1996).

To avoid complications related to excessive reddening, we excluded candidates lying within the dust lanes of their parent galaxies. Since this normally accounted for only a small fraction ($`5\%`$) of the total detected cluster population, this should have minimal effects on our global conclusions. However, we have attempted to correct this bias in the same manner as for our incomplete spatial coverage of the target galaxies (see below).

$`V`$ versus $`VI`$ CMDs for the candidate GCs in five of our program galaxies222For the CMDs and the photometry tables of the GCs in NGC 4565 and NGC 4594, we refer to Kissler-Patig et al. (1999) and Larsen, Forbes & Brodie (2001), respectively are shown in Fig. 2, together with that of the GC system of the Milky Way as a comparison. One rather unusual finding in this context is that the CMD of NGC 4517 reveals a relatively large number of GC candidates with $`0.3VI0.6`$, which is bluer than Galactic GCs. However, the sizes of all these candidates (which are located on the WF chips of WFPC2) are fully consistent with GCs at the distance of NGC 4517. Hence we choose to retain those sources as GC candidates. Spectroscopy will be needed to confirm the nature of these sources.

The photometry and positions of the 50 brightest GC candidates in the five galaxies for which GC photometry was not published before are listed in Tables 48 (in appendix B).

## 3 Results: Properties of the globular cluster systems

### 3.1 The total number of globular clusters

The method we used to estimate the total number of GCs in the target galaxies was to make corrections for (i) contamination by Galactic foreground stars and (ii) spatial coverage. The latter was done by comparing the target GC systems to that of the Milky Way, as detailed below.

#### 3.1.1 Correction for Galactic foreground stars

First, our estimates were corrected for foreground contamination, which was estimated from Galactic models (Bahcall & Soneira 1981). The predicted star counts at the Galactic latitude and longitude of each of the sample galaxies were subjected to the same colour and brightness cuts as the GC candidates and then subtracted from their respective bins. The total number of contaminating stars passing these criteria was small, generally $`35`$ stars per pointing.

#### 3.1.2 Correction for completeness

Artificial star experiments were performed by using the daophot-ii task addstar to estimate the completeness of the finding algorithm. These were carried out in the usual fashion, using artifical GCs generated from the PSFs constructed during the photometry measurements. These ‘fake’ GCs were added in groups of 100 in randomly placed positions on each chip. For the WF chips, the 50% completeness level was generally 1.0 – 1.5 mag beyond the turnover of the GC luminosity function (GCLF), ensuring accurate determination of the turnover point and dispersion of the GCLF. The WF-chip completeness functions are shown in Fig. 3. For each magnitude bin and each galaxy, a completeness fraction was calculated, and the number of objects in the bin was divided by this number to produce a completeness-corrected value.

#### 3.1.3 Extrapolation over the luminosity function

In preparation for the calculation of specific GC frequencies as defined by Harris & van den Bergh (1981), the total number of GCs is defined as twice the number of GCs brighter than the turnover magnitude of the GCLF, where the GCLF is assumed to be a Gaussian (in magnitude units). We therefore fitted the completeness-corrected GCLFs by a Gaussian with a standard deviation $`\sigma =1.1`$ mag (Harris 1996; Whitmore 1997) and derived the total number of clusters from the Gaussian fit.

#### 3.1.4 Correction for incomplete spatial coverage

To estimate the total number of GCs in each galaxy, we need to correct for the incomplete spatial coverage of our data. Since each pointing only covered part of the galaxy, an extrapolation over the entire galaxy was needed to estimate the total GC population. In the case of elliptical and S0 galaxies, this extrapolation is often done by evaluating the radial distribution of GCs within the observed radial extent of the galaxy. However, the low number statistics of GCs in most spiral galaxies does not allow this to be done accurately. Hence the correction for spatial coverage was done by means of a direct comparison to the positions of GCs in the Milky Way. Since this process has been described and illustrated in detail in Kissler-Patig et al. (1999), we will only summarize the process here.

We created a mask defined by our spatial coverage of each galaxy, and applied this mask to the Milky Way. By calculating the number of objects we would detect in the Milky Way mask, if placed at the distance of the target galaxy, we calculated the total number of GCs $`N_{GC}`$ as follows:

$$N_{GC}=N_{OBS}\frac{N_{MW}}{N_{MSK}}$$

where $`N_{MSK}`$ is the number of objects detected in the mask, and $`N_{MW}`$ is the total number of GCs in the Milky Way. It should be noted that there are four different orientations of the mask that would preserve the position angle of the target galaxy (these can be seen as reflections of the mask across its horizontal or vertical axes). Our final value for $`N_{MSK}`$ was the average of the values of all four possible orientations, and the standard deviation of the four values was taken into account in the uncertainty of $`N_{MSK}`$.

For this study, data regarding Milky Way GCs were taken from the McMaster catalog (Harris 1996), containing 141 GCs with appropriate ($`VI`$) colour information. However, we adopted $`N_{MW}=160\pm 20`$ (van den Bergh 1999) as the total number of GCs in the Milky Way. The undetected Milky Way GCs are assumed to be behind the Galactic bulge. Since the scale lengths and scale heights of the thin discs of the sample galaxies are very similar to those of the Milky Way (e.g., Bothun, Harris & Hesser 1992; Morrison, Boroson & Harding 1994; Baggett, Baggett & Anderson 1998; Gavazzi et al. 2000; Wu et al. 2002), the counterparts of such GCs would not be detected in our observations either. We therefore assume that the percentage of obscured GCs is the same in all sample galaxies.

We calculated a statistical uncertainty of the total number of GC candidates calculated this way by considering Poisson errors in the observed number of GCs, Poisson errors in the average number of Milky Way GCs in the mask (averaged over the four different orientations), and errors in the number of contaminating foreground objects. Our final numbers for the total population of GC candidates are listed in Table 3.

#### 3.1.5 Comparison with previous studies

A potential caveat of our method to correct for incomplete spatial coverage is the implicit assumption that the spatial distribution of the GCs in the sample galaxies is similar to that of the Milky Way globular cluster system. Hence a comparison with other methods is in order. Unfortunately, there have been only very few ground-based, wide-field studies of spiral galaxy GC systems, which is likely due to the complications discussed in the Introduction. Although this is slowly changing with the use of wide-field cameras on 4-m-class telescopes (Rhode & Zepf 2002), we can currently only make a useful comparison for one galaxy: NGC 7814, for which a previous (ground-based) study of the GC system was performed by Bothun et al. (1992). Bothun et al. found a much larger number of GC candidates than we did ($`N_{GC}=498\pm 164`$ versus our $`106\pm 28`$). From an inspection of our WFPC2 frames, we suspect that the different result is due to their counting method, which was based on object overdensities (relative to a background field) rather than colours and sizes. We find a large number of small galaxies in our WFPC2 field, which Bothun et al. (1992) likely counted along with ‘real’ GCs, while they were excluded from our GC candidate list. This once again illustrates HST’s unique power in discriminating GCs from small, nucleated galaxies.

A perhaps more useful comparison of our method with others can be done for the case of NGC 4594, for which Larsen et al. (2001) used the same HST data, but applied a different method: Larsen et al. constructed radial distribution functions of GCs to correct for incomplete spatial coverage. While that method might in principle be expected to produce more accurate results than the Milky Way mask method used here, it requires a large population of GCs to avoid complications related to small number statistics. (In our galaxy sample, this requirement is only met in the case of NGC 4594). In any case, a comparison of Larsen et al.’s results with ours is very encouraging: They find a total population of $`N_{GC}=1150\pm 575`$, compared to $`N_{GC}=1270\pm 308`$ (statistical error only) in this study. Especially heartening in this respect is that NGC 4594 (the Sombrero galaxy) is a typical Sa galaxy, which might perhaps not be expected to have a GC system particularly similar to that of the Milky Way, an Sbc galaxy.

### 3.2 Cluster specific frequencies

Traditionally, the specific frequency of GCs in galaxies, $`S_N`$, is defined as $`S_NN_{GC}\times 10^{0.4(M_V+15)}`$ (Harris & van den Bergh 1981), i.e., the total number of GCs per unit galaxy $`V`$-band luminosity (normalized to $`M_V=15`$). $`S_N`$ was introduced primarily for use in early-type (E and S0) galaxies, where there is little variation in stellar populations between individual galaxies. This is not necessarily the case for spirals, and indeed $`S_N`$ is generally known to increase along the Hubble sequence (going from late– to early-type). Averaging specific frequencies per Hubble type, the recent compilation of Ashman & Zepf (1998) supplemented by the study of two spiral galaxies by Kissler-Patig et al. (1999) yields $`S_N=0.4\pm 0.2`$ for Sc spirals to $`1.9\pm 0.5`$ for elliptical galaxies outside galaxy clusters (see Section 4). Zepf & Ashman (1993) attempted to account for differences in stellar mass-to-light ratios among galaxy Hubble types in a statistical sense by introducing a parameter $`T`$ (hereafter $`T_{\mathit{ZA93}}`$) to be the number of GCs per unit stellar mass (10$`{}_{}{}^{9}\mathrm{M}_{}^{}`$) of a galaxy. Conversion from luminosity to mass was achieved my assuming a characteristic $`M/L_V`$ value for each galaxy Hubble type.

#### 3.2.1 Total specific frequencies

The $`S_N`$ and $`T_{\mathit{ZA93}}`$ values of the target galaxies are listed in Table 3, calculated using the absolute magnitudes given in Table 1 and the number of GCs derived in Section 3.1. The quoted uncertainties were derived by taking into account the random errors discussed above as well as an uncertainty of 0.2 mag in the (total) absolute, dereddened $`V`$-band magnitudes of the galaxies. $`T_{\mathit{ZA93}}`$ values were calculated by converting the galaxy luminosities to masses following Zepf & Ashman (1993), using $`M/L_V`$ values of 5.4, 6.1, and 5.0 for Sa, Sab-Sb, and Sbc-Sc galaxies, respectively (cf. Faber & Gallagher 1979). We note that the uncertainty in the distance of a galaxy has a negligible effect on the $`S_N`$ and $`T_{\mathit{ZA93}}`$ values, since the change in the total number of GCs (due to a change in the physical area covered by the WFPC2 frames) is compensated by a similar change of the luminosity of the galaxy.

#### 3.2.2 Bulge specific frequencies

As already mentioned in the Introduction, a view is emerging that inner metal-rich GCs in spiral galaxies may be associated with their bulges rather than with their (thick) discs. Ample evidence for this is available for the Milky Way and M 31, where the metal-rich GC systems show kinematics, metallicities, and spatial distributions matching those of the underlying bulge stars (Minniti 1995; Côté 1999; Barmby et al. 2001; Perrett et al. 2002). Such a physical association is also likely present for elliptical galaxies, where the spatial distribution of metal-rich GCs typically follows that of the spheroidal galaxy light distribution very closely, whereas the metal-poor GC system typically has a more extended distribution (e.g., Ashman & Zepf 1998 and references therein). If indeed this association can be confirmed for a large sample of spirals as well, it would provide an important causal link between the formation of a spheroidal stellar system and that of metal-rich GCs. A suggestion that this might be the case was recently provided by Forbes et al. (2001) who compared the GC systems of the Milky Way (an Sbc galaxy), M 31 (Sb), and M 104 (Sa). They argued that the “bulge specific frequency (bulge $`S_N`$)” – which they defined as the number of metal-rich GCs within 2 bulge effective radii from the galaxy centres divided by the bulge luminosity (normalized to $`M_V=15`$) – was consistent among those three galaxies (bulge $`S_N0.5`$), and similar to typical values for field elliptical galaxies as well.

Structural bulge properties and bulge-to-total luminosity ratios (hereafter $`B/T`$ ratios) for the sample galaxies are listed in Table 1. For NGC 4517, NGC 4565, NGC 4594 and NGC 7814, these numbers were taken from the literature (Gavazzi et al. 2000; Wu et al. 2002; Baggett, Baggett & Anderson 1998; Bothun, Harris & Hesser 1992, respectively). The bulge/disc decompositions for the other sample galaxies are described in Appendix A. For NGC 4013 and IC 5176, we used our own WFPC2 images (in the F555W band). For NGC 3628, we used the WIYN $`V`$-band image published by Howk & Savage (1999) which was graciously made available to us. Among the sample galaxies, only the bulges of NGC 4013, NGC 4517, NGC 4594 and NGC 7814 were well fit by a de Vaucouleurs’ profile. The bulges of NGC 3628, NGC 4565 and IC 5176 were much better fit by an exponential profile. Half-light radii $`r_{1/2}`$ for the bulges of the latter galaxies were calculated as $`1.678r_0`$, where $`r_0`$ is the exponential scale length of the profile. The final half-light bulge radii $`r_{1/2}`$ of all sample galaxies are listed in the last row of Table 1.

To identify ‘metal-rich’ GCs in our sample galaxies, we considered a lower limit of \[Fe/H\] = $``$1, which is the location of the ‘dip’ in the metallicity distribution of the Milky Way GC system (Harris 1996). This is shown in Fig. 4, superposed onto the GC metallicity histograms of the sample galaxies. The metallicities were derived by transforming the $`VI`$ colours (corrected for Galactic extinction) into \[Fe/H\] values using the recent calibration of Kissler-Patig et al. (1998). The $`BI`$ colours of GCs in NGC 4565 were converted into \[Fe/H\] values using the relation for Milky Way GCs, given in Couture, Harris & Allwright (1990). Note that the \[Fe/H\] values for GC candidates in external spiral galaxies as derived from colours are likely overestimates given the possibility of reddening by dust, especially in the inner regions. Thus, the number of “metal-rich” clusters resulting from this exercise has to be regarded as an upper limit.

To allow a direct comparison with the ‘bulge $`S_N`$’ values defined in Forbes et al. (2001), we counted the total number of ‘metal-rich’ GC candidates as well as those with a projected galactocentric radius $`<2r_{1/2}`$ identified in our images. Correction for incomplete spatial coverage was done by dividing the observed number of GCs by the fraction of the total area within $`2r_{1/2}`$ sampled by our images. Since the surface density of GCs typically falls off with increasing galactocentric radius, this correction procedure was done in elliptical annuli, i.e., multiplying the number of GCs in each annulus by the ratio of total annulus area to that which was actually observed. We chose to skip this measurement for the case of NGC 3628 since the HST observations of NGC 3628 only covered a very small portion of the inner bulge, rendering the extrapolated result very uncertain. In the cases of NGC 4565 and IC 5176, we only counted GCs on one side of the disc for this purpose, namely the side that was visually not significantly impacted by dust extinction. The bulge $`S_N`$ values for the sample galaxies were then calculated using the $`B/T`$ ratios listed in Table 1. The results are listed in Table 3.

The bulge $`S_N`$ values within $`2r_{1/2}`$ for most of our galaxies are roughly consistent with a value of $``$ 0.1, with the exception of NGC 4594, the Sombrero galaxy, for which we find a value of $``$ 0.5. At face value, the low bulge $`S_N`$ values for most of our galaxies seem to be inconsistent with the prediction of Forbes et al. (2001) who suggested that bulge $`S_N`$ values within $`2r_{1/2}`$ are constant among spiral galaxies with a value around 0.5. However, systematic uncertainties of the inner bulge $`S_N`$ values as derived from our data can be substantial. E.g., we generally optimized our observing strategy to detect GCs in the halos of the target galaxies (cf. Fig. 1), and dust extinction effects are strong in the innermost bulge regions. Hence we defer further discussion on the universality of the bulge $`S_N`$ until high-resolution near-IR data are available. Our following discussion on trends of $`S_N`$ values along the Hubble sequence will therefore focus on the ‘total’ $`S_N`$ values.

## 4 Discussion

In the context of the ‘major merger’ (Ashman & Zepf 1992) or ‘multi-phase collapse’ (Forbes, Brodie & Grillmair 1997) scenarios for the formation of early-type galaxies, the higher $`S_N`$ values of ellipticals (with respect to spirals) are due to ‘extra’ GCs that were formed with high efficiency during the event that also formed the spherical stellar component (bulge). If this is indeed how bulges of spirals form in general, one would expect $`S_N`$ to vary systematically along the Hubble sequence of spirals (increasing from late-type to early-type spirals). If instead bulges were formed from disc stars by means of redistribution of angular momentum through minor perturbations and/or bar-like instabilities as in the scenario of Pfenniger & Norman (1990), one would not expect to see any significant variation along the Hubble sequence, since no GCs would be formed in this secular building of bulges.

Since our galaxy sample spans a large range of Hubble types (from Sa to Sc), we are in a good position to test this. We plot the $`S_N`$ and $`T_{\mathit{ZA93}}`$ values versus Hubble type and $`B/T`$ ratio in Figure 5. Values for the Milky Way and M31 (taken from Forbes et al. 2001 and references therein) are also included in that Figure. There are a number of remarkable aspects of these results, which we discuss in order below.

First of all, there does not seem to be any significant difference between trends derived from $`T_{\mathit{ZA93}}`$ and those from $`S_N`$. Hence, we will only discuss trends involving $`S_N`$, with the understanding that the trends involving $`T_{\mathit{ZA93}}`$ are consistent with them.

A general trend in Figure 5 is that $`S_N`$ stays basically constant at a value of $``$ 0.55 for spirals with $`B/T\begin{array}{c}<\\ \end{array}0.3`$ (roughly corresponding to Hubble types Sb and later). As these galaxies are clearly dominated in mass by their halos, this constant $`S_N`$ (and $`T_{\mathit{ZA93}}`$) suggests that the formation process of massive (late-type) spiral galaxies involves the creation of a quite constant number of GCs per unit galaxy luminosity (and mass). We suggest that this population of GCs represents a ‘universal’, old halo population that is present around each galaxy.

It is now well established that elliptical galaxies in environments similar to that of the spirals in our sample (i.e., outside galaxy clusters) have higher specific frequencies than these late-type spirals. A $`\chi ^2`$-weighted average $`S_N`$ value for ellipticals outside galaxy clusters published in the literature is $`S_N`$ = 1.9 $`\pm `$ 0.6 (using $`S_N`$ data from Ashman & Zepf 1998; Brown et al. 2000; Georgakakis, Forbes & Brodie 2001; Goudfrooij et al. 2001). If indeed the higher values of $`S_N`$ for elliptical galaxies are due to an ‘extra’ population of GCs which was formed during the formation process of spheroidal stellar systems (such as bulges) involving a rapid dissipative collapse in which the physical conditions are such that giant star clusters are formed very efficiently (i.e., as in the ‘merger’ and ‘multi-phase collapse’ scenarios mentioned above), then one would naively expect the $`S_N`$ values to scale linearly with $`B/T`$ ratio from $`S_N=0.55\pm 0.25`$ at $`B/T=0`$ to $`S_N=1.9\pm 0.6`$ at $`B/T=1`$. This area is hatched in the middle right panel of Fig. 5.

Interestingly, the $`S_N`$ values of most sample galaxies are indeed consistent with this simple prediction. There are however two galaxies for which this is clearly not the case: NGC 3628 (Sb) and NGC 7814 (Sab).

The result for NGC 3628 is quite remarkable. This Sb galaxy turns out to have a $`S_N`$ value that is significantly higher than all other Sb or Sbc galaxies with $`S_N`$ measurements to date. However, NGC 3628 is a peculiar galaxy. It is part of the Leo Triplet (Arp 317; Arp 1966). Large-area optical observations reveal a tidal plume extending $``$ 100 kpc towards the east, as well as a bridge between NGC 3628 and NGC 3627 (Burkhead & Hutter 1981; Chromey et al. 1998), both of which contain large amounts of H i gas ($`5.4\times 10^8\mathrm{M}_{}`$, which is $``$ 15% of the total H i mass observed in the main body of NGC 3628; Haynes, Giovanelli & Roberts 1979). Rots (1978) constructed three-body orbital models for the tidal interactions between NGC 3627 and NGC 3628 which reproduce the formation of the plume and bridge, yielding an age for the plume of $`8\times 10^8`$ yr since perigalacticon. This age was confirmed by photometric measurements of bright clumps in the plume by Chromey et al. (1998). It is therefore quite conceivable that the high $`S_N`$ value of NGC 3628 may be due to an ‘extra’ supply of star clusters which were formed during the gaseous interaction that also caused the tidal plume and bridge. This situation seems to have occurred in M82, another nearby starburst galaxy which has undergone a past starburst triggered by tidal interactions with a gas-rich neighbour, during which a population of GCs did indeed form (see, e.g., de Grijs, O’Connell & Gallagher 2001). As the gas from which such second-generation star clusters are formed is expected to be of relatively high metallicity (relative to halo GCs), this picture is consistent with the observation that the (total) bulge $`S_N`$ of NGC 3628 (i.e., the number of metal-rich GC candidates per unit bulge luminosity) is higher than those of the other Sb–Sc galaxies in our sample by a factor similar to the total $`S_N`$. This idea should be verified by performing wide-field imaging as well as spectroscopy of a significant number of GC candidates in NGC 3628.

As to the deviant result for NGC 7814, it is perhaps somewhat surprising to find that the two earliest-type spirals in this sample (NGC 4594 and NGC 7814) have very different specific frequencies. NGC 4594 (Sa) has $`S_N`$ values that are similar to those of field ellipticals (e.g., Ashman & Zepf 1998) and S0 galaxies (Kundu & Whitmore 2001), and thus seems to have formed GCs at a similarly high efficiency relative to its field stars. It is therefore highly unlikely that bulges like the one of NGC 4594 can be formed through the ‘secular evolution’ scenario of Pfenniger & Norman (1990). On the other hand, NGC 7814 has a much lower $`S_N`$, which is consistent with the values found for the bulk of the latest-type spirals. This difference in $`S_N`$ persists when considering the ‘bulge $`S_N`$’ values. At face value, this result suggests that the (star–) formation histories of galaxy bulges of early-type spirals can be significantly different from one galaxy to another. The formation process of some bulges (such as that of NGC 4594) may have involved dissipative collapse, being accompanied by the creation of a significant number of (metal-rich) GCs, while bulges in other spiral galaxies (like NGC 7814) may have been built from the redistribution of angular momentum through bar-like instabilities as in the ‘secular evolution’ scenario. In this context, it may well be relevant that NGC 7814 is the least luminous (and thus presumably least massive) spiral galaxy in our sample, being $``$ 5 times less luminous than NGC 4594. As Pfenniger (1991) showed, the accretion of small satellite galaxies (with a mass of $``$ 10% of that of the massive galaxy) can trigger the formation of a bar and produce a bulge through secular evolution. It is tempting to suggest that the difference in $`S_N`$ value between NGC 4594 and NGC 7814 is simply due to the fact that satellite galaxies at a satellite-to-giant galaxy mass ratio of $``$10% are more common at low (giant) galaxy luminosities than at high luminosities. Thus suggestion is certainly consistent with observed luminosity functions of galaxies in poor groups (Jerjen, Binggeli & Freeman 2000; Flint, Bolte & Mendes de Oliveira 2001).

Given the importance of this issue in the context of our understanding of the formation of the Hubble sequence of galaxies, we argue that this suggestion be tested with a statistically significant sample of (edge-on) Sa – Sab spirals which covers an appropriate range of luminosities. This can now be realized with only very modest amounts of observing time using the recent ACS camera aboard HST, which reaches some 2 magnitudes deeper than WFPC2 at a given exposure time. The current generation of large-area CCD mosaics on ground-based 4-m-class telescopes will be very valuable as well in terms of covering the whole GC populations around such galaxies in one shot.

## 5 Summary

We have studied the GC systems of 7 giant, edge-on spiral galaxies covering the Hubble types Sa to Sc, allowing us to study the variation of the properties of GC systems along the Hubble sequence. Our high-resolution HST/WFPC2 data (supplemented by archival WFPC2 data for two galaxies) reached $``$ 1.5 mag beyond the turn-over magnitude of the GC luminosity function for each galaxy.

Specific frequencies of GCs ($`S_N`$ values) were evaluated by comparing the numbers of GCs found in our WFPC2 pointings with the number of GCs in our Milky Way which would be detected in the same spatial region if placed at the distance of the target galaxies. Results from this method were compared with the more commonly used method of constructing radial distribution functions of GCs to correct for incomplete spatial coverage in the case of NGC 4594 (where both methods were possible), and found to be consistent with one another to within 1$`\sigma `$.

The $`S_N`$ values of spirals with $`B/T\begin{array}{c}<\\ \end{array}0.3`$ (i.e., spirals with a Hubble type later than about Sb) are consistent with a value of $`S_N=0.55\pm 0.25`$. We suggest that this population of GCs represents a ‘universal’, old halo population that is present around each galaxy.

Most galaxies in our sample have $`S_N`$ values that are consistent with a scenario in which GC systems are made up of (i) the aforementioned halo population plus (ii) a population that is associated with bulges of spirals, which grows linearly with the luminosity (and mass) of the bulge. Such scenarios include the ‘merger scenario’ for the formation of elliptical galaxies as well as the ‘multi-phase collapse’ scenario, but it seems inconsistent with the ‘secular evolution’ scenario of Pfenniger & Norman (1990), in which bulges are formed from disc stars by means of the redistribution of angular momentum through minor perturbations and/or bar instabilities.

On the other hand, the bulge-dominated spiral NGC 7814 shows a low $`S_N`$ value, consistent with those of the latest-type spirals. NGC 7814 is the least luminous galaxy in our sample. Based on observed luminosity functions of galaxies in poor groups, we suggest that the ‘secular evolution’ scenario to build bulges in early-type spirals is most viable for low-luminosity spirals.

Thus, our results suggest that the formation histories of galaxy bulges of early-type spirals can be significantly different from one galaxy to another. Given the importance of our understanding of the formation of the Hubble sequence of galaxies in the context of galaxy evolution, we argue that the GC systems of a statistically significant sample of luminous, edge-on Sa – Sab spirals be studied in the near future.

##### Acknowledgments.

This paper is based on observations obtained with the NASA/ESA Hubble Space Telescope, which is operated by AURA, Inc., under NASA contract NAS 5–26555. We are grateful to Chris Howk for making his WIYN images of our Northern targets available to us, and thank David Clements for his contribution during the proposal preparation stage. We thank the referee, Duncan Forbes, for a fast, thorough and constructive review. We have made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with the National Aeronautics and Space Administration. PG was affiliated with the Astrophysics Division of the Space Science Department of the European Space Agency during part of this project. DM is supported by FONDAP Center for Astrophysics 15010003. JS and LB would like to thank the Space Telescope Science Institute for financial support through its Summer Student Program.

## Appendix A Bulge/disc decompositions

Since bulge/disc decompositions are not the main topic of this paper, we will only provide a brief description of our method here. A more detailed description will be provided elsewhere. We first determined the scale heights of the discs at 4 independent positions along their major axes, outside the region where the bulge contribution dominates. In general, the positions of the galaxy planes were determined by assuming symmetrical light distributions with respect to the disc plane, and the vertical profiles on either side of the plane were averaged together. However, the latter was not done for the bulge fits on NGC 4013 and IC 5176, since the WFPC2 frames only covered a large enough spatial extent on one side of the galaxy centre.

The vertical brightness distributions of the discs were well fitted by an exponential profile

$$\mathrm{\Sigma }_\mathrm{D}(z)=\mathrm{\Sigma }_{\mathrm{D},0}\mathrm{exp}(z/z_0)\text{,}$$
(1)
where $`z_0`$ is the disc scale height. We found that as long as one stays outside the central areas where there is an obvious contribution from bulge light, the disc scale heights generally do not exhibit any significant variation with distance along the major axis. This is consistent with the findings of de Grijs, Peletier & van der Kruit (1997). The average disc scale heights are listed in Table 1.

To obtain structural parameters for the bulges of these galaxies, we used a least-square program that fits a two-component (disc+bulge) model to the minor-axis brightness profiles. The disc component was fit by an exponential profile and the bulge component was fitted by two different functions, namely (i) a de Vaucouleurs’ (1953) profile

$$\mathrm{\Sigma }_\mathrm{B}(z)=\mathrm{\Sigma }_{\mathrm{B},\mathrm{e}}\mathrm{exp}(7.688((z/z_\mathrm{e})^{1/4}1))\text{,}$$
(2)
where $`z_\mathrm{e}`$ is the effective radius along the minor axis and $`\mathrm{\Sigma }_{\mathrm{B},\mathrm{e}}`$ is the intensity at $`z=z_\mathrm{e}`$, and (ii) a (second) exponential profile. For each galaxy, the disc scale height $`z_0`$ as found above was fixed in the fit. Areas where dust extinction is prominent (such as the central dust lanes) were flagged and ignored during the fitting process. The choice between the de Vaucouleurs’ profile and the (second) exponential profile was dictated by the $`\chi ^2`$ value of the fit.

Among the sample galaxies, only the bulges of NGC 4013, NGC 4517, NGC 4594 and NGC 7814 were well fit by a de Vaucouleurs’ profile. The bulges of the other galaxies were better fit by an exponential profile. The effective radii (or exponential scale radii rexp) of the bulges are listed in Table 1, after converting them to equivalent radii $`r\sqrt{ab}=b/\sqrt{1ϵ}`$, where $`a`$ and $`b`$ are the semimajor and semiminor axes of the ellipse, and $`ϵ`$ its ellipticity. For NGC 3628, NGC 4013, NGC 4517 and IC 5176, the bulge ellipticities were determined by using the ellipse fitting program in the isophote package within IRAF. The dusty areas of the galaxies were flagged and ignored in the ellipse fitting process. For the galaxies with exponential bulge profiles, we also list the bulge half-light radii (= 1.678 rexp) for comparison with the effective radii of the bulges with de Vaucouleurs’ profiles.

## Appendix B Tables of globular cluster candidates

The photometry and astrometry of the 50 brightest GC candidates on the WFPC2 images of the galaxies for which GC photometry was not published before are given in Tables 48. For the photometry of GC candidates in NGC 4565 and NGC 4594 we refer to Kissler-Patig et al. (1999) and Larsen et al. (2001), respectively.