# Determination of the chiral condensate from QCD Dirac spectrum on the lattice

###### Abstract

We calculate the chiral condensate of QCD with 2, 2+1 and 3 flavors of sea quarks. Lattice QCD simulations are performed employing dynamical overlap fermions with up and down quark masses covering a range between 3 and 100 MeV. On –1.9 fm lattices at a lattice spacing 0.11 fm, we calculate the eigenvalue spectrum of the overlap-Dirac operator. By matching the lattice data with the analytical prediction from chiral perturbation theory at the next-to-leading order, the chiral condensate in the massless limit of up and down quarks is determined.

###### pacs:

11.15.Ha,11.30.Rd,12.38.Gc^{†}

^{†}preprint: OU-HET-689-2010

^{†}

^{†}preprint: UTHEP-617

^{†}

^{†}preprint: KEK-CP-249

JLQCD and TWQCD collaborations

## I Introduction

Chiral condensate in Quantum Chromodynamics (QCD) is not a directly accessible quantity in experiment, yet plays a crucial role in the low-energy dynamics of QCD as an order parameter of chiral symmetry breaking. When in the limit of massless quarks is nonzero, chiral symmetry is spontaneously broken and hadrons acquire a mass of order , the QCD scale. Only the pion remains massless as the Nambu-Goldstone boson; its dynamics is well described by an effective theory known as chiral perturbation theory (ChPT) Weinberg:1968de ; Gasser:1983yg . is one of two most fundamental parameters in ChPT (the other is the pion decay constant ) and appears in the predictions of various physical quantities.

Calculation of , the expectation value of the scalar density operator , from the first principles of QCD requires nonperturbative techniques. In this paper we report on a numerical calculation of in lattice QCD including the effects of up, down and strange sea quarks. We investigate the low-lying eigenvalue spectrum of the Dirac operator, which is related to through the Banks-Casher relation Banks:1979yr and its extension to the case of nonzero . Since only the low-lying eigenvalues are relevant, one can avoid contamination from ultraviolet divergence of the scalar density operator , which is of order at a finite quark mass and a lattice spacing .

The Banks-Casher relation is satisfied in the limit of massless quarks after taking the large volume limit (the thermodynamical limit), which is the meaning of the arrow in the relation . When the massless limit is taken at a finite volume, the vacuum expectation value of shows a critical fluctuation which leads to a divergent correlation length and vanishing chiral condensate. Taking the thermodynamical limit, on the other hand, is numerically expensive in practical lattice calculations. In this study, we use the low energy effective theory as a guidance of volume and quark mass scalings. Namely, once the scaling behavior predicted by the effective theory is confirmed by the lattice data, the infinite volume and chiral limit according to the scaling can be safely taken. The scaling we consider is that for varying the eigenvalue , volume and the quark mass .

We use the ChPT formula for the low-lying eigenvalue spectrum of the Dirac operator Damgaard:2008zs , which is valid in both the and regimes. The regime is a parameter region where the quark mass is so small that the Compton wavelength of the pion is longer than the extent of the finite-volume space-time. In this regime, the constant mode of the pion field has to be integrated out over the group manifold when the path integral is evaluated, which is a nonperturbative calculation in the sense that one does not use the expansion in small pion field. At the leading order of the so-called expansion, the system is equivalently described by the chiral Random Matrix Theory (ChRMT), with which a number of theoretical predictions for the low-lying Dirac spectrum have been derived Damgaard:1997ye ; Wilke:1997gf ; Akemann:1998ta ; Damgaard:2000ah . In the other regime ( regime), where the pion Compton wavelength fits in the volume, the conventional ChPT applies and the calculation of the Dirac operator spectrum is available as well Smilga:1993in ; Osborn:1998qb . The new method given in Damgaard:2008zs consistently combines the both results within a systematic expansion, and thus is valid in both regimes as well as in between. By lattice calculations, we produce the data at various sets of quark masses to firmly test this analytic expectation.

Application of the ChPT formula for the low-lying Dirac spectrum requires a good control of the chiral symmetry in the calculation. In this work, we use the overlap-Dirac operator Neuberger:1997fp ; Neuberger:1998wv which satisfies the Ginsparg-Wilson relation Ginsparg:1981bj and thus realizes a modified chiral symmetry on the lattice Luscher:1998pqa . Since the chiral effective theory is constructed only assuming the presence of chiral symmetry, the same construction as in the continuum theory can be applied to lattice QCD with overlap fermions. Although the overlap-Dirac eigenvalues lie on a circle on the complex plane, the physical imaginary part is uniquely identified up to discretization effects.

In our previous works Fukaya:2007fb ; Fukaya:2007yv we performed large-scale lattice simulations of two-flavor QCD using the overlap fermion formulation Aoki:2008tq and calculated low-lying Dirac eigenvalues. By matching the lowest eigenvalue spectrum with the ChRMT expectations, we extracted . Since the ChRMT corresponds to the leading order of the expansion in ChPT, this result is subject to NLO or corrections, which could be sizable on the lattice of size used in that study. Another limitation was that the quark mass must be very small to apply the regime formula, and the runs in the regime could not be used in the analysis.

Application of the new formula Damgaard:2008zs is attempted for the first time in our recent work Fukaya:2009fh ; Fukaya:2010ih in which next-to-leading order (NLO) corrections are included in the analysis. Using the 2+1-flavor QCD data generated with dynamical overlap fermions on lattices, the value of in the chiral limit of two light quarks is obtained. The present paper provides an extensive description of this work.

In this paper, we analyze the low-lying Dirac spectrum mainly on the 2+1-flavor QCD simulations, where the strange quark mass is fixed near its physical value. The mass of degenerate up and down quarks, , covers a range between and in the regime lattice ensembles, which are generated for calculations of various physical quantities including the pion mass and decay constant JLQCD:2009sk . We also generate an regime ensemble, where the up and down quarks are kept nearly massless while is fixed near the physical value. The NLO ChPT formula allows us to combine these data to obtain in the chiral limit of up and down quarks. We also extract the pion decay constant and one of the NLO low energy constants (LECs) from the correction terms.

As demonstrated in the following sections, the ChPT formula provides information on the shape of the low-lying Dirac spectrum, with which we can test the agreement between the formula and the lattice data in detail. The volume dependence gives a critical test, since it is essentially controlled by and . At some parameter points, we compare the data obtained on a larger () lattice to those on the smaller () volume to check if the ChPT prediction consistently describes the difference. The sea quark mass dependence of the chiral condensate is partly controlled by but we also expect a nonanalytic dependence due to the pion-loop effects that should be observed in the lattice data. These nontrivial consistency checks are performed to gain confidence in the final result for .

In addition to the main analysis in 2+1-flavor QCD, we also investigate two-flavor QCD and three-flavor QCD where and are degenerate. For the case of two-flavor QCD, the lattice data are the same as in our previous studies Fukaya:2007fb ; Fukaya:2007yv , but we use the ChPT formula valid at NLO in this new analysis. The degenerate three-flavor QCD configurations are newly generated for this study at two light quark masses. We thus obtain the chiral condensate for these variants of QCD.

This paper is organized as follows. In Section II, we briefly explain how the chiral condensate is determined from the Dirac eigenvalue density and how the finite volume effects are removed using ChPT. Some technical aspects of lattice QCD simulations are described in Section III and the numerical results are discussed in the following sections. First, we describe the strategy to extract the LECs from lattice data of the spectral density in Section IV. Second, numerical scaling tests of the NLO ChPT formula are presented in Section V. We then proceed to the determination of in the chiral limit in Section VI. Summary of this analysis and conclusions are given in Section VII.

## Ii Banks-Casher relation in a finite volume

Eigenvalues of the Dirac operator in four-dimensional continuum Euclidean space are pure imaginary (which we denote , , with real ’s). The spectral density at an eigenvalue is defined by

(1) |

where denotes an average over the gauge configuration space. Since a nonzero eigenvalue appears as a pair with its complex conjugate, , we only consider the region in the following.

The chiral condensate in the limit of massless quarks and infinite volume is an order parameter of the chiral symmetry breaking. Through the Banks-Casher relation Banks:1979yr , is related to as

(2) |

with which one can identify the spontaneous breaking of chiral symmetry by measuring instead of .

Even when the volume , sea quark masses and are all finite, a similar nonperturbative relation

(3) |

holds. Here, means that the quantity is evaluated with the valence quark mass analytically continued to a pure imaginary value . In this relation, the ultraviolet divergence in the definition of the operator cancels by taking its real part at an imaginary value (where the divergent part is pure imaginary), which is natural because the left hand side of the equation only refers to low-lying modes and is insensitive to the ultraviolet region of the dynamics.

We note that the relation (3) is valid even when the ensemble average is restricted to a given topological sector of in the gauge field configurations Fukaya:2006vs , that we denote . Namely, if we define the spectral density at a given topological sector as , it is obtained by computing .

Although the chiral condensate is different from at finite volumes, the difference can be described by the low-energy effective theory, provided that the energy scales of the theory are well below the QCD scale:

(4) |

It is, therefore, possible to directly compare the lattice QCD calculation of at finite , , and with the prediction of the effective theory. By taking the limit of after according to the effective theory and summing over , one can reproduce the physical , which in the limit of gives .

In this direction, studies in both lattice QCD DeGrand:2006nv ; Lang:2006ab ; Hasenfratz:2007yj ; DeGrand:2007tm ; Bar:2010zj ; Jansen:2009tt ; Hasenfratz:2008ce and the low-energy effective theory have been done. Smilga and Stern Smilga:1993in and Osborn et al. Osborn:1998qb calculated the Dirac eigenvalue spectrum in the conventional expansion of (partially quenched) ChPT to NLO. In the vicinity of , which corresponds to the limit of zero valence pion mass, a special treatment of the zero-momentum modes is needed because the correlation length exceeds the size of the volume. This special treatment is known as the expansion of ChPT, in which the zero-momentum modes are nonperturbatively integrated over the (or when the topological charge is fixed) manifold. At the leading order (LO), the finite size effect around was calculated in Damgaard:1997ye ; Wilke:1997gf ; Akemann:1998ta ; Damgaard:2000ah . Their results are expressed using the Bessel functions, which has a gap from zero, reflecting the fact that no spontaneous symmetry breaking occurs at finite volumes.

Recently, an interpolation between the and regimes was considered in Damgaard:2008zs . The recipe for the calculation is to keep the same counting rule as in the expansion but to integrate the zero-modes nonperturbatively like in the expansion. Partial quenching is performed with the so-called replica trick so that results at arbitrary nondegenerate set of quark masses can be compared to lattice QCD. Using this hybrid method, the results mentioned above (in the expansion Smilga:1993in ; Osborn:1998qb and in the expansion Damgaard:1997ye ; Wilke:1997gf ; Akemann:1998ta ; Damgaard:2000ah ) are smoothly connected. Comparison with the lattice data is, therefore, no longer limited in either the or regimes, and more precise determination of is possible.

Here we briefly reproduce the result of Damgaard:2008zs where we consider a general theory with flavors of sea quarks. The spectral density in a fixed topological sector of is given by

(5) |

where two terms and are given in the following. includes the leading finite quark mass correction to that modifies the overall normalization of the spectrum, and is therefore called the effective chiral condensate.

The spectrum of the near-zero modes () is mainly affected by the zero-momentum pion fields. The first term in (5) has the same functional form as the one at the leading order of the expansion Damgaard:1997ye ; Wilke:1997gf ; Akemann:1998ta ; Damgaard:2000ah , which is expressed in terms of dimensionless combinations and :

(6) |

with matrix and matrix defined by and , , , respectively (’s and ’s denote the (modified) Bessel functions.). The phase factor is 1 for and 3.

The second term in (5) is a logarithmic
NLO correction (chiral-logarithms) which is also partly seen
in the conventional expansion Osborn:1998qb .
With the meson mass , which is
made of either sea quark () or valence quark (), it is
given by^{1}^{1}1In this paper, we use simplified notations:
and correspond to and
in Damgaard:2008zs , respectively.

(7) |

where

(8) | |||||

(9) |

Here, the physical meson masses are given by the leading order relations , and . The scale (= 770 MeV in this work) is a subtraction scale. The function given by denotes a finite volume correction from nonzero momentum pion modes. In the expansion, it is expressed by the modified Bessel function Bernard:2001yj while in the expansion a polynomial expression is used Hasenfratz:1989pk . In this study, we need the both expressions:

(10) |

where denotes a four-vector with integer ’s and ’s are the shape coefficients defined in Hasenfratz:1989pk . Their formula and numerical values for the first several ’s are summarized in Appendix A. In our numerical study, we truncate the sum at and , which indeed shows a good convergence around the threshold . We note that both and are finite even in the limit of .

The effective chiral condensate in (5) is given by

(11) |

where (renormalized at ) is one of the low-energy constants at NLO Gasser:1983yg . From the sea quark mass dependence of , one can determine as well as and .

In the expression (5), dependence on the topological charge is encoded only in the first term and the second term does not depend on since it is a contribution from nonzero momentum modes. On the other hand, the chiral logarithm manifests itself in the both terms through . Since could also contain through , the spectral density shows a nonanalytic functional form.

The above ChPT results are subject to higher order corrections in the expansion, for which the expansion parameter is either or . Although the zero-mode contribution is treated nonperturbatively, there are two-loop contribution of nonzero momentum modes that could also couple to the zero-mode and introduce different types of group integrals. At the two-loop level, these contributions may have the order , or , whose coefficients are unknown. We therefore need to carefully check the convergence of the expansion at NLO for our parameter sets. In the following analysis, we test the NLO formula with various sets of quark masses, as well as different lattice volumes ( = 16, 24) for the runs and different topological sectors for the runs, in order to confirm the convergence.

## Iii Lattice QCD simulations

Numerical simulations of lattice QCD are performed with the Iwasaki gauge action Iwasaki:1985we at (except for the run of QCD at for which we choose ) including 2, 2+1 ( fixed), and 3 () flavors of dynamical quarks. For the quark action, we employ the overlap-Dirac operator Neuberger:1997fp

(12) |

where denotes the quark mass and is the Hermitian Wilson-Dirac operator with a large negative mass . We take throughout our simulations (here and in the following the parameters are given in the lattice unit.). For the details of numerical implementation of the overlap-Dirac operator, we refer to our previous paper Aoki:2008tq .

It is known that the numerical cost for the dynamical simulation of the overlap fermions becomes prohibitively large when has (near) zero-modes. To avoid this problem, we introduce extra Wilson fermions and associated twisted mass bosonic spinors to generate a weight

(13) |

in the functional integrals Izubuchi:2002pq ; Vranas:2006zk ; Fukaya:2006vs . Both of these fermions and ghosts are unphysical as their masses are of order of the lattice cutoff, and do not affect low-energy physics. The numerator suppresses the appearance of near-zero modes, while the denominator cancels unwanted effects from high modes. The twisted-mass parameter controls the value of threshold below which the eigenmodes are suppressed. In our numerical studies, we set = 0.2.

With the determinant (13) the index of the overlap-Dirac operator, or the topological charge in the continuum limit Hasenfratz:1998ri , never changes from its initial value during the molecular dynamics steps since its change always requires crossing zero eigenvalue of . In this work the simulations are mainly performed in the trivial topological sector, . In order to check the topological charge dependence, we also carry out independent simulations at , and at several sets of parameters.

(GeV) | |||||||||

2 | 2.35 | 1.776(38) | 0.002 | 0 | 4680 | 0.5 | 34(12) | ||

2.30 | 1.667(17) | 0.015 | 0 | 10000 | 0.5 | 48(21) | |||

0.025 | 0 | 10000 | 0.5 | 38(16) | |||||

0.035 | 0 | 10000 | 0.5 | 28(12) | |||||

0.050 | 0 | 10000 | 0.5 | 24(9) | |||||

0.050 | -2 | 5000 | 0.5 | 50(24) | |||||

0.050 | -4 | 5000 | 0.5 | 34(16) | |||||

0.070 | 0 | 10000 | 0.5 | 23(10) | |||||

0.100 | 0 | 10000 | 0.5 | 9(3) | |||||

2+1 | 2.30 | 1.759(10) | 0.002 | 0.080 | 0 | 5000 | 0.5 | 17(9) | |

0.015 | 0.080 | 0 | 2500 | 1.0 | 15(6) | ||||

0.015 | 0.080 | 1 | 1800 | 1.0 | 5(1) | ||||

0.025 | 0.080 | 0 | 2500 | 1.0 | 11(5) | ||||

0.035 | 0.080 | 0 | 2500 | 1.0 | 24(11) | ||||

0.050 | 0.080 | 0 | 2500 | 1.0 | 9(5) | ||||

0.015 | 0.100 | 0 | 2500 | 1.0 | 5(3) | ||||

0.025 | 0.100 | 0 | 2500 | 1.0 | 10(3) | ||||

0.035 | 0.100 | 0 | 2500 | 1.0 | 34(21) | ||||

0.050 | 0.100 | 0 | 2500 | 1.0 | 5(3) | ||||

2.30 | 1.759(10) | 0.015 | 0.080 | 0 | 2500 | 1.0 | 2(1) | ||

0.025 | 0.080 | 0 | 2500 | 1.0 | 3(2) | ||||

3 | 2.30 | 1.759(10) | 0.025 | 0.025 | 0 | 2500 | 1.0 | 4(1) | |

0.035 | 0.035 | 0 | 2500 | 1.0 | 5(1) | ||||

0.080 | 0.080 | 0 | 2500 | 1.0 | 20(12) | ||||

0.100 | 0.100 | 0 | 2500 | 1.0 | 24(15) |

Simulation parameters are summarized in Table 1. The lattice size is for the runs, while it is for the main and runs. In order to investigate the finite volume scaling, we also simulate on a lattice at the same lattice spacing with two choices of sea quark masses and . For the determination of the lattice spacing = 0.11–0.12 fm ( 1.7 – 1.8 GeV), we choose the -baryon mass as the input for the and 3 ensembles JLQCD:prep , while it is determined from the heavy quark potential with an input = 0.49 fm for the case. Our lattice size is then estimated as fm for the runs, and fm for the runs.

In the runs, seven different values of the up and down quark mass are taken. For the runs, we choose two different values of strange quark mass (= 0.080 and 0.100) and six (for the former) or five (for the latter) values of are chosen. For the degenerate flavor runs, we take four different values of . Note that the lightest up and down quark mass in the or runs roughly corresponds to 3 MeV in the physical unit, with which pions are in the regime while kaons still remain in the regime.

We compute 50–80 pairs of low-lying eigenvalues (that we denote ’s) of the massless overlap-Dirac operator at every 5 or 10 (depending on the parameters) trajectories. In the calculation of the eigenvalues, we employ the implicitly restarted Lanczos algorithm for the chirally projected operator (where ) of which eigenvalue corresponds to . From each eigenvalue of , the eigenvalue , as well as its complex conjugate, are extracted through the relation . In order to compare with ChPT, every complex eigenvalue is mapped onto the imaginary axis as . The difference between and is a discretization effect, which is negligible (within 1%) for 0.03. In the analysis, we consider positive only.

For each run, 1,800–10,000 (depending on the parameters) trajectories are accumulated using the hybrid Monte Carlo algorithm. The integrated auto-correlation time of the lowest in the unit of the trajectory length is also listed in Table 1. Because of its infra-red nature, the lowest Dirac eigenvalue is expected to be most difficult to decorrelate and thus has the longest auto-correlation time. The measurement is not stable and the statistical error is as large as 50% in some cases, but is typically or less. In the following analysis, the statistical error for the spectral density and other quantities is estimated by the jackknife method after binning the data in every 100 trajectories.

Details of configuration generation and other quantities will be reported in a separate paper.

## Iv Extraction of LECs at each set of sea quark masses

Although a global fit of the lattice data for spectral density to the NLO ChPT formula is possible in principle, we prefer to simplify the analysis, for better understanding of numerical sensitivity of the lattice data and the errors in the final result. We first consider the mode number below , or the integrated eigenvalue density,

(14) |

at each set of sea quark masses. The analytic ChPT result (5) is also integrated numerically from 0 to . In the second term of (5) we can replace by as their difference is a higher order effect. Then, there are two unknown parameters in the formula: and .

The data points of at two reference values of are hence sufficient to determine the parameters. As the reference points we take = 0.004 ( 7 MeV) and 0.017 ( 30 MeV) except for the case with = 0.002, for which we choose = 0.0125 and 0.017 () or 0.010 and 0.017 (). For the runs we take = 0.01 and 0.02. Effectively, the lower point determines , while the other point is more sensitive to the NLO effects that contain . We check that the resulting values of and are stable against the change of the reference points by varying them by a factor of 2 or 3 while keeping the higher point less than 0.030 to avoid possible higher order corrections. The reference points are different for the regime runs and for the runs, because the small eigenvalues are highly suppressed (the lowest eigenvalue is larger than 0.004) for these cases.

For the = 2+1 lattice data, we test both the and ChPT formulas. For the latter case, the strange quark is assumed to be decoupled from the theory, which we call reduced ChPT.

Numerical results are listed in Table 2. Before moving to further analysis of the results let us describe our observations for the spectral function.

Figures 1–8 show the spectral density and its integral (14) obtained in our lattice simulations. A typical example is that for QCD in the regime (Figure 1). The upper panel shows the histogram plot of the spectral function as a function of . The lattice data for each bin have a jackknife estimate of the statistical error. The solid (red) curve represents the NLO ChPT formula, while the dotted (blue) curve corresponds to the leading order result (in the expansion). Since the first reference point is 0.004, it probes the first peak, which corresponds to the lowest eigenvalue. Starting from the second peak, the effect of the NLO term appears, as clearly seen from the difference of the two curves. Therefore, by taking the second reference point at 0.017, the data have enough sensitivity to the NLO parameter . This observation is of course specific to the particular volume of our lattice; on larger lattices, the peaks move toward the origin and the impact of the NLO term would become less significant on the second or third peaks (see Figure 8).

The agreement with the formula can be seen more clearly by looking at , the mode number below (lower panel). The lattice data depart from the leading-order curve (dotted), which corresponds to the first term in (5), at around 0.005 (see the inset). Then, the data follow the nontrivial functional form of the NLO formula (solid), which comes from the chiral logarithm originating from pion loops. The NLO formula works precisely up to , and the deviation is still within two sigma at , which is about the half of the physical strange quark mass . This is a typical range where the NLO ChPT is valid, and the higher order corrections would become sizable above this value. In contrast to a recent work by Giusti and Lüscher Giusti:2008vb , where they take a wider range of (up to 95 MeV) into the analysis, we conservatively choose the reference points below 30 MeV, so that (partially quenched) ChPT with an imaginary valence quark mass can be safely applied.

In Figure 1, the difference between the and (dashed curve) formulas is not sizable below 0.03–0.04. This is natural because the strange quark (with ) decouples from the dynamics of the low-lying modes. As a result, the extraction of does not significantly depend on the formula we use ( or ).

The convergence of the chiral expansion is better in the regime as shown in Figure 2, in which the lattice data at = 0.002 and = 0.080 are plotted. In the plot of the mode number (lower panel), the LO and NLO curves coincide up to 0.025. Beyond this value, we observe some deviation, which is also seen in the histogram plot (upper panel).

The NLO ChPT correction in this work explains the disagreement of the lattice data with the expectation from the random matrix theory found in our previous work Fukaya:2007fb ; Fukaya:2007yv . Namely, if we adjust the parameter using the lowest eigenvalue distribution (the first peak of the histogram), then the second peak would not agree at the leading order. Indeed, the NLO contribution is responsible for this.

(lattice) | =3 ChPT | (reduced) =2 ChPT | comment | ||||

2 | 0.002 | – | – | 0.00218(19) | 0.059(65) | (=2.35) | |

0.015 | – | – | 0.00362(15) | 0.0527(20) | |||

0.025 | – | – | 0.00353(15) | 0.0664(90) | |||

0.035 | – | – | 0.00382(14) | 0.0681(64) | |||

0.050 | – | – | 0.00449(15) | 0.0644(20) | |||

0.050 | – | – | 0.00400(16) | 0.0728(60) | (=-2) | ||

0.050 | – | – | 0.00482(19) | 0.0636(19) | (=-4) | ||

0.070 | – | – | 0.00480(15) | 0.0707(23) | |||

0.100 | – | – | 0.00478(12) | 0.0862(72) | |||

2+1 | 0.002 | 0.080 | 0.00204(07) | 0.0469(102) | 0.00204(05) | 0.0425(49) | |

0.015 | 0.080 | 0.00314(18) | 0.0536(15) | 0.00305(17) | 0.0551(16) | ||

0.015 | 0.080 | 0.00354(48) | 0.0521(25) | 0.00319(58) | 0.0558(62) | () | |

0.015 | 0.080 | 0.00273(06) | 0.0520(25) | 0.00270(06) | 0.0545(26) | (=24) | |

0.025 | 0.080 | 0.00333(18) | 0.0624(20) | 0.00326(18) | 0.0647(20) | ||

0.025 | 0.080 | 0.00299(06) | 0.0600(23) | 0.00297(05) | 0.0629(24) | (=24) | |

0.035 | 0.080 | 0.00404(39) | 0.0636(17) | 0.00393(36) | 0.0666(16) | ||

0.050 | 0.080 | 0.00423(22) | 0.0696(16) | 0.00413(21) | 0.0738(16) | ||

0.015 | 0.100 | 0.00309(14) | 0.0564(19) | 0.00303(13) | 0.0578(19) | ||

0.025 | 0.100 | 0.00349(20) | 0.0622(17) | 0.00342(19) | 0.0642(17) | ||

0.035 | 0.100 | 0.00418(40) | 0.0647(14) | 0.00409(38) | 0.0673(14) | ||

0.050 | 0.100 | 0.00383(13) | 0.0713(16) | 0.00376(13) | 0.0747(16) | ||

3 | 0.025 | 0.025 | 0.00335(23) | 0.0531(10) | – | – | |

0.035 | 0.035 | 0.00334(21) | 0.0612(24) | – | – | ||

0.080 | 0.080 | 0.00453(23) | 0.0767(14) | – | – | ||

0.100 | 0.100 | 0.00520(22) | 0.0835(22) | – | – |

## V Scaling tests with various parameters

### v.1 Sea quark masses

Figures 3 and 4 compare the spectral density and its integral at various sea quark masses. The data are obtained for (Figure 3) and (Figure 4) with up and down quark masses in the regime ( = 0.050, 0.035, 0.015) and in the regime ( = 0.002). Note that the value of the regime run in QCD is slightly higher () than in other runs (). In the plot, we adjust the value of by a factor of 1.065 which corresponds to the ratio of lattice spacings between the two values.

In the plots, we clearly observe the sea quark mass dependence. For heavier sea quarks, the spectral density shows a higher peak near the lowest eigenvalue (around 0.004), and the peak height becomes lower by reducing the quark mass. As one enters the regime, the lowest eigenvalue is pushed up to 0.015. This is what we expect for the effect of the fermionic determinant . Namely, when the quark mass is reduced to the value around (or below) the lowest eigenvalue, those eigenvalues are suppressed.

The expectation from NLO ChPT precisely follows the lattice data (solid curves in Figures 3 and 4). Here the values of and are determined for each set of sea quark mass, so that the comparison is not parameter-free. But, still the precise agreement of the shape of the spectral density is encouraging. We investigate the mass dependence of in the next section.

In Figure 5, similar plots are also shown for the lattice data where three sea quarks have a degenerate mass (). We have four values of the quark mass (0.100, 0.080, 0.035 and 0.025) in the regime. In QCD, we observe a stronger dependence on the quark mass than in = 2 or , as suggested by the NLO formula for in (11).

### v.2 Topological charge

In the low eigenvalue region, the Dirac spectral density is known to be sensitive to the topological charge of the gauge fields, which is clearly seen in Figures 6 and 7. The solid curves in the plots show the expectation from the NLO ChPT with input parameters determined from the lattice data. Therefore, there is no free parameter to be adjusted in the curves for the cases. We observe that the dependence of the lattice data is qualitatively well described by ChPT.

Note, however, that the extracted values of and from data show a 2.2 difference for while they are consistent in the case. Since the topological charge dependence is a part of finite volume effects Brower:2003yx ; Aoki:2007ka ; Aoki:2009mx , which should be accounted for by the effective theory analysis, the deviation suggests possible higher order effects in the expansion of ChPT. Indeed, the runs are carried out with a shorter temporal extent than that of (), so that the lattice volume is % smaller in the physical unit. Higher order finite volume effects may therefore be enhanced for . In the final results of the data, we add this 11% deviation as an estimate of the systematic error due to the finite volume.

### v.3 Finite volume

The finite volume scaling can be tested more explicitly with the runs, for which the lattice data are available. Here we note that the comparison of obtained on the and lattices is not straightforward, because there is a finite volume effect encoded in in the definition of (11). It is still possible to analytically convert the values of in different volumes. The results for obtained on the lattice are converted to those at the smaller volume as = at and = at , which agrees with the values calculated on the lattice: = 0.00314(18) and 0.00333(18), respectively. In fact, as shown in Figure 8, we find that the same inputs of (converted) and describe the data at different volumes very well.

From these analysis, the systematic error due to finite volume is estimated as %, which is taken from the difference between the and (after the conversion) results.

## Vi Determination of the chiral condensate

The extracted values of and for each lattice ensemble are summarized in Table 2. Note that is extracted at the NLO accuracy, while the value of which first appears in the NLO term, might receive larger systematic corrections from the NNLO contributions.

As already noted above, since is fixed at a large value (0.080 or 0.100) in the ensembles, there is little difference between the reduced and ChPT analysis near the chiral limit of : and are almost equal within the statistical errors. The difference between = 0.080 and = 0.100 is also small (always less than 1), and therefore we concentrate on the data at in the following.

Now we analyze the sea quark mass dependence of . The NLO formula for (11) contains the low energy constants , and as parameters. The chiral condensate , in particular, appears at the leading order, and its determination through the quark mass dependence of is valid at the NLO accuracy, while the other parameters controlling the NLO correction can be determined at the leading-order.

In the fit of the lattice data, we attempt two procedures: a 3-parameter (, , ) fit without any additional inputs, and 2-parameter (, ) fits for several input values of . In the 2-parameter fits, the value of determined in our previous works both in the regime ( = 0.0474(30)) Noaki:2008iy and in the regime ( = 0.0524(34)) Fukaya:2007pn are used for the lattice ensembles. Any difference due to the different input values of suggests some systematic error (although it turns out to be only in the final result). For the = and 3 runs, we use a naive (linear) chiral limit of given in Table 2, of which values are for ChPT and for reduced ChPT, respectively.

The chiral extrapolation of in QCD is shown in Figure 9. From the plot, we can see the crucial role played by the data point at = 0.002 which is in the regime. Without this data point, one may naively expect that the data do not have enough sensitivity to probe the curvature due to the chiral logarithm and the chiral extrapolation favors a larger value of ( 0.0032). Taking the regime data into account, the presence of chiral logarithm is consistent with the negative curvature seen in the data.

fit range | =2 ChPT LECs | d.o.f. | ||
---|---|---|---|---|

2prm fit | (0.0474) | |||

0.002-0.025 | 0.00243(18) | [0.0474] | -0.00004(21) | 4.8 |

0.002-0.035 | 0.00246(15) | [0.0474] | -0.00009(13) | 2.5 |

0.002-0.050 | 0.00233(12) | [0.0474] | 0.00007(10) | 2.4 |

0.002-0.070 | 0.00233(09) | [0.0474] | 0.00006(07) | 1.8 |

0.002-0.100 | 0.00246(07) | [0.0474] | -0.00004(03) | 2.4 |

2prm fit | (0.0524) | |||

0.002-0.025 | 0.00250(19) | [0.0524] | 0.00012(30) | 5.2 |

0.002-0.035 | 0.00255(15) | [0.0524] | 0.00002(18) | 2.7 |

0.002-0.050 | 0.00243(12) | [0.0524] | 0.00021(14) | 2.4 |

0.002-0.070 | 0.00246(10) | [0.0524] | 0.00018(09) | 1.8 |

0.002-0.100 | 0.00262(08) | [0.0524] | 0.00001(04) | 2.8 |

3prm fit | ||||

0.002-0.035 | 0.00174(45) | 0.0290(71) | -0.00024(03) | 2.5 |

0.002-0.050 | 0.00246(68) | 0.0547(78) | 0.00038(90) | 3.5 |

0.002-0.070 | 0.00237(38) | 0.0489(15) | 0.00009(35) | 2.4 |

0.002-0.100 | 0.00206(26) | 0.0386(48) | -0.00009(02) | 2.3 |

The extracted values of the LECs in ChPT are summarized in Table 3. We attempt the 2-parameter fits of various number of data points, 3–7, taken from the lowest quark mass and the 3-parameter fits with 4–7 data points. In the table, the range of used in the fit is listed. For the 2-parameter fits, we take two input values of . The quality of the fits can be inferred from the value of d.o.f. also listed in the table.

As far as the heaviest point is discarded in the fits, the resulting value of is insensitive to the input value of , and it is consistent with the three-parameter fit as well. On the other hand, the determination of is unstable, but all the data suggest . We take and as the central values, which are from the 4-point fit with the input . For the final results of this paper, we take the deviation from the other input value of , as well as the other fitting ranges, as a systematic error due to the chiral extrapolation.

fit range | =3 ChPT LECs | d.o.f. | ||
---|---|---|---|---|

2prm fit | ||||

0.002-0.025 | 0.00186(09) | [0.0411] | -0.00013(09) | 1.0 |

0.002-0.035 | 0.00186(09) | [0.0411] | -0.00014(08) | 0.6 |

0.002-0.050 | 0.00186(09) | [0.0411] | -0.00014(07) | 0.5 |

0.002-0.080 | 0.00185(08) | [0.0411] | -0.00014(07) | 0.4 |

3prm fit | ||||

0.002-0.035 | 0.00185(10) | 0.0433(13) | -0.00023(53) | 1.3 |

0.002-0.050 | 0.00186(09) | 0.0406(05) | -0.00012(25) | 0.7 |

0.002-0.080 | 0.00186(08) | 0.0413(02) | -0.00015(09) | 0.5 |

fit range | reduced =2 ChPT LECs | d.o.f. | ||
---|---|---|---|---|

2prm fit | ||||

0.002-0.025 | 0.00199(06) | [0.0406] | -0.00044(10) | 0.9 |

0.002-0.035 | 0.00197(06) | [0.0406] | -0.00040(09) | 0.7 |

0.002-0.050 | 0.00197(05) | [0.0406] | -0.00039(06) | 0.4 |

0.002-0.080 | 0.00198(05) | [0.0406] | -0.00042(04) | 0.5 |

3prm fit | ||||

0.002-0.035 | 0.00197(09) | 0.0407(10) | -0.00042(51) | 1.3 |

0.002-0.050 | 0.00198(07) | 0.0416(05) | -0.00038(21) | 0.6 |

0.002-0.080 | 0.00196(07) | 0.0399(03) | -0.00044(07) | 0.5 |

For the lattice data, the chiral extrapolation is more stable because of the precise data point in the regime, which is simply due to higher statistics we accumulated. Figure 10 clearly shows the logarithmic curvature. Namely, a naive linear extrapolation of four data points in the regime ( = 0.015–0.050) would lead to 0.0028 in the chiral limit, while the regime point is lower ( 0.0020).

We analyze the lattice data listed in Table 2 using the 2- and 3-parameter fits. We consistently use the formula for the data obtained with ChPT (fourth column in Table 2). The same applies for the reduced ChPT analysis. The fit results are presented in Table 4 and in Table 5. The curves in Figure 10 represent the fits for various numbers of data points. We find that all the curves go through the data points except for the heaviest one. The chiral limit shown by a cross symbol is almost unchanged by taking different fit schemes (2-parameter or 3-parameter) and the number of data points included in the fit (4, 5, or 6). This is because the precise regime point works as an anchor near the chiral limit.

For ChPT, listed in Table 4 denotes the value of in the limit of and with a fixed strange quark mass . This corresponds to the physical value of the chiral condensate defined in the limit of vanishing up and down quark masses.

From Tables 4 and 5, we can see that results for obtained via and reduced ChPT formulas are consistent with each other. (We assume that in the theory corresponds to in reduced ChPT.) The result is also stable against the changes of the number of fitting parameters and the fitting range. For the other parameters, and , we find larger dependence on the choice of fit procedures, which is expected because they only appear in the NLO terms.

We take , and as the central values, which are obtained from the 5-point fit with three free parameters in ChPT. As mentioned above, the other results are used to estimate systematic errors.

fit range | =3 ChPT LECs | d.o.f. | |||
---|---|---|---|---|---|

2prm fit | |||||

0.025-0.100 | 0.00152(11) | 0.00204(04) | [0.0431] | 0.00010(10) | 3.1 |

0.025-0.100 | 0.00182(14) | 0.00242(06) | [0.0531] | 0.00031(17) | 2.8 |

3prm fit | (combined with +1 data) | ||||

0.002-0.100 | 0.00145(12) | 0.00191(08) | 0.0401(17) | 0.00003(7) | 1.7 |

For the degenerate () lattice results, the number of data points is limited to four. The 3-parameter fit, therefore, does not work and we restrict ourselves to the 2-parameter fit. We attempt the fit with two values of , 0.0431 and 0.0531, as the input. The former value is an chiral limit of in Table 2 taken with a linear function in the quark mass and the latter is the one at the lightest sea quark mass . Due to the lack of the regime data point, the chiral limit is not as stable as in the or data. In fact, the resulting value of strongly depends on the input value of . Between the two representative values of , changes about 20%. Since controls the NLO effects as seen in (11), this change suggests that the one-loop calculation is not sufficient to control the dependence.

We also attempt a combined fit of all the and data points using the ChPT formula. The total nine data points are simultaneously fitted with a reasonable ( 1.7). The result is given in Table 6 and plotted in Figure 11 (and 10) which indicate a strong effect of chiral logarithm in the data compared to the case of shown in Figure 9. The chiral limit is less than half of the value at the lowest regime point (). One should note, however, that the data points in the fit includes those with and 0.100 which are out of the typical convergence region of chiral expansion , so that the result may contain large systematic effects.

It is still remarkable that all the results suggest that the chiral condensate of QCD is smaller than of or of QCD. This is consistent with the view that the chiral condensate decreases and eventually disappears as the number of flavor increases and the asymptotic freedom is lost. We take , and determined from the combined fit as the central values of ChPT parameter, and will include a 26% deviation from the 2-parameter fit with as a systematic error in the final results.

So far we have treated the strange quark mass fixed at in the studies. In fact, the combined fit of the degenerate and 2+1 results implies that in the chiral limit of changes by less than 1% when varies between 0.060 and 0.100. We can, therefore, safely ignore the error due to a slight mismatch of the strange quark mass from its physical value. This weak sensitivity to supports the assumption that the strange quark at the physical mass is almost decoupled from the low energy theory and the use of the reduced ChPT formula to fit the lattice QCD data.

## Vii Summary and Conclusion

Before quoting the final results, let us discuss the possible systematic errors.

Our simulation uses an exactly chiral-symmetric Dirac operator and has reached almost the chiral limit: =0.002 for the and 2+1 runs, which corresponds to 3 MeV in the physical unit. The NLO ChPT formula (5) is valid in both the and regimes. As a result, the chiral extrapolation of in and 2+1 QCD is stable. As discussed in the previous section, by varying the fit range (and in the ChPT formula for the analysis), we estimate the systematic effects due to the chiral extrapolation as for of QCD and for of QCD, respectively. The upper and lower limits come from the variation of with the various fit schemes. We also found that the dependence is negligible in the range = 0.060–0.100, so that can be treated as the chiral condensate at the physical value of the strange quark mass.

On the other hand, for the degenerate lattice data, our estimate for the systematic error in the chiral fit is larger because of the bad convergence of ChPT and smaller number of data points. In fact, we observe that moves as large as 26% depending on the fit schemes, from which we estimate the systematic error from this source to be 26% for in the chiral limit.

The finite volume effects have also been discussed in the previous sections. For the lattice results, we expect that a possible higher order effect in the expansion beyond one-loop ChPT is partly reflected in the difference among different topological sectors. We thus estimate the systematic error from this source to be 11%. For the (and 3) case, we use more direct comparison with lattice results and the systematic error is estimated as . With a naive order counting, the leading finite volume effect is estimated as , which is the two-loop effect in the expansion. With = 71 MeV (see below), this gives a large value ( 0.52) as the size of the two-loop correction. In fact, including the numerical coefficient given in the Appendix A the one-loop correction is 0.06. If we assume that the numerical coefficient is also small ( 0.05) at the two-loop order, this naive estimate gives a 3% effect, which is in the same ball park as the estimate given above. The small numerical coefficient at the two-loop level is indeed obtained in a recent study Lehner:2010mv .

Since our lattice studies are done at only one value of , it is difficult to estimate the size of discretization effects. It should be partly reflected in the mismatch of the observables measured in different ways. For instance, the inverse lattice spacing determined from the -meson mass, (17) GeV, is 1% larger than the determination from the -baryon mass, (10) GeV JLQCD:prep . The latter corresponds to the Sommer scale = 0.51 fm, which is higher than the nominal value 0.49 fm or the recently favaored value 0.46–0.47 fm by about 4–10%. On the other hand, a naive order counting with MeV suggests a systematic effect of %, which is consistent with the above mismatch. We therefore add this naive estimate, , as the systematic error due to finite lattice spacing.

We convert the value of the condensate
to the definition in the standard renormalization scheme,
i.e. the scheme.
By using the nonperturbative renormalization technique
through the RI/MOM scheme we obtained the factor in our previous
work Noaki:2009xi as
for
and for and 3,
where the errors are statistical and systematic, respectively^{2}^{2}2
The value for is slightly changed from Noaki:2009xi
due to the different determination of the lattice scale,
that affects the renormalization group running of the factor.
.

ChPT | ChPT | |||
---|---|---|---|---|

renormalization | % | % | % | – |

chiral fit | % | % | % | % |

finite volume | % | % | % | % |

finite | % | % | % | % |

total | % | % | % | % |

Including all the systematic effects above, which are summarized in Table 7, we obtain the chiral condensate of up and down quarks in their massless limit as

(15) |

where the errors are again statistical and systematic, respectively. Here, the total systematic errors are obtained by adding each estimate by quadrature. Note that the result for is slightly changed from Fukaya:2009fh because a different input for the scale determination is used. We find a nontrivial dependence on the chiral condensate: as the strange quark mass goes down from ( QCD) to the chiral limit , the value of decreases.

The chiral condensate determined in this work (15) is consistent with those in our previous results obtained from the pseudoscalar meson mass Noaki:2008iy ; JLQCD:2009sk and from the topological susceptibility Aoki:2007pw ; Chiu:2008kt ; Chiu:2008jq . A similar work done using the Wilson fermion and the regime ChPT Giusti:2008vb quoted MeV, which is slightly higher than our result. A recent work Bernardoni:2010nf in a mixed action approach (overlap valence + Wilson sea) has also reported a larger value. More detailed study would be necessary to understand the source of the discrepancy, if it is significant.

From the NLO terms, we also obtain

(16) |

(or 100(4)(11) MeV) and

(17) |

where the systematic errors are estimated in a similar manner. For , their estimates are listed in Table 7 while for , the systematic error is dominated by the one from chiral extrapolation, as seen in Table 3 and 4. Although the accuracy for these quantities is not as good as that of the chiral condensate, they provide important consistency checks.

In this study, we have investigated the eigenvalue spectrum of the QCD Dirac operator, which is free from ultraviolet power divergences. The lattice QCD results show a good agreement with the ChPT calculation at NLO in the region of less than