# CERN/TH-99-197 CPT-99/PE.3856 FTUV/99-49IFIC/99-51 Finite-size scaling of the quark condensate in quenched lattice QCD

## Abstract

We confront the finite volume and small quark mass behaviour of the scalar condensate, determined numerically in quenched lattice QCD using Neuberger fermions, with predictions of quenched chiral perturbation theory. We find that quenched chiral perturbation theory describes the numerical data well, allowing us to extract the infinite volume, chiral limit scalar condensate, up to a multiplicative renormalization constant.

Introduction   

Chiral symmetry breaking plays a central role in our comprehension of low energy QCD and understanding it from first principle calculations is of great importance. One of the cleanest ways of determining a condensate associated with the breaking of a global symmetry is through a finite-size scaling analysis. This technique has proved very successful in the study of scalar $`\mathrm{O}(\mathrm{N})`$ models . For chiral symmetry breaking in QCD, this would correspond to placing the system in a box and studying the scaling of the scalar condensate as a function of the volume $`V`$ and of the quark mass $`m`$ as the limit of restoration of chiral symmetry is approached ($`m0`$, $`V`$ finite).

The very small quark mass limit of QCD is expected to be well described by the lowest orders of chiral perturbation theory ($`\chi `$PT), which predict how the restoration of chiral symmetry takes place in a finite volume, as a function of the quark mass . The only free parameter entering the leading order contribution in the chiral expansion is the infinite volume quark condensate. Thus, a comparison of the mass and volume dependence of the finite volume quark condensate with the predictions of $`\chi `$PT provides a very powerful test of the hypothesis of spontaneous chiral symmetry breaking and permits an extraction of the infinite volume scalar condensate $`\mathrm{\Sigma }`$. Such a study requires, however, a good control over the chiral properties of the theory, which is difficult to achieve with traditional formulations of fermions on the lattice.

The situation is different, however, when Dirac operators that satisfy the Ginsparg–Wilson (GW) relation  are considered. Actions constructed from such operators have been shown to have an exact lattice chiral symmetry . This symmetry ensures that the relations implied by chiral symmetry in the continuum, hold also on the lattice at finite lattice spacing $`a`$ . For a review of the GW relation and its implications, we refer to .

A particular realization of an operator satisfying the GW relation, which we will be using here, has been proposed by Neuberger :

$$D_\mathrm{N}[m+(1+s)(1\gamma _5Q(Q^2)^{1/2})],$$
(1)   
where $`Qc_0\gamma _5(1+sD_\mathrm{W})`$, $`D_\mathrm{W}`$ is the Wilson Dirac operator, and the factor $`c_0`$ is a convenient normalization to keep the spectrum of $`Q^2`$ bounded by 1. The parameter $`s`$ satisfies $`|s|<1`$ and $`m`$ is the bare quark mass.

$`D_\mathrm{N}`$ satisfies the GW relation at zero quark mass. In contrast to the standard Wilson formulation, the breaking of the chiral symmetry is soft, i.e. only due to the quark mass term. This opens the possibility to confront finite volume simulations with finite-size scaling predictions in the regime of restoration of the chiral symmetry.

The complexity of the operator $`D_\mathrm{N}`$ renders its numerical treatment very demanding. We therefore restrict to the quenched approximation. The predictions of $`\chi `$PT must then be modified to take into account the effect of quenching. The finite size scaling of the quark condensate has recently been worked out using the framework of quenched chiral perturbation theory (q$`\chi `$PT) . In particular, analytical expressions for this scaling have been obtained in sectors of fixed topology. Operators satisfying the GW relation also satisfy an index theorem . Thus, by computing the eigenvalues of $`D_\mathrm{N}`$ at zero quark mass and identifying the zero modes, a clean separation of different topological sectors can be achieved, which is not possible with other formulations of lattice fermions. As we will see below, using the q$`\chi `$PT results in fixed topological sectors to interpret our numerical data proves very useful. A preliminary account of this work was presented at Lattice 99.

Light quarks on a torus   
To study the volume dependence of the scalar condensate, we work on a four-dimensional torus of volume $`L^4`$. Under the assumption that chiral symmetry is spontaneously broken, a description of QCD in terms of a chiral Lagrangian should be a good approximation at momenta $`p4\pi F_\pi `$. To lowest order in $`p/4\pi F_\pi `$ and in the quark mass, this Lagrangian is given by

$`={\displaystyle \frac{F_\pi ^2}{4}}\mathrm{Tr}[_\mu U^{}(x)_\mu U(x)]\mathrm{\Sigma }\mathrm{Re}\mathrm{Tr}[Me^{i\theta /N_f}U]`$ (2)   
where $`U(x)=\mathrm{exp}[i2\mathrm{\Pi }(x)/F_\pi ]SU(N_f)`$, $`\mathrm{\Pi }(x)`$ being the pion fields; $`M`$ is the quark mass matrix, which we take to be proportional to the identity matrix (i.e. $`M=mI`$), and $`\mathrm{\Sigma }`$ is the infinite volume and zero quark mass scalar condensate. In eq. (2) we have included the expected $`\theta `$ angle dependence.

Let us now consider the regime

$`M_\pi 1/LF_\pi ,`$ (3)   
where $`M_\pi ^2=2m\mathrm{\Sigma }/F_\pi ^2`$ to leading order in $`\chi `$PT. In this regime, the partition function is dominated by the zero mode of the $`U(x)`$ field , since the action of the non-zero modes has a kinetic contribution that goes like $`F_\pi ^2L^21`$. The partition function then reduces, to leading order, to an integral over the $`SU(N_f)`$ group manifold:

$`Z={\displaystyle _{SU(N_f)}}𝑑U_0e^{V\mathrm{\Sigma }\mathrm{Re}\mathrm{Tr}[Me^{i\theta /N_f}U_0]},`$ (4)   
where $`U_0`$ is the global mode. We can also define the partition function restricted to fixed topology by Fourier transforming in $`\theta `$ :

$`Z_\nu `$ $`=`$ $`{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\theta }{2\pi }}{\displaystyle _{SU(N_f)}}𝑑U_0e^{i\theta \nu }\mathrm{exp}[V\mathrm{\Sigma }\mathrm{ReTr}[Me^{i\theta }U_0]]`$ (5)   
$`=`$ $`{\displaystyle _{U(N_f)}}𝑑U_0det(U_0)^\nu \mathrm{exp}[V\mathrm{\Sigma }\mathrm{ReTr}[MU_0]].`$   
These integrals and their derivatives with respect to the quark mass have been known for a long time. For details see .

In our case, however, we are interested in the quenched approximation. Recently a similar reasoning has been applied to quenched QCD. The main difference in the quenched case is that the chiral symmetry group is no longer $`SU(N_f)_L\times SU(N_f)_R\times U(1)`$, but a graded Lie group $`U(1|1)_L\times U(1|1)_R/U_A(1)`$. According to , the partition function for fixed topology is then given by

$`Z_\nu ={\displaystyle _{U(1|1)}}𝑑U_0sdet(U_0)^\nu \mathrm{exp}[V\mathrm{\Sigma }\mathrm{Re}s\mathrm{Tr}[MU_0]].`$ (6)   
This integral has been computed analytically in terms of Bessel functions . By differentiating its logarithm with respect to the quark mass, $`m`$, the quark condensate for fixed topology is found to be

$`\mathrm{\Sigma }_\nu =\mathrm{\Sigma }z[I_\nu (z)K_\nu (z)+I_{\nu +1}(z)K_{\nu 1}(z)]+\mathrm{\Sigma }{\displaystyle \frac{\nu }{z}},`$ (7)   
where $`zm\mathrm{\Sigma }V`$ and $`I_\nu (z),K_\nu (z)`$ are the modified Bessel functions.

This formula summarizes the scaling of the quark condensate in a periodic box with the volume and quark mass in the small $`m\mathrm{\Sigma }V`$ limit, as a function of only one non-perturbative parameter: the infinite volume condensate $`\mathrm{\Sigma }`$. For fixed volume, the limit as $`m0`$ is given by

$`\mathrm{\Sigma }_{\nu =0}`$ $`=`$ $`m\mathrm{\Sigma }^2V\left(1/2\gamma +\mathrm{log}2\mathrm{log}m\mathrm{\Sigma }V+𝒪\left(m\mathrm{\Sigma }V\mathrm{log}m\mathrm{\Sigma }V\right)^2\right)`$   
$`\mathrm{\Sigma }_{\nu =\pm 1}`$ $`=`$ $`{\displaystyle \frac{1}{mV}}+{\displaystyle \frac{1}{2}}m\mathrm{\Sigma }^2V\left(1+𝒪\left(m\mathrm{\Sigma }V\mathrm{log}m\mathrm{\Sigma }V\right)^2\right).`$ (8)   
where $`\gamma `$ is the Euler constant. These results have two interesting features that we wish to emphasize. First, there is a divergence $`1/m`$ in sectors with topology. From the point of view of the underlying theory, this is not surprising since it corresponds to the contribution of the fermionic zero modes. Note however that these terms do not contain information about the infinite volume condensate and vanish in the infinite volume limit, as expected. The second interesting feature is the appearance of a logarithmic enhancement in $`\mathrm{\Sigma }_{\nu =0}`$, which is also peculiar to the quenched approximation. This term contains information about the infinite volume condensate.

In principle, by fitting the dependence of the finite volume condensate in quark mass and volume to Monte Carlo data, we can extract the infinite volume condensate. However, the naive bare quark condensate that is measured on the lattice is UV-divergent. A simple dimensional analysis of the possible divergences shows that the bare scalar condensate has a leading cubic divergence. One important advantage of Neuberger’s operator is that the coefficient of this leading divergence is known analytically. It is $`6/(1+s)`$, for $`SU(3)`$. The cubic divergence can then be subtrated exactly. However, after this trivial subtraction, the condensate is still divergent and has the form:

$`\mathrm{\Sigma }_\nu ^{sub}(a)\overline{\mathrm{\Psi }}\mathrm{\Psi }_\nu ({\displaystyle \frac{6}{1+s}}){\displaystyle \frac{1}{a^3}}=C_2{\displaystyle \frac{m(a)}{a^2}}+C_1{\displaystyle \frac{m(a)^2}{a}}+\mathrm{\Sigma }_\nu ,`$ (9)   
where $`m(a)`$ is the bare lattice mass. The constants $`C_i`$ are not known a priori and have to be determined, preferably non-perturbatively. The linear divergence proportional to $`m(a)^2`$ is negligibly small for the values of the mass and the cutoff we consider in this work. However, the quadratic divergence is not and turns out to be very important numerically. The condensate extracted through a fit of the lattice data to eqs. (9) and (7), of course, still requires a multiplicative renormalization to eliminate a residual logarithmic UV divergence in $`\mathrm{\Sigma }_\nu `$.

After subtracting the unphysical contribution of fermionic zero modes to $`\mathrm{\Sigma }_\nu `$, the finite volume condensate $`\mathrm{\Sigma }_\nu ^{sub}`$ vanishes, as expected, in the limit of zero quark mass . Not surprisingly, the power divergences can be separated, in principle, from the physical contribution to the condensate, by a study of the volume dependence of $`\mathrm{\Sigma }_\nu ^{sub}`$, while keeping the quark mass small enough to stay in the region of validity of $`\chi `$PT.

Even with the cubic divergence already subtracted, separating the physical condensate from the remaining power divergences may not be easy, in practice, because the statistical errors in these divergences can hide the small physical contribution. Thus, $`\mathrm{\Sigma }_\nu ^{sub}`$ must be computed with very good accuracy. Clearly, the logarithmic enhacement of $`\mathrm{\Sigma }_{\nu =0}`$ in eq. (8) could be very helpful in this respect; however, as we will see, extracting the condensate from the logarithmic term at zero topology requires much larger statistics than available to us at this time. We will concentrate instead on the study of the condensate in the topological charge one (or minus one) sector.

Numerical results   
For our numerical simulations we work on hypercubic lattices of size $`L^4`$ with periodic boundary conditions for both the gauge and the fermion fields. We work in the quenched approximation and use standard methods to obtain decorrelated gauge field configurations.

In selecting the value of $`\beta =6/g_0^2`$, some care has to be taken. On the one hand, the quadratic divergence $`1/a^2`$ should not hide the physical effect. On the other hand, in choosing too small values of $`\beta `$ there is the risk that Neuberger’s operator falls into a different universality class . Indeed, by computing the low-lying eigenvalues of Neuberger’s operator at $`\beta =5.7`$ and $`s=0`$, we only found eigenvalues $`\mathrm{O}(1)`$ and hence no light physical modes. A scan of the lowest eigenvalue of $`Q^2`$ as a function of $`s`$ showed that $`\lambda _{\mathrm{min}}(Q^2)`$ decreased with increasing $`s`$, contrary to what is expected (and found) at larger values of $`\beta `$.

The situation at $`\beta =5.85`$ appeared to be different, however. The values of $`\lambda _{\mathrm{min}}(Q^2)`$ reach a maximum around $`s=0.6`$, where the localization properties should also be optimal . In accordance, the eigenvalues of Neuberger’s operator became very small, so that there is little doubt that light, physical modes are present. Since the lattice spacing at $`\beta =5.85`$ is $`a^11.5\mathrm{GeV}`$  we estimated that the quadratic divergence term would not hide the physical signal, at least for reasonable values of the physical condensate.

A technical challenge is the numerical treatment of the square root appearing in Neuberger’s operator. We have chosen a Chebyshev approximation for this task, which allows us to reach a well controlled accuracy. In order to avoid any systematic effects in the values of physical observables, we demand that

$$XQ^2P_{n,ϵ}(Q^2)^2X^2/(2X)^2<10^{16}.$$
(10)   
In eq. (10) $`X`$ denotes a random vector and $`P_{n,ϵ}`$ denotes a standard Chebyshev approximation of the function $`1/\sqrt{x}`$ in the range $`ϵx1`$. $`P_{n,ϵ}`$ is a matrix-valued polynomial of degree $`n`$, which is constructed through numerically stable recursion relations . We require tantamount accuracies for all inversions. We note in passing that with the requirement of eq. (10) also the GW relation itself is satisfied to a similar accuracy for zero mass.

In order to decrease the degree of the polynomial employed, we have computed the 11 lowest eigenvalues of $`Q^2`$ and their corresponding eigenvectors and have set $`ϵ`$ to be the value of the largest. The contributions of these lowest lying eigenvectors are then treated exactly and projected out of the operator $`Q^2`$. Through this procedure, near-zero modes of $`Q^2`$ are taken into account automatically. All eigenvalue computations performed in our work are based on minimizing the Ritz functional .

As pointed out in , it is advantageous for the computation of the eigenvalues of Neuberger’s operator, and for its inversion as well, to stay in a given chiral subspace. This is possible because $`D_N^{}D_N`$ commutes with $`\gamma _5`$.

We computed the scalar condensate at several values of the quark mass using a multiple mass solver  on lattices of size $`8^4`$, $`10^4`$ and $`12^4`$. We checked through the calculation of the two lowest eigenvalues of Neuberger’s operator to which topological sector each gauge field configuration belonged. We then obtained $`\mathrm{\Sigma }_\nu ^{sub}`$ by computing

$$\mathrm{\Sigma }_\nu ^{sub}=\frac{1}{V}\mathrm{Tr}^{}\left\{\frac{1}{D_N}+\frac{1}{D_N^{}}\frac{a}{1+s}\right\}_\nu ,$$
(11)   
where the trace was performed in the chiral sector opposite to that with the zero modes and the gauge average was done in a sector of fixed topology $`\nu `$. With this definition, we take into account the contribution of all the non-zero eigenvalues of $`D_N`$ to the condensate <sup>4</sup><sup>4</sup>4With this definition the real eigenvalues at the cut-off level, $`m+2/a`$, are doubly counted. Although it is a completely negligible effect, we took it into account.. In this way, the term $`1/m`$ in eq. (8) is absent. Three gaussian sources and standard inverters were used to compute the trace in eq. (11). Topological charge zero configurations are very rare at larger volumes. For this reason we did not compute the condensate in this sector on the $`10^4`$ and $`12^4`$ lattices since the statistics we have gathered are too small.

We show in fig. 1 our results for $`a^3\mathrm{\Sigma }_{\nu =\pm 1}^{sub}/am`$ on our lattice volumes as a function of bare quark mass. We have 15, 10 and 7 gauge configurations on our $`8^4`$, $`10^4`$ and $`12^4`$ lattices, respectively. The solid lines are a fit of the data for all volumes and masses to eqs. (9) and (7). This fit has only two parameters, namely the infinite volume, zero quark mass, scalar condensate $`\mathrm{\Sigma }`$ and the coefficient of the quadratic divergence. We find $`a^3\mathrm{\Sigma }=0.0032(4)`$ and $`C_2=0.914(8)`$.

Clearly, the formulae derived in $`\chi `$PT give a very good description of the numerical data. The infinite volume condensate that we extract from this fit in physical units is $`\mathrm{\Sigma }(\mu 1.5\text{ GeV})=(221_9^{+8}\text{ MeV})^3`$, up to a multiplicative renormalization constant, which has not been computed yet for Neuberger’s operator. We stress that the quoted error on the condensate is purely statistical. It does not include, for instance, the expected systematic errors from finite lattice spacing effects, nor the possible contributions from higher orders in chiral perturbation theory of $`O(F_\pi L)^2`$. An additional cautoniary remark is that the statistics for the largest volume, $`12^4`$, is rather small as indicated by the large statistical error. We plan to increase the statistics in the future and include in the analysis also higher topologies, which are more frequent at larger volumes.

The condensate we obtain is quite close to that reported in ref.  using Wilson fermions and a different method. However, a meaningful comparison can only be made when our systematic errors are quantified and the multiplicative renormalization included.

Recently the authors of  also studied the quark condensate as a function of quark mass and volume, using Neuberger’s operator. However, a comparison with the predictions of q$`\chi `$PT in fixed topological sectors was not attempted and no definite conclusion on the existence or value of the infinite volume condensate was reached.

According to Random Matrix Theory (RMT), the value of $`\mathrm{\Sigma }`$ may also be extracted from the distribution of the lowest non-zero eigenvalue $`\lambda _{\mathrm{min}}(D_\mathrm{N})`$, defined by the square root of the lowest non-zero eigenvalue of $`D_N^{}D_N`$ at zero quark mass. We have only gathered a reasonable statistics for the smaller lattice of size $`8^4`$. In the topological charge $`\nu =0,1`$ sectors, the corresponding distributions are given by :

$$P_{\nu =0}(z)=\frac{z}{2}e^{\frac{1}{4}z^2},P_{\nu =\pm 1}(z)=\frac{z}{2}I_2(z)e^{\frac{1}{4}z^2},$$
(12)   
where $`z\lambda _{\mathrm{min}}(D_\mathrm{N})\mathrm{\Sigma }V`$. Recently, the authors of ref.  found very good agreement with these distributions in very small lattices. Inserting our value of $`a^3\mathrm{\Sigma }=0.0032(4)`$ in the distribution for zero topology of eq. (12), we get for the expectation value of this eigenvalue $`\lambda _{\mathrm{min}}(D_\mathrm{N})=0.135(15)`$ (where the error comes from the statistical error in the condensate), while from our data on the $`8^4`$ lattice (with 41 topology zero configurations) we obtain $`\lambda _{\mathrm{min}}(D_\mathrm{N})=0.170(12)`$. In topology one sectors, the expected value for the $`\lambda _{\mathrm{min}}(D_\mathrm{N})=0.237(17)`$. From the data, again obtained on the $`8^4`$ lattice ,with an accummulated statistics of 29 configurations, we obtain $`\lambda _{\mathrm{min}}(D_\mathrm{N})=0.218(13)`$. In addition, we generated a sample of eigenvalues according to the distribution of eq. (12) for $`\nu =\pm 1`$ with $`a^3\mathrm{\Sigma }=0.0032(4)`$ and a statistics identical to that of the corresponding simulation. The resulting mean value of $`\lambda _{\mathrm{min}}(D_\mathrm{N})`$ and error are fully compatible with those given by our simulation, indicating that our data do not suffer from autocorrelation effects. This provides a nice cross check on the value of $`\mathrm{\Sigma }`$ obtained from our finite-size scaling analysis.

We finally briefly comment on the difficulty in measuring the condensate in the topology zero sector. This is due to the logarithmic enhancement, which can be shown to originate from the contribution of the single lowest eigenvalue, $`\lambda _{\mathrm{min}}(D_N)`$ to the condensate, if the distribution of this eigenvalue is that given by RMT in eq. (12). As is clear from eq. (12), the distribution of the lowest eigenvalue does not have a gap and it is easy to check that the contribution of this single eigenvalue to the condensate, $`\mathrm{\Sigma }_{\mathrm{min}}`$, has a logarithmic IR divergence in the sector of zero topology:

$`\mathrm{\Sigma }_{\mathrm{min}}={\displaystyle \frac{1}{V}}{\displaystyle 𝑑zP_{\nu =0}(z)\frac{2m}{m^2+(z/\mathrm{\Sigma }V)^2}}=m\mathrm{\Sigma }^2V\mathrm{log}m\mathrm{\Sigma }V+O(m\mathrm{\Sigma }^2V),`$ (13)   
which reproduces exactly the logarithmic dependence in eq. (8). We have generated a sample of eigenvalues according to eq. (12) with $`a^3\mathrm{\Sigma }=0.0032`$. In doing so we find that reconstructing the logarithmic behaviour of the $`\nu =0`$ condensate requires a statistics much larger than the one available to us. This is reflected in the results from our actual simulation where we find that the $`\nu =0`$ condensate on the $`8^4`$ lattice displays very large statistical errors. For this reason, we have not included these data in our determination of the scalar condensate.

Conclusion

Chiral perturbation theory assumes that chiral symmetry is spontaneously broken in QCD. Under this assumption it provides, for small enough values of the quark mass and large enough volumes, the mass and volume behaviour of the scalar condensate. This behaviour is determined to lowest order by only one free parameter, namely $`\mathrm{\Sigma }`$, the scalar condensate in infinite volume and for zero quark mass.

Recent developments in lattice QCD have revealed that, contrary to a long-standing belief, chiral symmetry can be realized on the lattice. This theoretical advance is connected to the Ginsparg–Wilson relation. Neuberger proposed a particular operator that satisfies this relation and we have used this operator in our numerical work.

Although Neuberger’s operator is very difficult to treat numerically, it can be used in practice. In this work we computed the scalar condensate on lattices of various sizes and for a number of quark masses, in the regime of chiral symmetry restoration. The results of this numerical computation are shown in fig. 1, where we confront our numerical data with the finite volume and mass behaviour predicted by quenched chiral perturbation theory. Obviously, chiral perturbation theory describes the numerical data well, providing evidence for the spontaneous breaking of chiral symmetry.

Although our results are very encouraging, a number of cautionary remarks have to be made. The lattices we used are rather small and it would be desirable to probe the system further on larger lattices. In addition, it would be important to repeat the calculation at a larger value of $`\beta `$ to estimate the lattice spacing effects. Here we were only able to determine a value of the scalar condensate up to a multiplicative renormalization constant, which would clearly be needed for quoting a physical value. Finally, all our results are obtained in the quenched approximation but, given the complexity of Neuberger’s operator, it would be very difficult to go beyond this approximation.

During the completion of this work, a paper  using Neuberger’s operator to compute the scalar condensate appeared. The data presented in this paper are, however, taken in the strong coupling regime on only one lattice (with small L/$`a`$=4) and can hence not directly be compared to our work.

Acknowldegments We thank Martin Lüscher, Massimo Testa and Peter Weisz for many useful and stimulating discussions and T. Wettig for useful correspondance about RMT. We acknowledge the computer centre at NIC, (Jülich) and CIEMAT (Madrid) for providing computer time and technical support. L.L. acknowledges support from the EEC through the TMR network EEC-CT98-00169.