We generated configurations with the parametrized fixed-point Dirac operator on a box at a lattice spacing . We compare the distributions of the three lowest eigenvalues in the topological sectors with that of the Random Matrix Theory predictions. The ratios of expectation values of the lowest eigenvalues and the cumulative eigenvalue distributions are studied for all combinations of and . After including the finite size correction from one-loop chiral perturbation theory we obtained for the chiral condensate in the scheme , where the error is statistical only.

2+1 Flavor QCD simulated in the -regime in different topological sectors

P. Hasenfratz, D.  Hierl, V. Maillart, F. Niedermayer, A. Schäfer, C. Weiermann and M. Weingart

Institute for Theoretical Physics

University of Bern

Sidlerstrasse 5, CH-3012 Bern, Switzerland

Institute for Theoretical Physics, University of Regensburg, D-93040 Regensburg, Germany

Spontaneous chiral symmetry breaking and the related existence of light Goldstone bosons is a basic feature of QCD. Chiral Perturbation Theory (ChPT) provides a systematic description of this physics in terms of a set of low energy constants which encode the related non-perturbative features of QCD. The method of low energy effective Lagrangians simplifies the calculations significantly [1, 2, 3, 4] and over the years ChPT became a refined powerful technique. Gasser and Leutwyler recognized very early, decades before the numerical calculations could attack such problems, that these constants can also be fixed using physical quantities which, presumably, will never be measured in real experiments. They can be studied, however, in lattice QCD.

The - regime [5, 6, 7, 8] describes physics close to the chiral limit in a box whose size is larger than the QCD scale. On the other hand, the size of the box relative to the Goldstone boson correlation length must be small. Under these conditions the Goldstone bosons, as opposed to other excitations, feel the effect of boundaries strongly. ChPT provides a powerful and systematic way to calculate the finite size corrections. A nice additional feature of the -regime is that Random Matrix Theory (RMT) [9] makes precise predictions for microscopic observables. RMT relates in particular the distribution of low-lying eigenvalues of the Dirac operator in different topological sectors to the chiral condensate.

The pioneering numerical works [10, 11, 12, 13, 14, 15, 16] in the -regime in quenched QCD suggested that this regime can be an excellent tool to study low-energy physics in QCD. The special problems of the -regime called for new numerical procedures [14, 15, 16, 17] which were first tested also in this approximation. Combining these numerical developments with the renormalization properties of the spectral density the work [18] present a state of the art analysis for the quark condensate and the first Leutwyler-Smilga sum rule in quenched QCD. The quenched approximation, however, tends to be singular in the chiral limit and it is expected that quenching is even more problematic here than in other cases [19].

Unfortunately, full QCD simulations are expensive. In addition, in the regime the Dirac operator should have excellent chiral properties. The standard choice is the overlap Dirac operator [20] with a hybrid Monte Carlo algorithm.

The results in [21] reflect already the basic physics features of the -regime. The distribution of the low-lying eigenmodes could be fitted quite well to the RMT predictions and a reasonable value for the chiral condensate was obtained. These results are promising given the fact that the box was small (), the quark mass was rather large () and the lattice was coarse (). The results in a larger box of and at smaller quark mass obtained in [22] were consistent with those of [21].

The first serious simulation results in the -regime have been presented by the JLQCD group recently [23]. In that work QCD was simulated with overlap fermions (created with the Wilson kernel) on a lattice of size with a resolution obtained from . Different bare quark masses () were considered in the range . The smallest quark mass (corresponding to ) was certainly small enough to reach the -regime. Hybrid Monte Carlo with overlap fermions has problems when the topological charge changes. To avoid this an action was used which prevented topology change and the whole run stayed in the sector. The authors observed an overall good agreement with RMT. The fermion condensate, which is the only parameter to be fitted, was found to be in the scheme.

Our work is based on the parametrized fixed-point (FP) action which has been tested in detail in quenched QCD [24]. The exact FP Dirac operator satisfies the Ginsparg-Wilson relation


where is a local operator and is trivial in Dirac space. The parametrized FP action has many gauge paths and involves a special smearing with projection to the gauge group . As a consequence, hybrid Monte Carlo algorithms can not be used. In the work [25] we advised a partially global update with three nested accept/reject steps which can reach small quark masses even on coarse lattices. Several steps of this algorithm were developed in [26] using some earlier suggestions [27, 28, 29, 30]. All the three steps are preconditioned. Pieces of the quark determinant are switched on gradually in the order of their computational expenses. The largest and most fluctuating part of the determinant is coming from the ultraviolet modes. We reduce these fluctuations by calculating the trace of [31, 32, 26]. The lowest-lying modes are calculated and subtracted. The determinant of the reduced and subtracted is calculated stochastically [33, 34, 35, 36]. For the determinant breakup we generalized the mass shifting method of [31]. For the strange quark we performed a root operation. Since the strange quark mass is not very small, we used a polynomial expansion to approximate this root operator. Recently we replaced the polynomial expansion by a rational approximation [37] which brought a performance gain in this part.

Due to the partial global updating procedure the algorithm, beyond a certain volume, scales with which constrains the size of lattices which can be considered. Due to the large number of subtracted low-lying modes, small quark mass is not a barrier.

We generated configurations in the Markov chain on an lattice using the partially global algorithm discussed above. The distance between every second configurations in this chain is similar to that of two gauge configurations separated by a typical Metropolis gauge update sweep. We considered every tenth configurations from the Markov chain and performed the measurements on the remaining configurations. When calculating statistical errors we formed 20 bins and used the jackknife method.

We fixed the lattice spacing from the Sommer parameter [38] and found . Our box has the size .

The Dirac operator has no exact chiral symmetry due to parametrization errors. As a consequence, the quark masses have an additive mass renormalization. The degenerate and the quark masses in the code are and . The additive renormalization was measured following the steps described in [39] and has the value , as shown in Fig. 1. Subtracting the additive renormalization we get the bare masses and which corresponds to and , respectively.

Figure 1: The AWI mass for the three combinations of quark masses. The linear extrapolation to gives an additive mass renormalization .

As we mentioned above, a simulation with lighter quark masses would not be more expensive. In our case, however, a smaller were not useful either. RMT predicts the probability distribution for the -th low-lying eigenvalue () of the Dirac operator in the topological sector . Denoting the corresponding eigenvalues of the (continuum) Dirac operator by , the variable is related to the bare chiral condensate as . Here V is the volume of the box. Fig. 2 shows the prediction of RMT for the cumulative distributions for and . The distributions depend on where are the quark masses. Decreasing the mass by more than a factor 10 (at fixed ) the cumulative distributions practically remain unchanged. On the other hand sending the strange quark mass to infinity (at fixed ) we land on a flavor theory with a visibly different distribution. The strange quark has a (modest) effect on our observables.

Figure 2: Random Matrix Theory prediction for the cumulative distribution of , where is the k-th eigenvalue of the continuum Dirac operator in a gauge background with topological charge , and . Here the results are shown, but the picture is similar for .

While the eigenvalues of the continuum Dirac operator lie on the imaginary axis the spectrum of lattice Dirac operators is more complicated. As it is well known, the spectrum of the Dirac operator satisfying the Ginsparg-Wilson relation with lies on the circle . In our case where is different from 1, it is convenient to introduce the rescaled operator , for which the factor is eliminated in the GW relation, and the spectrum lies on the circle. (In fact, our operator for the low-lying modes is effectively a constant close to 1 within a few percent.) To relate the eigenvalues on the GW circle to those appearing in the RMT (or in general, in continuum expressions) it is natural to use the stereographic projection


As Fig. 3 shows, the low-lying eigenvalues of (which we obtain during our simulation) are close but not exactly on the GW circle, hence we project them first horizontally onto the GW circle, and make the stereographic projection in the next step. This is justified by the observation that making a systematic GW improvement of towards the imaginary part of the eigenvalues stays practically constant – they move horizontally to the GW circle. The values of obtained by this procedure are then compared with the RMT predictions.

Figure 3: The histogram of the deviation from the circle, , for the eigenvalues of on and lattices for the complex eigenvalues with , using the present mass parameters.

Having the probability distributions from RMT we can calculate the ratios , where the factor cancels. These predictions are compared with the measured ratios in Fig. 4, where all the ratios in the topological sectors for the first three lowest eigenvalues are shown.

Figure 4: Ratios of expectation values obtained from simulations compared to the results of RMT, where is indexing the -th lowest eigenvalue in the topological sector . The different symbols refer to the denominator while of the numerator is indicated at the data points. For example, the highest ratio with value in the figure refers to


At this point we have to discuss the way we identify the topological sectors. According to the index theorem [40] one can identify the topological charge from the zero modes of a Ginsparg-Wilson Dirac operator. However, our has no exact chiral symmetry and has modes on the real axis. Occasionally, the real eigenvalue might be even far from zero in which case the topological interpretation is uncertain. We note here that using an exact Ginsparg-Wilson operator (with kernel ) for measuring the topological charge is not a good solution to this problem. The overlap just projects most of the real eigenvalues to the point . Following the intuitive picture that topology is related to extended objects we investigated the correlation between the inverse participation ratio (IPR) of the eigenvector having a real eigenvalue and the value of . Here IPR is given by for a normalized eigenvector . The IPR is inversely proportional to the size of the effective support of the eigenfunction. As Fig. 5 shows, there is a strong correlation indeed. As is moving away from zero, the IPR increases indicating that the wave function becomes more and more localized. Fig. 6 demonstrates that the real eigenvalues are strongly concentrated in a small region close to zero. The histogram of has a large, narrow peak followed by a long tail.

Figure 5: The inverse participation ratio for real eigenvalues vs. the eigenvalue. Eigenvectors with larger real part have smaller extension in agreement with the expectations.
Figure 6: Histogram of for the real eigenvalues. In determining the topological charge we considered only the real eigenvalues with .
Figure 7: Histogram of inverse participation ratio for the real eigenmodes of .

Fig. 7 shows the histogram of IPR for the real eigenmodes of . There is a strong, narrow peak which corresponds to extended real modes and a long tail corresponding to wave functions with smaller support.

For the questions we study in this paper it is not necessary to use a procedure which leaves the relative weights of the different topological sectors intact. (Due to the long autocorrelation time of the topological charge it would be anyhow difficult to control these weights.) We introduce therefore a cut above which the real eigenvalues are discarded. This cut must be sufficiently small to allow for a good identification of the topological charge, but large enough to allow for a reasonable number of configurations in each sectors we want to investigate. The cut should be smaller than the spectral gap (for the spectrum see Fig. 8), otherwise the effect of topology will be strongly reduced. We introduced a cut of size 0.03 which is approximately equal to the minimal gap in the sector. We accepted a real eigenvalue as a sign of topology if . This cut splits our sample of 414 configurations into 328, 51 and 35 configurations for the sectors, respectively. This cut corresponds also to the end of the peak of the IPR distribution shown in Fig. 7.

Figure 8: The spectrum of for 50 configurations on and lattices, at the same bare mass parameters and gauge coupling.

Fig. 9 compares the cumulative distributions of the first three eigenvalues in the topological sector obtained using the overlap improved fixed-point Dirac operator vs. in the measurement. The distributions are very similar showing that the two operators are close to each other. The distribution from has a tail at small eigenvalues (i.e. has a smaller gap than ) demonstrating a systematic error one obtains when using different Dirac operators for the generation of the configurations and for the measurements: the small eigenvalues are less effectively suppressed.

In Figs. 10, 11, 12 the cumulative distributions of the operator are compared with RMT for the topological sectors. The only matching parameter is the bare condensate which enters the RMT predictions through and in . We obtained the result for this bare quantity . Both errors are statistical. The first error comes from the statistical errors of the measured distributions, the second one is due to the error in the lattice scale .

The result above receives a small correction due to the fact that, for simplifying the presentation, we suppressed a technical complication. As mentioned before the exact fixed-point operator satisfies a Ginsparg-Wilson relation with a local operator . For this reason the quark mass enters in a simple additive way in the form , where . Effectively the operator behaves in the infrared like a constant close to 1. Its expectation value for the lowest eigenvectors is within 1%. Using the spectrum of (which is in fact identical to that of ) the matching with RMT gives the slightly changed result .

In the language of a ferromagnetic model one should interpret as the absolute value of the magnetization in the finite volume . This differs from the value defined in the infinite volume by a finite-size correction which can be calculated in ChPT. In the presence of the magnetic field the orientation of magnetization is controlled by the Boltzmann factor where is the total magnetization. This gives


where is related to the modified Bessel functions. Comparing this expression with the result of ChPT [7] we get for (corresponding to two flavors) , where


Here the shape coefficient takes the value 0.14046 for a symmetric box [7]. The one-loop finite-size correction to the order parameter (magnetization in the example above) has been calculated for the symmetry group in [5, 6], for the group in [41] and up to the two-loop level in [7].

Neglecting the strange quark contribution in the finite size effects we use the two-flavor result to correct the measured bare by the one-loop finite size correction to obtain .

In the previous version of this paper we quoted only this bare value. Since then we have calculated the renormalization factor connecting the bare lattice result to the scheme. The main points of this calculation are summarized in the Appendix. With the conversion factor we obtain

Figure 9: The cumulative distributions in the topological sector for the three lowest complex eigenvalues as measured with and with the overlap operator using the fixed-point operator as kernel.
Figure 10: The cumulative distribution of for , . Here is obtained from by stereographic projection described in the text, and is rescaled using . The continuous line is the RMT prediction at and .
Figure 11: The cumulative distribution of for , .
Figure 12: The cumulative distribution of for , .

Acknowledgements We thank Gilberto Colangelo, Stephan Dürr, Jürg Gasser, Anna Hasenfratz and the members of the BGR Collaboration for valuable discussions. We also thank the LRZ in Munich and CSCS in Manno for support. The analysis was done on the PC clusters of ITP in Bern. This work was supported by the Schweizerischer Nationalfonds.


The bare lattice scalar density is not universal, it depends on the lattice action and it needs renormalization. The conventional way is to express the result in scheme at some given scale, e.g. .

To relate the bare lattice result to scheme one generally uses the Rome-Southhampton RI/MOM method [42]. In this method one introduces an intermediate scheme (or 111These schemes differ in their definition of the quark field renormalization factor .) where some Green’s functions of the actual operator (calculated in a fixed gauge) are equated to their corresponding Born terms at a given scale .

In the first step one relates the bare lattice operator non-perturbatively to its renormalized counterpart in the scheme


The matching scale is restricted by two conditions: It should be sufficiently large to avoid non-perturbative effects and small enough to avoid large cut-off effects. For the scalar density, the non-perturbative matching was done in [43], using a technique which removes a large part of cut-off effects.

In the second step one relates the renormalized operators in the and schemes in perturbation theory


Obviously, the value of must lie in the perturbative regime. The factor has been calculated up to NNLO [44] and NNNLO [45]. Combining (A-1) and (A-2), we obtain the matching factor connecting the bare lattice and the renormalized results


As the expansion parameter in the perturbative series is the running coupling of QCD, we need to determine it for the scale where the matching is done. The running of is described by the differential equation


where the 4-loop -function in scheme is given in [48]. As initial condition we use from [49] 222This result is obtained by using the PDG value [50] and running it from 5-flavors down to 3-flavors, across the and thresholds.The procedure described above is carried out at several matching scales .

In the third step Renormalization Group technique is used to obtain the desired conversion factor for the scalar density at ,


Here, due to the relation , the anomalous dimension is related to the mass anomalous dimension by . The latter has been calculated to four loops in [46, 47].

In Fig. 13, we plot our result at vs. the matching scale for the intermediate schemes and . If the non-perturbative, cut-off and higher order perturbative effects were negligible, the two curves would coincide and would be independent of the choice of . Taking the average values for and taking into account the slight difference between the curves, we get

Figure 13: The renormalization factor for the scalar density, connecting the and lattice schemes, at , obtained using and schemes at the intermediate step vs. the matching scale .


  • [1] S. Weinberg, Physica A96 (1979) 327.
  • [2] J. Gasser and H. Leutwyler, Phys. Lett. B 125 (1983) 321, 325.
  • [3] J. Gasser and H. Leutwyler, Ann. Phys. (N.Y.). 158 (1984) 142.
  • [4] J. Gasser and H. Leutwyler, Nucl. Phys. B 250 (1985) 465.
  • [5] J. Gasser and H. Leutwyler, Phys. Lett. B 184 (1987) 83; 188 (1987) 477.
  • [6] J. Gasser and H. Leutwyler, Nucl. Phys. B 307 (1988) 763.
  • [7] P. Hasenfratz and H. Leutwyler, Nucl. Phys. B 343 (1990) 241.
  • [8] F. C. Hansen, Nucl. Phys. B 345 (1990) 685.
  • [9] E. V. Shuryak and J. J. M. Verbaarschot, Nucl. Phys. A 560 (1993) 306 [arXiv:hep-th/9212088],
    J. J. M. Verbaarschot and T. Wettig, Ann. Rev. Nucl. Part. Sci 50 (2000) 343 [arXiv:hep-ph/0003017],
    P. H. Damgaard and S. M. Nishigaki, Phys. Rev. 63 (2001) 045012 [arXiv:hep-th/0006111]
  • [10] R. G. Edwards, U. M. Heller, J. E. Kiskis and R. Narayanan, Phys. Rev. Lett. 82 (1999) 4188.
  • [11] W. Bietenholz, K. Jansen and S. Shcheredin, JHEP 0307 (2003) 033.
  • [12] W. Bietenholz, Th. Chiarappa, K. Jansen Kei-ichi Nagai and S. Shcheredin, JHEP 0402 (2004) 023.
  • [13] L. Giusti, M. Lüscher, P. Weisz and H. Wittig, JHEP 0311 (2003) 23 [arXiv:hep-lat/0309189].
  • [14] L. Giusti, P. Hernandez, M. Laine, P. Weisz and H. Wittig, JHEP 0404 (2004) 013 [arXiv:hep-lat/0402002].
  • [15] R. G. Edwards, Nucl. Phys. Proc. Suppl. 106 (2002) 38 [arXiv:hep-lat/0111009].
  • [16] T. DeGrand and S. Schaefer, Comput. Phys. Commun. 159 (2004) 185 [arXiv:hep-lat/0401011].
  • [17] L. Giusti, C. Hoelbling, M. Lüscher and H. Wittig, Comput. Phys. Commun. 153 (2003) 31 [arXiv:hep-lat/0212012].
  • [18] L. Giusti and S. Necco, JHEP 0704 (2007) 090 [arXiv:hep-lat/0702013].
  • [19] P. H. Damgaard, Nucl. Phys. Proc. Suppl. 128 (2004) 47 [arXiv:hep-lat/0310037].
  • [20] H. Neuberger, Phys. Lett. B 417 (1998) 141 [arXiv:hep-lat/9707022]; 427 (1998) 353 [arXiv:hep-lat/9801031].
  • [21] T. DeGrand and S. Schaefer, Phys. Rev. D 71 (2005) 034507 [arXiv:hep-lat/0412005]; 72 (2005) 054503 [arXiv:hep-lat/0506021]; PoS LAT2005 (2006) 140 [arXiv:hep-lat/0508025].
  • [22] Th. DeGrand, Zh. Liu, S. Schaefer, Phys. Rev. D 74 (2006) 094504, Erratum-ibid. D 74 (2006) 099904 [arXiv:hep-lat/0608019];
  • [23] H. Fukaya et al. [JLQCD Collaboration], Phys. Rev. Lett. 98, 172001 (2007) [arXiv:hep-lat/0702003].
  • [24] P. Hasenfratz, S. Hauswirth,  K. Holland, T. Jorg, F. Niedermayer and U. Wenger, Int. J. Mod. Phys. C 12 (2001) 691 [arXiv:hep-lat/0003013],
    P. Hasenfratz, S. Hauswirth, T. Jorg, F. Niedermayer and K. Holland, Nucl. Phys.B 643 (2002) 280 [arXiv:hep-lat/0205010],
    C. Gattringer et al. Bern-Graz-Regensburg Collaboration, Nucl.Phys.B 677 (2004) 3 [arXiv:hep-lat/0307013],
    P. Hasenfratz, K. J. Juge and F. Niedermayer, Bern-Graz-Regensburg Collaboration, JHEP 0412:030,2004 [arXiv:hep-lat/0411034].
  • [25] A. Hasenfratz, P. Hasenfratz and F. Niedermayer, Phys.Rev. D72 (2005) 114508 [arXiv:hep-lat/0506024].
  • [26] A. Hasenfratz and F. Knechtli, Comput. Phys. Commun. 148 (2002) 81 [arXiv:hep-lat/0203010].
  • [27] M. Grady, Phys. Rev. D 32 (1985) 1496 [arXiv:hep-lat/0203010].
  • [28] M. Creutz, Algorithms for Simulating Fermions (1992), in Quantum Fields on the Computer, World Scientific Publishing.
  • [29] A. Alexandru and A. Hasenfratz, Phys. Rev. D 65 (2002) 114506 [arXiv:hep-lat/0203026]; 66 (2002) 094502 [arXiv:hep-lat/0207014].
  • [30] F. Knechtli and U. Wolff, (Alpha Coll.) Nucl. Phys. B 663 (2003) 3 [arXiv:hep-lat/0303001].
  • [31] M. Hasenbusch, Phys. Lett. B 519 (2001) 177 [arXiv:hep-lat/0107019].
  • [32] P. de Forcrand, Nucl. Phys.proc. Suppl. 119 (1999) [arXiv:hep-lat/9809145].
  • [33] A. D. Kennedy and J. Kuti, Phys. Rev. Lett. 54 (1985) 2473.
  • [34] A. D. Kennedy, J. Kuti, S. Meyer and B. J. Pendleton, Phys. Rev. D 38 (1988) 627.
  • [35] B. Joo, I. Horvath and K, F. Liu, Phys. Rev. D 67 (2003) 074505 [arXiv:hep-lat/0112033].
  • [36] I. Montvay and E. Scholz, Phys. Lett.  B 623 (2005) 73 [arXiv:hep-lat/0506006].
  • [37] M. Weingart, ’Stochastic Estimator of the Quark Determinant in Full QCD Simulation’, Diploma work, Uni. Bern, 2007.
  • [38] R. Sommer, Nucl. Phys. B 411 839 (1994) [arXiv:hep-lat/9310022].
  • [39] A. Hasenfratz, P. Hasenfratz, D. Hierl, F. Niedermayer and A. Schäfer, PoS LAT2006 (2006) 178 [arXiv:hep-lat/0610096].
  • [40] P. Hasenfratz,V. Laliena and F. Niedermayer, Phys. Lett. B 427 (1998) 125 [arXiv:hep-lat/9801021].
  • [41] H. Neuberger, Phys. Rev. Lett, 60 (1988) 880; Nucl. Phys. B 300 [FS22] (1988) 180.
  • [42] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, Nucl. Phys. B 445 (1995) 81 [arXiv:hep-lat/9411010].
  • [43] V. Maillart and F. Niedermayer, arXiv:0807.0030 [hep-lat].
  • [44] E. Franco and V. Lubicz, Nucl. Phys. B 531 (1998) 641 [arXiv:hep-ph/9803491].
  • [45] K. G. Chetyrkin and A. Retey, Nucl. Phys. B 583 (2000) 3 [arXiv:hep-ph/9910332].
  • [46] K. G. Chetyrkin, Phys. Lett. B 404 (1997) 161 [arXiv:hep-ph/9703278].
  • [47] J. A. M. Vermaseren, S. A. Larin and T. van Ritbergen, Phys. Lett. B 405 (1997) 327 [arXiv:hep-ph/9703284].
  • [48] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin, Phys. Lett. B 400 (1997) 379 [arXiv:hep-ph/9701390].
  • [49] Y. Aoki et al., Phys. Rev. D 78 (2008) 054510, arXiv:0712.1061 [hep-lat].
  • [50] W. M. Yao et al. [Particle Data Group], J. Phys. G 33 (2006) 1.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description