format: text

text:
# Lemma 5.1

CERN-TH/2002-350

DESY 02-212

CPT-2002/P.4458

Numerical techniques for lattice QCD in the $`ϵ`$–regime

L. Giusti1,2, C. Hoelbling2, M. Lüscher1, H. Wittig3

1 CERN, Theory Division, CH-1211 Geneva 23, Switzerland

2 Centre Physique Théorique, CNRS, Case 907, Luminy, F-13288 Marseille, France

3 DESY, Theory Group, Notkestrasse 85, D-22603 Hamburg, Germany

Abstract

In lattice QCD it is possible, in principle, to determine the parameters in the effective chiral lagrangian (including weak interaction couplings) by performing numerical simulations in the $`ϵ`$–regime, i.e. at quark masses where the physical extent of the lattice is much smaller than the Compton wave length of the pion. The use of a formulation of the lattice theory that preserves chiral symmetry is attractive in this context, but the numerical implementation of any such approach requires special care in this kinematical situation due to the presence of some very low eigenvalues of the Dirac operator. We discuss a set of techniques (low-mode preconditioning and adapted-precision algorithms in particular) that make such computations numerically safe and more efficient by a large factor.

1. Introduction

At low energies the physics of the light pseudo-scalar mesons can be described by an effective chiral $`\sigma `$–model with $`\mathrm{SU}(3)_{\mathrm{left}}\times \mathrm{SU}(3)_{\mathrm{right}}`$ symmetry group. Purely strong interaction phenomena and electroweak transitions (such as the non-leptonic kaon decays) are both covered by the effective theory if the appropriate interaction terms are included in the lagrangian. The challenge is then to compute the associated coupling constants from the underlying field theory of the strong and the electroweak interactions.

In the effective chiral theory the physical amplitudes are obtained in the form of an asymptotic expansion in powers of the quark masses and the meson momenta. The coupling constants appear in the coefficients of this expansion, and to determine their values it is therefore essential to choose a computational strategy where the limit of small masses and momenta can be safely reached.

The so-called $`ϵ`$–regime of QCD \[?–?\], combined with a formulation of lattice QCD in which chiral symmetry is exactly preserved \[?–?\], provides a framework that may prove to be particularly suitable in this context. Some encouraging results have in fact already been obtained along these lines in the case of the quark condensate \[?–?\] which is one of the free parameters in the chiral lagrangian. In general numerical simulations in the $`ϵ`$–regime are technically demanding, however, because the lattice Dirac operator tends to be ill-conditioned. The problem is linked to the spontaneous breaking of chiral symmetry \[?,?–?\] and is hence present independently of how precisely the lattice theory is defined.

In this paper we discuss a set of fast numerical methods for the computation of the quark propagator and the zero-modes of the Dirac operator. Being efficient is clearly very important here, but we wish to emphasize that an ill-conditioned system also gives rise to questions of numerical stability and accuracy \[?\] that must be properly dealt with.

Although some of our techniques are more widely applicable, we focus on the case of Neuberger’s lattice Dirac operator \[?\] in this paper. After a brief review of the $`ϵ`$–regime in sect. 2, we then first discuss a particular numerical implementation of this operator, following previous work on the subject \[?,?–?\], which is uniformly accurate to a specified level of precision (sects. 3–5). The ability to rigorously control the approximation error is seen to be very useful in sect. 6, where we describe an efficient way to calculate the index of the Dirac operator, and it also provides the basis for the adapted-precision inversion algorithm discussed in sect. 9. Before this, in sects. 7 and 8, we show how the zero-mode contribution to the quark propagator can be safely separated and introduce a method, referred to as low-mode preconditioning, that takes care (to some extent) of the potentially very large condition number of the Dirac operator in the subspace orthogonal to the zero-modes.

2. The $`ϵ`$–regime of QCD

In numerical simulations of lattice QCD the adjustable parameters are the gauge coupling, the quark masses and the lattice size. Usually a euclidean $`T\times L^3`$ lattice is considered, with periodic boundary conditions, where $`T`$ and $`L`$ refer to the time and the space directions respectively. We shall also assume in this paper that the lattice Dirac operator satisfies the Ginsparg–Wilson relation so that chiral symmetry is exactly preserved on the lattice \[?\]. The $`ϵ`$–regime is then reached at fixed gauge coupling by scaling $`T`$ and $`L`$ to values much larger than the confinement radius and simultaneously the light quark masses $`m=m_u,m_d,m_s`$ to zero in such a way that the combination

$$\begin{array}{c}\text{ }x=m\mathrm{\Sigma }V,VTL^3,\text{(2.1)}\text{ }\hfill \end{array}$$

remains fixed. In this equation the parameter

$$\begin{array}{c}\text{ }\mathrm{\Sigma }=\underset{m0}{lim}\underset{V\mathrm{}}{lim}\left|\overline{u}u\right|\text{(2.2)}\text{ }\hfill \end{array}$$

is the (unrenormalized) $`u`$–quark condensate in infinite volume.$`\mathrm{}`$ For simplicity we do not consider the case of quenched or partially quenched QCD in this section. Moreover we stick to a continuum notation and drop all factors that formally converge to $`1`$ in the continuum limit.

2.1 Effective chiral theory

To leading order the purely strong interaction part of the euclidean action density of the effective theory reads

$$\begin{array}{c}\text{ }=\frac{1}{4}F^2\mathrm{tr}\left\{_\mu U^{}_\mu U\right\}\frac{1}{2}\mathrm{\Sigma }\mathrm{tr}\left\{UM+M^{}U^{}\right\},\text{(2.3)}\text{ }\hfill \end{array}$$

where $`U`$ denotes the SU(3)–valued chiral field and $`M=\text{diag}(m_u,m_d,m_s)`$ the quark mass matrix. The normalizations chosen here are such that the parameters $`F`$ and $`\mathrm{\Sigma }`$ coincide with the pion decay constant and the chiral condensate at tree-level of the chiral perturbation expansion.

In the chiral limit the effective theory is expected to provide a correct description of the low-energy properties of lattice QCD, up to lattice spacing effects, if $`F`$ and $`\mathrm{\Sigma }`$ are set to the proper values in units of the lattice spacing. These parameters can thus be determined by calculating any suitable correlation functions on the lattice at small masses and momenta and by comparing the results with the predictions of the effective theory. In particular, such computations may be performed in the $`ϵ`$–regime, where the quark mass term breaks chiral symmetry only very weakly and where the zero-momentum modes $`U(x)=V`$ dominate the partition function of the effective theory \[?\]. There are many different ways to study the matching between the effective theory and lattice QCD in this regime, most of which still need to be explored.

An example that we wish to quote for illustration is the two-point function of the left-handed current

$$\begin{array}{c}\text{ }J_\mu ^a=\overline{\psi }\lambda ^a\gamma _\mu P_{}\psi ,P_\pm \frac{1}{2}(1\pm \gamma _5),\text{(2.4)}\text{ }\hfill \end{array}$$

which is represented by

$$\begin{array}{c}\text{ }\widehat{J}_\mu ^a=\frac{1}{2}F^2\mathrm{tr}\left\{\lambda ^aU_\mu U^{}\right\}\text{(2.5)}\text{ }\hfill \end{array}$$

in the effective theory. Standard notations are being used here, where $`\psi `$ denotes the quark field with flavour components $`(u,d,s,\mathrm{})`$ and $`\lambda ^a`$ an SU(3) generator that acts on the flavour indices of the fields. At non-zero time separations and vanishing spatial momentum, this correlation function turns out to be proportional to $`F^2/T`$ \[?,?,?\] and it is thus expected to be particularly suitable for the calculation of the parameter $`F`$.

2.2 Weak-interaction couplings

The effective couplings associated with hadronic weak transitions can perhaps be computed in a similar way by matching correlation functions in the $`ϵ`$–regime. In the case of the CP conserving $`K\pi \pi `$ decays, for example, the weak interaction density that is obtained after integrating out the $`W`$–bosons is, to a good approximation, equal to a linear combination of \[?–?\]

$$\begin{array}{c}\text{ }𝒪^\pm =\left\{(\overline{s}\gamma _\mu P_{}u)(\overline{u}\gamma _\mu P_{}d)\pm (\overline{s}\gamma _\mu P_{}d)(\overline{u}\gamma _\mu P_{}u)\right\}\left(uc\right),\text{(2.6)}\text{ }\hfill \\ \\ \text{ }𝒪_m=\left(m_c^2m_u^2\right)\left\{m_d\left(\overline{s}P_+d\right)+m_s\left(\overline{s}P_{}d\right)\right\}.\text{(2.7)}\text{ }\hfill \end{array}$$

The operator $`𝒪_m`$ is usually omitted from this list, because it does not contribute in transition matrix elements where the total energy-momentum is conserved. Mixings of $`𝒪^\pm `$ and $`𝒪_m`$ with other operators are, incidentally, excluded by the exact symmetries of the lattice theory.

In the effective chiral theory, the operators $`𝒪^\pm `$ and $`𝒪_m`$ are expected to be represented by certain linear combinations of the fields \[?\]

$$\begin{array}{c}\text{ }\widehat{𝒪}_1=\frac{1}{8}F^4\{(U_\mu U^{})_{13}(U_\mu U^{})_{21}+\text{ }\hfill \\ \\ \text{ }(U_\mu U^{})_{23}(U_\mu U^{})_{11}+\frac{1}{5}(_\mu U_\mu U^{})_{23}\},\text{(2.8)}\text{ }\hfill \\ \\ \text{ }\widehat{𝒪}_2=\frac{1}{8}F^4(_\mu U_\mu U^{})_{23},\text{(2.9)}\text{ }\hfill \\ \\ \text{ }\widehat{𝒪}_3=\frac{1}{2}F^2\mathrm{\Sigma }(UM+M^{}U^{})_{23},\text{(2.10)}\text{ }\hfill \end{array}$$

and our task is then to determine the dimensionless coefficients in these linear combinations. An obvious possibility in the $`ϵ`$–regime is to match three-point correlation functions such as

$$\begin{array}{c}\text{ }C^{ab}(x_0,y_0)=_0^L\mathrm{d}^3x\mathrm{d}^3yJ_\mu ^a(x)𝒪^\pm (0)J_\mu ^b(y)\text{(2.11)}\text{ }\hfill \end{array}$$

plus a similar correlation function with only one left-handed current (so as to be able to disentangle the contributions of the octet operators $`\widehat{𝒪}_2`$ and $`\widehat{𝒪}_3`$). Explicit computations in chiral perturbation theory suggest that these correlation functions are indeed well-behaved and that, in principle, the desired coefficients can be calculated along these lines \[?\].

2.3 Spectrum of the Dirac operator

The last ten years or so have seen a remarkable development, involving both random matrix and chiral perturbation theory, that culminated in the analytical determination of the distributions of the low-lying eigenvalues of the Dirac operator in the $`ϵ`$–regime \[?,?–?\] (for a review see ref. \[?\]). An important insight gained in the course of this work is that the eigenvalues scale like $`(\mathrm{\Sigma }V)^1`$ at large volumes. In particular, the spectral gap in the subspace orthogonal to the chiral zero-modes is, on average, roughly of this size.

For the calculation of the quark propagator in the $`ϵ`$–regime, the presence of some very low eigenvalues (an effect that has been confirmed numerically on small lattices \[?,?,?\]) is a source of difficulty. To make this a bit more concrete, let us consider a $`32^4`$ lattice with spacing $`a=0.1\mathrm{fm}`$ and let us assume that $`\mathrm{\Sigma }`$ is around $`(250\mathrm{MeV})^3`$ on this lattice. From the discussion above we then deduce that the expectation value of the spectral gap is about $`5\times 10^4`$ in lattice units. Note that the quark masses are scaled proportionally to $`(\mathrm{\Sigma }V)^1`$ in the $`ϵ`$–regime and thus do not provide a strong infrared cutoff on the spectrum.

The computation of the quark propagator, using the conjugate gradient algorithm for example, is perhaps still possible under these conditions, but may require a very large number of iterations. Moreover, the accuracy of the solution vector that can be reached on any given machine may be unsatisfactory, especially if the propagator is calculated via the so-called normal equations (as is the case if the conjugate gradient algorithm is used) \[?\]. Effectively these technical problems set a limit on the space-time volumes where numerical simulations can be performed and thus on how close to the chiral point one can get in practice.

3. Lattice Dirac operator

We set up the lattice theory as usual on a finite four-dimensional lattice with spacing $`a`$ and periodic boundary conditions. Our notational conventions are summarized in appendix A. Following refs. \[?–?\] we consider a formulation of lattice QCD where the lattice Dirac operator $`D`$ satisfies the Ginsparg–Wilson relation

$$\begin{array}{c}\text{ }\gamma _5D+D\gamma _5=\overline{a}D\gamma _5D,D^{}=\gamma _5D\gamma _5,\text{(3.1)}\text{ }\hfill \end{array}$$

for some positive constant $`\overline{a}`$ proportional to $`a`$. In the presence of a given background gauge field, the propagator of a quark with mass $`m`$ is then given by

$$\begin{array}{c}\text{ }\psi (x)\overline{\psi }(y)=\left\{D_m^1\right\}(x,y),\text{(3.2)}\text{ }\hfill \end{array}$$

where $`\{\mathrm{}\}(x,y)`$ stands for the kernel in position space of the operator in the curly bracket and

$$\begin{array}{c}\text{ }D_m=(1\frac{1}{2}\overline{a}m)D+m\text{(3.3)}\text{ }\hfill \end{array}$$

is the massive lattice Dirac operator. We shall always assume that $`0\overline{a}m2`$.

As already mentioned in sect. 1, we concentrate on the case of the Neuberger–Dirac operator in this paper, even though many of the techniques that will be described are more generally applicable. Explicitly this operator is given by \[?\]

$$\begin{array}{c}\text{ }D=\frac{1}{\overline{a}}\left\{1+\gamma _5\text{sign}(Q)\right\},\text{(3.4)}\text{ }\hfill \\ \\ \text{ }Q=\gamma _5\left(aD_\mathrm{w}1s\right),\overline{a}=\frac{a}{1+s},\text{(3.5)}\text{ }\hfill \end{array}$$

where $`s`$ is an adjustable parameter in the range $`|s|<1`$ and $`D_\mathrm{w}`$ denotes the standard Wilson–Dirac operator (appendix A).

4. Minmax polynomial approximation

For the numerical implementation of the Neuberger–Dirac operator we need to find an algorithm that evaluates $`\text{sign}(Q)\eta `$ on any given fermion field $`\eta `$ to a specified precision. In the present section we take a first step in this direction by constructing an optimal polynomial approximation to $`\text{sign}(x)`$ in a range that excludes a small interval around the origin. More precisely we are looking for a polynomial $`P(y)`$ of degree $`n`$ that minimizes the error

$$\begin{array}{c}\text{ }\delta =\underset{ϵy1}{\mathrm{max}}\left|h(y)\right|,h(y)1\sqrt{y}P(y),\text{(4.1)}\text{ }\hfill \end{array}$$

for specified $`ϵ>0`$. In the range $`\sqrt{ϵ}|x|1`$ the function $`xP(x^2)`$ then approximates $`\text{sign}(x)`$ uniformly with a maximal deviation $`\delta `$.

4.1 Existence and uniqueness

Polynomial approximations that minimize the maximal relative error are referred to as minmax polynomials. They can be shown to exist and to be uniquely determined provided the function that is to be approximated does not vanish. Moreover, practical methods have been devised, based on the characteristic properties of the associated error function, that allow one to compute them reliably (see ref. \[?\], for example, or any other good book on approximation theory).

For reasons of numerical stability, it is advantageous to represent the polynomial in the form of a series

$$\begin{array}{c}\text{ }P(y)=\underset{k=0}{\overset{n}{}}c_kT_k(z),z=(2y1ϵ)/(1ϵ),\text{(4.2)}\text{ }\hfill \end{array}$$

of Chebyshev polynomials (appendix A). Our task is then to adjust the coefficients $`c_k`$ so that the deviation (4.1) is minimized.

4.2 Error bounds

The characteristic property alluded to above of the minmax polynomial is that the error function $`h(y)`$ has (at least) $`n+2`$ extrema of equal height and alternating sign, i.e. $`h(y)`$ is oscillating around zero with constant amplitude (fig. 1). Evidently it is not possible to know in advance where these extrema are located.

Let us now choose an arbitrary sequence

$$\begin{array}{c}\text{ }ϵy_0<y_1<\mathrm{}<y_{n+1}1\text{(4.3)}\text{ }\hfill \end{array}$$

of $`n+2`$ values. We can then find a polynomial such that

$$\begin{array}{c}\text{ }h(y_l)=(1)^lu\text{for all}l=0,\mathrm{},n+1.\text{(4.4)}\text{ }\hfill \end{array}$$

In fact these are just $`n+2`$ linear equations in the unknowns $`c_0,\mathrm{},c_n,u`$. This polynomial has the alternating property on the chosen set of points, but not necessarily so on the approximation interval. It can be shown, however, that

$$\begin{array}{c}\text{ }\left|u\right|\underset{ϵy1}{\mathrm{max}}\left|h_{}(y)\right|\underset{ϵy1}{\mathrm{max}}\left|h(y)\right|\text{(4.5)}\text{ }\hfill \end{array}$$

where $`h_{}(y)`$ denotes the error function of the minmax polynomial. The solution of the linear equations (4.4) thus yields a nearly optimal polynomial if these bounds are tight.

Fig. 1. Minmax polynomial approximation of sign$`(x)`$ for $`ϵ=0.0025`$ and $`n=22`$ (full line). In the range $`\sqrt{ϵ}|x|1`$ the relative error of this approximation is $`5\%`$, while for an error of $`0.1\%`$ a polynomial of degree $`n=57`$ is required (dashed line).

4.3 Maehly’s exchange algorithm

The basic idea of this algorithm is to start from some well-chosen initial set of points (4.3) and to adjust them subsequently so as to tighten the bounds (4.5). This can be achieved, for example, by choosing the new points to be the nearby extremal points of the current error function.

A good initial set of points for this recursion are the extremal points of the Chebyshev polynomial $`T_{n+1}(z)`$. Only few iterations are then required until a polynomial is obtained that is very nearly optimal. There is in fact little to be gained by performing further iterations once the lower and the upper bound in (4.5) differ by no more than (say) $`0.01\left|u\right|`$. Evidently this procedure is extremely robust, because the quality of the current polynomial (with respect to the best possible approximation) can be rigorously controlled.

4.4 Implementation details

Inspection shows that the linear system (4.4) is well-conditioned and it can thus be solved straightforwardly by $`QR`$ factorization, using Householder transformations, followed by backward substitution (see ref. \[?\] for example).

To localize the extrema of $`h(y)`$, a few bisection steps are sufficient in general, since there is no need to be very precise at this point. An important detail here is that many digits are lost when calculating the error function, even if the Clenshaw recursion is used to evaluate the Chebyshev series (4.2) \[?\]. This becomes a problem when $`\delta /12n`$ approaches the machine precision, and with standard 64 bit arithmetic the approximation accuracy $`\delta `$ that can be reached is thus limited to values greater than $`n\times 10^{15}`$ or so.

4.5 Efficiency of the approximation

The precision of the approximation to the sign function provided by the minmax polynomial scales approximately like

$$\begin{array}{c}\text{ }\delta =A\mathrm{e}^{bn\sqrt{ϵ}}\text{(4.6)}\text{ }\hfill \end{array}$$

where $`A=0.41`$ and $`b=2.1`$. At large values of $`n\sqrt{ϵ}`$, this empirical formula is fairly accurate, but it slightly over-estimates the actual value of $`\delta `$ when the degree of the polynomial is less than 30 or so. With respect to the Chebyshev approximation that was used in ref. \[?\], the precision is improved by about a factor of $`2`$.

We finally mention that, on a current PC processor and for degrees $`n400`$, only a few seconds are required to obtain the coefficients $`c_k`$ of the minmax polynomial. It is hence practical to recompute these polynomials whenever they are needed.

5. Uniform approximation of the Neuberger–Dirac operator

The minmax polynomial $`P_{n,ϵ}(y)`$ of degree $`n`$ on the interval $`[ϵ,1]`$ yields an approximation of the operator $`\text{sign}(Q)`$ through

$$\begin{array}{c}\text{ }\text{sign}(Q)XP_{n,ϵ}(X^2),XQ/Q.\text{(5.1)}\text{ }\hfill \end{array}$$

If $`ϵ`$ is chosen so that $`Q^2ϵQ^2`$, the error in this formula is an operator with norm less than or equal to $`\delta `$. In other words, the approximation error is always bounded by $`\delta \eta `$, uniformly in the field $`\eta `$ to which the operator is applied.

This method is straightforward but actually not recommendable if $`Q^2`$ has some exceptionally low eigenvalues (as is often the case in practice \[?\]). It is far more efficient under these conditions to first separate the few lowest modes and to treat them exactly. Evidently this should be done in such a way that the total approximation error remains under control.

5.1 Low-mode projectors

The spectrum of $`Q`$ in the vicinity of the origin can be reliably determined by minimizing the Ritz functional of $`Q^2`$ \[?,?\]. This technique also yields an approximation to the associated eigenvectors, and we now need to discuss by how much these vectors deviate from the true eigenvectors.

So let us assume that a specified number $`l`$ of approximate eigenvectors has been computed. We denote by $`V`$ the linear space spanned by these vectors and by $`P`$ the corresponding orthonormal projector. It is then trivial to calculate the parameter

$$\begin{array}{c}\text{ }\varrho =\underset{vV,v=1}{\mathrm{max}}(1P)Qv\text{(5.2)}\text{ }\hfill \end{array}$$

which measures the deviation of $`V`$ from being an exact eigenspace of $`Q`$. The importance of this parameter is clarified by

###### Lemma 5.1

Let $`\nu _1,\mathrm{},\nu _l`$ be the eigenvalues of $`PQP`$ in the subspace $`V`$. Then there are $`l`$ linearly independent eigenvectors of $`Q`$ with eigenvalues $`\lambda _1,\mathrm{},\lambda _l`$ such that $`|\nu _k\lambda _k|\varrho `$ for all $`k=1,\mathrm{},l`$.

The proof of the lemma (which holds independently of how $`V`$ was obtained) is given in appendix B.

In the following we shall take it for granted that the eigenvalues $`\nu _k`$ are separated from zero and from the rest of the spectrum of $`Q`$ (all eigenvalues $`\lambda `$ except $`\lambda _1,\mathrm{},\lambda _l`$) by a distance greater than $`\varrho `$. So far we have in fact never encountered a situation where this condition could not be satisfied, and we thus omit any further discussion of this point. The presence of a spectral gap around zero implies that the subsets of positive and negative eigenvalues can be identified without any numerical ambiguity. If we introduce the eigenvectors

$$\begin{array}{c}\text{ }u_kV,PQu_k=\nu _ku_k,(u_k,u_j)=\delta _{kj},\text{(5.3)}\text{ }\hfill \end{array}$$

the associated orthonormal projectors are given by

$$\begin{array}{c}\text{ }P_+=\underset{\nu _k>0}{}u_k(u_k)^{}\text{and}P_{}=\underset{\nu _k<0}{}u_k(u_k)^{}\text{(5.4)}\text{ }\hfill \end{array}$$

respectively. In the same way we may define the projectors $`(P_\pm )_{\mathrm{exact}}`$ to the subspaces spanned by the corresponding exact eigenvectors of $`Q`$ (whose existence is guaranteed by lemma 5.1) with positive and negative eigenvalues $`\lambda _k`$.

We may now ask how accurately the computed projectors $`P_\pm `$ approximate the exact projectors $`(P_\pm )_{\mathrm{exact}}`$. The answer to this question depends on the size of the residues

$$\begin{array}{c}\text{ }\varrho _k=(Q\nu _k)u_k\text{(5.5)}\text{ }\hfill \end{array}$$

and also on the distances between the eigenvalues of $`Q`$. A small value of $`\varrho _k`$ alone does, in fact, not exclude sizeable mixings of $`u_k`$ with several eigenvectors of $`Q`$ if these have eigenvalues that are within a distance $`\varrho _k`$ of $`\nu _k`$.

Eventually we are interested in estimating the deviation of the projectors rather than that of the individual eigenvectors, and what counts in this case is the distance $`d_k`$ of $`\nu _k`$ from the exact spectrum of $`Q`$ in the subspace that is orthogonal to the range of $`(P_+)_{\mathrm{exact}}`$ if $`\nu _k>0`$ or $`(P_{})_{\mathrm{exact}}`$ if $`\nu _k<0`$. The quality of the approximation is then controlled by the parameters

$$\begin{array}{c}\text{ }\kappa _\pm ^2=\underset{\pm \nu _k>0}{}\varrho _k^2/d_k^2,\kappa _\pm >0,\text{(5.6)}\text{ }\hfill \end{array}$$

as the following lemma shows.

###### Lemma 5.2

Provided the spectral conditions listed above are satisfied, the bound

$$\begin{array}{c}\text{ }P_\pm (P_\pm )_{\mathrm{exact}}\frac{\kappa _\pm (1+2\kappa _\pm )}{12\kappa _\pm (1+2\kappa _\pm )}\text{(5.7)}\text{ }\hfill \end{array}$$

holds where it is assumed that $`2(l+1)\kappa _\pm (1+2\kappa _\pm )<1`$.

The proof of the lemma is given in appendix C. Here we only note that the parameters $`\kappa _\pm `$ are very small in the cases of interest. The right-hand side of eq. (5.7) is then practically equal to $`\kappa _\pm `$ and the condition on the last line is trivially fulfilled.

5.2 Approximation formula and error bound

It is now straightforward to write a program that computes the projectors to a specified accuracy. In this program the number of low modes to be included in the projectors should be determined dynamically in such a way that the spectral distance from the other modes is not accidentally very small. The parameters $`\kappa _\pm `$ can then be estimated without difficulty and the minimization of the Ritz functional is stopped when the desired level of precision is reached.

Once the projectors are computed, eq. (5.1) may be replaced by

$$\begin{array}{c}\text{ }\text{sign}(Q)P_+P_{}+\left(1P_+P_{}\right)XP_{n,ϵ}(X^2),XQ/Q,\text{(5.8)}\text{ }\hfill \end{array}$$

where $`ϵ`$ is set to a value equal to (or perhaps slightly less than) the smallest eigenvalue of $`Q^2/Q^2`$ in the sector orthogonal to the low modes. The properties of the minmax polynomial, lemma 5.2 and simple triangle inequalities then imply that the associated approximation $`\stackrel{~}{D}_m`$ to the massive Neuberger–Dirac operator satisfies

$$\begin{array}{c}\text{ }\stackrel{~}{D}_mD_m\frac{1}{a}(1+s\frac{1}{2}am)\left\{2(\kappa _++\kappa _{})+\delta \right\}\text{(5.9)}\text{ }\hfill \end{array}$$

up to terms proportional to $`\kappa _\pm \delta `$ and $`\kappa _\pm ^2`$ (from now on these will be neglected).

5.3 Miscellaneous remarks

In the tests that we have performed we never had any difficulty in reaching accuracy levels $`\kappa _\pm `$ of $`10^8`$ or even $`10^{10}`$. Moreover, the actual approximation error is significantly smaller than suggested by lemma 5.2, by a factor $`10`$ at least. Since the calculation of the projectors consumes a relatively small amount of computer time, they will normally be obtained to a precision that is estimated to be sufficient, by a wide margin, for the applications that one has in mind. The approximation error (5.9) is then essentially determined by the degree $`n`$ of the minmax polynomial.

A numerically safe method to evaluate the polynomial in eq. (5.8) is provided by the Clenshaw recursion \[?\]. The accumulated rounding errors are then on the order of $`n`$ times the machine precision. In particular, if standard 32 bit floating-point arithmetic is used (which can be profitable at intermediate stages of the algorithms described later), a rounding error of about $`5n\times 10^7`$ should be added to the theoretical approximation error $`\delta `$.

We finally note that the degree of the minmax polynomial will usually be significantly larger than the number of the separated low modes. The application of the projectors in eq. (5.8) then requires only a small fraction of the total execution time and the associated rounding errors can also be safely neglected.

6. Calculation of the index of the Dirac operator

6.1 Preliminaries

The operator $`D^{}D`$ commutes with $`\gamma _5`$ and thus leaves the subspaces of fermion fields with definite chirality invariant. Using the Ginsparg–Wilson relation, it can be shown that the action of the operator in these subspaces is given by the hermitian operators

$$\begin{array}{c}\text{ }D^\pm =P_\pm DP_\pm \text{(6.1)}\text{ }\hfill \end{array}$$

up to a proportionality factor equal to $`2/\overline{a}`$. Moreover, the spectrum of the low-lying eigenvalues of $`D^+`$ and $`D^{}`$ is exactly the same, including degeneracies, apart from the fact that there can be a surplus of zero-modes in one of the sectors.

It happens only accidentally (with probability zero strictly speaking) that there are zero-modes with both chiralities, and we shall thus exclude this exceptional case in the following. The index $`\nu `$ of the Dirac operator is then given by

$$\begin{array}{c}\text{ }\nu =\sigma n_0\text{(6.2)}\text{ }\hfill \end{array}$$

where $`\sigma `$ denotes the chirality of the sector that contains the zero-modes and $`n_0`$ the number of these modes.

6.2 Zero-mode counting

The low-lying eigenvalues of $`D^+`$ and $`D^{}`$ can, in principle, be found by minimizing the associated Ritz functionals \[?,?\]. In terms of computer time such calculations tend to be very expensive, however, since each application of the Neuberger–Dirac operator to a given fermion field involves from say 50 to several 100 applications of $`Q^2`$. Any numerical method that requires a precise computation of the eigenvalues should hence be avoided.

The strategy that we propose is to run the minimization program in both chirality sectors simultaneously in a controlled way, where the precision is improved in steps by a specified factor. After every iteration the current estimate of the lowest eigenvalue in each sector provides a rigorous upper bound on the gap above zero in the spectrum of $`D^+`$ and $`D^{}`$ respectively. The program then proceeds to lower the larger bound and continues in this manner until the estimated relative error on any one of the eigenvalues in the two sectors is determined to be less than $`10\%`$. At this point we know that the zero-modes (if any) must be in the other chirality sector.

To count the zero-modes the minimization program must now be run once more in that sector, improving the precision on the calculated eigenvalue in steps as before. As soon as the eigenvalue drops below the spectral gap in the other sector (which was determined to an accuracy of $`10\%`$ in the previous calculation), there has to be at least one zero-mode, independently of the current level of precision, because the Ritz functional always provides an upper bound on the true eigenvalue. We can then restart the minimization program, requiring a second eigenvalue to be computed, and so on.

This procedure terminates when an eigenvalue is found that is definitely above zero, which will be the case if the relative error is estimated to be less than $`10\%`$ for example. It should be emphasized that in this way we never require any eigenvalues to be calculated to high precision and yet the zero-mode counting is rigorously correct.

6.3 Using reduced-precision operators

In the discussion above we have implicitly assumed that the approximation error (5.9) of the numerical representation $`\stackrel{~}{D}^\pm `$ of the operators $`D^\pm `$ can be neglected. An important point to note is that the associated error on the eigenvalues is at most

$$\begin{array}{c}\text{ }\omega =\frac{1}{a}(1+s)\left\{2(\kappa _++\kappa _{})+\delta \right\},\text{(6.3)}\text{ }\hfill \end{array}$$

and the same is true for the gradients of the Ritz functionals of $`D^\pm `$. We may, therefore, start the minimization of the Ritz functionals using a relatively poor approximation of the Dirac operator, and then gradually decrease $`\delta `$ so that the error $`\omega `$ is always smaller than the current magnitude of the gradients.

Apart from being in control of the numerical errors, an obvious advantage of this procedure is that the degree of the minmax polynomial is not larger than absolutely necessary, at all stages of the calculation. The total computational effort is thus minimized. Moreover, a further acceleration may be achieved by using single-precision arithmetic in the programs for the Neuberger–Dirac operator as long as $`\delta `$ is not too small (cf. subsect. 5.3). On current PC processors this saves almost a factor of $`2`$ in execution time \[?\].

6.4 Error estimation

In the algorithm described in subsect. 6.2, the error on the calculated eigenvalues must be estimated in a reliable way. For any given eigenvalue, the error is a sum of two contributions, one coming from the numerical approximation of the lattice Dirac operator and the other from the fact that, in each step of the procedure, the minimum of the relevant Ritz functional is only found to some specified accuracy.

A reliable but often poor estimate of the latter error is provided by the gradient of the Ritz functional \[?,?\]. However, for reasons of numerical stability it is in any case advisable to compute an additional eigenvalue in each sector, with possibly reduced precision \[?\], and the Temple inequality may then be used to obtain a more accurate estimate of the error on the other eigenvalues \[?\].

7. The quark propagator: general strategies

The presence of chiral zero-modes (which is the normal case in large volumes) complicates the calculation of the quark propagator, and it is mainly this issue that we wish to address in this section. We assume in the following that the quark mass $`m`$ is strictly positive and that a gauge field configuration is being considered where the zero-modes (if any) have chirality $`\sigma =+1`$.

7.1 Negative chirality components

The computation of the quark propagator amounts to solving the linear equation

$$\begin{array}{c}\text{ }D_m\psi =\eta \text{(7.1)}\text{ }\hfill \end{array}$$

for a set of source fields $`\eta `$. We first write the solution in the form

$$\begin{array}{c}\text{ }\psi =(D_mD_m)^1D_m\eta \text{(7.2)}\text{ }\hfill \end{array}$$

and note that the operator in brackets commutes with $`\gamma _5`$. The components of the solution vector with negative chirality are thus given by

$$\begin{array}{c}\text{ }P_{}\psi =(D_mD_m)^1P_{}D_m\eta .\text{(7.3)}\text{ }\hfill \end{array}$$

In this equation the inversion of $`D_mD_m`$ takes place in the chirality sector that does not contain the zero-modes, and the only potential difficulty is then that there can be some very low non-zero eigenvalues of $`D^{}D`$ in this subspace (cf. subsect. 2.4).

We shall come back to this problem in the next two sections, where we discuss the inversion of $`D_mD_m`$ in the negative chirality sector. For the time being, we merely assume that an efficient algorithm is available which allows us to invert the operator in this sector, to a specified precision and at arbitrarily small quark masses.

7.2 Positive chirality components

Once the negative chirality component $`P_{}\psi `$ has been computed, it can be inserted on the right-hand side of the equation

$$\begin{array}{c}\text{ }D_mP_+\psi =\eta D_mP_{}\psi \text{(7.4)}\text{ }\hfill \end{array}$$

for the positive chirality components $`P_+\psi `$. A possible expression for the solution of this equation is

$$\begin{array}{c}\text{ }P_+\psi =(P_+D_mP_+)^1\left\{P_+\eta P_+D_mP_{}\psi \right\}\text{(7.5)}\text{ }\hfill \end{array}$$

which involves an inversion of $`P_+D_mP_+`$ in the sector that contains the zero-modes.

It is now very important to note that

$$\begin{array}{c}\text{ }P_+D_mP_+=P_+\left\{\frac{1}{2}\overline{a}\left(1\frac{1}{2}\overline{a}m\right)D^{}D+m\right\}P_+\text{(7.6)}\text{ }\hfill \end{array}$$

is a positive operator in this sector whose condition number is at most $`2/\overline{a}m`$ and thus comparatively small. If the quark mass in lattice units is not too close to zero, the inversion in eq. (7.5) is hence well-conditioned and the solution can be computed straightforwardly, using the conjugate gradient algorithm for example.

7.3 Separation of the zero-modes

Arbitrarily small quark masses can also be handled but require a subtraction of the zero-mode component from the source field. Equation (7.5) then becomes

$$\begin{array}{c}\text{ }P_+\psi =\frac{1}{m}P_0P_+\eta +(P_+D_mP_+)^1\left\{(1P_0)P_+\eta P_+D_mP_{}\psi \right\},\text{(7.7)}\text{ }\hfill \end{array}$$

where $`P_0`$ denotes the projector to the subspace spanned by the zero-modes. Note that the field in the curly bracket is orthogonal to this subspace, and the operator inversion in eq. (7.7) thus remains well-defined even in the limit of a vanishing quark mass. The situation is then practically the same as in the negative chirality sector discussed in subsect. 7.1. In particular, the acceleration techniques described in the next two sections can be applied here too.

At this point we still need to explain how to compute the zero-modes (and thus the projector $`P_0`$). First recall that a basis of approximate zero-modes is obtained in the course of the calculation of the index of the Dirac operator (cf. subsect. 6.2). This basis may not be very accurate, but we can now easily improve on the accuracy by noting that

$$\begin{array}{c}\text{ }P_0\chi =\left(1\frac{1}{2}\overline{a}D\right)P_+\chi D(D^{}D)^1P_{}D^{}P_+\chi \text{(7.8)}\text{ }\hfill \end{array}$$

for any field $`\chi `$. The operator inversion in this formula is in the sector with negative chirality (where there are no zero-modes) which we have assumed to be feasible.

It is our experience that this procedure is numerically stable and yields accurate results when applied to the approximate zero-modes. More precisely, the error on the calculated zero-modes \[after the application of eq. (7.8)\] is estimated to be at most $`r+2\omega /|\lambda _1|`$ in this case, where $`\omega `$ is the accuracy of the approximation to the Neuberger–Dirac operator used in the numerator, $`\lambda _1`$ the first non-zero eigenvalue of $`D`$ and $`r`$ the relative residuum associated with the numerical inversion of $`D^{}D`$.$`\mathrm{}`$ An argumentation similar to the one at the beginning of appendix C shows that the deviation $`(1P_0)\chi `$ of any field $`\chi `$ from being an exact zero-mode is less than or equal to $`D\chi /|\lambda _1|`$. The approximation error can thus be rigorously controlled.

7.4 Left-handed propagator

The correlation functions that we have proposed in sect. 2 to probe the weak-interaction operators involve only left-handed quark and antiquark fields. To evaluate these correlation functions, it thus suffices to compute the purely left-handed components of the quark propagator, which are given by

$$\begin{array}{c}\text{ }\left\{P_{}(D_mD_m)^1P_{}D_mP_+\right\}(x,y).\text{(7.9)}\text{ }\hfill \end{array}$$

From this expression it is immediate that the complications arising from the zero-modes are absent in this case. In fact we only need to invert $`D_mD_m`$ in the chirality sector where there are no zero-modes. Note that if the zero-modes had the opposite chirality, we would write the propagator in the form

$$\begin{array}{c}\text{ }\left\{P_{}D_mP_+(D_mD_m)^1P_+\right\}(x,y)\text{(7.10)}\text{ }\hfill \end{array}$$

so that the inversion of $`D_mD_m`$ is again in the “good” sector.

Such correlation functions thus appear to be particularly attractive from the numerical point of view. Moreover, their behaviour in the chiral limit tends to be less singular than is generally the case \[?,?,?\].

8. Low-mode preconditioning

We now proceed to discuss the linear system

$$\begin{array}{c}\text{ }A\psi =\eta ,AD_mD_m,\text{(8.1)}\text{ }\hfill \end{array}$$

in the negative chirality sector, assuming as before that the zero-modes (if any) are in the other sector.

As already emphasized in sect. 2, the presence of some extremely low eigenvalues of $`A`$ forbids a straightforward application of the established inversion algorithms to solve the equation. To reduce the condition number of the system, an obvious possibility is to calculate the few lowest eigenvalues of $`A`$ and to compute the components of the solution vector along the associated eigenspace and its orthogonal complement separately. However, on large lattices this proposition is not practical, because the computation of the eigenvectors to the required level of precision is far too expensive in terms of computer time.

In this section we describe a version of low-mode preconditioning that works out even if the eigenvectors cannot be calculated very accurately. Effectively the algorithm compensates for the approximation errors through a simple block diagonalization that can be implemented exactly.

8.1 Approximate eigenvectors

So let us assume that $`e_1,\mathrm{},e_n`$ is a set of orthonormal Dirac fields with negative chirality that satisfy

$$\begin{array}{c}\text{ }Ae_k=\alpha _ke_k+r_k,(e_l,r_k)=0,r_k\omega _k\alpha _k,\text{(8.2)}\text{ }\hfill \end{array}$$

for all $`k,l`$ and some positive numbers $`\alpha _k`$ and $`\omega _k`$. Evidently, if the bounds $`\omega _k`$ on the residues $`r_k`$ are small, these fields can be regarded as approximate eigenvectors of $`A`$ with eigenvalues $`\alpha _k`$. In the following we shall also take it for granted that

$$\begin{array}{c}\text{ }(v,Av)\gamma v^2,\gamma \underset{k}{\mathrm{max}}\alpha _k,\text{(8.3)}\text{ }\hfill \end{array}$$

for all negative chirality vectors $`v`$ in the orthogonal complement of the vector space $`E`$ spanned by $`e_1,\mathrm{},e_n`$. This guarantees that $`e_1,\mathrm{},e_n`$ approximate the lowest $`n`$ eigenvectors of $`A`$ rather than an arbitrary set of eigenvectors.

It should be emphasized that no further assumptions need to be made and that the error bounds $`\omega _k`$, in particular, do not need to be extremely small (in most cases values of $`0.1`$ or so will do). In practice the vectors $`e_1,\mathrm{},e_n`$ can be computed by minimization of the Ritz functional of $`A`$ and stopping the program when the estimated relative errors of the calculated eigenvalues (as determined from the gradient of the Ritz functional \[?,?\]) drop below the specified values of $`\omega _k`$.$`\mathrm{}`$ Since the eigenvectors of $`A`$ are mass-independent, we may set $`m=0`$ in this calculation and use the same vectors $`e_1,\mathrm{},e_n`$ at all values of $`m`$. Only the eigenvalues $`\alpha _k`$ and the residual vectors $`r_k`$ then need to be recalculated if the quark mass changes.

8.2 Preconditioned system

We now introduce the projector $`P`$ to the subspace $`E`$ and note that the linear system (8.1) is equivalent to

$$\begin{array}{c}\text{ }(1P)\left\{AAP(PAP)^1PA\right\}(1P)\psi \text{ }\hfill \\ \\ \text{ }=(1P)\eta (1P)AP(PAP)^1P\eta ,\text{(8.4)}\text{ }\hfill \\ \\ \text{ }PAP\psi =P\eta PA(1P)\psi .\text{(8.5)}\text{ }\hfill \end{array}$$

The first of these equations determines the component $`(1P)\psi `$ of the solution in the subspace orthogonal to $`E`$. Once it has been calculated, the parallel component, $`P\psi `$, is obtained from the second equation.

To make this a little more explicit, it is helpful to introduce the field

$$\begin{array}{c}\text{ }\varphi (1P)\psi =\psi \underset{k=1}{\overset{n}{}}e_k(e_k,\psi ).\text{(8.6)}\text{ }\hfill \end{array}$$

Equation (8.4) may then be written as

$$\begin{array}{c}\text{ }(1P)A\varphi \underset{k=1}{\overset{n}{}}r_k\frac{1}{\alpha _k}(r_k,\varphi )=(1P)\eta \underset{k=1}{\overset{n}{}}r_k\frac{1}{\alpha _k}(e_k,\eta ).\text{(8.7)}\text{ }\hfill \end{array}$$

Together with the constraint $`P\varphi =0`$, this is a well-defined system whose solution immediately yields the complete field through

$$\begin{array}{c}\text{ }\psi =\varphi +\underset{k=1}{\overset{n}{}}e_k\frac{1}{\alpha _k}\left\{(e_k,\eta )(r_k,\varphi )\right\}.\text{(8.8)}\text{ }\hfill \end{array}$$

The scalar products in these formulae give rise to a computational overhead that is completely negligible compared to the effort required for the application of the operator $`A`$. This assumes, of course, that the decomposition (8.2) has been worked out before the conjugate gradient program is started and that it is possible to keep the vectors $`e_k`$ and $`r_k`$ in memory during the whole calculation.

8.3 Condition numbers

The operator on the left-hand side of eq. (8.7) is the difference of two operators,

$$\begin{array}{c}\text{ }(1P)A(1P)\underset{k=1}{\overset{n}{}}\frac{1}{\alpha _k}r_k(r_k)^{},\text{(8.9)}\text{ }\hfill \end{array}$$

that act in the orthogonal complement of $`E`$ and that are both positive. The lowest eigenvalue of the first operator is greater or equal to $`\gamma `$, while the norm of the second is bounded by $`_{k=1}^n(\omega _k)^2\alpha _k`$. This implies that the condition number of the preconditioned system (8.7) is at most

$$\begin{array}{c}\text{ }(4/\overline{a}^2)\times \left\{\gamma _{k=1}^n(\omega _k)^2\alpha _k\right\}^1,\text{(8.10)}\text{ }\hfill \end{array}$$

which is practically equal to $`4/\overline{a}^2\gamma `$ if the error bounds $`\omega _k`$ are smaller than $`0.2\times n^{1/2}`$ for example. Since the condition number of the original system is about $`4/\overline{a}^2`$ divided by the lowest eigenvalue of $`A`$, a reduction in the condition number by roughly a factor of $`\mathrm{max}_k\alpha _k/\mathrm{min}_k\alpha _k`$ is thus achieved.

Evidently this ratio depends on the gauge field configuration and on the number of modes that are being subtracted. With $`n`$ as small as $`4`$ we have found factors of $`30`$ and more. An interesting remark in this context is that the improvement factor is roughly independent of the lattice parameters, because the low-lying eigenvalues of $`A`$ scale proportionally to $`(\mathrm{\Sigma }V)^2`$ (cf. subsect. 2.3). The probability distribution of the ratio of the $`n`$th to the first eigenvalue is in fact universal in random matrix theory and analytically computable.

8.4 The low-mode contribution --- a precision issue

The program that minimizes the Ritz functional of $`A`$ initially yields a set of approximate eigenvectors $`\stackrel{~}{e}_1,\mathrm{},\stackrel{~}{e}_n`$ that are orthonormal to machine precision but that are not guaranteed to satisfy eq. (8.2) very accurately. We can, however, easily correct for this deficit by defining the vectors $`e_k`$ to be the exact eigenvectors of $`A`$ in the subspace spanned by the given vectors. The numbers $`\alpha _k`$ are then the eigenvalues of the matrix

$$\begin{array}{c}\text{ }M_{kl}=(\stackrel{~}{e}_k,A\stackrel{~}{e}_l)\text{(8.11)}\text{ }\hfill \end{array}$$

and the residues $`r_k`$ are obtained through $`r_k=(1P)Ae_k`$.

While the subspace spanned by the eigenvectors is by definition given to machine precision, the calculation of the eigenvalues $`\alpha _k`$ is a possible source of numerical inaccuracy. Note that the error in these numbers is amplified in eq. (8.8) and may thus easily affect the precision of the solution vector.

The matrix $`M_{kl}`$ should therefore be calculated using a high-precision approximation to the Neuberger–Dirac operator. It is also important to apply the operator $`A`$ in its original form, eq. (8.1), so that the eigenvalues $`\alpha _k`$ of the matrix are obtained with an error equal to the approximation error on $`D_m`$ times $`(\alpha _k)^{1/2}`$ (rather than just the approximation error).

9. Adapted-precision inversion algorithm

Iterative improvement combined with the use of low-precision arithmetic is mentioned in many books on numerical mathematics as an effective method to speed up the computation of the solution of linear systems (see refs. \[?,?\] for example). In the present context an additional acceleration can be achieved by adjusting the accuracy of the approximation to the Neuberger–Dirac operator in the course of the inversion algorithm. This has first been described in ref. \[?\], and we now wish to present a second scheme that results in an even larger acceleration factor.

We again consider eq. (8.1) in the chirality sector that does not contain the zero-modes. The adapted-precision algorithm that will be discussed should then be applied to the preconditioned system (8.7), but in order to keep the presentation as simple as possible we shall stick to the original system in this section.

9.1 Conjugate gradient algorithm

For the basic inversion algorithm we choose the standard conjugate gradient method. To fix our notations we now first write down the defining equations of this algorithm even though these are very well known \[?,?\].

The algorithm generates a sequence $`\psi _i`$, $`i=1,2,3,\mathrm{}`$, of approximate solutions recursively together with an accompanying sequence $`p_i`$ of so-called search directions. If we introduce the residues

$$\begin{array}{c}\text{ }g_i=\eta A\psi _i,\text{(9.1)}\text{ }\hfill \end{array}$$

the recursion is given by

$$\begin{array}{c}\text{ }\psi _{i+1}=\psi _i+a_ip_i,a_i(g_i,g_i)/(p_i,Ap_i),\text{(9.2)}\text{ }\hfill \\ \\ \text{ }g_{i+1}=g_ia_iAp_i,\text{(9.3)}\text{ }\hfill \\ \\ \text{ }p_{i+1}=g_{i+1}+b_ip_i,b_i(g_{i+1},g_{i+1})/(g_i,g_i),\text{(9.4)}\text{ }\hfill \end{array}$$

and is usually started by setting

$$\begin{array}{c}\text{ }\psi _1=0,g_1=\eta ,p_1=g_1.\text{(9.5)}\text{ }\hfill \end{array}$$

Further information (on orthogonality and convergence properties in particular) can be found in the books quoted above.

9.2 Iterative improvement

Once an approximate solution $`\chi _1`$ of the linear system is found, by performing a certain number of conjugate gradient iterations for example, it can be improved in the following way. First note that the norm of the residue

$$\begin{array}{c}\text{ }\eta _1=\eta A\chi _1\text{(9.6)}\text{ }\hfill \end{array}$$

is much smaller than $`\eta `$ as otherwise we would not consider $`\chi _1`$ to be an approximate solution. Secondly it is obvious that the solution $`\chi `$ of the subtracted system

$$\begin{array}{c}\text{ }A\chi =\eta _1\text{(9.7)}\text{ }\hfill \end{array}$$

yields the solution of the original system through

$$\begin{array}{c}\text{ }\psi =\chi _1+\chi .\text{(9.8)}\text{ }\hfill \end{array}$$

The important point to note is that $`\chi `$ will in general be a small correction to $`\chi _1`$ since $`\eta _1`$ is small. In particular, if we solve the subtracted system to a relative precision of a few decimal places, the accuracy of the total field $`\psi `$ is improved by these many digits.

We can now iterate this procedure and obtain the solution in the form

$$\begin{array}{c}\text{ }\psi =\underset{k=1}{\overset{\mathrm{}}{}}\chi _k,\text{(9.9)}\text{ }\hfill \\ \\ \text{ }\eta _0=\eta ,\eta _k=\eta _{k1}A\chi _k(k=1,2,3,\mathrm{}),\text{(9.10)}\text{ }\hfill \end{array}$$

where the fields $`\chi _k`$ are assumed to be approximate solutions of the subtracted linear systems $`A\chi =\eta _{k1}`$. Depending on how accurate the solutions are, their magnitude $`\chi _k`$ will decrease more or less rapidly, but the algorithm is in any case always exact since the inaccuracies are eventually corrected for.

9.3 Computation of $`\chi _k`$

So far no acceleration has been achieved and we have merely reorganized the calculation in a fairly complicated way. The idea is now to use a numerical representation of $`A`$ with low precision when calculating $`\chi _k`$. In doing so we should of course be very careful to choose a sufficiently accurate approximation as otherwise it is possible that the conjugate gradient algorithm becomes inefficient or even unstable.

Before addressing this question, we recall that $`A`$ acts in the subspace of negative chirality where

$$\begin{array}{c}\text{ }A=\frac{2}{\overline{a}}P_{}D_m^{}P_{},m^{}=\frac{1}{2}\overline{a}m^2.\text{(9.11)}\text{ }\hfill \end{array}$$

So if we use this representation of the operator and if a minmax polynomial is chosen that achieves a certain relative precision $`\delta `$, the total approximation error is

$$\begin{array}{c}\text{ }\stackrel{~}{A}A\tau \frac{2}{a^2}\left\{(1+s)^2\frac{1}{4}(am)^2\right\}\left\{2(\kappa _++\kappa _{})+\delta \right\}.\text{(9.12)}\text{ }\hfill \end{array}$$

Note incidentally that the original form (8.1) of the operator offers no numerical advantage at this point. The representation (9.11) is in fact more efficient by nearly a factor of $`2`$ for a given uniform approximation error which is what counts here.

Now when the conjugate gradient algorithm is applied to solve the linear system $`A\chi =\eta _{k1}`$, a sequence $`\chi _{k,i}`$ of approximate solutions is generated (cf. subsect. 9.1). The decrease of the associated residues

$$\begin{array}{c}\text{ }g_{k,i}=\eta _{k1}A\chi _{k,i}\text{(9.13)}\text{ }\hfill \end{array}$$

can then be observed and the recursion is stopped as soon as $`g_{k,i}`$ falls below a certain level or if the number of iterations reaches a specified maximal value. In the course of this calculation a minimal requirement is that the residues are obtained to some precision at least, i.e. the chosen approximation for the operator $`A`$ should be such that the error $`\tau \chi _{k,i}`$ is always smaller than $`g_{k,i}`$.

Since we cannot know in advance how large the approximate solutions $`\chi _{k,i}`$ are going to be, the accuracy $`\tau `$ should be adjusted dynamically by increasing the degree of the minmax polynomial when needed. A relatively low precision, say $`\delta =10^3`$, is certainly good enough for the first few iterations, but as soon as $`\tau \chi _{k,i}`$ approaches $`g_{k,i}`$, the program should switch to a better approximation. In general this may happen several times before the recursion is halted and the calculated approximate solution $`\chi _k`$ is returned to the calling program.

9.4 Computation of $`\eta _k`$ and global stopping criterion

The residues $`\eta _k`$ can in principle be calculated via eq. (9.10), but it is safer to obtain them directly from the current total solution vector through

$$\begin{array}{c}\text{ }\eta _k=\eta A\underset{l=1}{\overset{k}{}}\chi _l.\text{(9.14)}\text{ }\hfill \end{array}$$

In this way $`\eta _k`$ is guaranteed to be the true residue of the current approximation to the solution $`\psi `$ of the linear system that we wish to solve eventually. In particular, if the global stopping criterion is based on the norm of the so computed residue, we need not worry about a possible accumulation of approximation errors in the course of the calculation. The algorithm simply stops when and only when an algorithm-independent criterion is fulfilled.

It is not difficult to figure out how good the numerical approximation of $`A`$ must be in eq. (9.14). Basically the approximation error $`\tau `$ times the norm of the current total solution vector should be small compared to $`\eta _k`$ so that the residue is correctly obtained to a few decimal places. The stopping criterion for the next subtracted system, $`A\chi =\eta _k`$, then has to be such that the accuracy of the solution is not driven to a point, where the random digits in the source $`\eta _k`$ become important.

9.5 Fine tuning of the algorithm

To achieve the maximal possible acceleration, the degree of the minmax polynomial should never become very large in the course of the solution of the subtracted systems. It is therefore important that the conjugate gradient algorithm stops before this happens, which will be the case if the limits on the residue and the maximal number of iterations are properly chosen.

A small additional improvement is obtained if the search direction is passed from the current subtracted system to the next as if the whole algorithm was a single sequence of conjugate gradient iterations. This is a correct procedure up to rounding and approximation errors that can partly be corrected by adjusting the initial search directions $`p_{k,1}`$ such that the orthogonality relations

$$\begin{array}{c}\text{ }((p_{k,1}g_{k,1}),g_{k,1})=((p_{k,1}g_{k,1}),Ap_{k,1})=0\text{(9.15)}\text{ }\hfill \end{array}$$

are restored to machine precision.

If the stopping criteria are well chosen, this algorithm in general leads to significant savings in computer time (factors of $`3`$ to $`4`$ are not uncommon). As in the case of the calculation of the index of the Dirac operator discussed in sect. 6, the use of 32 bit arithmetic can give another factor of $`2`$ or so on current PC processors \[?\]. Single-precision arithmetic is in fact completely adequate for the approximate solution of the subtracted systems \[but will usually be insufficient for the accumulation of the total solution vector and the calculation of the residues (9.14)\].

10. Concluding remarks

Numerical simulations in the $`ϵ`$–regime of QCD are difficult, not only because they tend to require large amounts of computer time but also because special care needs to be taken to remain in control of the approximation errors. In this paper we have presented an algorithmic framework that is fully satisfactory in this respect and that makes the calculations significantly faster.

Although we focused on a particular formulation of lattice QCD, it is clear that many of the methods that we have discussed are more widely applicable. An obvious case in this connection are the variants of lattice Dirac operators and gauge field actions that have recently been studied in refs. \[?,?–?\]. It is also conceivable that, under certain conditions, even ordinary formulations of lattice QCD can profit from mixed-precision algorithms and low-mode preconditioning.

A disadvantage of our algorithms is perhaps that the propagators for different quark masses cannot be obtained simultaneously. In the $`ϵ`$–regime the quark masses are very small, however, and apart from the zero-mode contribution (which must be calculated only once), the quark propagator should not vary a lot from one value of the quark mass to another. By taking linear combinations of previously computed propagators as a first approximation, the computer time needed for the calculation of the propagator at the next quark mass should thus be significantly reduced.

The computation of the quark propagator in the $`ϵ`$–regime, using just the standard conjugate gradient algorithm, requires a computational effort that grows roughly like the square of the lattice volume at large volumes. Once low-mode preconditioning is switched on, and if the number of subtracted modes is scaled proportionally to some fractional power of the volume, the situation becomes less transparent, but it is possible that a more favourable asymptotic behaviour will result. In any case the total acceleration factor achieved by the methods described in this paper is quite impressive, and we hope that physically interesting computations in the $`ϵ`$–regime will now be feasible on much larger lattices than previously considered.

We are indebted to Pilar Hernández, Mikko Laine, Laurent Lellouch and Peter Weisz for helpful discussions on chiral perturbation theory and on the distribution of the low-lying eigenvalues of the Dirac operator. We also wish to thank DESY for computer time and the staff of the computer centre for their support. L.G. was supported in part by the EU under contract HPRN-CT-2000-00145. C.H. acknowledges support from the EU under contracts HPMF-CT-2001-01468 and HPRN-CT-2000-00145.

Appendix A. Notations

A.1 Index conventions and Dirac matrices

Lorentz indices are taken from the middle of the Greek alphabet and run from $`0`$ to $`3`$. Where appropriate repeated indices are automatically summed over. The Dirac matrices are assumed to satisfy

$$\begin{array}{c}\text{ }(\gamma _\mu )^{}=\gamma _\mu ,\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu },\text{(}\text{A}\text{.1)}\text{ }\hfill \end{array}$$

but are otherwise left unspecified. It is often advantageous, however, to choose a chiral representation where

$$\begin{array}{c}\text{ }\gamma _5=\gamma _0\gamma _1\gamma _2\gamma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).\text{(}\text{A}\text{.2)}\text{ }\hfill \end{array}$$

In particular, pairs of Weyl fields (such as $`e_k`$ and $`r_k`$ in sect. 8) can then conveniently be stored in the memory space allocated for a Dirac field.

A.2 Wilson-Dirac operator

In terms of the gauge-covariant forward and backward difference operators

$$\begin{array}{c}\text{ }_\mu \psi (x)=\frac{1}{a}\left\{U(x,\mu )\psi (x+a\widehat{\mu })\psi (x)\right\},\text{(}\text{A}\text{.3)}\text{ }\hfill \\ \\ \text{ }_\mu \psi (x)=\frac{1}{a}\left\{\psi (x)U(xa\widehat{\mu },\mu )^1\psi (xa\widehat{\mu })\right\},\text{(}\text{A}\text{.4)}\text{ }\hfill \end{array}$$

the Wilson-Dirac operator is given by

$$\begin{array}{c}\text{ }D_\mathrm{w}=\frac{1}{2}\left\{\gamma _\mu (_\mu +_\mu )a_\mu _\mu \right\},\text{(}\text{A}\text{.5)}\text{ }\hfill \end{array}$$

where $`a`$ denotes the lattice spacing, $`U(x,\mu )\mathrm{SU}(3)`$ the link variables and $`\widehat{\mu }`$ the unit vector in direction $`\mu `$.

A.3 Chebyshev polynomials

The Chebyshev polynomials $`T_n(z)`$, $`n=0,1,2,\mathrm{}`$, are defined through

$$\begin{array}{c}\text{ }T_n(\mathrm{cos}\theta )=\mathrm{cos}(n\theta ).\text{(}\text{A}\text{.6)}\text{ }\hfill \end{array}$$

In particular, $`T_0(z)=1`$ and $`T_1(z)=z`$ while for larger degrees the polynomials may be obtained algebraically using

$$\begin{array}{c}\text{ }T_{n+1}(z)2zT_n(z)+T_{n1}(z)=0,n1.\text{(}\text{A}\text{.7)}\text{ }\hfill \end{array}$$

The Clenshaw recursion \[?\] is a direct consequence of these equations.

Appendix B. Proof of lemma 5.1

We first note that

$$\begin{array}{c}\text{ }(1P)QP=PQ(1P)=\varrho .\text{(}\text{B}\text{.1)}\text{ }\hfill \end{array}$$

This is a direct consequence of the definition (5.2) of $`\varrho `$ and of the hermiticity of the operators $`Q`$ and $`P`$. If we define the projected operator

$$\begin{array}{c}\text{ }\stackrel{~}{Q}=PQP+(1P)Q(1P),\text{(}\text{B}\text{.2)}\text{ }\hfill \end{array}$$

it follows that

$$\begin{array}{c}\text{ }Q\stackrel{~}{Q}=(1P)QP+PQ(1P)=\varrho \text{(}\text{B}\text{.3)}\text{ }\hfill \end{array}$$

since $`P`$ and $`1P`$ project to orthogonal subspaces. The min-max principle (ref. \[?\], Theorem XIII.1) now implies that the eigenvalues of $`Q`$ and $`\stackrel{~}{Q}`$, ordered in ascending order and counting multiplicities, are by no more than a distance $`\varrho `$ apart.

Appendix C. Proof of lemma 5.2

We focus on the case of $`P_+`$ since the argument is exactly the same for the other projector. So let us consider an eigenvalue $`\nu _k>0`$ and let $`u_kV`$ be the associated normalized eigenvector of $`PQP`$ \[eq. (5.3)\]. From the definition (5.5) of $`\varrho _k`$ we then infer that

$$\begin{array}{c}\text{ }(Q\nu _k)u_k=w_k,w_k=\varrho _k.\text{(}\text{C}\text{.1)}\text{ }\hfill \end{array}$$

Next we write down the equations

$$\begin{array}{c}\text{ }\left[1(P_+)_{\mathrm{exact}}\right]u_k^2=(u_k,\left[1(P_+)_{\mathrm{exact}}\right]u_k)\text{ }\hfill \\ \\ \text{ }=((Q\lambda )u_k,\frac{1(P_+)_{\mathrm{exact}}}{(Q\lambda )^2}(Q\lambda )u_k)\text{(}\text{C}\text{.2)}\text{ }\hfill \end{array}$$

that are trivially true as long as $`\lambda `$ is not in the spectrum of $`Q`$. We can actually safely set $`\lambda =\nu _k`$ because the distance $`d_k`$ of $`\nu _k`$ from the spectrum in the range of the projector $`1(P_+)_{\mathrm{exact}}`$ does not vanish. The norm of the operator ratio in eq. (C.2) is then at most $`1/d_k^2`$ and this leads to the inequality

$$\begin{array}{c}\text{ }\left[1(P_+)_{\mathrm{exact}}\right]u_k\varrho _k/d_k\text{(}\text{C}\text{.3)}\text{ }\hfill \end{array}$$

when eq. (C.1) is taken into account.

Together with the representation (5.4) of the projector $`P_+`$ in terms of the eigenvectors $`u_k`$, this bound implies

$$\begin{array}{c}\text{ }\left[1(P_+)_{\mathrm{exact}}\right]P_+v\underset{\nu _k>0}{}(\varrho _k/d_k)\left|(u_k,v)\right|\text{(}\text{C}\text{.4)}\text{ }\hfill \end{array}$$

for any fermion field $`v`$, from which we deduce that

$$\begin{array}{c}\text{ }P_+\left[1(P_+)_{\mathrm{exact}}\right]=\left[1(P_+)_{\mathrm{exact}}\right]P_+\kappa _+.\text{(}\text{C}\text{.5)}\text{ }\hfill \end{array}$$

This is not quite what the lemma says, but the inequality is in fact more restrictive than it seems, because $`(P_+)_{\mathrm{exact}}`$ and $`P_+`$ are, by definition, projectors of the same rank. In particular, the trace of the difference $`HP_+(P_+)_{\mathrm{exact}}`$ vanishes, and a few lines of algebra show that

$$\begin{array}{c}\text{ }H^2+H=2P_+\left[1(P_+)_{\mathrm{exact}}\right]P_+\text{ }\hfill \\ \\ \text{ }+P_+\left[1(P_+)_{\mathrm{exact}}\right](1P_+)+(1P_+)\left[1(P_+)_{\mathrm{exact}}\right]P_+.\text{(}\text{C}\text{.6)}\text{ }\hfill \end{array}$$

The operator thus satisfies

$$\begin{array}{c}\text{ }H^{}=H,\mathrm{Tr}H=0,H^2+H\kappa _+\left(1+2\kappa _+\right)\text{(}\text{C}\text{.7)}\text{ }\hfill \end{array}$$

and consequently has to be of order $`\kappa _+`$.

Fig. 2. Plot of $`f(h)=h^2+h`$ and graphical determination of the ranges of $`h`$ (hashed bars) where the bounds (C.8) are satisfied.

To make this explicit, we set $`\kappa =\kappa _+\left(1+2\kappa _+\right)`$ and note that the eigenvalues $`h`$ of $`H`$ have to be such that

$$\begin{array}{c}\text{ }\kappa h^2+h\kappa \text{(}\text{C}\text{.8)}\text{ }\hfill \end{array}$$

(see fig. 2). In terms of $`\kappa `$ the condition stated on the last line of the lemma reads

$$\begin{array}{c}\text{ }2(l+1)\kappa <1\text{(}\text{C}\text{.9)}\text{ }\hfill \end{array}$$

so that $`\kappa <\frac{1}{4}`$ and

$$\begin{array}{c}\text{ }|h|\frac{\kappa }{12\kappa }\text{or}|h+1|\frac{\kappa }{12\kappa }.\text{(}\text{C}\text{.10)}\text{ }\hfill \end{array}$$

Now since $`H`$ acts in a subspace of dimension $`2l`$ at most, and since it has vanishing trace, the second range in (C.10) is excluded if

$$\begin{array}{c}\text{ }(2l1)\frac{\kappa }{12\kappa }<1\frac{\kappa }{12\kappa }\text{(}\text{C}\text{.11)}\text{ }\hfill \end{array}$$

which is the case in view of (C.9). All eigenvalues of $`H`$ thus satisfy the first inequality in (C.10) so that

$$\begin{array}{c}\text{ }H\frac{\kappa }{12\kappa }.\text{(}\text{C}\text{.12)}\text{ }\hfill \end{array}$$

Recalling the definition of $`H`$ and $`\kappa `$, this proves the lemma.

References

 J. Gasser, H. Leutwyler, Phys. Lett. B188 (1987) 477; Nucl. Phys. B307 (1988) 763

 H. Neuberger, Phys. Rev. Lett. 60 (1988) 889; Nucl. Phys. B300 (1988) 180

 P. Hasenfratz, H. Leutwyler, Nucl. Phys. B343 (1990) 241

 F. C. Hansen, Nucl. Phys. B345 (1990) 685

 F. C. Hansen, H. Leutwyler, Nucl. Phys. B350 (1991) 201

 H. Leutwyler, A. Smilga, Phys. Rev. D46 (1992) 5607

 P. H. Ginsparg, K. G. Wilson, Phys. Rev. D25 (1982) 2649

 D. B. Kaplan, Phys. Lett. B288 (1992) 342; Nucl. Phys. B (Proc. Suppl.) 30 (1993) 597

 Y. Shamir, Nucl. Phys. B406 (1993) 90

 V. Furman, Y. Shamir, Nucl. Phys. B439 (1995) 54

 P. Hasenfratz, Nucl. Phys. B (Proc. Suppl.) 63 (1998) 53; Nucl. Phys. B525 (1998) 401

 P. Hasenfratz, V. Laliena, F. Niedermayer, Phys. Lett. B427 (1998) 125

 H. Neuberger, Phys. Lett. B417 (1998) 141; ibid. B427 (1998) 353

 M. Lüscher, Phys. Lett. B428 (1998) 342

 P. Hernández, K. Jansen, M. Lüscher, Nucl. Phys. B552 (1999) 363

 R. G. Edwards, U. M. Heller, R. Narayanan, Phys. Rev. D59 (1999) 094510

 P. Hernández, K. Jansen, L. Lellouch, Phys. Lett. B469 (1999) 198

 T. DeGrand (MILC collab.), Phys. Rev. D63 (2001) 034503

 P. Hasenfratz, S. Hauswirth, T. Jörg, F. Niedermayer, K. Holland, Nucl. Phys. B643 (2002) 280

 E. V. Shuryak, J. J. Verbaarschot, Nucl. Phys. A560 (1993) 306

 J. J. Verbaarschot, I. Zahed, Phys. Rev. Lett. 70 (1993) 3852

 J. J. Verbaarschot, Phys. Rev. Lett. 72 (1994) 2531

 P. H. Damgaard, Phys. Lett. B424 (1998) 322

 S. M. Nishigaki, P. H. Damgaard, T. Wettig, Phys. Rev. D58 (1998) 087704

 J. C. Osborn, D. Toublan, J. J. Verbaarschot, Nucl. Phys. B540 (1999) 317

 P. H. Damgaard, J. C. Osborn, D. Toublan, J. J. Verbaarschot, Nucl. Phys. B547 (1999) 305

 P. H. Damgaard, S. M. Nishigaki, Phys. Rev. D63 (2001) 045012

 G. H. Golub, C. F. van Loan, Matrix computations, 2nd ed. (Johns Hopkins University Press, Baltimore, 1989)

 B. Bunk, Nucl. Phys. B (Proc. Suppl.) 63 (1998) 952

 H. Neuberger, Phys. Rev. Lett. 81 (1998) 4060; Int. J. Mod. Phys. C 10 (1999) 1051

 R. G. Edwards, U. M. Heller, R. Narayanan, Nucl. Phys. B540 (1999) 457

 J. van den Eshof, A. Frommer, T. Lippert, K. Schilling, H. A. van der Vorst, Comput. Phys. Commun. 146 (2002) 203

 P. H. Damgaard, P. Hernández, K. Jansen, M. Laine, L. Lellouch, Finite-Size Scaling of Vector and Axial Current Correlators, hep-lat/0211020

 S. Weinberg, Phys. Rev. D8 (1973) 605 and 4482

 M. K. Gaillard, B. W. Lee, Phys. Rev. Lett. 33 (1974) 108

 G. Altarelli, L. Maiani, Phys. Lett. B52 (1974) 351

 L. Maiani, G. Martinelli, G. C. Rossi, M. Testa, Nucl. Phys. B289 (1987) 505

 S. Capitani, L. Giusti, Phys. Rev. D64 (2001) 014506

 C. W. Bernard, T. Draper, A. Soni, H. D. Politzer, M. B. Wise, Phys. Rev. D32 (1985) 2343

 P. Hernández, M. Laine, in preparation

 P. H. Damgaard, Nucl. Phys. B (Proc. Suppl.) 106 (2002) 29

 R. G. Edwards, U. M. Heller, J. E. Kiskis, R. Narayanan, Phys. Rev. Lett. 82 (1999) 4188; Phys. Rev. D61 (2000) 074504

 P. H. Damgaard, R. G. Edwards, U. M. Heller, R. Narayanan, Phys. Rev. D61 (2000) 094503

 D. C. Handscomb (ed.), Methods of numerical approximation (Pergamon Press, Oxford, 1966)

 W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in FORTRAN, 2nd ed. (Cambridge University Press, Cambridge, 1992)

 B. Bunk, K. Jansen, M. Lüscher, H. Simma, Conjugate gradient algorithm to compute the low-lying eigenvalues of the Dirac operator in lattice QCD, notes (September 1994)

 T. Kalkreuter, H. Simma, Comput. Phys. Commun. 93 (1996) 33

 M. Lüscher, Nucl. Phys. B (Proc. Suppl.) 106 (2002) 21

 M. Reed, B. Simon, Methods of modern mathematical physics, vol. IV (Academic Press, New York, 1978)

 P. H. Damgaard, M. C. Diamantini, P. Hernández, K. Jansen, Nucl. Phys. B629 (2002) 445

 L. Giusti, C. Hoelbling, C. Rebbi, Phys. Rev. D64 (2001) 114508 \[E: ibid. D65 (2002) 079903\]; Nucl. Phys. B (Proc. Suppl.) 106 (2002) 739

 P. Hasenfratz, S. Hauswirth, K. Holland, T. Jörg, F. Niedermayer, U. Wenger, Nucl. Phys. B (Proc. Suppl.) 94 (2001) 627; Int. J. Mod. Phys. C12 (2001) 691

 C. Gattringer, Phys. Rev. D63 (2001) 114501

 C. Gattringer, I. Hip, C. B. Lang, Nucl. Phys. B597 (2001) 451

 W. Bietenholz, Nucl. Phys. B644 (2002) 223

 T. DeGrand, A. Hasenfratz, T. G. Kovacs, hep-lat/0211006


text:
# Parton Percolation in Nuclear Collisions

## 1 Initial State Conditions

Statistical QCD predicts that with increasing temperature, strongly interacting matter will undergo a transition from a hadronic phase to a plasma of deconfined quarks and gluons. In the hadronic state, the chiral symmetry of the QCD Lagrangian (for massless quarks) is spontaneously broken; in the quark-gluon plasma, it is restored. These predictions are the result of finite temperature lattice QCD studies, and for the calculations it is crucial that they deal with a thermal medium, i.e., with equilibrium thermodynamics.

The initial state of two colliding nuclei is clearly a non-equilibrium configuration. The canonical view of its evolution is schematically illustrated in Fig. 1. After the collision, there is a short pre-equilibrium stage, in which the primary partons of the colliding nuclei interact, multiply and then thermalize to form a quark-gluon plasma. This QGP then expands, cools and hadronizes. A prerequisite for the equilibration process is evidently that the partons originating from different nucleons form a large-scale interconnected system. If there is no “cross talk” between partons from different nucleons, thermalization is not possible. The problem of color connection has been studied in hadron production through $`W^+W^{}`$ decay at LEP. The $`W`$’s are produced essentially at rest in the annihilation of an energetic $`e^+e^{}`$ pair (see Fig. 2), and it is possible to compare the reaction in which both $`W`$’s produce hadronic jets to that in which one decays leptonically. If there is cross talk between the decay quarks of one $`W`$ with those from the other, the multiplicity of the four-jet decay is predicted to be less than twice that in the two-jet decay . The data show no such reduction, suggesting that the decay quarks from different $`W`$’s don’t communicate .

It is therefore necessary to determine under what conditions the initial state parton configurations can lead to color connection, and more specifically, if variations of the initial state can lead to a transition from disconnected to connected partonic clusters. The results of such a study of the pre-equilibrium state in nuclear collisions do not depend on the subsequent evolution and thus in particular not require any kind of thermalization.

The structural problem underlying the transition from disconnected to connected systems of many components is a very general one, ranging from clustering in spin systems to the formation of galaxies. The formalism is given by percolation theory, which describes geometric critical behavior . We shall return to the basic idea a little later on.

## 2 Partons in Nuclei

If you look at a fast nucleon coming at you, what do you see? The answer depends on who’s looking. Another nucleon or a pion sees a disc of radius $`r`$ 1 fm and a certain greyness. A hard photon, with a resolution scale $`Q^11`$ fm, sees a swarm of partons. How many there are depends on the resolution scale: given a finer scale, you can see smaller partons, and there are more the harder you look (Fig. 3). The partons in a nucleon have a transverse size $`r_T`$ determined by their intrinsic transverse momentum $`k_T`$, with $`r_T1/k_T`$. The scale $`Q^1`$ specifies the minimum $`k_T^1`$ resolved, so the probing photon sees all partons in the range $`0k_TQ`$.

The momentum $`p`$ of the incoming nucleon is distributed among the partons; a parton of momentum $`k`$ carries the fraction $`x=k/p`$. In deep inelastic scattering experiments, the distribution of partons in a nucleon is determined for given $`x`$ and $`Q`$. Denoting the gluon content by $`g(x,Q)`$, that of quarks and antiquarks by $`q(x,Q)`$ and $`\overline{q}(x,Q)`$, respectively, we have the overall momentum conservation sum rule

$$_0^1𝑑xx\{g(x,Q)+\underset{i}{}[q_i(x,Q)+\overline{q}_i(x,Q)]\}=1,$$
(1)
where $`i`$ counts the number of quark flavors.

The number of partons in a nucleon at rapidity $`y`$, as seen by a photon of scale $`Q`$, is thus given by

$$\frac{dN}{dy}=x\{g(x,Q)+\underset{i}{}[q_i(x,Q)+\overline{q}_i(x,Q)]\},$$
(2)
where $`x=(k_0+k_L)/(p_0+p_L)`$ in terms of the parton and nucleon energies $`k_0,p_0`$ and longitudinal momenta $`k_L,p_L`$. Since the scale $`Q`$ specifies the maximum $`k_T`$ resolved, $`dN/dy`$ in Eq. (2) gives us the total number of partons in the range $`0k_TQ`$.

In a minimum bias nucleon-nucleon collision, the transverse parton size itself determines the resolution: it sets the scale at which partons ‘probe each other’ in the colliding nucleons, so that here the highest relevant $`k_T`$ fixes $`Q`$. Since at $`y=0`$, the fractional momentum is $`x=k_T/\sqrt{s}`$, Eq. (2) provides at given $`\sqrt{s}`$ the total number of partons of transverse momenta up to $`Q`$.

As mentioned, the quark and gluon distributions in a nucleon are determined from deep inelastic scattering data. In Fig. 4, we show the resulting variation of $`(dN/dy)_{y=0}`$ as function of $`Q^2<k_T^2>`$ for two values of the c.m.s. energy, $`\sqrt{s}=20`$ GeV (SPS) and 200 GeV (RHIC), using the GRV94DIS parametrization , which is particularly suitable for our kinematic range. It is evident that with increasing $`Q^2`$, more and more partons of decreasing size come into play, so that $`(dN/dy)_{y=0}`$ increases strongly with $`Q`$. It is also clear that at higher collision energy, there are more partons – at RHIC two times more than at SPS. The slow decrease at large $`Q^2`$ is due to kinematic constraints: at fixed $`\sqrt{s}`$, increasing $`Q^2`$ means increasing $`x`$ and hence decreasing gluon or sea quark density.

## 3 Partons in Nuclear Collisions

Consider now the collision of two heavy nuclei at high energy, as seen in the overall center of mass. The Lorentz-contraction in the longitudinal direction makes it a collision of two thin discs, so that in the transverse plane, the parton density increases with $`A`$. The partons from different nucleon begin to overlap and form clusters: see Fig. 5. How does the cluster size grow with parton density, and when does it reach the dimension of the total transverse collision area? These are precisely the questions addressed by percolation theory, so that here we make a small interlude.

### 3.1 Percolation Theory

Consider placing $`N`$ small circular discs (‘partons’) of radius $`r`$ onto a large circular manifold (‘the transverse nuclear plane’) of radius $`Rr`$; the small discs may overlap. With increasing parton density $`nN/\pi R^2`$, this overlap will lead to more and larger connected partonic clusters. The striking feature of this phenomenon is that the average cluster size $`S_{\mathrm{cl}}`$ does not grow as some power of $`n`$; instead, it increases very suddenly from very small to very large values (see Fig. 6). This suggests some kind of geometric critical behavior, and in fact in the ‘thermodynamic limit’ $`N\mathrm{},R\mathrm{}`$, the cluster size diverges at a critical threshold value $`n_c`$ of the density $`n`$,

$$S_{\mathrm{cl}}(n_cn)^\gamma ,$$
(3)
as $`nn_c`$ from below. This appearence of infinite clusters at $`n=n_c`$ is defined as percolation: the size of the cluster reaches the size of the system. The divergence is governed by the critical exponent $`\gamma =43/18`$, determined analytically, while the threshold

$$n_c=\frac{1.128}{\pi r^2}$$
(4)
is obtained numerically or through analytical approximation . Hence we obtain in the limit of large systems

$$\frac{N}{\pi R^2}=\frac{1.128}{\pi r^2}$$
(5)
as percolation condition. Note that because of parton overlap, the manifold is at percolation not totally covered by discs, even though the overall disc area slightly exceeds that of the manifold: $`N\pi r^2=1.128\pi R^2`$. In fact, one can show that at $`n=n_c`$, the fraction

$$1\mathrm{exp}\{1.128\}0.68$$
(6)
of the area $`\pi R^2`$ is covered by partonic discs.

Before we return to the study of nuclear collisions, we want to comment briefly on the relation between percolation and thermal phase transitions . Some thermal critical behavior, such as the magnetization transition for ferromagnetic spin systems, can be equivalently formulated as percolation. However, percolation seems to be a more general phenomenon and in particular can occur even when the partition function is analytic, i.e., when there is no thermal critical behavior. A specific example of this is the Ising model in a non-vanishing external field, which has a percolation transition even though there is no magnetization transition.

### 3.2 Parton Percolation

The results of the previous subsection tell us that in nuclear collisions there is indeed, as function of parton density, a sudden onset of large-scale color connection. There is a critical density at which the partons form one large cluster, losing their independent existence and their relation to the parent nucleons. Parton percolation  is thus the onset of color deconfinement and although it is a necessary prerequisite for any subsequent QGP formation, it does not require or imply any kind of parton thermalization.

To obtain quantitative predictions, we have to specify the relevant scales. The partonic size is through the uncertainty relation determined by its average transverse momentum,

$$\pi r^2\frac{\pi }{<k_T^2>},$$
(7)
and as mentioned, for a given resolution scale, $`<k_T^2>Q^2`$. Given the parton density in a nucleon, we now have to specify the density in a nucleus-nucleus collision. At SPS energy ($`\sqrt{s}20`$ GeV), the wounded nucleon model appears to work quite well, so that we have in a central $`AA`$ collision

$$\left(\frac{dN}{dy}\right)_{y=0}^{AA}2A\left(\frac{dN}{dy}\right)_{y=0}.$$
(8)
It is clear, however, that at higher energies, collision dependent contributions will come into play .

For central $`AA`$ collisions, we thus obtain the percolation condition

$$\frac{2A}{\pi A^{2/3}}\left(\frac{dN}{dy}\right)_{y=0}=\frac{1.128}{\pi Q^2}$$
(9)
in terms of $`A`$, the resolution scale $`Q`$ and the nucleonic parton density obtained in deep inelastic scattering. Let us separate the basic parton contributions from the nuclear dependence and rewrite eq. (9) as

$$\frac{1}{Q^2}\left(\frac{dN}{dy}\right)_{y=0}=\frac{1.128}{2A^{1/3}}.$$
(10)
When the l.h.s. of this equation, determined by P.D.F.’s, becomes equal to the r.h.s., fixed by nuclear size, we have the onset of percolation. In Fig. 7, the results are shown for typical SPS and RHIC energies. We thus find that for $`\sqrt{s}=20`$ GeV, there is percolation for $`A>60`$, while for ($`\sqrt{s}=200`$ GeV) it sets in somewhat earlier, for $`A>40`$.

In the case of non-central $`AA`$ collisions, the manifold disc size $`\pi R^2`$ has to be replaced by the actual transverse overlap area at the given impact parameter. This overlap area can be determined in a Glauber study, using Woods-Saxon nuclear profiles , and the resulting counterpart of eq. (10) then leads to Figs. 8. Here the number $`N_{\mathrm{part}}`$ of wounded or participant nucleons is used to specify the centrality of the collision, since this quantity is directly measurable. For $`PbPb`$ collisions at $`\sqrt{s}=20`$ GeV, this leads to an onset of percolation at $`N_{\mathrm{part}}150`$ (corresponding to an impact parameter $`b6`$ fm), while for $`AuAu`$ at $`\sqrt{s}=200`$ GeV, $`N_{\mathrm{part}}80`$ (with $`b10`$ fm) is the threshold.

Beyond the percolation point, we then have a parton condensate, containing interacting and hence color-connected partons of all scales $`k_TQ`$. The percolation point thus specifies the onset of color deconfinement; it says nothing about any subsequent thermalization. If there is eventual thermalization, the partonic momentum $`k_T`$ will be related to the temperature $`T`$; hence the resolution scale $`Q`$, which determines the range of $`k_T`$, is in some sense a precursor of $`T`$. It is thus of interest to check how the percolation value $`Q_s`$ is related to the effective nuclear size and to $`\sqrt{s}`$. This is illustrated in Fig. 9 for the $`A`$-dependence at two energies studied here. In the same figure, we also show the corresponding results as function of centrality at $`A=200`$. It is seen that bigger $`\sqrt{s}`$, larger $`A`$ or more central collisions lead to a ‘hotter’ parton condensate, in the mentioned pre-thermal sense.

## 4 Observable Consequences

We have seen that in nuclear collisions, the parton structure of nucleons leads to critical behavior in the form of parton percolation in the transverse collision plane. This critical behavior is independent of any subsequent thermalization; it is determined by the initial collision conditions in the pre-equilibrium stage.

The parton condensate which is formed through percolation is closely related to the color glass condensate  studied for nuclear collisions in the limit of large $`A`$ and/or $`\sqrt{s}`$; the different approaches simply focus on different aspects. In percolation studies, the central topic is the onset of parton condensation in terms of geometric critical behavior. In contrast, the color glass condensate describes the features of the high density limit for the pre-equilibrium partonic medium, in particular in terms of classical fields.

The occurrence of geometric critical behavior in the pre-equilibrium stage of nuclear collisions can lead to observable consequences even if there never is any subsequent thermalization. We note here three possible effects of this kind: charmonium suppression, strangeness enhancement, and energy-independent hadronic source radii.

### 4.1 Charmonium Suppression

Charmonium states are formed very early in nuclear collisions, with a typical $`J/\psi `$ formation time of some 0.2 - 0.3,fm obtained from binding energy or radius. This is also the time needed for the formation of the parton condensate, as determined by $`Q_s^1`$. The $`J/\psi `$ thus finds itself in the non-equilibrium medium provided by the parton condensate. The typical scale of the charmonium state thus has to be compared to the intrinsic scale $`Q_s`$ of the parton condensate. If the latter is indeed the precursor of temperature, it is also a precursor form of the screening mass. We will therefore assume a charmonium state $`i`$ of radius $`r_i`$ to be dissociated when $`Q_s>r_i`$. With

$$r_{J/\psi }(0.9\mathrm{GeV})^1,r_\chi (0.6\mathrm{GeV})^1r_\psi ^{}(0.45\mathrm{GeV}^1,$$
(11)
for the different subthreshold charmonium states, we then obtain the suppression values of $`Q_s`$ shown in Fig. 10.

To see what these different suppression points imply for $`J/\psi `$ production in nuclear collisions, we recall that in nucleon-nucleon collisions, $`J/\psi `$’s are produced in part through feed-down. Only about 60% of the observed $`J/\psi `$’s are directly produced $`1S`$ states; of the remainder, about 30% come from $`\chi _c`$ and about 10% from $`\psi ^{}`$ decay. Since these decays occur very late in the collision evolution, the parton condensate ‘sees’ and hence suppresses the different charmonium states in the given fractions. From Fig. 10 we thus expect that the survival probability of $`J/\psi `$ s in $`PbPb`$ collisions at the SPS will show a first anomalous suppression step at about $`N_{\mathrm{part}}150`$, since at this point the production fraction from $`\chi _c`$ and $`\psi ^{}`$ decays is removed. A second drop would be expected around $`N_{\mathrm{part}}250`$, since now the directly produced $`J/\psi `$’s are dissociated.

The crucial consequence is the predicted two-step suppression pattern; to obtain reliable numerical values clearly requires the inclusion of more details. In particular, the nuclear geometry, resulting in percolating and non-percolating (‘hot’ vs. ‘cold’) regions in the transverse plane, has to be taken into account correctly in its effect on the suppression. A two-step suppression pattern was first obtained for $`J/\psi `$ production in a quark-gluon plasma, since the different charmonium radii lead to dissociation for different screening radii. It now appears that in fact already pre-equilibrium parton condensation leads to a similar result, so that such a pattern, if observed, implies deconfinement but not the formation of a thermalized quark-gluon plasma.

In this connection we note also that the onset points for charmonium dissociation through parton condensation agree fairly well with the ‘steps’ seen in the measured $`J/\psi `$ survival probability (see Fig. 11). In view of the mentioned theoretical uncertainties and also because of possible kinematic suppression effects for very central collisions , this agreement should so far be considered in a more qualitative way. However, the observed threshold value for the onset of anomalous suppression measured in terms of the energy density is by more than a factor two above the energy density at deconfinement. This seems to support pre-equilibrium deconfinement through parton condensation as underlying mechanism.

### 4.2 Strangeness Enhancement

The relative abundance of the hadrons produced in high energy collisions, from nucleon-nucleon to heavy ion interactions, is well accounted for by an ideal gas of hadronic resonances  The temperature of this gas converges at high energies to $`T_H150180`$ MeV, i.e., to the confinement temperature obtained in lattice QCD for low baryon density (see Fig. 12 for a compilation from nucleus-nucleus collisions ). The baryochemical potential depends on the baryon number content of the initial state, decreasing from values around $`\mu _B0.5`$ GeV in heavy ion collisions at the AGS to near zero for $`pp/p\overline{p}`$ and RHIC heavy ion data.

The only ‘flaw’ in this picture is that in the elementary $`pp/p\overline{p}`$ collisions one observes reduced strangeness production: the relative abundance of strange hadrons is reduced by a factor dependent only on its content of strange quarks/antiquarks. This strangeness suppression disappears for high energy heavy ion collisions.

The dependence of strangeness abundance on the density of interacting hadrons can be accounted for if strangeness conservation is assumed to hold locally . In a thermal medium of temperature $`T`$ and overall volume $`V`$, the relative phase space weight (Boltzmann factor) for the presence of a strange particle of mass $`m`$ is given by $`\mathrm{exp}\{m/T\}`$. If, however, the medium contains only one such strange particle, and if strangeness conservation is taken to hold locally within some correlation volume $`V_0V`$, then the correct Boltzmann factor should be $`\mathrm{exp}\{2m/T\}`$, since the simultaneous presence of the strange particle and its antiparticle requires the expenditure of energy $`2m`$. The crucial point here is the assumed local nature of strangeness conservation: the factor $`\mathrm{exp}\{m/T\}`$ would be correct if the single strange particle could be compensated by a ‘far away’ antiparticle somewhere else in the overall volume $`V`$. It is only the requirement of zero strangeness within $`V_0`$ that leads to the enhanced suppression. Since a given strange particle and its antiparticle in a really ideal gas will eventually separate beyond $`V_0`$, the requirement of local strangeness conservation appearently implies the introduction of a dynamical (i.e., non-ideal) correlation.

To illustrate the effect of this phenomenon, we consider the abundance of kaons in a hadron gas. Given the ideal gas density of kaons, $`n_K(T)=m_K^2K_2(m_K/T)`$, local strangeness conservation leads to the suppressed form

$$n_K(T,x)=n_K(T)\left\{\frac{I_1(x)}{I_0(x)}\right\},$$
(12)
where $`x=V_0n_K(T)`$ specifies the number of kaons inside the correlation volume $`V_0`$; $`I_0,I_1`$ and $`K_2`$ are the corresponding Bessel and Hankel functions of imaginary argument. In the limit of high density or large correlation volume, $`x\mathrm{}`$, $`I_1(x)/I_0(x)1`$, and the ideal gas abundances are correct. For low density or small correlation volume, i.e., for $`x0`$, $`I_1(x)/I_0(x)1`$, thus resulting in an effective strangeness suppression as compared to the ideal gas abundance of kaons (see Fig. 13).

We are thus confronted with the puzzling question of how the transition from suppressed to normal strangeness production occurs. What specifies the correlation volume $`V_0`$? A rather natural solution to this puzzle is that the abundance of strange hadrons is already determined by the initial state partonic cluster size in the transverse plane: the extension of this cluster specifies $`V_0`$. This means that parton percolation and the corresponding sudden increase of $`V_0`$ result in $`I_1/I_01`$ and thus trigger the transition to genuine ideal gas abundances of strange particles. More detailed work on this is in progress; an interesting aspect is obviously the resulting correspondence between charmonium suppression and strangeness enhancement.

### 4.3 Thermal or Statistical?

This point is undoubtedly the most speculative of the present considerations. At various stages of the high energy heavy ion program, it has been asked whether a single collision event already produces a thermal medium (‘matter’), or whether each individual event is something like one member of a Gibbs ensemble in thermodynamics, so that only an average over many events results in a thermal pattern.

If we take the extreme point of view that the observable phenomena are determined fully determined by the initial state parton configuration of the colliding nuclei, without any subsequent thermalization, then a single collision does not lead to a thermodynamic medium. Everything is specified by the given nuclear collision configuration, and there will in particular not be any kind of ‘expanding and cooling matter’. This suggests that the source size as determined by HBT interferometry should essentially measure the initial nuclear collision configuration; it should show no dependence on the collision energy and not lead to any ‘medium life-time’, with always $`R_{\mathrm{out}}R_{\mathrm{side}}`$. Moreover, the transverse expansion velocity should also become constant with increasing collison energy. These features are indeed observed, as shown in Figs. 12 and 14, contrary to all earlier predictions.

Perhaps there exist models which nevertheless can reproduce the entire set of present observations in thermal terms. However, it does seem worthwhile to pursue further the possibility that initial state conditions alone already determine them rather naturally, and to search for some crucial additional observable which would unambiguously indicate event-by-event thermalization.

Acknowledgements

It is a pleasure to thank S. Digal, S. Fortunato, F. Karsch, C. Lourenço, M. Nardi, P. Petreczky and K. Redlich for many stimulating discussions and helpful comments.