# Feasibility study of using the overlap-Dirac operator for hadron spectroscopy. Presented by C. McNeile.

## 1 INTRODUCTION

Two critical systematic errors in the calculation of the $`f_B`$ decay constant from lattice QCD are the chiral extrapolations and the unquenching errors . The only way to reduce these errors is to simulate QCD with lighter quark masses. Unfortunately, because of exceptional configurations, it is difficult to further reduce the masses of the light quarks in quenched simulations with the clover operator. Progress in reducing the sea quark masses in dynamical fermion simulations with Wilson like quarks is slow .

It seems plausible that the difficulty of simulating with light quark masses with the clover operator is due to explicit chiral symmetry breaking in the action. Neuberger has derived  a fermion operator, called the overlap-Dirac operator, that has a lattice chiral symmetry .

Our goal is to simulate the overlap-Dirac operator in a mass region: ($`M_{PS}/M_V=0.30.5`$), which is inaccessible to clover quarks (but not staggered quarks ). Most of the techniques developed in the quenched theory can be used for full QCD simulations .

## 2 THE OVERLAP-DIRAC OPERATOR

The massive overlap-Dirac operator  is

$$D^N=\frac{1}{2}(1+\mu +(1\mu )\gamma _5\frac{H(m)}{\sqrt{H(m)H(m)}})$$
(1)   
where $`H(m)`$ is the hermitian Wilson fermion operator with negative mass, defined by

$$H(m)=\gamma _5(D^Wm)$$
(2)   
where $`D^W`$ is the standard Wilson fermion operator. The parameter $`\mu `$ is related to the physical quark mass and lies in the range $`0`$ to $`1`$. The $`m`$ parameter is a regulating mass, in the range between a critical value and 2. The physics should be independent of the mass $`m`$, but the value of $`m`$ effects the locality of the operator and the number of iterations required in some of the algorithms used to compute the overlap-Dirac operator.

## 3 NUMERICAL TECHNIQUES

Quark propagators are calculated using a sparse matrix inversion algorithm. The inner step of the inverter is the application of the fermion matrix to a vector. For computations that use the overlap-Dirac, the step function

$$ϵ(H)\underset{¯}{b}=\frac{H}{\sqrt{HH}}\underset{¯}{b}$$
(3)   
must be computed using some sparse matrix algorithm. The nested nature of the algorithm required to calculate the quark propagators for the overlap-Dirac operator makes the simulations considerably more expensive than those that use traditional fermion operators.

Practical calculations of the overlap operator are necessarily approximate. To judge the accuracy of our approximate calculation we used the (GW) Ginsparg-Wilson error: 

$$\gamma _5D^N\underset{¯}{x}+D^N\gamma _5\underset{¯}{x}2D^N\gamma _5D^N\underset{¯}{x}\frac{1}{\underset{¯}{x}}$$
(4)   
which just checks that the matrix obeys the Ginsparg-Wilson relation .

Our numerical simulations were done using $`\beta =5.9`$ quenched gauge configurations, with a volume of $`16^332`$. The quark propagators were generated from point sources. For all the algorithms we investigated, we used $`m`$ equal to $`1.5`$.

## 4 LANCZOS BASED METHOD

Borici  has developed a method to calculate the action of the overlap-Dirac operator on a vector, using the Lanczos algorithm. In exact arithmetic, the Lanczos algorithm generates an orthonormal set of vectors that tridiagonalises the matrix. 

$$HQ_n=Q_nT_n$$
(5)   
where $`T_n`$ is a tridiagonal matrix. The columns of $`Q_n`$ contain the Lanczos vectors.

The “trick”, to evaluate the step function (Eq. 3), is to set the target vector $`\underset{¯}{b}`$, as the first vector in the Lanczos sequence. An arbitrary function $`f`$ of the matrix $`H`$ acting on a vector is constructed using

$`(f(H)b)_i`$ $`=`$ $`{\displaystyle \underset{j}{}}(Q_nf(T_n)Q_n^{})_{ij}b_j`$ (6)   
$`=`$ $`b(Q_nf(T_n))_{i\mathrm{\hspace{0.33em}1}}`$ (7)   
where the orthogonality of the Lanczos vectors has been used. The $`f(T_n)`$ matrix is computed using standard dense linear algebra routines. For the step function the eigenvalues of $`T_n`$ are replaced by their moduli. Eq. 7 is linear in the Lanczos vectors and thus can be computed in two passes.

The major problem with the Lanczos procedure is the loss of the orthogonality of the sequence of vectors due to rounding errors. It is not clear how this lack of orthogonality effects the final results. Some theoretical analysis has been done on this method . It is claimed that the lack of orthogonality is not important for some classes of functions.

In Fig. 1, we plot the eigenvalues of the Ginsparg-Wilson operator, as a function of the number of Lanczos steps, for a $`2^4`$ hot $`SU(3)`$ configuration. As the number of Lanczos steps increases, the eigenvalue spectrum moves closer to a circle (the correct result). Even after 50 iterations of the Lanczos algorithm, there are still small deviations from the circle.

Unfortunately, it is much harder to look at the eigenvalue spectrum using a production gauge configuration, so we computed the GW error instead. The GW error was: $`\mathrm{5\hspace{0.33em}10}^3`$ (50 iterations), $`\mathrm{6\hspace{0.33em}10}^4`$ (100) iterations, and $`\mathrm{3\hspace{0.33em}10}^4`$ (300 iterations) on a single gauge configuration.

Fig. 2, is an effective mass plot of the pion, for two choices of mass and number of Lanczos sweeps. The “plateau” in the pion effective mass plot for approximate operator that used 50 Lanczos iterations is higher than the lowest pion mass that can be reached with non-perturbatively improved clover. It is not clear what causes the ”shoulder” in Fig. 2. We would like to compute the eigenvalues of the overlap-Dirac operator on the bigger gauge configuration, to check how accurately we are computing the overlap-Dirac operator.

## 5 RATIONAL APPROXIMATION

The step function can be be approximated by a rational approximation.

$$ϵ(H)H(c_0+\underset{k=1}{\overset{N}{}}\frac{c_k}{H^2+d_k})$$
(8)   
The rational approximation typically approximates the step function, between two values. The eigenvalues of the matrix $`H`$ should lie in the region where the approximation is good. The coefficients $`c_k`$ and $`d_k`$ can be obtained from the Remez algorithm . The number of iterations required in the inverter is controlled by the smallest $`d_k`$ coefficient, which acts like a mass. We have not yet implemented the technique of projecting out some of the low lying eigenmodes .

On one configuration we obtained GW errors of: $`\mathrm{1\hspace{0.33em}10}^4`$, and $`\mathrm{4\hspace{0.33em}10}^5`$, for the $`N=6`$, and $`N=8`$, optimal rational approximations . The multiplicative scaling factor was tuned to obtain the best results. Unfortunately, the above results required up to 600 iterations for the smallest $`d_k`$, which was too large to use as the inner step of a quark propagator inverter.

One feature of the optimal rational approximation , is that the lowest $`d_k`$ factor is smaller than the square of the validity of the approximation, which means that the condition number of the inversion is that of the matrix $`H^2`$. We experimented with a hybrid quadrature and series approximation to Robertson’s integral representation of the step function.

$`ϵ(H)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{2H}{\pi (t^2+H^2)}}𝑑t`$   
$``$ $`{\displaystyle _0^{\theta _S}}{\displaystyle \frac{2H}{\pi (t^2+H^2)}}𝑑t+{\displaystyle _{\theta _S}^{\theta _L}}{\displaystyle \frac{2H}{\pi (t^2+H^2)}}𝑑t`$   
In Eq. 5, the first integral was approximated using an open quadrature rule and the second integral was approximated by a Chebyshev series. The step size in the quadrature formulae reduces the condition number of the required inversion. In our preliminary tests of the algorithm, the hybrid method produced a substantial reduction in the number of iterations required over the optimal rational approximation. However the computed solution was less accurate than that produced by the optimal rational approximation, because the Ginsparg-Wilson error was $`\mathrm{6\hspace{0.33em}10}^4`$.

Clearly more work is required on the algorithms that calculate the step function, before the overlap-Dirac operator can be used in the quark mass region we are interested in.

This work is supported by PPARC. The computations were carried out on the T3D and T3E at EPCC in Edinburgh.