format: text

text:
# Contents

## Abstract and Summary

Quantum Chromodynamics (QCD), the theory of the strong interaction, is one of the most prominent examples for a beautiful and successful physical theory. At large distance, or equivalently at low energy, perturbative expansions in the coupling constant—the standard tool to treat quantum field theories analytically—break down, and a non-perturbative formulation is required to calculate physical quantities. In this thesis, we construct the Fixed-Point fermion action for lattice QCD, which is a highly improved discretization of the continuum theory that preserves the chiral symmetry inherent in the original formulation. We perform studies in quenched light hadron spectroscopy to examine the properties of this action and investigate in detail the chiral limit of pseudoscalar mesons, which is inaccessible to non-chiral lattice formulations.

To start with, Chapter 1 provides a brief introduction to the field of elementary particle physics, to Quantum Chromodynamics and the lattice as a tool to probe the non-perturbative regime of the strong interaction, and motivates the construction of improved transcriptions of the theory to discrete space-time. A long standing problem, namely the formulation of chiral symmetric lattice fermions, is addressed in Chapter 2. An elegant solution has been found using Renormalization Group methods, leading to the classically perfect Fixed-Point actions. Chapter 3 describes the parametrization and construction of the Fixed-Point fermion action for lattice QCD and presents some elementary properties of the resulting Dirac operator. A different possibility to obtain chiral lattice fermions is the overlap construction. We combine the Fixed-Point and the overlap approach in Chapter 4 to remove the residual chiral symmetry breaking of our parametrized Dirac operator, getting a fermion action which inherits the advantages of both formulations at a higher computational cost. The chirality and locality properties of this overlap-improved Dirac operator are then tested in the artificial framework of smooth instanton gauge configurations.

Next, we turn to one of the most fundamental applications of lattice QCD, namely the calculation of hadron masses. Chapter 5 gives an introduction to the technical details of how the light hadron mass spectrum is extracted from lattice simulations. With chiral symmetric fermion actions, it is possible to perform lattice simulations at quark masses very close to or even at the physical mass of up and down quarks, thus allowing to study the chiral limit, which is complicated by non-analytic terms in the quenched approximation to QCD. At such small quark masses, additional quenching effects appear in a finite lattice volume which contaminate in particular the pseudoscalar meson channel and are related to the zero modes of the Dirac operator. We devote Chapter 6 to the study of these topological finite-volume effects and examine possible solutions for the problem of extracting reliable pseudoscalar meson masses at small volumes and quark masses.

In Chapter 7, we present the results of a spectroscopy simulation with the Fixed-Point fermion and gluon lattice actions. This study is the one of the first hadron spectroscopy calculations with a chiral symmetric action including checks for cut-off and finite-volume effects. After estimating the magnitude of the topological quenching effects, we closely examine the chiral limit of the pseudoscalar meson and extract the coefficient of the quenched chiral logarithm in two different ways. We also consider the chiral extrapolations for vector mesons and baryons and present part of the light hadron spectrum at finite lattice spacing. Then we study the dependence of the hadron masses on the physical volume and the lattice spacing for the parametrized Fixed-Point Dirac operator. The scaling properties of the vector meson mass is compared to other formulations of lattice fermions. Finally, we investigate how well the continuum energy-momentum hadron dispersion relation is preserved by our lattice action, and examine the effect of overlap-improvement on the spectrum and dispersion relation. The final chapter contains our conclusions and prospects for the future.

The work covered in this thesis is part of an ongoing project of parametrizing, testing and applying Fixed-Point fermions in lattice QCD, carried out in collaboration with Thomas Jörg, Peter Hasenfratz, Ferenc Niedermayer and Kieran Holland. The simulations in the last chapter were performed in the framework of the BGR collaboration. Part of the results presented here have already been published in papers  and conference proceedings . While the focus of this thesis is on simulations of the light hadron spectrum, we will recapitulate some of the basic issues discussed in the PhD thesis of Thomas Jörg  which are relevant for understanding the applications and results in the later chapters in order to keep this work as self-contained as possible.

## Chapter 1 Introduction

This introductory chapter provides some background information for the work covered in the body of the thesis. We start at the very beginning and give a short overview of the history and evolution of the field of elementary particle physics. Then we briefly present in Section 1.2 the foundations of Quantum Chromodynamics, the theory of the strong nuclear force, and introduce the important concepts of symmetries and asymptotic freedom. In order to calculate physical quantities in a quantum field theory, it is necessary to introduce a regularization. The lattice, described in Section 1.3, provides a regularization that allows to probe the non-perturbative regime of strong coupling, where phenomena related to the hadronic world can be examined. We define the most simple lattice actions and the basic tools needed to carry out lattice computations. Finally, in Section 1.4 we present arguments why it is worthwhile to search for improved formulations of lattice QCD. This motivates the construction and application of the Fixed-Point Dirac operator that we perform in this thesis.

### 1.1 The Search for the Fundamental Properties of Nature

Understanding nature is the ultimate goal of every physicist. The basic questions lying at the foundations of a work like this are: How does nature work? Can we explain the phenomena we see? Can we make predictions about what can be seen? From the beginnings of history people have witnessed the phenomena of nature and tried to explain them. Starting at observations accessible to everyday life experience, the interest has moved to objects beyond human perception. At the end of this journey towards finding the fundamental laws of nature, there are two areas: the very small and the very large. The world of the very large is studied in cosmology, where one tries to understand the origin, evolution and fate of the universe as a whole. At the other end of the spectrum one asks what the basic building blocks of the universe are and how they interact. These questions are addressed by the field that is today called elementary particle physics, and it is there where this work tries to add an almost infinitely small fraction to scientific knowledge.

##### The World beneath the Atom

For most people, including those working in sciences like biology and classical chemistry, the smallest structures of interest are atoms or even molecules, and the subatomic world is not considered relevant. This is justified if one is dealing with objects large compared to the atom, but if our interest lies in how nature works at the fundamental level, the fact that the atom is not undividable, as its Greek name implies, can no longer be ignored and the subatomic structure of matter needs to be examined. Thanks to Rutherford’s experiments it has been known for more than 100 years that atoms are built from a tiny nucleus and a surrounding cloud of electrons. Rutherford concluded that the nucleus is made of positively charged particles which he called protons, and for a certain time in the early 20th century, it seemed like with protons and electrons and Einstein’s photon the basic constituents of matter were found. Paul Dirac’s formulation of Quantum Electrodynamics (QED) in 1926 explained beautifully how electrons interact by exchange of photons. However, Dirac’s equation implied the existence of an electron with exactly the same properties, but opposite charge. This looked first as if the theory would be wrong, since such a particle had never been seen before. As a theoretical physicist however, Dirac trusted the beauty of his theory more than the experimental possibilities at that time and drew the conclusion that this antiparticle—the so-called positron—had to exist. Dirac’s prediction turned true when in 1932 the existence of the positron was confirmed in experiments. The observation that our universe is mainly made of matter, and not of antimatter like positrons and antiprotons, is related to a small asymmetry known as CP-violation and is a subject of present research.

There were also a number of other problems which implied that protons, electrons, photons and the electromagnetic force alone were not sufficient to explain the structure of matter. Among them was the unsolved question why the atomic nucleus is stable: Protons are positively charged, so there should be a strong electromagnetic repulsion between the protons in the nucleus, which drives them apart. The newly discovered neutron could not help in solving this problem, as it is not electrically charged and therefore not able to hold the nucleus together. Obviously there had to be some other force which would explain why atomic nuclei didn’t fall into pieces. Another problem was the anomalous magnetic moment of the proton. While for the electron the measurements for this quantity were in perfect agreement with the theoretical prediction of QED, there was almost a factor of 3 difference for the proton, which was a sign that the proton has some non-trivial internal structure and is not an elementary particle. Again, Quantum Electrodynamics alone was not able to explain this phenomenon. Yet another problem was found in the nuclear beta decay, where in an unstable atomic nucleus a proton decays into a neutron and a positron. Here the energy of the positron leaving the nucleus was found to be considerably smaller than the energy difference between the proton and the neutron, and it was not clear where the missing energy was lost. To solve this problem, Wolfgang Pauli postulated in 1931 the existence of the neutrino, an uncharged particle which carries the remaining energy in the beta decay. This particle would be very difficult to observe, as its interactions with other matter are very limited, and in fact the neutrino was experimentally found only in 1956. Altogether, it became clear that while for some time it seemed as if the world of elementary particles was almost fully explained, the theory was obviously not complete and there had to be other, yet unknown mechanisms responsible for these phenomena. The situation changed dramatically with the discovery of a wealth of new particles in cosmic ray observations and in experiments with the newly invented particle accelerators.

##### Handling Elementary Particles

The way to get experimental information on subatomic particles is to collide two particles with as much energy as possible and then to observe what happens. In general new particles are created, and one just needs to detect them and check their properties. In the early 20th century, the only way to observe such high-energy collisions was to wait for cosmic particles to crash into the atmosphere. These particles are emitted in cosmic events like supernovae and therefore carry a lot more energy than what was possible to reach on earth at that time. When such a fast-moving particle hits a nitrogen or oxygen atom of the earth’s atmosphere, the collision products can be examined in suitable detectors. It was in cosmic ray experiments where in 1937 the muon and ten years later the pion and kaon particles were found. Unfortunately almost all cosmic radiation is absorbed in the outmost layers of the atmosphere. Hence for many interesting experiments with cosmic radiation it is necessary to equip a balloon or an airplane with the appropriate instruments and send them into the stratosphere. Furthermore, it is not possible to design a cosmic ray experiment at own will, as the properties of the incoming and the target particles can not be set up freely.

These drawbacks were overcome by the development of particle accelerators. With such a device one takes a particle, accelerates it to very high energies and lets it collide with a target. It is then possible to measure all interesting quantities of the collision products. The accelerated particles, which can be charged particles like electrons or protons, move in a ring-like structure, where they are kept by strong magnetic fields. The larger the diameter of the ring and the stronger the magnetic field, the faster the particles can move and the more energy is set free in the collision. As an example, the LEP collider at CERN which was running until 2001 has a diameter of 27 km and reaches a total energy of 100 GeV in electron-proton collisions. The Large Hadron Collider (LHC) which is under construction at CERN will collide protons and antiprotons at energies of 14 TeV.

Reaching high energies in a collider experiment is crucial because the total energy provides a threshold for the mass of the created particle. If the rest energy of a particle is larger than the total energy of the collided particles, it can not be created in the collision process. Thus for example to create a $`\rho `$ meson, a total energy of 770 MeV, which corresponds to its mass, is required. The problem is that often the particles predicted by theorists have masses too large to be created in current colliders, and therefore larger and larger colliders have to be constructed in order to confirm or falsify the theoretical predictions.

##### Bringing Order into the Chaos

The availability of particle accelerators lead to an enormous growth in the number of newly found particles in the 1950s and 1960s, and there was a definite need for a theory which explained why all these particles were there. All one could do at that time was to bring some order into the wealth of particles and to classify them according to their properties. While most of the particles were very short-lived and had life-times on the order of $`10^{24}`$ seconds, a few of them decayed only after a much longer time of about $`10^{10}`$ s. These particles were called “strange” due to this unexplained property by Murray Gell-Mann in 1953. Gell-Mann found that this whole wealth of particles could be explained in a systematic way when assuming an underlying structure, namely a small number of constituents which, when grouped in different combinations, form the experimentally found particles. These constituents, introduced by Gell-Mann and Zweig, were called quarks, a name taken from James Joyce’s novel “Finnegans Wake”. The quark model could not only explain the known particles, but also predict new ones, which were needed in order to fill the gaps in the tables of possible combinations of quarks. The problem with the quark model was just that no one had ever seen a quark as a separate object in an experiment. All the detected collision products were made out of two or three quarks. In 1973, work of t’Hooft, Politzer, Gross and Wilczek explained this puzzle with the concept of asymptotic freedom, implying that the strong force between two quarks increases when the quarks are pulled apart. In particular, a state with a single quark is not allowed, as it would need infinite energy to separate it from the others. Moreover, when the force between quarks pulled apart reaches a certain threshold, new quarks can be created out of the vacuum, and what remains are again bound states of two or three quarks. Taking the quarks as fundamental building blocks and the color force introduced by Gell-Mann, Fritzsch and Leutwyler as an interaction between the quarks, the quantum theory of the strong nuclear force, Quantum Chromodynamics (QCD), was born. Finally there was a tool to describe the strong interaction, and all the different particles that were found could be explained from common grounds with only a few basic elements and from underlying symmetry principles.

At about the same time, Glashow, Weinberg and Salam developed a quantum theory for the weak interaction, which is responsible for the nuclear beta decay mentioned before. They postulated the existence of the $`W`$ and $`Z`$ bosons as mediators of the weak interaction, and these particles were indeed found at CERN in 1983. Furthermore, the theory of Glashow, Weinberg and Salam allowed unifying the electromagnetic and weak interactions into the so-called electroweak theory. Today, QCD as a theory of the strong interactions and Glashow, Weinberg and Salam’s electroweak theory form the Standard Model (SM) of elementary particle physics, which has been very successful up to date in explaining what nature does at a very small scale. The Standard Model does not include gravitation, which is the last of the four fundamental forces listed in Table 1.1. At the subatomic level however, the gravitational force is negligibly small, and thus it is ignored in SM particle physics. The constituents of matter appearing in the Standard Model are on one hand the six quarks listed in Table 1.2 and the six leptons $`e`$, $`\nu _e`$, $`\mu `$, $`\nu _\mu `$, $`\tau `$, $`\nu _\tau `$, which all are fermions and thus follow the Pauli exclusion principle, and on the other hand the photon, gluon and the $`W^\pm `$ and $`Z`$ particles which are bosons and carry the electromagnetic, strong and weak interactions between the fermions. Finally, the SM predicts the existence of a Higgs particle, which gives a non-vanishing mass to the weak bosons. The existence of the Higgs boson is not yet confirmed by experiment, but it is expected that the particle will be found as soon as the next generation of colliders start operation.

It is obvious that the Standard Model is not yet the ultimate theory of nature, not only because it does not contain gravity, but also because quite a large number of unknown input parameters are needed. Therefore many theoretical physicists work on finding candidates for an even more fundamental theory that unifies all the four interactions. These attempts lead to exciting discoveries like superstring theories living in 10-dimensional space-time, and more recently 11-dimensional $`M`$-Theory. While from the theoretical point of view these theories are very attractive, from what we know today it is extremely difficult to connect them to phenomenological information and thus to test their predictions, as the typical energy scales involved are far beyond reach of any foreseeable experiment.

In the following, we will stay within the bounds of the Standard Model. We concentrate on the strong interaction and the particles participating therein, the quarks and gluons. Many fundamental questions in particle physics are related to the strong force, hence the study of Quantum Chromodynamics is a highly rewarding task, both from the phenomenological and the theoretical point of view.

### 1.2 Quantum Chromodynamics

The strong interactions between elementary particles are described by Quantum Chromodynamics (QCD), the quantum theory of the color force. The basic degrees of freedom of the theory are the quark and gluon fields. Like all quantum field theories in the Standard Model, QCD is a local gauge theory. The gluons, which are the gauge fields of the theory, are introduced to ensure local gauge invariance and thus generate the interaction among the particles. The gauge group has to be chosen as an external input when constructing the theory. From particle phenomenology follows that the quarks appear in three different colors, and that in nature the gauge group of the color force is the special unitary group $`SU(3)`$.<sup>1</sup><sup>1</sup>1 As a theoretical generalization, the theory can also be set up with the gauge group $`SU(N_c)`$ for an arbitrary number of colors $`N_c`$. The beauty and strong predictive power of QCD lies in the fact that only a small number of parameters need to be fixed to define the theory and to get physical predictions.

#### 1.2.1 The QCD Lagrangian

The fermions from which QCD is constructed are the $`n_f`$ flavors of quark fields $`q^k(x)\{u,d,s,c,t,b\}`$, $`k=1,\mathrm{},n_f`$, which are Grassmann-valued Dirac spinors and $`SU(3)`$ triplets in color space. Thus, under a local gauge transformation $`U(x)SU(3)`$ the quark and antiquark fields $`\overline{q}^k=(q^k)^{}\gamma _0`$ transform like

$`q^k(x)`$ $``$ $`U(x)q^k(x),`$ (1.1)
$`\overline{q}^k(x)`$ $``$ $`\overline{q}^k(x)U^{}(x),`$ (1.2)
The gauge bosons are the $`N_c^21`$ gluon fields $`A_\mu ^a(x)SU(N_c)`$.

A field theory is defined by its Lagrangian density $``$, from which the equations of motion and thus the dynamics of the theory can be derived. The QCD Lagrangian

$$_{QCD}(x)=_F(x)+_G(x),$$
(1.3)
can be split into the fermionic (quark) part

$$_F(x)=\underset{k=1}{\overset{n_f}{}}\overline{q}^k(x)(i\gamma ^\mu D_\mu m)q^k(x),$$
(1.4)
and the purely gluonic part

$$_G(x)=\frac{1}{4}F_{\mu \nu }^a(x)F^{\mu \nu a}(x),$$
(1.5)
which in itself defines a non-trivial Yang-Mills theory and describes the kinematics of the gluons. The sum over the repeated color index runs from $`a=1,\mathrm{},N_c^21`$. The gluon field strength tensor appearing in the Lagrangian $`_G`$ is defined by

$$F_{\mu \nu }^a(x)=_\mu A_\nu ^a(x)_\nu A_\mu ^a(x)g_sf_{abc}A_\mu ^b(x)A_\nu ^c(x),$$
(1.6)
where $`g_s`$ is the strong coupling constant and $`f_{abc}`$ are the structure constants of the gauge group $`SU(N_c)`$. To ensure local gauge invariance, in the fermionic Langrangian $`_q(x)`$ the covariant derivative

$$D_\mu (x)=_\mu ig_sA_\mu (x),$$
(1.7)
has to be taken, with the gauge field $`A_\mu (x)`$ being an element of the gauge group,

$$A_\mu (x)=A_\mu ^a(x)\frac{\lambda ^a}{2},$$
(1.8)
where the group generators $`\lambda ^a/2`$ follow the commutation relation

$$[\frac{\lambda ^a}{2},\frac{\lambda ^b}{2}]=if_{abc}\frac{\lambda ^c}{2}.$$
(1.9)
Requiring the Lagrangian (1.3) to be gauge invariant, the transformation rules for the gauge field and the field strength tensor are found to be

$`A_\mu (x)`$ $``$ $`U(x)A_\mu (x)U^{}(x){\displaystyle \frac{1}{g_s}}_\mu U(x)U^{}(x),`$ (1.10)
$`F_{\mu \nu }(x)`$ $``$ $`U(x)F_{\mu \nu }(x)U^{}(x).`$ (1.11)
Having specified the QCD Lagrangian, the theory is defined, and it remains to prescribe how to extract physical quantities. This is done most elegantly in the Feynman path integral formalism, thus promoting the classical field theory to a quantum theory. Let us now switch to Euclidean space (see Appendix E.2), which will be natural for setting up a lattice formulation. Expectation values for physical observables $`𝒪`$, that can be arbitrary operators built from quark and gluon fields, are defined by the path integral

$$𝒪=\frac{1}{Z}D\overline{q}DqDA𝒪e^{S^E[\overline{q},q,A]},$$
(1.12)
where the normalization in the denominator is given by the partition function

$$Z=D\overline{q}DqDAe^{S^E[\overline{q},q,A]}.$$
(1.13)
Eq. (1.13) shows that a quantum field theory in imaginary time formally resembles a system in classical statistical mechanics, where the probability of a state is proportional to the Boltzmann factor $`\mathrm{exp}(E/kT)`$. In QCD, the Euclidean action

$$S^E[\overline{q},q,A]=d^4x_{QCD}^E(x),$$
(1.14)
where $`_{QCD}^E`$ is the QCD Lagrangian transformed to Euclidean space, appears in the exponent of the Boltzmann factor.

#### 1.2.2 Global Vector and Axial Symmetries

In the limit of $`n_f`$ massless quarks, the QCD Lagrangian (1.3)–(1.5) exhibits a global symmetry

$$U_V(1)\times SU_V(n_f)\times U_A(1)\times SU_A(n_f),$$
(1.15)
acting on the flavor and spin degrees of freedom. Writing the $`n_f`$ quark fields as a vector, the corresponding symmetry transformations are

$$q(x)e^{i\varphi (T_DT_F)}q(x),$$
(1.16)
where $`T_D\{\mathrm{𝟏},\gamma _5\}`$ acts on the Dirac structure and generates the vector ($`V`$) or axial vector ($`A`$) transformations, and $`T_F`$ works in flavor space to create the $`U(1)`$ (for $`T_F=\mathrm{𝟏}`$) or $`SU(n_f)`$ transformations. The conserved currents related to these global symmetries through the Noether theorem are

$$j^\mu (x)=\overline{q}(x)\gamma ^\mu (T_DT_F)q(x).$$
(1.17)
The vector $`U_V(1)`$ symmetry is unbroken even for finite quark mass and gives rise to baryon number conservation. The $`SU_V(n_f)`$ leads to the multiplet structure of the hadrons. The axial $`U_A(1)`$ is explicitly broken on the quantum level by instanton contributions, leading to the Adler-Bell-Jackiw (ABJ) anomaly  of the flavor-singlet axial current and the massiveness of the $`\eta ^{}`$ meson. The $`SU_A(n_f)`$ is believed to be spontaneously broken by a non-zero vacuum expectation value of the quark condensate $`\overline{q}q`$, and the associated $`(n_f^21)`$ massless Goldstone bosons for $`n_f=2`$ are the pions. In the real world, the global symmetry (1.15) arises from the smallness of the light quark masses (see Table 1.2), where setting $`m_u=m_d=0`$ and in some cases even $`m_s=0`$ is a good approximation. For non-zero, but small quark masses, the pions are no longer real Goldstone bosons, but quasi-Goldstone particles that acquire a small mass. This would explain why the experimentally observed pion masses $`m_{\pi ^0}=135`$ MeV and $`m_{\pi ^\pm }=140`$ MeV are so small compared to the masses of other hadrons.

#### 1.2.3 Asymptotic Freedom

The coupling constant $`g_s`$ of the strong interaction is actually not a constant, but depends on the momentum transfer $`Q`$ of a given process through quantum corrections, leading to the emergence of a generic scale $`\mathrm{\Lambda }`$. Often not the coupling constant $`g_s`$ itself is used, but the fine-structure constant $`\alpha _s=g_s^2/4\pi `$, which is to leading order given by

$$\alpha _s(Q^2)=\frac{\alpha _s(\mathrm{\Lambda })}{1+\alpha _s(\mathrm{\Lambda })\frac{332n_f}{12\pi }\mathrm{ln}(\frac{Q^2}{\mathrm{\Lambda }^2})}.$$
(1.18)
In the running coupling (1.18), asymptotic freedom of QCD shows up in the fact that $`\alpha _s`$ gets small at large momenta $`Q`$. On the other hand, the coupling increases with larger momenta or equivalently smaller distance, leading to confinement of quarks<sup>2</sup><sup>2</sup>2Quark confinement is a non-perturbative phenomenon which does not follow from perturbation theory.. At the mass of the $`Z`$-boson $`m_Z=91`$ GeV, measurements of the coupling constant give a value of $`\alpha _s(m_Z)=0.118`$, which is a reasonably small value that a perturbative expansion in $`\alpha _s`$ around the free theory makes sense. For deep inelastic scattering processes studied in collider experiments, the momentum transfer is of this order, so in this region QCD can be treated perturbatively. At scales around 1 GeV however that are typical for the hadronic world, $`\alpha _s`$ is on the order of 1 and thus no longer a small parameter in which an expansion is possible. Perturbation theory therefore breaks down when small momenta or large distances are involved. In this non-perturbative region of Quantum Chromodynamics, where one would like to investigate issues like the hadron spectrum, hadronic matrix elements of operators, spontaneous chiral symmetry breaking, confinement or the topological structure of the vacuum, it is necessary to use another approach to perform calculations. This is where lattice QCD comes into play.

### 1.3 QCD on the Lattice

The lattice formulation of Quantum Chromodynamics in Euclidean space, originally proposed in 1974 by Wilson , was designed as a tool to calculate observables in the non-perturbative region of QCD from first principles. Lattice QCD is at present the only method which allows to compute low-energy hadronic quantities in terms of the fundamental quark and gluon degrees of freedom without having to tune additional parameters. The only input parameters are the bare quark masses and the bare coupling constant, and from these all other quantities like the masses of the hundreds of experimentally observed hadrons can be calculated. Hence the lattice is a very powerful tool in checking that QCD is the correct theory for the strong interactions and in making predictions for the dependence of hadronic quantities on the input parameters. Formulating QCD on a discrete space-time lattice opens the possibility to treat these problems on computers, using methods analogous to those in Statistical Mechanics. However, due to the large number of degrees of freedom involved, lattice QCD simulations are computationally very demanding, and it is still necessary to use a number of tricks and approximations in order to cope with these demands. It is then also important to examine whether the effects introduced by these approximations are under control. Since the first numerical measurements in a lattice gauge theory by Creutz, Jacobs and Rebbi in 1979 , the progress in computer technology and the theoretical developments in the field have allowed to get closer to examining in a systematical manner the deep questions which Lattice QCD is able to answer. We present in the following a brief introduction to the basics of lattice QCD that is necessary to follow the rest of the work. For a more detailed discussion, we refer to the standard textbooks  or recent introductory articles . An extensive overview of the status of current research in lattice QCD can be found in the proceedings of the annual lattice conference .

#### 1.3.1 The Lattice Regularization

Quantum field theories have to be regularized in order to give the path integrals in Eq. (1.12) that define physical observables a meaning. In perturbation theory, a convenient way to do this is by dimensional regularization, where the space-time dimension $`d`$ is modified by a small parameter $`ϵ`$ to $`d=4ϵ`$, or by introducing a momentum cut-off $`\mathrm{\Lambda }`$. At the end of a calculation, the regularization has to be removed by taking the limit $`ϵ0`$ or $`\mathrm{\Lambda }\mathrm{}`$. The lattice is nothing else than such a regulator for the theory. In the lattice regularization, the continuum Euclidean space-time variable $`x_\mu `$ is replaced by a discrete hypercubic space-time lattice,

$$x_\mu n_\mu a,n_\mu 𝐙,$$
(1.19)
with lattice spacing $`a`$. This introduces an ultraviolet cut-off by restricting the momenta to lie within the Brioullin zone $`|p_\mu |<\pi /a`$, removing the ultraviolet divergent behaviour of integrals. Restricting the space-time extent to a finite lattice $`n_\mu <N_\mu `$, $`n_\mu 𝐍_\mathrm{𝟎}`$, the momenta take discrete values $`p_\mu =k_\mu \pi /N_\mu a`$, with $`|k_\mu |<N_\mu `$ and $`k_\mu 𝐙`$. Every quantity that is calculated in the lattice regularized theory is finite, since the integrals are transformed into finite sums.

The continuum quark and gluon fields are replaced by lattice fields living on the sites and links of the lattice, respectively, and also the derivatives in the QCD Lagrangian (1.3)–(1.5) have to be discretized in some way. It is obvious that in the process of discretization some of the original symmetries of the Lagrangian are partially or fully lost. As an example, the Poincaré symmetry in continuum space-time is replaced by a cubic symmetry on the lattice. As mentioned before, it is necessary to remove the regularization at the end of the calculation to get physical results, and for the lattice regularization this means that the continuum limit $`a0`$ has to be taken. In this process one expects the lost symmetries to be restored. However, one requires that the most important symmetries like gauge invariance, which lies at the foundations of QCD, are also present at finite lattice spacing. The lattice formulation of the QCD Lagrangian should therefore respect these symmetries.

The lattice quark and antiquark field are Grassmann variables $`\mathrm{\Psi }(n)`$, $`\overline{\mathrm{\Psi }}(n)`$ defined at every lattice site $`n`$. In natural units physical quantities can be expressed in units of powers of length or inverse mass. For numerical applications, it is convenient to work with dimensionless quantities. This can be done by absorbing the dimension through appropriate powers of the lattice spacing $`a`$. For the quark fields, the transcription on the lattice is then given by

$$q(x)a^{3/2}\mathrm{\Psi }(n).$$
(1.20)
The lattice gauge fields $`U_\mu (n)`$ are defined by the path-ordered Schwinger line integral

$$U_\mu (n)=𝒫\mathrm{exp}\left(ig_{na}^{(n+\widehat{\mu })a}𝑑xA_\mu (x)\right),$$
(1.21)
acting as parallel transporters of color between neighbouring lattice points. They are thus elements of the gauge group $`SU(3)`$ defined on the links between lattice sites. To lowest order in the lattice spacing, (1.21) reduces to

$$U_\mu (n)e^{iagA_\mu (x)}.$$
(1.22)
Under a gauge transformation $`G(n)`$, the lattice quark and gluon fields transform like

$`\mathrm{\Psi }(n)`$ $``$ $`G(n)\mathrm{\Psi }(n),`$ (1.23)
$`\overline{\mathrm{\Psi }}(n)`$ $``$ $`\overline{\mathrm{\Psi }}(n)G^{}(n),`$ (1.24)
$`U_\mu (n)`$ $``$ $`G(n)U_\mu (n)G^{}(n+\widehat{\mu }).`$ (1.25)
The construction of the lattice gauge fields is done such in order to ensure gauge invariance of non-local quark operators. With these definitions, there are two different types of gauge invariant objects: color traces of closed loops of gauge links like the Wilson plaquette

$$U_{\mu \nu }(n)=U_\mu (n)U_\nu (n+\widehat{\mu })U_\mu ^{}(n+\widehat{\nu })U_\nu ^{}(n),$$
(1.26)
and quark bilinears like $`\overline{\mathrm{\Psi }}(n)U_\mu (n)\mathrm{\Psi }(n+\mu )`$, where gauge links connect the quark and the antiquark field along an arbitrary path in space-time.

Integrals over continuum space-time variables, as they appear in the action (1.14), are replaced by sums over the lattice sites,

$$d^4xf[\overline{q}(x),q(x),A_\mu (x)]a^4\underset{n}{}\widehat{f}[\overline{\mathrm{\Psi }}(n),\mathrm{\Psi }(n),U_\mu (n)],$$
(1.27)
where $`\widehat{f}`$ is a discretized version of the function $`f`$. The path integral over the quark and gluon fields in the expression for the expectation value of observables (1.12) and in the partition function (1.13) is transformed into a product of ordinary integrals over the fields at all lattice sites,

$$D\overline{q}DqDA\underset{\alpha ,l}{}𝑑\overline{\mathrm{\Psi }}_\alpha (l)\underset{\beta ,m}{}𝑑\mathrm{\Psi }_\beta (m)\underset{\rho ,n}{}𝑑U_\rho (n),$$
(1.28)
which yields finite expressions on a finite lattice and can be evaluated numerically.<sup>3</sup><sup>3</sup>3 The integration over the gluon fields is an integration over the gauge group SU(3).

#### 1.3.2 Simple Lattice Actions

Although the discretization of the continuum QCD Lagrangian (1.3)–(1.5) might appear trivial at first sight, there are some complications. The lattice action given by Wilson  is the most simple working version and is still widely used in simulations, although it is not free of problems, as we will see later. The action

$$S[\overline{\mathrm{\Psi }},\mathrm{\Psi },U]=S_G[U]+S_F[\overline{\mathrm{\Psi }},\mathrm{\Psi },U],$$
(1.29)
can again be split in separate gauge and fermion parts. The Wilson gauge action $`S_G^{(\mathrm{W})}`$ is constructed from the plaquette $`U_{\mu \nu }`$ in (1.26) by

$$S_G^{(\mathrm{W})}[U]=\beta \underset{n}{}\underset{\mu <\nu }{}\left(1\frac{1}{N_c}\mathrm{Re}\mathrm{Tr}U_{\mu \nu }(n)\right),$$
(1.30)
which in the limit $`a0`$ goes over to the continuum form up to $`𝒪(a^2)`$ errors. The parameter $`\beta =2N_c/g_s^2`$ takes over the role of the bare coupling constant. The fermionic lattice action

$$S_F[\overline{\mathrm{\Psi }},\mathrm{\Psi },U]=\underset{n,n^{}}{}\overline{\mathrm{\Psi }}(n)D(n,n^{})\mathrm{\Psi }(n^{}),$$
(1.31)
is bilinear in the quark fields. The Wilson fermion action $`S_F^{(\mathrm{W})}`$ is defined by setting $`D=D^\mathrm{W}`$, with the Wilson Dirac operator

$$\begin{array}{c}D^\mathrm{W}(n,n^{})=\frac{1}{2a}\underset{\mu }{}\left[(\gamma _\mu r)\delta _{n^{},n+\widehat{\mu }}U_\mu (n)(\gamma _\mu +r)\delta _{n,n^{}+\widehat{\mu }}U_\mu ^{}(n\widehat{\mu })\right]\hfill \\ \hfill +\left(m+\frac{4r}{a}\right)\delta _{nn^{}},\end{array}$$
(1.32)
where the bare quark mass $`m`$ is another parameter of the theory. The Wilson Dirac operator differs from the naive discretization of the Euclidean continuum Dirac operator $`\gamma _\mu D_\mu +m`$ by a dimension $`d=5`$ term proportional to the unphysical parameter $`r`$,

$$S_F^{(\mathrm{W})}=S_F^{(\mathrm{naive})}a\frac{r}{2}\underset{n}{}\overline{\mathrm{\Psi }}(n)\mathrm{}\mathrm{\Psi }(n),$$
(1.33)
which is called the Wilson or doubler term and is needed to remove unphysical particles that appear through poles at the corners of the Brioullin zone from the spectrum. In the continuum limit $`a0`$, the Wilson term vanishes as required. However, as we show later, the cost for introducing the Wilson term is the explicit breaking of chiral symmetry, leading to many theoretical and practical problems and limitations.

#### 1.3.3 Monte Carlo Integration

The integral over the fermion fields, which are anticommuting Grassmann variables, in the lattice version of the partition function (1.13) can be performed analytically. For a bilinear fermion action (1.31), the integration over quark and antiquark fields gives the determinant of the fermion matrix, and the partition function on the lattice then reads

$$Z=\underset{\mu ,n}{}𝑑U_\mu (n)detDe^{S_G[U]},$$
(1.34)
where $`D`$ is the lattice Dirac operator and $`S_G[U]`$ is a lattice version of the gluon action. Expectation values for observables $`𝒪`$ are calculated from

$$𝒪=\frac{1}{Z}\underset{\mu ,n}{}𝑑U_\mu (n)𝒪detDe^{S_G[U]}.$$
(1.35)
At this point, it is obvious that the theory is ready to be put on a computer, since in Eq. (1.35) only an integration over the $`SU(3)`$ gauge fields is left. However, for standard numerical integration, the number of degrees of freedom is still far too large, therefore one has to resort to statistical methods. The way the gauge field integral is usually handled is by Monte Carlo integration: A finite number $`N`$ of gauge configurations $`U^{(i)}`$, $`(i=1,\mathrm{},N)`$, are statistically sampled with the probability distribution given by the fermion determinant $`detD`$ times the Boltzmann factor $`\mathrm{exp}(S_G[U^{(i)}])`$. Observables are then estimated from the sample mean

$$𝒪\overline{𝒪}=\frac{1}{N}\underset{i=1}{\overset{N}{}}𝒪[U^{(i)}].$$
(1.36)
In practice this amounts to generating a set of $`N`$ independent gauge configurations with a Markov chain algorithm that respects the required probability distribution and measuring the observable on the resulting set of gauge configurations. As we have seen, after integrating out the fermions in the partition function (1.34), their influence on the weighting is given by the determinant. It turns out that the calculation of the fermion determinant is by far the most time-consuming part in a lattice simulation. This is why most lattice QCD calculations up to date have been done in the quenched approximation, which is explained in the next section.

#### 1.3.4 The Quenched Approximation

The determinant of the Dirac operator is a non-local quantity. Even for moderate lattice sizes, its exact calculation is not feasible on today’s computers. Various algorithms have been developed to tackle this problem, but keeping dynamical fermion loops in the simulation is still a very demanding task. The easiest way out is to consider quenched QCD, as done in this work, where the fermion determinant in Eq. (1.35) is set to

$$detD=1.$$
(1.37)
This approximation is equivalent to making the virtual quarks infinitely heavy, leading to the complete suppression of internal quark loops. Neglecting the fermion determinant simplifies the technical treatment enormously, as then in the generation of the gauge configurations the probability distribution is given by the Boltzmann factor alone, whereas in unquenched QCD the determinant has to be calculated in every Monte Carlo update step. It is however obvious that quenched QCD is not the correct theory to describe nature, as for example two quarks can be pulled apart to an arbitrary distance<sup>4</sup><sup>4</sup>4The energy needed increases linearly with distance. in quenched QCD, while in nature at a certain point two additional quarks are created and string breaking occurs. Quenched QCD is not even mathematically clean, as it is not unitary. The only reason for using the quenched approximation is that the computation of the fermion determinant is extremely demanding, and by setting $`detD=1`$ a factor of several orders of magnitude in time is gained.

The reason why the investigation of quenched lattice QCD is nevertheless interesting is that since neglecting of the determinant amounts just to a different weighting of the gauge configurations, the quenched theory still shows the crucial properties of QCD like asymptotic freedom and spontaneous symmetry breaking. Therefore it is possible to examine many non-trivial questions first in the quenched theory, giving qualitative hints what in full, unquenched QCD might occur. A phenomenological argument why ignoring virtual quark loops is not a completely useless approximation is given by the Zweig rule, which states that processes where the constituent quarks do not survive are suppressed. Many years of lattice simulations have shown that the errors due to quenching in physical observables are in most cases only on the order of 10%, allowing to make also quantitative predictions. However, it is very important that quenching effects are well investigated.

#### 1.3.5 Continuum Limit, Renormalization and Scaling

The parameters $`g_s`$ (or $`\beta `$) and $`m`$ which are put into a lattice simulation are bare quantities, and when taking the continuum limit $`a0`$, physical quantities and not the bare parameters have to be kept fixed. Consider the physical observable $`𝒪_{\mathrm{phys}}`$ with mass dimension $`d_𝒪`$ and its dimensionless lattice counterpart $`𝒪_{\mathrm{lat}}`$, which depends on the lattice spacing through the coupling $`g_s(a)`$ and the quark mass $`m(a)`$. The continuum limit

$$𝒪_{\mathrm{phys}}=\underset{a0}{lim}a^{d_𝒪}𝒪_{\mathrm{lat}}(g_s(a),m(a)),$$
(1.38)
is taken by measuring $`𝒪_{\mathrm{lat}}`$ at different values of $`g_s`$. At sufficiently small $`a`$, the dependence of the coupling constant on the lattice spacing $`g_s(a)`$ should be a universal function, independent of the observable under consideration. The same holds for the function $`m(a)`$ for the quark mass in unquenched QCD. This property is called scaling, and the range of lattice spacings or gauge couplings the hypothesis is valid is called the scaling window. Lattice simulations have to be performed within this scaling window in order to extract reasonable continuum results.

The scale dependence can be removed by forming ratios of particle masses

$$\frac{am_1(a)}{am_2(a)}=\frac{m_1(0)}{m_2(0)}+𝒪(m_1a)+𝒪((m_1a)^2)+\mathrm{}.$$
(1.39)
The $`a`$-dependent terms on the right hand side are artifacts from discretization errors and depend on the choice of lattice action. If these terms are small, a controlled continuum extrapolation is possible, and we speak of scaling of the quantity under consideration. It is certainly desirable that a lattice action shows good scaling, that is small scaling violations.

When results in physical units are wanted, the lattice results, which are always dimensionless, have to be converted to physical units by matching the result of one observable with the experimental data. This observable might be the mass of a hadron like $`m_\rho 770`$ MeV or $`m_N940`$ MeV or a decay constant like $`f_\pi 93`$ MeV. In the quenched approximation, the lattice spacing does not depend on the quark mass, as there are only external quarks. It is then also possible to fix the scale from a purely gluonic quantity like the string tension $`\sqrt{\sigma }420`$ MeV. More reliable than the string tension are Sommer-type scales , which are also related to the quark-antiquark potential.

### 1.4 Why Improved Formulations of Lattice QCD?

In principle the choice of how to discretize the continuum QCD Lagrangian (1.3)–(1.5) is free, as long as the correct continuum limit is reached. However, there are a number of reasons why it is worthwhile to search for improved lattice formulations of the Lagrangian. One reason is the reduction of discretization errors: Working at finite lattice spacing $`a`$ introduces discretization errors which affect simulation results and have to be identified and removed. Simulations are normally carried out at lattice spacings between 0.05 fm and 0.2 fm, and the typical size of hadrons is on the order of 1 fm. As explained above, the continuum limit (1.38) is taken by measuring observables at several lattice spacings and extrapolating the results to $`a0`$. While the discretization errors should disappear in the continuum limit, it is advantageous to have a lattice formulation of the theory with small discretization errors. First of all it is often not clear how controlled and safe the continuum extrapolation is, therefore having results which show smaller $`a`$-dependence leads to a more reliable extrapolation. On the other hand, working at small lattice spacings is numerically very demanding, as the computational effort grows roughly like $`a^6`$ for quenched simulations and like $`a^{10}`$ for the unquenched case , thus it would obviously be helpful to have a formulation of the theory which gives results of the same quality at larger lattice spacing. Hence, a lattice QCD action with small lattice artifacts allows either for a more reliable extrapolation or to work at larger lattice spacings and to save computation time. The usual way to reduce discretization errors is to improve the lattice action and operators systematically in orders of the lattice spacing by adding irrelevant terms which remove the artifacts order by order. This improvement program, proposed by Symanzik , has been applied to various cases, the best known of which is the $`𝒪(a)`$-improved Wilson clover fermion action .

Improved actions are also expected to show better behaviour in restoring rotational and internal symmetries. Most prominent is the $`U_A(1)SU_A(n_f)`$ chiral symmetry in (1.15), which is explicitly broken for Wilson-type fermions by the Wilson term in Eq. (1.33). Even in the continuum limit, this explicit breaking of chiral symmetry leads to unwanted effects like an additive renormalization of the quark mass, which means that in lattice simulations the bare quark mass is a parameter which needs to be tuned. Another consequence of explicit chiral symmetry breaking in the quenched theory is the appearance of exceptional configurations, for which the quark propagator diverges although the bare quark mass is still far from the critical value. This makes it impossible to simulate quarks much lighter than the strange quark, and therefore a long and unreliable chiral extrapolation from the simulated quark masses to the physical masses of the up and down quarks is needed. The partially conserved axial vector current also needs to be renormalized. Furthermore, mixing between operators of different chiral representations occurs, leading to technical difficulties in calculations of weak matrix elements. There also exists a close connection between chiral zero modes of the Dirac operator and the topological structure of the gauge fields, and with standard formulations of lattice QCD neither topology nor chiral fermion zero modes are well-defined notions. For a long time, it has been believed that chiral symmetry can not be preserved on the lattice. Only after the resurrection of the Ginsparg-Wilson relation , it has been realized that it is possible to retain an exact, slightly modified chiral symmetry  on the lattice.

A radical approach to improvement is the classically perfect Fixed-Point action , which is defined at the fixed point of Renormalization Group transformations. The Fixed-Point action gives exact continuum results for classical predictions even at non-zero lattice spacing. Thus it allows for scale invariant instanton solutions, satisfies the fermionic index theorem and preserves chiral symmetry . Even for quantum results, discretization errors are expected to be considerably reduced. In this work we will construct and apply a parametrization of the Fixed-Point fermion action. Since we will also make use of a recently constructed Fixed-Point action  for the gluons, the results in this thesis serve as a first extensive test for Fixed-Point actions in QCD.

## Chapter 2 Chiral Fermions and Perfect Actions

Only very recently, it has become possible to simulate chiral fermions on the lattice. This exciting discovery lead to growing activity in applying and testing chiral lattice actions. In this chapter we recapitulate the problems with formulating chiral lattice fermions and different solutions, which all obey the ubiquitous Ginsparg-Wilson relation. Fixed-Point (FP) fermions are not only chiral, but also classically perfect. We present the conceptual basics of perfect actions and the application to free fermions.

### 2.1 Chiral Symmetry on the Lattice

In Section 1.2.2 we have presented the global flavour symmetries inherent in the continuum QCD Lagrangian. In this section we describe the problems arising when the theory is transcribed onto the lattice, and how it is possible to retain chiral symmetry in lattice QCD. Consider the global flavour-singlet $`U_V(1)`$ vector transformation

$`\mathrm{\Psi }(n)`$ $``$ $`e^{i\varphi }\mathrm{\Psi }(n),`$ (2.1a)
$`\overline{\mathrm{\Psi }}(n)`$ $``$ $`\overline{\mathrm{\Psi }}(n)e^{i\varphi },`$ (2.1b)
and the $`U_A(1)`$ axial vector transformation

$`\mathrm{\Psi }(n)`$ $``$ $`e^{i\varphi \gamma _5}\mathrm{\Psi }(n),`$ (2.2a)
$`\overline{\mathrm{\Psi }}(n)`$ $``$ $`\overline{\mathrm{\Psi }}(n)e^{i\varphi \gamma _5},`$ (2.2b)
acting on the lattice fermion fields. It is obvious that the fermion lattice action $`_{n,n^{}}\overline{\mathrm{\Psi }}(n)(D(n,n^{})+m)\mathrm{\Psi }(n^{})`$ satisfies the $`U_V(1)`$ symmetry for all quark masses $`m`$. Setting $`m=0`$, the chiral $`U_A(1)`$ symmetry is present only if the Dirac operator anticommutes with $`\gamma _5`$:

$$\{D,\gamma _5\}=0.$$
(2.3)
The major obstacle for formulating lattice fermions respecting chiral symmetry is the Nielsen-Ninomiya no-go theorem , which states that it is not possible to have a lattice Dirac operator which is local, has the correct continuum limit, is free of doublers and satisfies Eq. (2.3). If the continuum fermion action is discretized naively, chiral symmetry is preserved, but instead of one fermion there appear 16 massless particles. To remove these doublers, in the Wilson action (1.33) a term is added to the action which gives the doublers a mass, but breaks chiral symmetry explicitly by violating the anticommutation relation (2.3). It is clear that all the other properties in the Nielsen-Ninomiya theorem have to be conserved in order to obtain a reasonable lattice Dirac operator, and therefore breaking chiral symmetry seems to be the only way to get around the theorem. However, instead of the hard breaking by the Wilson term, a better approach is to slightly modify Eq. (2.3) to the so-called Ginsparg-Wilson relation

$$\{D,\gamma _5\}=aD\gamma _52RD,$$
(2.4)
where the newly introduced term on the right hand side vanishes in the continuum limit. It is useful to express Eq. (2.4) in terms of the quark propagator:

$$\{D^1,\gamma _5\}=a\gamma _52R.$$
(2.5)
The term $`2R`$, which is denoted like this for historical reasons, is a local operator. From this requirement follows that Eq. (2.5) is a highly non-trivial condition, since the quark propagator $`D^1`$ on the left-hand side is a non-local object.

It has been shown by Lüscher  that for a Dirac operator fulfilling the Ginsparg-Wilson relation (2.4), it is possible to define an exact lattice chiral symmetry, which is a modified version of the continuum $`U_A(1)`$ symmetry. When the transformation (2.2) is replaced by

$`\mathrm{\Psi }(n)`$ $``$ $`e^{i\varphi \gamma _5(1aD/2)}\mathrm{\Psi }(n),`$ (2.6a)
$`\overline{\mathrm{\Psi }}(n)`$ $``$ $`\overline{\mathrm{\Psi }}(n)e^{i\varphi (1aD/2)\gamma _5},`$ (2.6b)
the fermion action is invariant. An analogous statement holds for the flavour non-singlet axial transformation. At this point, it might seem that there is more symmetry than expected, because due to the ABJ anomaly the $`U_A(1)`$ symmetry should be broken at the quantum level. The solution comes from the observation that the fermionic integration measure is not invariant under the modified transformation (2.6), but transforms like

$$D\overline{\mathrm{\Psi }}D\mathrm{\Psi }\mathrm{exp}(2N_f\times \mathrm{index}(D))D\overline{\mathrm{\Psi }}D\mathrm{\Psi },$$
(2.7)
thus creating the expected anomaly for topologically non-trivial gauge configurations. The fermionic index in (2.7),

$$\mathrm{index}(D)n_{}n_+,$$
(2.8)
is the difference between the number of zero eigenmodes of the Dirac operator with positive and negative chirality, and is related to the topological charge $`Q_{\mathrm{top}}`$ through the Atiyah-Singer index theorem

$$\mathrm{index}(D)=Q_{\mathrm{top}}=\frac{1}{32\pi ^2}d^4xϵ_{\mu \nu \rho \sigma }\mathrm{tr}(F_{\mu \nu }F_{\rho \sigma }).$$
(2.9)
### 2.2 Fermions with Exact or Approximate Chiral Symmetry

While the Ginsparg-Wilson relation (2.4) has been known for a long time, no solution was found until recently, when three different formulations were independently discovered which all fulfill the Ginsparg-Wilson relation and thus retain exact chiral symmetry on the lattice. These solutions are the domain wall  and overlap fermions , which were originally proposed to formulate chiral gauge theories, and the Fixed-Point fermions .

Domain wall fermions are defined by extending Wilson fermions into a non-physical fifth dimension with lattice spacing $`a_s`$, lattice size $`N_s`$ and a negative mass. The different chiralities are then located on the two opposite domain walls, with the mixing exponentially suppressed by the size of the fifth dimension $`N_s`$. In the limit $`N_s\mathrm{}`$, exact chiral symmetry is obtained.

For overlap fermions, there exists an explicit construction with exact chiral symmetry: Defining the kernel

$$A=1aD^\mathrm{W},$$
(2.10)
with the Wilson operator $`D^\mathrm{W}`$, the overlap Dirac operator $`D^{(\mathrm{ov})}`$ is given by

$$D^{(\mathrm{ov})}=\frac{1}{a}\left(1\frac{A}{\sqrt{A^{}A}}\right).$$
(2.11)
Domain wall fermions with infinite fifth dimension $`N_s`$ are in fact equivalent to overlap fermions, when a different kernel $`A`$  is put into (2.11). It is not obvious that the overlap construction generates a local operator, which requires that the couplings decrease exponentially with distance. Losing locality would render the whole formulation useless. However, it has been shown that both the overlap operator with Wilson kernel (2.11) and the 4-d effective formulation for the domain wall operator are local .

Fixed-Point fermions are defined through Renormalization Group transformations, as discussed in detail in Section 2.3. It has been first shown for FP fermions that the index theorem on the lattice remains valid . For FP fermions, there is no explicit expression, except for the non-interacting case. They have to be constructed in an iterative procedure, which we present in Section 3.3. An important difference to domain wall and overlap fermions is that FP fermions do not only respect chiral symmetry, but are classically perfect and therefore are expected to have small cut-off effects.

Having presented lattice fermion formulations with exact chiral symmetry, it is important to show where approximations have to be taken which introduce again some residual chiral symmetry breaking. For domain wall fermions, it is obvious that the extension of the fifth dimension $`N_s`$ has to be finite in actual simulations. Mixing of the two chiralities is then still possible, and in recent simulations by the RBC  and CP-PACS  collaborations and in  the effects of this residual chiral symmetry breaking have been investigated closely. The exponential decay was found to be surprisingly slow, and although rather large extensions $`N_s`$ of $`𝒪`$(32–64) have been used in the simulations, getting close to the limit of exact chiral symmetry does not seem to be easily possible. This might be an effect arising from the small eigenvalues of the hermitean Wilson operator . Although there are several proposals how to cope with this problem , it is not obvious why one should work with domain wall fermions, given the equivalence to overlap fermions, where the chiral symmetry breaking effects are much better under control and can even be eliminated completely.

Also for FP fermions, approximations have to be taken. First of all, the Dirac operator has to be restricted to finite extension, in our case to the hypercube. All couplings outside the hypercube are truncated. Second, only a limited set of gauge paths are considered for the couplings in the hypercube. It is therefore clear that only approximate chiral symmetry is present for the parametrized FP operator. The question then is whether the residual chiral symmetry breaking is negligibly small for the task under consideration. We took advantage of the freedom in the choice of Dirac operator for the overlap kernel (2.10) and used the overlap construction (2.11) with the parametrized FP operator as an input kernel to remove the residual chiral symmetry breaking of our parametrization for some applications.

The difficulty in simulations with overlap fermions arises from the inverse square root in Eq. (2.11), which is hard to calculate numerically and requires again an iterative procedure. Using tricks like the exact treatment of small eigenvalues of $`A^{}A`$, it is however possible to make the calculation of the inverse square root up to machine precision feasible within a few hundred iterative steps, rendering the chiral symmetry exact.

From the above considerations, one can quantify the computational demand of the different fermion formulations compared to the standard case of Wilson fermions: The simulation of domain wall fermions requires a factor of $`N_s`$ more computer time due to the additional fifth dimension. For overlap fermions, the factor is given by the number of iterative steps that has to be taken in order to compute the inverse square root. For typical applications, this factor is of order $`𝒪(200)`$. The computational cost of FP fermions depends on the parametrization that is chosen. We will come back to this point in Section 3.3.

Another way to get an approximately chiral symmetric fermion action is to optimize a parametrization of a general Dirac operator for solving the Ginsparg-Wilson relation . By truncating the expansion of the general operator in terms of the number of gauge paths and couplings and putting the truncated operator into the Ginsparg-Wilson relation, the free parameters can be fixed. The resulting operator approximates the Ginsparg-Wilson relation to a precision which depends on the truncation.

### 2.3 Perfect Actions from Renormalization Group Transformations

The perfect action approach to improving lattice actions followed to construct the Fixed-Point Dirac operator is inspired by the Renormalization Group flow of asymptotically-free theories . A Renormalization Group (RG) transformation  reduces the number of degrees of freedom by integrating some of them out in the path integral, taking into account their effect on the remaining variables exactly. This allows to get rid of short-distance fluctuations without changing the physical content of the theory. Consider a lattice action which contains all possible interactions. The RG transformation is defined by some blocking function which averages over the fields to produce a new action on a coarser lattice with fewer fields. The new action generally has different couplings from the original action, thus we can imagine the blocking step as a flow in the coupling parameter space. Repeated RG steps generate a trajectory in this space. In Figure 2.1, we show the RG trajectory for QCD with massless quarks. The fixed point has the property that the couplings are reproduced after a blocking step. For asymptotically free theories, the fixed point is on the surface of vanishing coupling $`g=0`$.

Starting on this surface, the RG trajectory flows to the fixed point. If one starts close to this surface at some small coupling $`g`$, or equivalently small lattice spacing $`a`$, the RG trajectory flows quickly towards the fixed point and then flows away from it. Let us assume we have an action with couplings lying on the RG trajectory at an arbitrarily small $`a`$, thus having arbitrarily small lattice artifacts. From there one can reach any point on the RG trajectory by making sufficiently many blocking steps, and all actions on the RG trajectory describe the same physics. The physical observables of the continuum quantum theory are thus identical to those of any lattice quantum theory on the RG trajectory, independently of the lattice spacing. Such lattice actions are called quantum perfect. The Fixed-Point action is an approximation to the RG trajectory for small couplings $`g`$ and is classically perfect, which means it completely describes the continuum classical theory without discretization errors .

Fixed-Point actions have many desirable features. By closely approximating the RG trajectory, they are expected to have largely reduced quantum lattice artifacts. They can be optimized for locality. The FP Dirac operator satisfies the Ginsparg-Wilson relation and so has good chiral behavior. The FP QCD action has well-defined topology and satisfies the index theorem on the lattice. The properties of FP actions have first been tested in models like the two-dimensional non-linear $`\sigma `$-model  and the $`CP^3`$-model . The approach has then been extended to $`SU(2)`$ and $`SU(3)`$ Yang-Mills theories and fermions in 2 and 4 dimensions , and first steps towards applications in partial differential equations have been taken . Recently, a new parametrization for $`SU(3)`$ Yang-Mills was constructed, showing reduced scaling violations in glueball masses and finite temperature measurements . We use this gluon action together with our parametrization of the FP Dirac operator for the simulations in Chapters 6 and 7. An extension of this FP gluon action to anisotropic lattices was constructed and tested in . For a pedagogical introduction to perfect actions, consult .

### 2.4 Free Fixed-Point Fermions

For the case of free fermions without mass, the Renormalization Group construction is relatively easy. Because the fermionic action is quadratic in the fermion fields, the Renormalization Group step for the fermion fields amounts to Gaussian integration, which can be done exactly. On the lattice, a RG transformation relates an action on a fine lattice with spacing $`a`$ to a different action on a coarser lattice with spacing $`2a`$. The blocking step thus connects the Dirac operators $`D_\mathrm{f}`$ on the fine and $`D_\mathrm{c}`$ on the coarse lattice by

$$D_\mathrm{c}^1=\frac{\mathrm{𝟏}}{\kappa }+\omega D_\mathrm{f}^1\omega ^{},$$
(2.12)
provided $`D_\mathrm{f}`$ has no zero modes, where $`\kappa `$ is an optimizable free parameter of the blocking and $`\omega `$ is the blocking function that relates the fine fields to the coarse fields. The Fixed-Point Dirac operator $`D^{\mathrm{FP}}`$ is reproduced under the blocking step,

$$(D^{\mathrm{FP}})^1=\frac{\mathrm{𝟏}}{\kappa }+\omega (D^{\mathrm{FP}})^1\omega ^{},$$
(2.13)
and depends on the choice of the blocking function $`\omega `$. For free fermions, that is in the absence of gauge fields, this equation can be solved analytically. The FP Dirac operator is local, and the rate of fall off for the couplings can be maximized by varying the parameter $`\kappa `$. However, $`D^{\mathrm{FP}}`$ contains infinitely many couplings. For practicality, the FP Dirac operator is approximated with an ultralocal operator, for which each point is only coupled to its neighbors on the hypercube. The effect of this truncation can be examined for the energy-momentum dispersion relation, which is equivalent to the continuum for the exact FP Dirac operator. The truncated operator deviates from the exact result, but shows still considerably smaller discretization errors than for example the Wilson operator .

In fact, the Renormalization Group procedure does not only generate a FP Dirac operator, but also a FP $`R`$ operator appearing on the right hand side of the Ginsparg-Wilson relation (2.4). Combining (2.4) with the blocking transformation (2.12) that connects the propagators $`D^1`$ on the coarse and fine lattices, we get the Renormalization Group relation

$$R_\mathrm{c}=\frac{\mathrm{𝟏}}{\kappa }+\omega R_\mathrm{f}\omega ^{},$$
(2.14)
for the $`R`$ operator, and at the fixed point, $`R_\mathrm{c}=R_\mathrm{f}=R^{\mathrm{FP}}`$. For free fermions, this equation can also be solved analytically. Choosing a symmetric overlapping block transformation $`\omega `$ with a scale factor 2 that averages over hypercubes , the exact $`R^{\mathrm{FP}}`$ has only hypercubic couplings, and therefore a truncation like for the Dirac operator is not required. With other methods to build a Dirac operator satisfying the Ginsparg-Wilson relation, for example the overlap construction (2.11), $`R`$ is unconstrained and typically $`R=1/2`$ is taken for simplicity.

## Chapter 3 The Parametrized Fixed-Point Dirac Operator

The concept of perfect actions is theoretically very attractive. The main obstacle for its application to a theory of general interest like QCD is to find a parametrization that is rich enough to capture all the beautiful properties of perfect actions, but is still feasible to use in numerical simulations. While the Fixed-Point Dirac operator is local, which means its couplings decrease exponentially, an ultralocal parametrization will introduce a truncation. This truncation does of course disturb the FP properties, and it is a non-trivial task to find a parametrized form whose properties do not deviate strongly from those of the FP operator. A parametrization of the FP Dirac operator will be more costly to simulate in terms of computer time than comparably simple Dirac operators like the Wilson or clover operators, as there are more couplings between lattice sites than just those to the nearest neighbor involved, and also the Clifford structure can be richer. However, one expects that the rewards compensate the additional cost of a more complicated action. In the case of the Dirac operator, a strong argument is certainly that chiral symmetry is preserved, in contrast to the standard actions. Additionally, the scaling violations are expected to be reduced for a parametrized FP operator. Smaller scaling violations allow to simulate at larger lattice spacings, while the results are still of unchanged quality. Since the computer time for a quenched lattice QCD simulation increases like $`a^6`$$`a^7`$, being able to simulate at lattice spacing $`2a`$ instead of $`a`$ brings a factor of $`𝒪(100)`$ in computational savings. Even more pronounced is the situation in the unquenched case, where the cost increases like $`a^8`$$`a^{10}`$ with the lattice spacing, so that the expected gain can even be of $`𝒪(1000)`$.

In this chapter we first derive the structure of general lattice Dirac operators respecting the appropriate discrete symmetries, and show how complicated operators with a large number of couplings and gauge paths can be calculated efficiently. This has been examined in detail in , and we present here the key concepts of the paper. Then we explain our procedure of fitting the parameters of the general Dirac operator to the FP operator, using the Renormalization Group recursion relations. Finally we show as a test for the properties of the resulting parametrization the eigenvalue spectrum, which measures how well the Ginsparg-Wilson relation is obeyed.

### 3.1 General Lattice Dirac Operators

The starting point for constructing any lattice Dirac operator is the question what general structure is allowed if the basic lattice symmetries have to be respected. These symmetries are discrete translation invariance, gauge invariance, $`\gamma _5`$-hermiticity, charge conjugation, permutation and reflection symmetry. Let us summarize the constraints which the discrete symmetries impose on any lattice Dirac operator $`D(n,n^{};U)`$.

#### 3.1.1 Discrete Symmetries and Gauge Invariance

From translation invariance follows that $`D(n,n+r)`$ depends on the lattice variable $`n`$ only through the $`n`$-dependence of the gauge fields. In particular, the coefficients in front of the different gauge paths which enter the Dirac operator do not depend on $`n`$. The hermiticity properties of the lattice operator are required to be the same as in the continuum,

$$D(n,n^{})=\gamma _5D(n^{},n)^{}\gamma _5,$$
(3.1)
where the hermitean conjugation acts in color and Dirac space. Permutations of the coordinate axes are defined in a straightforward way, as just the Dirac indices appearing in $`D`$ can be permuted.<sup>1</sup><sup>1</sup>1Note that cubic rotations on the lattice can be replaced by reflections and permutations of the coordinate axes. In Appendix C.1 and C.2, we derive the conditions from charge conjugation,

$$D(U_\mu )=C^1D(U_\mu ^{})^TC,$$
(3.2)
where $`C\gamma _\mu ^TC^1=\gamma _\mu `$, and reflections of a coordinate axis $`\eta `$,

$$D(n,n^{};U_\mu (n))=P_\eta ^1D(\stackrel{~}{n},\stackrel{~}{n}^{};U_\mu ^{𝒫_\eta }(\stackrel{~}{n}))P_\eta ,$$
(3.3)
where in our convention $`P_\eta =\gamma _\eta \gamma _5`$, and $`\stackrel{~}{n}`$ is the reflected space-time variable defined in (C.3).

It remains to ensure gauge invariance, which is the most crucial ingredient. If the Dirac operator transforms under the gauge transformation $`G(n)SU(N)`$ in (1.23)–(1.25) like

$$D(n,n^{};U)G(n)^{}D(n,n^{};U^g)G(n^{}),$$
(3.4)
where $`U^g`$ is the gauge transformed background field, the fermion action stays gauge invariant. This can be achieved by connecting the lattice sites $`n`$ and $`n^{}`$ along an arbitrary path

$$l=[l_1,l_2,\mathrm{},l_k],$$
(3.5)
of length $`k`$, where $`|l_i|=1,\mathrm{},4`$ is the direction of the path at step $`i`$, with a parallel transporter

$$U^l(n)=U_{l_1}(n)U_{l_2}(n+\widehat{l}_1)\mathrm{}U_{l_k}(n+\widehat{l}_1+\mathrm{}+\widehat{l}_{k1}),$$
(3.6)
made from products of gauge links. For every offset $`r=n^{}n`$ appearing in the Dirac operator, one or several paths $`l`$ connecting both end points have to be chosen to make $`D`$ gauge covariant.

#### 3.1.2 General Construction

The symmetry conditions (3.1)–(3.3) prescribe in which combinations the different permutations and reflections of the gauge paths (3.6) have to enter the Dirac operator.

To be more specific, a general gauge covariant lattice operator with color, space and Dirac indices can be written as

$$D(n,n^{})=\underset{i=1}{\overset{16}{}}\mathrm{\Gamma }_i\underset{l}{}c(\mathrm{\Gamma }_i,l)U^l(n),$$
(3.7)
where $`\mathrm{\Gamma }_i`$ are elements of the Clifford algebra and $`c(\mathrm{\Gamma }_i,l)`$ is the coupling for the given path $`U^l(n)`$ and Clifford algebra element. The basic symmetries of the Dirac operator impose the following restrictions on Eq. (3.7):

Translation invariance requires that the couplings $`c(\mathrm{\Gamma }_i,l)`$ are constants in space-time or gauge invariant functions of gauge fields, respecting locality and invariance under the symmetry transformations. Charge conjugation and $`\gamma _5`$-hermiticity together imply that the couplings $`c(\mathrm{\Gamma }_i,l)`$ are real for our choice of the Clifford algebra basis. From hermiticity it follows that the path $`l`$ and the opposite path $`\overline{l}=[l_k,\mathrm{},l_1]`$, or equivalently $`U^l(n)`$ and $`U^l(n)^{}`$, should enter in the combination

$$\mathrm{\Gamma }\left(U^l(n)+ϵ_\mathrm{\Gamma }U^l(n)^{}\right),$$
(3.8)
where the sign $`ϵ_\mathrm{\Gamma }`$ is defined by $`\gamma _5\mathrm{\Gamma }^{}\gamma _5=ϵ_\mathrm{\Gamma }\mathrm{\Gamma }`$. Permutations and reflections of the coordinate axes (hypercubic rotations) imply that for a given reference path $`l_0`$ a whole class of paths belongs to the Dirac operator. These paths are related to $`l_0`$ by all the $`16\times 24=384`$ reflections and permutations of the coordinate axes. Under such a symmetry transformation $`\alpha =1,\mathrm{},384`$ the Clifford algebra element $`\mathrm{\Gamma }_0\{\mathrm{𝟏},\gamma _\mu ,\sigma _{\mu \nu },\gamma _5,\gamma _\mu \gamma _5\}`$ associated with $`l_0`$ generally transforms<sup>2</sup><sup>2</sup>2The transformed Clifford algebra element $`\mathrm{\Gamma }^{(\alpha )}`$ is of the same type (S, V, P, T, A) as $`\mathrm{\Gamma }_0`$.  into $`\mathrm{\Gamma }^{(\alpha )}`$ and the parallel transporter $`U^l(n)`$ transforms to $`U^{l^{(\alpha )}}(n)`$. Furthermore the sign of the couplings may change, whereas their absolute value remains unchanged.

A Dirac operator satisfying all the basic symmetries can be written as

$$D(n,n^{})=\underset{\mathrm{\Gamma }_0,l_0}{}c(\mathrm{\Gamma }_0,l_0)\underset{\alpha }{}O^\alpha (n),$$
(3.9)
where the sum runs over a set of reference paths defined by $`\mathrm{\Gamma }_0`$ and $`l_0`$ as well as over all the symmetry transformations $`\alpha `$ defined by the group of permutations and reflections of the coordinate axes. The operators $`O^\alpha `$ are defined by

$$O^\alpha (n)=\mathrm{\Gamma }^{(\alpha )}\left(U^{l^{(\alpha )}}(n)+ϵ_\mathrm{\Gamma }U^{l^{(\alpha )}}(n)^{}\right).$$
(3.10)
To make the above construction clear, consider the example of the Wilson Dirac operator (1.32), which contains only the elements $`\mathrm{𝟏}`$ and $`\gamma _\mu `$ of the Clifford algebra and extends to nearest neighbors as sketched in Fig. 3.1. For the scalar element $`\mathrm{\Gamma }_0=\mathrm{𝟏}`$, the reference paths are $`l_0=[]`$, which amounts to the contact term, and the nearest neighbor coupling $`l_0=[1]`$, while for the vector element $`\mathrm{\Gamma }_0=\gamma _1`$ there is only the nearest neighbor $`l_0=[1]`$. The sum over these reference paths and Clifford algebra elements in Eq. (3.9) gives then the Wilson Dirac operator (1.32), when the coefficients

$`c(\mathrm{𝟏},[])`$ $`=`$ $`m+4r,`$
$`c(\mathrm{𝟏},[1])`$ $`=`$ $`r/2,`$
$`c(\mathrm{\Gamma }_1,[1])`$ $`=`$ $`1/2,`$ (3.11)
are taken. It is quite natural to include the full Clifford algebra to parametrize the FP Dirac operator. The scalar and vector elements are already present in the continuum. The tensor element $`\sigma _{\mu \nu }`$ appears in the $`𝒪(a)`$-improved Sheikoleslami-Wohlert clover operator and thus already in lowest order of the Symanzik improvement program. Without the pseudoscalar element $`\gamma _5`$, the Ginsparg-Wilson relation (2.4) could not be fulfilled, so it is crucial for the chiral properties of the Dirac operator. Moreover, the topological charge is proportional to $`\mathrm{Tr}(\gamma _5D)`$, which would be zero if $`D`$ does not contain the pseudoscalar element. Finally, as the Renormalization Group procedure which leads to the Fixed-Point operator generates all the elements of the Clifford algebra, also the axial vector element $`\gamma _\mu \gamma _5`$ should be included.

### 3.2 Efficient Implementation of General Dirac Operators

At first sight one might think that it is not feasible to calculate such a general Dirac operator with many different couplings, where every coupling might contain as many as $`768`$ paths. But one has to keep in mind that in applications like hadron spectroscopy, the calculation of propagators for small quark masses needs several hundreds of conjugate gradient steps and therefore one can afford to spend some time to precalculate and store the whole Dirac operator before starting to calculate the propagator. The preparation of the Dirac operator then only needs a small fraction of the overall time for a calculation. On top of this there are two reasons why the calculation of the gauge paths for general operators can be done in a very efficient way: First, there are a lot of paths which are invariant under certain subgroups of the reflections and permutations. This reduces the number of terms significantly and in some case even leads to a cancellation of certain terms because they have opposite signs. A less trivial fact is that the sum of paths for many couplings can be factorized in an efficient way, which means that large sums of many paths can be written as a product of smaller sums of fewer paths.

As an example, consider the nearest neighbor coupling with $`\mathrm{\Gamma }_0=\gamma _5`$ and reference path $`l_0=[2,1,2,3,4,3,4]`$. All the paths of this coupling can be written in the following compact way,

$$\gamma _5\underset{\mu \nu \rho \sigma }{}ϵ_{\mu \nu \rho \sigma }\left(S_{\mu \nu }P_{\rho \sigma }+P_{\rho \sigma }S_{\mu \nu }+\text{h.c.}\right),$$
(3.12)
where the color matrices $`S_{\mu \nu }`$ and $`P_{\rho \sigma }`$ are certain combinations of staples or plaquettes, respectively. When all the plaquettes and staples and the most frequent combinations like $`P_{\rho \sigma }`$ are precalculated, operators like the one in Eq. (3.12) can be calculated very quickly. As an illustration of this we consider the parametrization of the Fixed-Point Dirac operator used in our spectroscopy simulations, which has $`241`$ couplings in total, at least one per offset on the hypercube and per type of Clifford algebra element. Building this operator on a workstation takes only $`𝒪(20)`$ times as long as the multiplication of the operator with a vector and therefore it is a very small fraction of the time used to perform a calculation of one propagator. On a supercomputer, this relation gets worse, because the construction of the gauge paths is not as well-suited for optimization or vectorization like the matrix-vector multiplication used in the inversion of the Dirac operator. The measurements from our spectroscopy runs in Appendix B.2 show that on the Hitachi SR8000, the time to construct $`D`$ is on the order of 10–20% of the calculation time for 12 quark propagators that are needed to construct hadron correlators, which is not negligible anymore. However, when using the parametrized operator in an overlap construction, the build-up time is negligible again, because in that case the matrix-vector multiplication gets more expensive by a factor which is given by the expansion order of the inverse square root in the overlap operator.

Another question is how fast manipulations with the Dirac operator can be executed after it has been constructed. The basic operation required to calculate propagators or eigenvalues is the multiplication of the Dirac operator on a vector. For a Dirac operator with 81 hypercubic fermion offsets which contains the complete Clifford algebra, the matrix-vector multiplication requires $`9\times 4=36`$ times more $`c`$-number multiplications than for the Wilson Dirac operator, which only connects 9 offsets and whose Dirac structure can be treated trivially. The actual performance on a specific computer however depends a lot on the architectural and implementational details, and as we do not have an optimized code for the Wilson operator, we can not confirm this number from performance measurements. There are however additional arguments why in actual simulations with the parametrized FP operator the factor is considerably smaller than 36: The most striking one is that for small quark masses, the Wilson operator runs into problems with exceptional configurations, where the inversion converges very slowly or not at all. We did not encounter such problems for the parametrized FP Dirac operator at the comparably small quark masses covered in our spectroscopy simulations.

### 3.3 Parametrization of the Fixed-Point Dirac Operator in QCD

Finding a good parametrization of a FP action is a non-trivial task. In the last years, some attempts were taken to parametrize the FP Dirac operator , but these were limited to moderate generalizations of the Wilson operator, including only a few additional couplings and part of the Clifford structure. These parametrizations were rather thought to be taken as a starting point for the overlap construction than to be used on their own. We took a different approach: Our goal was to get a parametrization as close as possible to the massless FP Dirac operator, that can be directly used in QCD simulations near the chiral limit. We therefore use a general Dirac operator as defined in Eq. (3.9) with all the couplings of the hypercube (see Fig. 3.1) and all elements of the Clifford algebra. In order to get an operator which is close to the fixed point for a range of values of gauge couplings $`\beta `$, the coefficients $`c(\mathrm{\Gamma }_0,l_0)`$ in Eq. (3.9) were chosen not to be constants, but gauge invariant polynomials in local fluctuations of the gauge fields. Furthermore we apply a RG-inspired smearing procedure  to the gauge configurations and project them back to SU(3), that is we are using so-called fat links. For the gluon sector, we use the recent parametrization of the FP gauge action , which also makes use of fat gauge links.

Like in the case of free fermions, the QCD FP action is quadratic in the fermion fields and the Renormalization Group step can again be done analytically. The QCD FP equation is the generalization of the free case given by including gauge fields. In case $`D_f`$ has zero modes, it is most conveniently written as

$$D_c(V)=\kappa \mathbf{\hspace{0.17em}1}\kappa ^2b^2\omega (U)[D_f(U)+\kappa b^2\omega ^{}(U)\omega (U)]^1\omega ^{}(U),$$
(3.13)
where $`\kappa `$ is an optimizable free parameter and $`b`$ is a scale factor of the blocking transformation, while $`U`$ and $`V`$ are the gauge fields on the fine and coarse lattice, respectively. They are related through the FP equation of the pure $`SU(3)`$ gauge theory

$$S_G^{\mathrm{FP}}(V)=\underset{\{U\}}{\mathrm{min}}\left(S_G^{\mathrm{FP}}(U)+T(U,V)\right),$$
(3.14)
where $`S_G^{\mathrm{FP}}`$ is the FP action of the pure $`SU(3)`$ gauge theory and $`T(U,V)`$ is the kernel of the blocking transformation. An important fact for the parametrization of the FP Dirac operator is that Eq. (3.13) can also be given in terms of the propagators,

$$D_c^1(V)=\frac{\mathrm{𝟏}}{\kappa }+b^2\omega (U)D_f^1(U)\omega ^{}(U),$$
(3.15)
as long as $`D_f`$ has no zero modes. In contrast to Eq. (3.13) the equation for the propagator gives much more weight to the small physical modes of the Dirac operator and can therefore be used to improve the properties of the small modes of the parametrized FP Dirac operator.

#### 3.3.1 Fitting the Parameters

The parametrization is an iterative procedure, as sketched in Fig. 3.3. We start at a large value of $`\beta =100`$, generate a few thermal coarse gauge configurations $`V^{(i)}`$, $`i=1,\mathrm{},n_{\mathrm{conf}}`$, with the FP gauge action and determine the corresponding fine configurations $`U^{(i)}`$ through the minimization in Eq. (3.14). As a Dirac operator on the right-hand side of the FP equations (3.13) and (3.15), we choose the free FP Dirac operator. This is justified from the fact that the minimizing configurations $`U^{(i)}`$ have very small fluctuations at such a large value of $`\beta `$ and are therefore very close to the free field case. For each configuration, we take two sets of vectors $`\{x_k^{(i)}\}`$ and $`\{y_k^{(i)}\}`$, $`k=1,\mathrm{},n_{\mathrm{vec}}`$, of dimension $`12N_{\mathrm{site}}`$ on the coarse lattice. The $`x_k^{(i)}`$ can either be random vectors or approximate small eigenmodes of $`D_c`$, and the $`y_k^{(i)}`$ are just random vectors. From the right-hand sides of Eqs. (3.13) and (3.15), the vectors

$`\xi _k^{(i)}`$ $`=`$ $`D_cx_k^{(i)},`$ (3.16)
$`\eta _k^{(i)}`$ $`=`$ $`D_c^1y_k^{(i)},`$ (3.17)
are calculated. The couplings of the parametrized Dirac operator $`D_{\mathrm{par}}`$ are then determined by minimizing

$$\chi ^2=\underset{i,k}{}D_{\mathrm{par}}(V^{(i)})x_k^{(i)}\xi _k^{(i)}^2+\lambda \underset{i,k}{}D_{\mathrm{par}}^1(V^{(i)})y_k^{(i)}\eta _k^{(i)}^2,$$
(3.18)
where the sum runs over the different configurations and the set of vectors per configuration and $`\lambda `$ is a weighting factor that has to be appropriately chosen. During this procedure we keep leading terms in the naive continuum limit fixed such that the tree level mass is zero, the $`O(a)`$ Symanzik condition is fulfilled, the dispersion relation starts with slope $`1`$ and the normalization of the topological charge is correct  <sup>3</sup><sup>3</sup>3In Eq. (31) of , the sign for the condition on the Dirac operator which fixes the topological charge is incorrect. This influenced the first steps in our parametrization procedure. In the last steps this condition was not used. The forcing of the topological charge to the wrong sign in the earlier part of the iterative procedure does not affect the end result significantly, because the overlap reparametrization in the very last step straightened out this error. In fact, even during the phase when the wrong condition was applied, the linear terms in the fluctuation polynomials, which are not affected by this condition, helped keeping the FP properties in the pseudoscalar sector present (see also ). . Furthermore we fix the free field limit such that the truncated free FP operator is recovered on the trivial gauge configuration $`U=1`$. The minimization of the $`\chi ^2`$-function (3.18) yields a parametrized FP Dirac operator $`D_{\mathrm{par}}`$ which has good chiral properties over a larger range of gauge couplings than the initial truncated free FP Dirac operator. The fitted parametrized operator $`D_{\mathrm{par}}`$ is now used on fine configurations $`U^{}`$, determined via minimization from coarse configurations $`V^{}`$ generated thermally at a smaller value of $`\beta `$. Minimizing the $`\chi ^2`$-function (3.18) again gives $`D_{\mathrm{par}}(V^{})`$, which performs well on an even larger range of gauge couplings. The whole procedure is repeated until we reach $`\beta 3.0`$, corresponding to $`a0.16`$ fm. In the final phase of this iterative procedure, some of the naive continuum limit constraints on the parameters of the Dirac operator are released. Furthermore, in the last blocking step a low-order overlap expansion is applied to the parametrization which is put into the right hand side in order to reduce the remaining fluctuations of the small eigenvalues even further .

Let us make a few comments on this procedure of iteratively finding the parameters for the FP Dirac operator at lattice spacings typically used in simulations. First, the use of vectors for the calculation of a $`\chi ^2`$-function for $`D_{\mathrm{par}}`$ is mandatory because the definition of $`D_c`$ requires a matrix inversion which we can only afford for a limited number of vectors. Even then, on our workstations we were restricted to lattices of maximum size $`5^4`$ for the coarse and $`10^4`$ for the fine configurations. We usually worked with sets of $`5n_{\mathrm{conf}}15`$ different configurations and $`n_{\mathrm{vec}}=5`$ vectors per configuration. The use of large enough lattices is important because when the lattice is too small, there are essentially no small eigenvalues of $`D`$, which are however crucial for fitting the chiral properties of the operator and become particularly important when going to smaller values of $`\beta `$. If these small eigenvalues are missing, the fit only captures the FP properties in the region of large eigenvalues of $`D`$ well, and the resulting parametrization then suffers for example from large additive mass renormalization. For the same reason it is important to include the RG relation (3.15) for the propagator in the fit. At the largest value of $`\beta `$ however, the procedure was not so sensitive for the presence of small eigenvalues, and we worked on smaller lattices of size $`3^4`$ and $`6^4`$, respectively.

Second, we checked that the whole procedure does not strongly depend on the choice of input Dirac operator for the right-hand side of the RG relation at the largest value of $`\beta `$. When instead of the free-field FP operator the Wilson operator is used in the first step, the set of parameters after a few iterative steps agrees well with the one resulting from the free FP operator as an input. This observation confirms that at $`\beta =100`$, the minimized gauge configurations $`U`$ have such small fluctuations that essentially any Dirac operator can be taken as an input for the blocking without changing the result. The fluctuations can be quantified by measuring the average value of the plaquette

$$u=\frac{1}{6N_{\mathrm{site}}}\underset{n}{}\underset{\mu <\nu }{}\mathrm{Re}\mathrm{Tr}U_{\mu \nu }(n),$$
(3.19)
normalized to $`u=3`$ in the free field limit, which is listed for the fine and coarse gauge configurations at different values of the coupling in Table 3.1. The measurements show that the minimizing fine configurations $`U`$ have very small fluctuations already at $`\beta =3`$, and are very close to the free field limit at the largest $`\beta `$. We have to remark that the average value of the plaquette hides the somewhat different distribution for thermal and minimizing configurations, which one also has to take into account.

The last set of coarse and fine gauge configurations was chosen to cover a range of values of the gauge coupling in the interval $`\beta [2.7,3.4]`$, corresponding to lattice spacings in the range $`0.1`$ fm $`<a<0.22`$ fm. These values were chosen such that the final parametrization can be used at somewhat larger lattice spacings than typically used for simulations of unimproved actions. That the fit captures the FP properties equally well on both ends of this range can be seen in the plot on the left of Fig. 3.4, where we show the correlation between the values in the propagator part of the $`\chi ^2`$-function (3.18) from two gauge configurations at $`\beta =2.9`$ and $`\beta =3.4`$, respectively. Configurations with $`\beta <2.9`$ were only used in the fitting of the FP relation (3.13) for $`D`$, that is for the first summand in the $`\chi ^2`$-function.

We also tested whether adding another intermediate step with fine and coarse configurations at $`\beta =5`$ changes the outcome of the parametrization significantly, but this was not the case. Furthermore we performed a simple check whether it makes sense to explicitly minimize the breaking of the Ginsparg-Wilson relation (2.5) by including it in the fit. This is a redundant constraint, as the FP relation itself makes sure that the operator fulfills the Ginsparg-Wilson relation. The right-hand side of (2.5) should be zero outside the hypercube, therefore we measured the Ginsparg-Wilson breaking by calculating

$$\chi _{\mathrm{GW}}^2=\underset{x}{}x\left[\{D_{\mathrm{par}}^1,\gamma _5\}2\gamma _5\right]y,$$
(3.20)
where the vector $`y`$ only has non-zero entries at the lattice origin (0,0,0,0) and the sum runs over vectors $`x`$ with non-zero entries at one single lattice site outside the hypercube around the origin. This is computationally quite expensive due to the inversion of the Dirac operator which is needed and therefore increases the time for an iteration step in the minimization of $`\chi ^2`$ considerably. The plot on the right-hand side of Fig. 3.4 shows that $`\chi _{\mathrm{GW}}^2`$ from the Ginsparg-Wilson breaking seems to correlate highly with the $`\chi ^2`$ from the propagator FP relation, so we did not pursue this path further.

In order to parametrize the operator $`R`$ in Eq. (2.14) we proceed in a similar way as for the Dirac operator. We also use a general operator with fat links and fluctuation polynomials. The parametrization of $`R`$ is however simpler as it is trivial in Dirac space and therefore contains a smaller number of operators. In contrast to the equation for the blocking transformation of $`D^{\mathrm{FP}}`$ the corresponding equation for $`R`$ (2.14) contains no inversion and therefore the $`\chi ^2`$-function which we minimize in order to find the optimal parametrization of $`R_{\mathrm{par}}`$ can be defined as

$$\chi ^2=R_{\mathrm{par}}R_\mathrm{C}^2,$$
(3.21)
where the norm here is the matrix norm $`A^2=_{i,j}a_{ij}^2`$.

A final remark on the hypercubic truncation: The Fixed-Point $`R`$ operator is hypercubic, hence no truncation is needed. For the free FP Dirac operator, which is known analytically, the couplings outside the hypercube are very small, and thus the truncation does not distort the FP properties too much. As our results will show, also in the interacting case the FP operator can be well described by a hypercubic parametrization, although some couplings which lie outside the hypercube and are therefore left out tend to grow larger than for the free case. We did however not consider extending the parametrization beyond the hypercube, as then the computational demand both for the construction of $`D`$ and for the matrix-vector multiplication would grow very rapidly.

We have to mention another approximation that was taken in this work when parametrizing the FP Dirac operator: We only constructed the FP operator for zero quark mass, whereas in principle for every mass value a different parametrization would be necessary. At larger masses, our parametrization is therefore expected to deviate from the fixed point.

### 3.4 Eigenvalue Spectrum

An easily accessible observable that quantifies the quality of a given Dirac operator in terms of fulfilling the Ginsparg-Wilson relation is the eigenvalue spectrum. Consider the case of a non-trivial $`R`$, and define a rescaled Dirac operator $`\stackrel{~}{D}=\sqrt{2R}D\sqrt{2R}`$. Setting the lattice spacing $`a=1`$, the Ginsparg-Wilson relation (2.4) can be written as

$$\stackrel{~}{D}+\stackrel{~}{D}^{}=\stackrel{~}{D}^{}\stackrel{~}{D},$$
(3.22)
implying that the eigenvalues of $`\stackrel{~}{D}`$ lie on a circle of radius 1 and center $`(1,0)`$. Using $`R^{\mathrm{FP}}`$ and the hypercubic approximation of $`D^{\mathrm{FP}}`$, we compare the eigenvalue spectrum of $`\stackrel{~}{D}`$ on a gauge configuration at $`\beta =3.0`$ with the spectrum of the Wilson operator in Figs. 3.5 and 3.6. While the eigenvalues of the Wilson operator spread over a large region in the complex plane, they lie almost exactly on the circle for the parametrized FP Dirac operator, indicating that the hypercubic truncation and the restriction to a finite set of gauge paths has only slightly affected the Ginsparg-Wilson property. Another very striking observation is that the additive mass renormalization, which is given by the point where the smallest eigenvalues intersect with the real axis, is of $`𝒪(1)`$ and thus very large for the Wilson operator, while for the parametrized FP Dirac operator no additive mass renormalization is seen in the eigenvalue spectrum at this value of the gauge coupling. It has to be stressed that this property is by no means enforced in the parametrization procedure by shifting or constraining the parameters. It originates only from the fact that the parametrization describes the FP properties well.

Checking the eigenvalue spectrum on a set of different gauge configurations at given $`\beta `$, the fluctuation of the smallest or near-zero eigenvalues can be measured, which is an important quantity for simulating light quarks. Large fluctuations lead to the appearance of exceptional configurations, making it impossible to invert the Dirac operator at small quark mass. While the fluctuations are very large for the Wilson operator, they are on the order of $`10^2`$ for the parametrized FP operator at $`\beta =3.0`$. Due to this crucial property it is possible to perform lattice calculations at pseudoscalar to vector meson mass ratios as small as $`m_{\mathrm{PS}}/m_\mathrm{V}0.27`$ for lattice spacings of $`a=0.16`$, as we show in the results of our spectroscopy simulations in Chapter 7.

## Chapter 4 The Overlap-Improved Fixed-Point Dirac Operator

The overlap construction proposed by Neuberger  allows to formulate a lattice Dirac operator with exact chiral symmetry. In the last few years, a lot of activity has been going on testing and applying overlap fermions, but while extensive calculations of quenched QCD spectroscopy in the chiral limit have been done with the approximately chiral domain-wall fermions , studies of lattice spectroscopy at small quark mass with overlap fermions are only very recently emerging .

Most groups working with overlap fermions use the Wilson Dirac operator as a starting point for the overlap, which might not be an optimal choice. Wilson fermions show large scaling violations, and in the overlap construction only the $`𝒪(a)`$ artifacts are removed. The $`𝒪(a^2)`$ effects however will remain present in the resulting Dirac operator and might even get enhanced. Furthermore, the ultralocality of the Wilson operator is lost in the overlap, and while the Wilson overlap operator is still local, the couplings decrease only with a rather small exponent . Due to the strong chiral symmetry breaking of the Wilson operator, also the numerical calculation of $`1/\sqrt{A^{}A}`$ is not easy because the condition number of the matrix $`A^{}A`$ is large. Better kernels for the overlap have been considered only by a few groups  up to now. The clover action as a kernel seems to perform even worse than the Wilson operator .

The FP Dirac operator, being already an exact Ginsparg-Wilson operator, remains unchanged under the overlap construction. Because our parametrization of the FP Dirac operator is approximating the Ginsparg-Wilson relation very well, it will be changed only to a small extent by the overlap. The good scaling properties of the FP operator are then expected to be preserved in the end result. Also only relatively few iterative steps in the overlap expansion are necessary to ensure exact chirality. This property has been used to calculate the finite-volume scaling of the chiral condensate , which is a task where chiral symmetry is required to be present to a very high level. For the spectroscopy calculations in this work, we follow a slightly different strategy. We use an expansion to very low order in the inverse square root of the overlap, thus removing the already small chiral symmetry breaking effects introduced by the parametrization to a large extent, but not to machine precision. This strategy is from a computational point of view more than competitive to standard Wilson overlap simulations, and we expect to have the additional advantages of better localization of the resulting operator and improved scaling.

### 4.1 Implementation of the Overlap

Since we have parametrized the FP Dirac operator with a non-trivial $`R`$ operator in the Ginsparg-Wilson relation (2.4), the massive overlap-improved FP Dirac operator $`D(m)`$ which we use in our simulations has to be specified with the corresponding covariant scalar density  by

$$D(m)\left(1\frac{m}{2}\right)D+\frac{m}{2R},$$
(4.1)
where the massless overlap Dirac operator $`D`$ for non-trivial Ginsparg-Wilson $`R`$ is

$$D\frac{1}{\sqrt{2R}}\left(1\frac{A}{\sqrt{A^{}A}}\right)\frac{1}{\sqrt{2R}}.$$
(4.2)
In the kernel $`A`$ of the overlap expansion we use the parametrized FP Dirac operator $`D^{\mathrm{FP}}`$ from Chapter 3:

$$A1\sqrt{2R}D^{\mathrm{FP}}\sqrt{2R}.$$
(4.3)
The inverse square root in Eq. (4.2) is approximated using a Legendre expansion. The convergence of this expansion can be tremendously improved when treating the smallest eigenmodes of $`A^{}A`$ exactly, as then the condition number of the matrix becomes much smaller. We show in Fig. 4.1 the smallest 100 eigenvalues on 80 gauge configurations from our spectroscopy simulation. There are some isolated eigenvalues very close to zero and a rapidly increasing density of eigenvalues closer to 1. The threshold where the eigenvalues get dense decreases with increasing lattice spacing and volume. We showed in  that this threshold lies much higher for the FP kernel (4.3) than for the standard kernel with the Wilson operator, implying that the overlap expansion for the FP kernel needs orders of magnitudes fewer iterations. In the simulations in Chapter 7, we treat on all gauge configurations the smallest 100 eigenmodes exactly. The eigenvalue of the largest exactly treated mode ($`i=100`$ in Fig. 4.1) lies in this case between $`0.33<\lambda <0.42`$ for all considered gauge configurations. The largest eigenvalue of $`A^{}A`$ is typically close to $`\lambda =1.7`$, and therefore the subspace where the iterative solver works is very well-behaved.

Due to the good chiral properties of the starting FP operator, we can restrict the Legendre expansion to rather low order $`n10`$. In Fig. 4.2 we plot the residual breaking of the Ginsparg-Wilson relation determined by $`B=|(\stackrel{~}{D}+\stackrel{~}{D}^{}\stackrel{~}{D}^{}\stackrel{~}{D})v|^2`$ where $`v`$ is a random vector normalized to 1 and $`\stackrel{~}{D}=\sqrt{2R}D\sqrt{2R}`$ is the rescaled Dirac operator, as a function of the overlap order. We see that the chiral symmetry breaking decreases very quickly. While $`n10`$ has been used for measurements of the chiral condensate , for the spectroscopy calculations in this work we take $`n=3`$, which gives a Dirac operator with improved, but not exactly chiral properties relative to the parametrized FP Dirac operator. We will therefore call this the overlap-improved FP Dirac operator. As an order $`n`$ Legendre expansion of the inverse square root of $`A^{}A`$ requires $`2n`$ multiplications of the Dirac operator on a vector, the computational requirements of the overlap-improved Dirac operator from calculating $`A/\sqrt{A^{}A}`$ are by a factor of 7 larger than for the parametrized FP Dirac operator.<sup>1</sup><sup>1</sup>1The time for multiplications with $`R`$ is neglected here.

In a multi-mass inverter, the massive overlap Dirac operator can not be used in the form (4.1) due to the presence of the $`R`$ operator. Instead we actually invert the operator $`D2R`$ by writing

$$D(m)^1=\frac{1}{1m/2}2R\left(D2R+\frac{m}{1m/2}\right)^1,$$
(4.4)
where the term in brackets defines a shifted linear system that can be solved by the inverter. As one can see, every matrix-vector multiplication in the algorithm requires a multiplication of both the $`D`$ and the $`R`$ operator on the vector. Because the $`R`$ operator is trivial in Dirac space, this leads however only to a comparably small computational overhead.

### 4.2 Locality of Couplings

Every reasonable lattice Dirac operator has to be local, with exponentially decreasing couplings. If the Dirac operator is restricted to a finite number of lattice sites, like the Wilson or the parametrized FP Dirac operator, it is called ultralocal. The overlap construction takes an ultralocal operator as an input, but the result is no longer ultralocal. For the Wilson kernel, locality of the resulting overlap operator has been shown, but the exponential decay rate of the couplings is quite small. In Fig. 4.3 we compare the locality, measured by

$$f(r)=\underset{𝑦}{\mathrm{max}}\{|\mathrm{D}\mathrm{v}|;\mathrm{y}\mathrm{x}=\mathrm{r}\},$$
(4.5)
where $`v`$ is vector with a point source at $`x`$, for the overlap Dirac operator with the Wilson and the FP kernel. Together with other tests for the locality of the Wilson overlap operator , it follows that the locality is significantly improved when the FP kernel is used.

### 4.3 Locality of Instanton Zero Modes

The properties of approximately or exactly chiral Dirac operators can be tested on smooth instanton gauge configurations. Phenomenological models suggest that instantons might be responsible for dynamical chiral symmetry breaking . Single instantons, which are gauge configurations with topological charge $`Q=1`$, produce a zero mode of the Dirac operator through the index theorem (2.9). A pair of an instanton and an anti-instanton produces two complex-conjugate near-zero eigenmodes. The chiral condensate, which is the order parameter for chiral symmetry breaking, is related to the density of eigenvalues of the Dirac operator near the origin through the Banks-Casher relation. In the infinite volume limit, the density of exact zero modes becomes negligible compared to the density of near-zero modes. Only with the contribution of near-zero modes, the eigenvalue density at zero therefore does not vanish at infinite volume. Initiated by , many groups have recently studied whether the local chirality properties of these near-zero modes are consistent with the instanton model of chiral symmetry breaking .

In this section, we analyze eigenvalues and eigenvectors of the overlap-improved FP Dirac operator on gauge configurations describing a discretized exact instanton, as done in  for the chirally improved Dirac operator. For the above mentioned studies on local chirality, it is helpful to know how well a given Dirac operator reproduces the continuum zero mode of an instanton. We show that the Wilson overlap Dirac operator is in this respect not optimal, which might be due to its comparably bad locality properties. The gauge configurations are constructed from discretized $`SU(2)`$-instantons trivially embedded as $`2\times 2`$-submatrices in $`SU(3)`$ . We work on lattices of size $`12^4`$ and apply antiperiodic boundary conditions in time for the fermions. Eigenvalues and eigenvectors are calculated with the implicitly restarted Arnoldi algorithm .

First, we compare the flow of the zero eigenvalue with instanton size for various Dirac operators in Fig. 4.4. For an exactly chiral operator (and also for the overlap-improved FP operator on the scale of this plot), the eigenvalue is exactly zero. Due to the residual chiral symmetry breaking, it can move away from zero for approximately chiral Dirac operators. The eigenvalue is however restricted to the real axis, as long as the Dirac operator respects $`\gamma _5`$-hermiticity. For varying instanton radius $`\rho `$, the position $`x`$ of the zero eigenvalue on the real axis can be monitored, providing a measure for the chiral properties of the Dirac operator.

For the Wilson operator, the eigenvalue quickly flows away from zero towards the center of the circle with decreasing instanton radius $`\rho `$. Calculating only the few smallest eigenmodes, we lost track of it already at $`a\rho <2`$. The parametrized FP Dirac operator shows a much better behavior, with the eigenvalue staying close to zero. For $`a\rho <1.25`$, it even moves back and takes a negative value at $`a\rho =0.5`$. This can be interpreted as an ’overimprovement’ caused by the fairly large coefficients for the fluctuation polynomials  in our parametrization: If at some locations the gauge fields fluctuate very strongly, as it is the case for such an artificial small instanton, the fluctuation polynomials shift the couplings in the Dirac operator away from reasonable values. The histogram of plaquette values in Fig. 4.5 illuminates the qualitative difference of the fluctuations in these instanton and in thermal Monte Carlo gauge configurations. While on a Monte Carlo configuration the plaquette distribution is smooth, most of the plaquettes of an instanton configuration are very close to $`u=3`$, but there is also a peak at $`u1.6`$ from the plaquettes at the center of the instanton. As a consequence, the Dirac operator is strongly affected there through the fluctuation polynomials which are proportional to terms like $`3u`$, causing this unusual behavior of the real eigenvalue. We never observed such real eigenvalues shifted towards negative values in Monte Carlo gauge configurations used in actual lattice simulations, which were also used to parametrize the FP operator. Only in such an artificial environment as given by the discretized instantons the parametrized FP Dirac operator shows this effect.

We also investigate the effect of smearing for the FP operator. Obviously the RG inspired smearing, which we use together with the FP Dirac operator, does not change the flow of the real eigenvalue much. Only for the smallest instanton $`a\rho =0.5`$ the eigenvalue is pushed back slightly towards zero when the gauge configurations are smeared. Another approximately chiral operator that behaves very well on these instanton configurations is the chirally improved Dirac operator by Gattringer et al. , which is optimized for fulfilling the Ginsparg-Wilson relation.

As a second check we investigate the locality properties of the corresponding zero eigenmode $`\varphi _0`$, which is centered at the place where the instanton sits. With the gauge invariant density

$$p(x)\underset{\alpha ,a}{}\varphi _0^{\alpha a}(x)^{}\varphi _0^{\alpha a}(x)\stackrel{\mathrm{cont}.}{=}\frac{2\rho ^2}{\pi ^2(\rho ^2+x^2)^3},$$
(4.6)
where $`d^4xp(x)=1`$ for normalized eigenvectors, we can define as a measure for the localization of the zero mode $`\varphi _0`$ the inverse participation ratio

$$Id^4xp(x)^2\stackrel{\mathrm{cont}.}{=}\frac{1}{5\pi ^2\rho ^4},$$
(4.7)
which is plotted in Fig. 4.6 and compared to the continuum value. All Dirac operators agree for the largest instanton radii $`\rho 4`$, where the lattice results deviate from the continuum due to finite-volume effects. The Wilson overlap operator however deviates strongly from the continuum values also for smaller instantons, seemingly having problems to reproduce the continuum zero mode. This can be interpreted as a consequence of the relatively bad localization properties of the overlap with Wilson kernel. Varying the optimizing parameter $`s`$ in the Wilson overlap does not lead to a significant change in these result . The situation is different for the ultralocal Dirac operators like the parametrized FP and the Wilson operator. The parametrized FP operator reproduces the continuum value well even for very small instantons, while the Wilson operator performs somewhat worse, but still much better than the Wilson overlap.

The main result from this study is that the third order overlap-improved FP Dirac operator almost shows no change compared to the parametrized FP operator and captures the localization of the zero mode very well. The same observation has been made for an overlap operator with a chirally improved kernel, giving evidence that the Wilson kernel for the overlap construction, which misses the localization properties of the instanton zero mode, is not the best choice. Another observation we make when comparing the parametrized FP operator on both smeared and unsmeared gauge configurations is that the RG-inspired smearing of the gauge fields does not change the results for the inverse participation ratio. The smearing therefore does not seem to lead to a significant modification of the locality properties of the Dirac operator.

## Chapter 5 Hadron Spectroscopy in Lattice QCD

Among the most basic quantities that are calculable in lattice QCD are the masses and decay constants of the various bound states of quarks and gluons. Reproducing the experimentally observed spectrum of hadronic particles is one of the strongest tests that QCD is the correct theory to describe nature at the corresponding energy scale. In this chapter we present how masses of light hadrons, which are made from up, down and strange quarks<sup>1</sup><sup>1</sup>1Because the masses of the charm and bottom quarks are on the order of the lattice cutoff, lattice systems with heavy quarks have to be treated either in a non-relativistic approach , in the static approximation  or on anisotropic lattices., are extracted from quantities accessible on the lattice, and we discuss some refinements to improve the quality of these measurements.

First we derive in Section 5.1 the basic observation that the exponential decay of hadronic correlation functions is related to the hadron mass. The correlators can be calculated by inverting the lattice Dirac operator and contracting the resulting quark propagators together with the Dirac matrices corresponding to a particular hadronic state. Then we show in Section 5.2 how the overlap of the creation and annihilation operators with the ground state of the hadron can be increased by using extended wave functions at the source and sink locations. In the last section, we concentrate on the technical issues related to the fitting of lattice data to the predicted functional forms and show how the statistical error of the resulting hadron masses is estimated.

### 5.1 Fermionic Observables from Correlation Functions

The typical quantity to measure in lattice hadron spectroscopy is the two-point function

$$C(x)0|T\{𝒪_f(x)𝒪_i^{}(0)\}|0,$$
(5.1)
describing the space-time propagation of a particle created at the origin by the operator $`𝒪_i^{}`$ and annihilated at space-time coordinate $`x=(\stackrel{}{x},t)`$ by the operator $`𝒪_f`$. To create a meson, a quark bilinear operator of the form

$$𝒪^\mathrm{M}(x)\overline{\psi }_\mu ^{a,f_1}(x)\mathrm{\Gamma }_{\mu \nu }\psi _\nu ^{a,f_2}(x),$$
(5.2)
is used. The Clifford algebra element $`\mathrm{\Gamma }`$ determines the quantum numbers of the desired quark-antiquark state, and $`f_1`$, $`f_2`$ denote the flavors $`u`$,$`d`$,$`s`$ of the quark constituents. A baryon is created by the three-quark operator

$$𝒪^\mathrm{B}(x)ϵ_{abc}\chi _{\mu \nu \rho }\psi _\mu ^{a,f_1}(x)\psi _\nu ^{b,f_2}(x)\psi _\rho ^{c,f_3}(x),$$
(5.3)
with some appropriate spin function $`\chi _{\mu \nu \rho }`$. In the simulations in Chapter 7, we will work with the baryon operators used by the MILC collaboration , which create an equal mixture of a forward propagating baryon and a backward propagating antibaryon on a periodic lattice. In Tables 5.1 and 5.2 we list the spin content of our meson and baryon operators.

Let us demonstrate how to extract physical quantities from the correlation function (5.1) on a lattice of infinite volume. In order to single out particles with defined momenta, consider the spatial Fourier transform

$$C(\stackrel{}{p},t)=\underset{\stackrel{}{x}}{}e^{i\stackrel{}{p}\stackrel{}{x}}0|T\{𝒪_f(\stackrel{}{x},t)𝒪_i^{}(0,0)\}|0.$$
(5.4)
Inserting a complete set of eigenstates $`|n,\stackrel{}{q}`$ with spatial momentum $`\stackrel{}{q}`$, we get

$$C(\stackrel{}{p},t)=\underset{\stackrel{}{x}}{}e^{i\stackrel{}{p}\stackrel{}{x}}\underset{n}{}\frac{d^3\stackrel{}{q}}{(2\pi )^32E_n(\stackrel{}{q})}0|𝒪_f(\stackrel{}{x},t)|n,\stackrel{}{q}n,\stackrel{}{q}|𝒪_i^{}(0,0)|0,$$
(5.5)
where $`E_n(\stackrel{}{q})`$ is the energy of the intermediate state $`|n,\stackrel{}{q}`$. Applying the space-time translation $`𝒪_f(x)=e^{i𝒫x}𝒪_f(0)e^{i𝒫x}`$ with four-momentum $`𝒫=(,\stackrel{}{𝒫})`$ to the annihilation operator, the correlation function (5.5) can be written as an exponentially weighted sum over all intermediate states,

$`C(\stackrel{}{p},t)`$ $`={\displaystyle \underset{n}{}}{\displaystyle \frac{d^3\stackrel{}{q}}{(2\pi )^32E_n(\stackrel{}{q})}e^{iE_n(\stackrel{}{q})t}0|𝒪_f|n,\stackrel{}{q}n,\stackrel{}{q}|𝒪_i^{}|0\underset{\stackrel{}{x}}{}e^{i(\stackrel{}{p}\stackrel{}{q})\stackrel{}{x}}}`$
$`={\displaystyle \underset{n}{}}{\displaystyle \frac{0|𝒪_f|n,\stackrel{}{p}n,\stackrel{}{p}|𝒪_i^{}|0}{2E_n(\stackrel{}{p})}}e^{iE_n(\stackrel{}{p})t},`$ (5.6)
where $`_\stackrel{}{x}e^{i\stackrel{}{p}\stackrel{}{x}}=(2\pi )^3\delta (\stackrel{}{p})`$ has been used to get rid of one momentum variable, and $`𝒪_{i,f}𝒪_{i,f}(0,0)`$. For large Euclidean time $`\tau =it`$, only the state with lowest energy contributes, therefore the asymptotic form of the correlation function becomes

$$C(\stackrel{}{p},\tau )\stackrel{\tau \mathrm{}}{}\frac{0|𝒪_f|1,\stackrel{}{p}1,\stackrel{}{p}|𝒪_i^{}|0}{2E_1(\stackrel{}{p})}e^{E_1(\stackrel{}{p})\tau }.$$
(5.7)
Considering only intermediate states with zero momentum by setting $`\stackrel{}{p}=0`$, the mass of the lightest state can be extracted from the exponential decay of the Euclidean time correlation function

$$C(\tau )\stackrel{\tau \mathrm{}}{}\frac{A_fA_i^{}}{2m_1}e^{m_1\tau },$$
(5.8)
where the amplitudes $`A_{i,f}=0|𝒪_{i,f}|1`$ represent the overlap of the operators $`𝒪_{i,f}`$ with the lightest particle state having the quantum numbers of $`𝒪_{i,f}`$.

In Fig. 5.1, we plot typical examples of meson correlation functions on the lattice at varying quark masses, showing clearly the exponential decay with Euclidean time.

#### 5.1.1 Lattice Quark Propagators

In order to measure the correlation function (5.1) in lattice QCD, it is expressed in terms of the Euclidean quark propagator

$$G_{\mu \nu }^{ab}(x,y)0|T\{\psi _\mu ^a(x)\overline{\psi }_\nu ^b(y)\}|0,$$
(5.9)
using Wick contractions of the quark fields. Consider for example a flavor non-singlet meson operator (5.2), where $`f_1f_2`$. The two-point function is then

$`𝒪^\mathrm{M}(\stackrel{}{x},t)𝒪^\mathrm{M}(0,0)`$ $`0|T\{𝒪^\mathrm{M}(\stackrel{}{x},t)𝒪^\mathrm{M}(0,0)\}|0`$
$`=0|T\{\overline{\psi }_\sigma ^{a,f_1}(x)\mathrm{\Gamma }_{\sigma \mu }\psi _\mu ^{a,f_2}(x)\overline{\psi }_\nu ^{b,f_2}(0)\mathrm{\Gamma }_{\nu \rho }\psi _\rho ^{b,f_1}(0)\}|0`$
$`=G_{\mu \nu }^{ab,f_2}(x,0)\mathrm{\Gamma }_{\nu \rho }G_{\rho \sigma }^{ba,f_1}(0,x)\mathrm{\Gamma }_{\sigma \mu }`$
$`=\mathrm{Tr}G^{f_2}(x,0)\mathrm{\Gamma }G^{f_1}(0,x)\mathrm{\Gamma },`$ (5.10)
where the trace is taken over spin and color indices. Making use of the $`\gamma ^5`$-hermiticity of the quark propagator $`G(x,y)=\gamma ^5G(y,x)^{}\gamma ^5`$, where the hermitean conjugation also acts in spin and color space, Eq. (5.1.1) becomes

$$𝒪^\mathrm{M}(\stackrel{}{x},t)𝒪^\mathrm{M}(0,0)=\mathrm{Tr}G^{f_2}(x,0)\mathrm{\Gamma }\gamma _5G^{f_1}(x,0)^{}\gamma _5\mathrm{\Gamma }.$$
(5.11)
The importance of this last step comes from the fact that in lattice QCD calculations, the quark propagator $`G(x,y)=D^1(x,y)`$ is determined by a matrix inversion of the lattice Dirac operator, which is by far the most expensive part of quenched simulations. A complete inversion amounts to solving the linear system of $`12V`$ equations

$$D_{\mu \nu }^{ab}(x,y)G_{\nu \rho }^{bc}(y,z)=\delta _{ac}\delta _{\mu \rho }\delta (xz),$$
(5.12)
for $`G`$, which is in most cases not feasible. For the calculation of the right hand side of (5.11), it is however only necessary to know the quark propagator $`G(x,0)`$ from a fixed source point at the origin, where the particle is created, to all points on the lattice. Hence it is sufficient to evaluate 12 columns of the inverted matrix (one per spin and color) by solving

$$D_{\mu \nu }^{ab}(x,y)G_{\nu \rho }^{bc}(y,0)=\delta _{ac}\delta _{\mu \rho }\delta (x),$$
(5.13)
reducing the numerical size of the problem by a factor of $`V`$. This trick does not work in the case of flavor singlet mesons, where the disconnected contribution in the quark line graph requires the knowledge of the full quark propagator, which is the reason why their treatment on the lattice is much more demanding.

Finally, to get the zero momentum correlator (5.8), it is sufficient to sum Eq. (5.11) over all sink locations on a given timeslice,

$$C(\tau )=\underset{\stackrel{}{x}}{}\mathrm{Tr}G^{f_2}(x,0)\mathrm{\Gamma }\gamma _5G^{f_1}(x,0)^{}\gamma _5\mathrm{\Gamma },$$
(5.14)
which projects out the $`\stackrel{}{p}=0`$ contribution in (5.4).

### 5.2 Extended Source and Sink Operators

In order to get a good signal-to-noise ratio for the measured hadron correlators, the operators $`𝒪_i`$ and $`𝒪_f`$ in Eq. (5.1) should have a large overlap with the desired state. Local operators like (5.2) and (5.3) are not expected to fulfill this criterion well, as they do not take into account the spatial extension of the hadron. Especially for small quark masses or small lattice spacings, neglecting the hadron extension becomes a problem, as light hadrons typically have a size of $`𝒪`$(1 fm), and lattice spacings in current simulations are mostly between 0.05 fm and 0.2 fm. Maximum overlap for a meson would be reached for a non-local operator

$$𝒪_{\mathrm{ext}}^\mathrm{M}(t)=\underset{\stackrel{}{x},\stackrel{}{y}}{}\phi (\stackrel{}{x},\stackrel{}{y})\overline{\psi }_\mu ^{a,f_1}(\stackrel{}{x},t)\mathrm{\Gamma }_{\mu \nu }\psi _\nu ^{a,f_2}(\stackrel{}{y},t),$$
(5.15)
where $`\phi (\stackrel{}{x},\stackrel{}{y})`$ is the wave function of the meson. For a delta function $`\phi (\stackrel{}{x},\stackrel{}{y})=\delta (\stackrel{}{x}\stackrel{}{y})`$, the local operator (5.2) is reproduced.

Spatially extended operators like (5.15) are often referred to as smeared operators. They are in this form not gauge invariant quantities, and therefore their average over gauge configurations would vanish due to Elitzur’s theorem . To prevent that, one either has to include the parallel transporters in the operator or to work in a fixed gauge background. While gauge fixing is technically easy and imposes no restrictions on the wave function $`\phi (\stackrel{}{x},\stackrel{}{y})`$, it introduces a possible source of errors due to the Gribov copy problem (see Appendix A). To avoid this problem, various kinds of gauge invariant operators like Jacobi-smeared  or Wuppertal sources , have been constructed. In the following however, we will concentrate on gauge non-invariant operators and measure them on gauge-fixed configurations, which is a common procedure adopted in many large-scale simulations .

The wave function $`\phi (\stackrel{}{x},\stackrel{}{y})`$ of the simulated particle is a priori not known, so one has to make a more or less reasonable guess. A convenient, but not very physical choice is the totally factorized shell-model wave function

$$\phi (\stackrel{}{x},\stackrel{}{y})=\varphi (\stackrel{}{x})\varphi (\stackrel{}{y}),$$
(5.16)
where the quark and antiquark move independently inside a region given by the function $`\varphi (\stackrel{}{x})`$, which might for example chosen to be a wall , a hard sphere or cube , a Gaussian , or a radial exponential . General experience from lattice simulations has shown that using these kinds of extended operators, it is significantly easier to extract a reliable estimate for hadron masses from a fit to the correlation function (5.1), because due to the larger overlap with the ground state the contributions from higher states in (5.1) vanish at much smaller time separation $`\tau `$ than for local hadron operators <sup>2</sup><sup>2</sup>2To account for the higher state contributions, refined strategies are to perform a double exponential fit to both the lowest and the first excited state or to make Bayesian fits ..

The advantage of the shell-model wave function (5.16) over more physical functions which depend on the relative coordinates between the quark and antiquark is that due to the factorization into separate quark and antiquark parts, the quark propagator can be calculated like for a local operator with only one inversion per spin and color. In the following, we suppose that $`\varphi (\stackrel{}{x})`$ is a real, radial symmetric function,

$$\varphi _\stackrel{}{r}(\stackrel{}{x})=\varphi (\stackrel{}{r}\stackrel{}{x})$$
(5.17)
around a center $`\stackrel{}{r}`$ of the source or sink. Consider first a meson correlation function with a shell-model source operator (5.15)–(5.17)

$$𝒪_{\mathrm{sm}}^\mathrm{M}(\stackrel{}{r},t)=\underset{\stackrel{}{x},\stackrel{}{y}}{}\varphi _\stackrel{}{r}(\stackrel{}{x})\overline{\psi }_\mu ^{a,f_1}(\stackrel{}{x},t)\mathrm{\Gamma }_{\mu \nu }\varphi _\stackrel{}{r}(\stackrel{}{y})\psi _\nu ^{a,f_2}(\stackrel{}{y},t),$$
(5.18)
at timeslice $`t=0`$ and centered at $`\stackrel{}{r}=0`$, and a local sink operator at time $`t`$:

$`𝒪^\mathrm{M}(\stackrel{}{x},t)𝒪_{\mathrm{sm}}^\mathrm{M}`$ $`(0,0)=0|\overline{\psi }^{f_1}(\stackrel{}{x},t)\mathrm{\Gamma }\psi ^{f_2}(\stackrel{}{x},t)`$
$`\times \left({\displaystyle \underset{\stackrel{}{x}_0,\stackrel{}{y}_0}{}}\varphi _0(\stackrel{}{x}_0)\overline{\psi }^{f_2}(\stackrel{}{x}_0,0)\mathrm{\Gamma }\varphi _0(\stackrel{}{y}_0)\psi ^{f_1}(\stackrel{}{y}_0,0)\right)^{}|\mathrm{\hspace{0.17em}0}`$
$`=\mathrm{Tr}{\displaystyle \underset{\stackrel{}{x}_0}{}}G^{f_2}(\stackrel{}{x},t;\stackrel{}{x}_0,0)\varphi _0(\stackrel{}{x}_0)\mathrm{\Gamma }{\displaystyle \underset{\stackrel{}{y}_0}{}}G^{f_1}(\stackrel{}{y}_0,0;\stackrel{}{x},t)\varphi _0(\stackrel{}{y}_0)\mathrm{\Gamma }`$
$`=\mathrm{Tr}G_{\mathrm{sm}}^{f_2}(\stackrel{}{x},t;0,0)\mathrm{\Gamma }\gamma _5G_{\mathrm{sm}}^{f_1}(\stackrel{}{x},t;0,0)^{}\gamma _5\mathrm{\Gamma }.`$ (5.19)
The spatial distribution of the source is taken into account when inverting the Dirac operator on the vector $`\varphi _0(\stackrel{}{x})\delta (t)`$ instead of $`\delta (x)`$, defining the smeared source quark propagator $`G_{\mathrm{sm}}`$ in Eq. (5.2) as the solution of

$$D(x,y)G_{\mathrm{sm}}(y,0)=\varphi _0(\stackrel{}{x})\delta (t).$$
(5.20)
While introducing an extended source amounts to inverting the Dirac operator on a different source vector, using an extended sink leads to a weighting of the quark propagator at different lattice sites on a given timeslice, as can be seen from the smeared-source, smeared-sink meson correlator

$$\begin{array}{c}𝒪_{\mathrm{sm}}^\mathrm{M}(\stackrel{}{r},t)𝒪_{\mathrm{sm}}^\mathrm{M}(0,0)=\mathrm{Tr}\underset{\stackrel{}{x}_t}{}\varphi _\stackrel{}{r}(\stackrel{}{x}_t)G_{\mathrm{sm}}^{f_2}(\stackrel{}{x}_t,t;0,0)\mathrm{\Gamma }\gamma _5\hfill \\ \hfill \times \underset{\stackrel{}{y}_t}{}\varphi _\stackrel{}{r}(\stackrel{}{y}_t)G_{\mathrm{sm}}^{f_1}(\stackrel{}{y}_t,t;0,0)^{}\gamma _5\mathrm{\Gamma },\end{array}$$
(5.21)
where the smearing of the source again is absorbed in the quark propagator.

As shown above, the hadron mass is extracted from the Fourier transform of the correlation function at zero momentum, which implies summing Eq. (5.21) over all sink locations $`\stackrel{}{r}`$. For a smeared sink, the numerical effort can get quite large, as there are then three sums over all lattice points on a given time slice. A technical trick to accelerate the calculation of smeared-sink meson correlators is to rewrite Eq. (5.21) in Fourier space, which allows to make use of efficient Fast Fourier Transform (FFT) algorithms to speed up the calculation. With the discrete Fourier representations of both the sink wave function

$$\varphi _\stackrel{}{r}(\stackrel{}{x})=\frac{1}{N_s^{3/2}}\underset{\stackrel{}{k}}{}e^{i\frac{2\pi }{N}\stackrel{}{k}(\stackrel{}{r}\stackrel{}{x})}\varphi (\stackrel{}{k}),$$
(5.22)
where $`k_i=0,\mathrm{},N_s1`$ for $`i=1,\mathrm{},3`$, and the quark propagator

$$G(\stackrel{}{x},t;0,0)=\frac{1}{N_s^{3/2}}\underset{\stackrel{}{k}}{}e^{i\frac{2\pi }{N}\stackrel{}{k}\stackrel{}{x}}G(\stackrel{}{k},t;0,0),$$
(5.23)
the convolutions in the meson correlator (5.21) can be expressed as

$$\underset{\stackrel{}{x}}{}\varphi _\stackrel{}{r}(\stackrel{}{x}_t)G(\stackrel{}{x}_t,t;0,0)=\underset{\stackrel{}{k}}{}e^{i\frac{2\pi }{N}\stackrel{}{k}\stackrel{}{r}}\varphi (\stackrel{}{k})G(\stackrel{}{k},t;0,0),$$
(5.24)
leading to the zero-momentum smeared-source, smeared-sink meson correlator

$$\begin{array}{c}\underset{\stackrel{}{r}}{}𝒪_{\mathrm{sm}}^\mathrm{M}(\stackrel{}{r},t)𝒪_{\mathrm{sm}}^\mathrm{M}(0,0)=N_s^3\underset{\stackrel{}{k}}{}\mathrm{Tr}\varphi (\stackrel{}{k})G_{\mathrm{sm}}^{f_2}(\stackrel{}{k},t;0,0)\mathrm{\Gamma }\gamma _5\hfill \\ \hfill \times \varphi (\stackrel{}{k})G_{\mathrm{sm}}^{f_1}(\stackrel{}{k},t;0,0)^{}\gamma _5\mathrm{\Gamma }.\end{array}$$
(5.25)
With this trick, the smeared-sink meson correlator at zero momentum is calculated in the same manner as the point-sink case (5.14) after replacing the 3-space fields by their Fourier transforms.

The effect of improving the overlap of the interpolating operators with the desired hadron state can be seen in the effective mass plots<sup>3</sup><sup>3</sup>3Effective masses are explained in the next section. in Fig. 5.2. Compared are effective masses of pseudoscalar and vector mesons at large quark mass $`am=0.23`$ and lattice spacing $`a=0.16`$ fm for meson correlators made with point and Gaussian smeared sources, respectively. In both cases point sinks were chosen, and 70 gauge configurations of lattice size $`9^3\times 24`$ were used. The Gaussian smearing is defined by the shell-model wave function

$$\varphi _\stackrel{}{r}(\stackrel{}{x})=e^{\gamma (\stackrel{}{x}\stackrel{}{r})^2},$$
(5.26)
and the source center $`r_i=N_i/2`$ with $`i=1,\mathrm{},3`$ is located at the center of the time slice. The extension of the source used in the plots is given by $`\gamma =0.15`$. The plateau in the effective mass starts much earlier in most cases for smeared sources, providing a larger time interval over which the correlators can be fitted. The lattice spacing in the plots is rather large, and the effect will be enhanced at smaller $`a`$, where the point source correlators might not even reach a plateau within the given number of temporal lattice points. Fig. 5.3 shows plots for the effective energy of the pseudoscalar and vector meson at the lowest non-zero momentum $`|ap|=2\pi /9`$. There the signal is worse, and the length of the plateau for smeared correlators is only slightly increased, as the signal starts to deviate at large $`t`$. However, again the plateau region starts much earlier when using smeared sources, thus raising the confidence in that really the asymptotic behavior of the correlator is reached. Also, the statistical error of the energy resulting from the fit is significantly smaller for the smeared source correlators.

### 5.3 Fitting Hadron Propagators

As we have shown, hadron masses are extracted from the exponential fall-off of Euclidean time correlation functions at zero momentum (5.8), which are expressed in terms of the quark propagator and can thus be evaluated on the lattice. Consider the case of flavor non-singlet mesons Eq. (5.14), of which the correlation function $`C(t)`$ is measured on a lattice of temporal size $`T`$ for all $`0t<T`$.<sup>4</sup><sup>4</sup>4Staying from now on in Euclidean space, we will use the notation $`t`$ also for Euclidean time. Using periodic boundary conditions, the data points can be fitted against the asymptotic form

$$\begin{array}{cc}\hfill C_{\mathrm{}}(t)& =Z\left(e^{mt}+e^{m(Tt)}\right),\hfill \\ & =2Ze^{mT/2}\mathrm{cosh}\left(m(T/2t)\right),\hfill \end{array}$$
(5.27)
of the meson correlator in Euclidean time $`t`$ in order to determine the mass $`m`$ and the coefficient

$$Z=\frac{A_fA_i^{}}{2m},$$
(5.28)
containing the information on the matrix elements of the chosen operators between the vacuum and the hadron state.

In order to find out at what time $`t_0`$ the asymptotic range of the hadron propagator is reached, it is helpful to plot the effective mass $`m_{\mathrm{eff}}(t+1/2)`$, which is determined from a zero parameter fit of the asymptotic function $`C_{\mathrm{}}(t)`$ to the measured propagators $`C(t)`$ and $`C(t+1)`$ at two consecutive time slices. When $`t`$ gets large enough that the higher lying states have disappeared, $`m_{\mathrm{eff}}`$ starts to show a plateau. In general, only the measured correlators in a time interval $`t[t_0,t_1]`$ are then used to fit the parameters. The upper bound $`t_1`$ can be set to the point where the signal disappears in the statistical noise, which happens for all particles but pseudoscalar mesons after a certain temporal range.

In Figs. 5.55.8, we show examples of effective masses, fitted masses and the quality of the fit given by the value of $`\chi ^2/df`$ for several particles at intermediate quark mass. In all cases, $`t_1`$ was set to $`T/2`$. For all hadrons, $`\chi ^2/df`$ quickly decreases and stays then at a value of order 1, and the optimal fit interval starts in the range $`3t_06`$.

#### 5.3.1 Correlated Fits

Suppose we have measured the hadron propagators<sup>5</sup><sup>5</sup>5To avoid confusion between correlation functions and correlation matrices, we will denote hadronic correlation functions as hadron propagators where necessary. $`C^{(i)}(t)`$ for $`t=0,\mathrm{},T`$ on $`N`$ independent, importance sampled gauge configurations $`U^{(i)}`$ with $`i=1,\mathrm{},N`$. Trying to fit the gauge average $`\overline{C}(t)=1/N_iC^{(i)}(t)`$, we are faced with the problem that while the data is uncorrelated in Monte Carlo time $`i`$, it is strongly correlated in the temporal direction. The goal is to find the optimal parameters $`m`$ and $`Z`$ in the asymptotic form (5.27), taking into account the time correlations. This is done by minimizing the $`\chi ^2`$-function

$$\chi ^2=\underset{t,t^{}}{}\left\{\overline{C}(t)C_{\mathrm{}}(t;m,Z)\right\}(\mathrm{Cov})^1(t,t^{})\left\{\overline{C}(t^{})C_{\mathrm{}}(t^{};m,Z)\right\},$$
(5.29)
where the time correlations are encoded in the symmetric, positive definite covariance matrix

$$(\mathrm{Cov})(t,t^{})=\frac{1}{N(N1)}\underset{i=1}{\overset{N}{}}\left(C^{(i)}(t)\overline{C}(t)\right)\left(C^{(i)}(t^{})\overline{C}(t^{})\right).$$
(5.30)
As an illustration of the typical size of time correlations we plot a row of the normalized correlation matrix

$$\mathrm{\Sigma }(t,t^{})=\frac{(\mathrm{Cov})(t,t^{})}{\sqrt{(\mathrm{Cov})(t,t)(\mathrm{Cov})(t^{},t^{})}},$$
(5.31)
for different smeared-source hadron propagators in Fig. 5.4. The data is strongly correlated for all hadrons under consideration. For the pseudoscalar meson, the time correlations do not die out at large $`t`$, as the signal-to-noise ratio remains constant, whereas for the other hadrons what mostly remains is uncorrelated statistical noise. The time correlations of pseudoscalar propagators are even stronger for point sources , therefore it is in any case mandatory to perform correlated fits by including the covariance matrix.

#### 5.3.2 Resampling Methods for Error Estimates

The results of the above described procedure are the parameters of the fit function, that is the hadron mass $`m`$ and the amplitude $`Z`$. Since these fit parameters are statistical estimates from a Monte Carlo integration, it is necessary to provide an estimate of their statistical errors in order to judge their reliability. The measured observables are the hadron correlators, and what is needed is a tool to estimate the errors of quantities which depend in a complicated way on these observables. For the sample mean of the correlators, the standard deviation can be calculated in the usual way,

$$\sigma _{\overline{C}}=\left(\frac{1}{N(N1)}\underset{i=1}{\overset{N}{}}(C^{(i)}\overline{C})^2\right)^{1/2},$$
(5.32)
but for less trivial functions of the observables, there is no such function to estimate the error. With the advent of powerful computers, robust statistical methods have been developed which allow to estimate errors of arbitrarily complicated functions of observables with unknown probability distributions in a straightforward way. These widely used methods are known under the names jackknife and bootstrap, and they are both based on resampling of the measured data. We briefly present the jackknife and the bootstrap resampling techniques in the following.

Suppose we have a set of $`N`$ independent and identically distributed measurements

$$X_1,\mathrm{},X_N\stackrel{\mathrm{iid}}{}F,$$
(5.33)
following an unknown distribution function $`F`$. From these data points, an arbitrarily complicated secondary quantity $`\widehat{\theta }(X_1,\mathrm{},X_N)`$ is calculated. What we aim at is an expression for the standard deviation $`\sigma _{\widehat{\theta }}(F)`$ of the estimator $`\widehat{\theta }`$. Jackknife resampling requires to calculate the estimator

$$\widehat{\theta }_{(i)}=\widehat{\theta }(X_1,\mathrm{},X_{i1},X_{i+1},\mathrm{},X_N),$$
(5.34)
on the sample where the data point $`X_i`$ has been dropped, for all $`i=1,\mathrm{},N`$. The jackknife estimate for the error of $`\widehat{\theta }`$ is then defined by

$$\sigma _{\widehat{\theta }}^{\mathrm{jack}}=\left(\frac{N1}{N}\underset{i=1}{\overset{N}{}}(\widehat{\theta }_{(i)}\widehat{\theta }_{()})^2\right)^{1/2},$$
(5.35)
where $`\widehat{\theta }_{()}=_{i=1}^N\widehat{\theta }_{(i)}/N`$ is the average of the estimator over all jackknife samples.

The jackknife utilizes only $`N`$ of the $`2^N1`$ non-empty subsets of the data set. The error estimate might thus be improved when more of the subsets are used. This lead to the development of the bootstrap : A bootstrap sample

$$X_1^{},\mathrm{},X_N^{}\stackrel{\mathrm{iid}}{}\widehat{F},$$
(5.36)
is a random sample drawn with replacement from the observed values $`X_1,\mathrm{},X_N`$ and follows the empirical probability distribution $`\widehat{F}`$ of the data. The bootstrap error estimate for $`\widehat{\theta }`$ is defined from the estimator $`\widehat{\theta }^{}=\widehat{\theta }(X_1^{},\mathrm{},X_N^{})`$ calculated on $`B`$ bootstrap samples by

$$\sigma _{\widehat{\theta }}^{\mathrm{boot}}=\left(\frac{1}{B1}\underset{b=1}{\overset{B}{}}(\widehat{\theta }_{(b)}^{}\widehat{\theta }_{()}^{})^2\right)^{1/2},$$
(5.37)
where $`\widehat{\theta }_{()}^{}=_{b=1}^B\widehat{\theta }_{(b)}^{}/B`$ denotes the average over the bootstrap samples. In the limit $`B\mathrm{}`$, the bootstrap error is exactly the standard deviation of the estimator as a function of the empirical probability distribution $`\widehat{F}`$,

$$\sigma _{\widehat{\theta }}^{\mathrm{boot}}\stackrel{B\mathrm{}}{}\sigma _{\widehat{\theta }}(\widehat{F})\sigma _{\widehat{\theta }}(F).$$
(5.38)
In practice, the number of bootstrap samples $`B`$ is finite, and one has to make sure that it is large enough by varying $`B`$ and checking whether the error changes significantly. But since $`\widehat{F}`$ is only an estimate for the unknown probability distribution $`F`$, taking too many bootstrap samples does not help improving the error estimate for $`\sigma _{\widehat{\theta }}(F)`$. For most cases, it is considered safe to work with values in the range $`100<B<500`$. In our spectroscopy simulations, the computational cost for the bootstrap is negligible, and therefore we always calculate $`B=1000`$ bootstrap samples.

The jackknife and bootstrap procedures also provide an estimate for the bias $`\widehat{b}_{\widehat{\theta }}`$ of the estimator $`\widehat{\theta }`$. The jackknife estimate of bias is given by

$$\widehat{b}_{\widehat{\theta }}^{\mathrm{jack}}=(N1)(\widehat{\theta }_{()}\widehat{\theta }),$$
(5.39)
while for the bootstrap it is just the difference of the value of the estimator on the original sample and its mean value on the bootstrap samples,

$$\widehat{b}_{\widehat{\theta }}^{\mathrm{boot}}=\widehat{\theta }_{()}^{}\widehat{\theta }.$$
(5.40)
The bias-corrected estimator $`\stackrel{~}{\theta }=\widehat{\theta }\widehat{b}`$ might then be used instead of $`\widehat{\theta }`$ for a better estimate of the true value $`\theta `$.

## Chapter 6 Topological Finite-Volume Artifacts in Pion Propagators

In the chiral limit of quenched QCD, pion<sup>1</sup><sup>1</sup>1We denote in this chapter the pseudoscalar meson as a pion also for unphysical quark masses. propagators suffer from unphysical quenching effects which make a thorough examination of light pions difficult. These quenching artifacts are caused by zero eigenmodes of the Dirac operator and lead to unphysical divergences of the pion propagators at $`m0`$. In unquenched QCD, gauge configurations with zero modes are suppressed by the fermion determinant in the effective action, therefore these divergences are absent. Because the determinant is set to unity in the quenched theory, the suppression falls away, and at small quark mass the propagators are dominated by the zero mode effects. This can be clearly seen in Fig. 6.1, where the chiral limit of $`m_{\mathrm{PS}}^2(m_q)`$ is shown in dependence of the topology of the gauge configurations. While the pseudoscalar meson mass goes to zero in the chiral limit on configurations with trivial topology, it deviates as soon as also configurations with non-zero topological charge are considered in the Monte Carlo average.

For the study of light pseudoscalar mesons with masses around the physical mass of the pion, it is therefore unavoidable to get rid of these artificial effects. In particular for the investigation of quenched chiral logarithms in the pseudoscalar mass, which we undertake in Chapter 7, the effects from zero modes have to be properly disentangled from the chiral logarithm, which produces a measurable signal only at small quark masses. Due to the explicit breaking of chiral symmetry, traditional formulations of lattice fermions like Wilson or Symanzik-improved clover fermions do not allow to identify topological zero modes unambiguously. Only with the development of chiral symmetric lattice Dirac operators, it has become possible to identify the zero modes responsible for these unwanted effects.

In a study with the Wilson overlap operator , the authors have reported a change in the behavior of the pseudoscalar correlator at large time, suggesting that the zero modes only contaminate the small $`t`$ range. As a possible solution, they proposed to fit the mass from the large $`t`$ tail of the correlator. Fig. 6.2 shows the reported kink in the correlator as seen in our data, but even a fit to the flatter region does not give a pion mass which goes to zero in the chiral limit. Furthermore, we could not clearly identify such a kink in all our simulations.

In this chapter we derive and examine two other solutions of the problem. One solution is based on explicit identification and subtraction of the zero modes in the quark propagator. The other solution, originally proposed in , makes use of the fact that zero mode effects enter pseudoscalar and scalar meson propagators equivalently. The zero mode effects can then be subtracted in the meson propagators. We study these two solutions on a very small lattice of spatial physical extension $`L_s1`$ fm, where the zero mode effects are large.

### 6.1 Zero Mode Subtraction of the Quark Propagator

The most straightforward way to get rid of zero mode contributions is to subtract them directly from the quark propagator. An exact subtraction however is only possible for an exactly chiral Dirac operator. We derive in the following the subtraction of the zero modes for the propagator of the overlap-improved FP Dirac operator described in Chapter 4, which is not trivial due to the Fixed-Point $`R`$ appearing in the Ginsparg-Wilson relation. From the subtracted quark propagators, meson and baryon correlators can be constructed in the standard way, and measurements of hadronic quantities derived from these correlators should then be free of quenched topological finite-volume artifacts.

#### 6.1.1 Spectral Decomposition of the Massless Normal Dirac Operator

The overlap-improved FP Dirac operator is a solution of the Ginsparg-Wilson relation $`\{D,\gamma _5\}=D\gamma _52RD`$, where a non-trivial $`R`$ appears on the right-hand side. Defining the operator

$$\stackrel{~}{D}\sqrt{2R}D\sqrt{2R},$$
(6.1)
the Ginsparg-Wilson relation reduces to the simpler case $`\{\stackrel{~}{D},\gamma _5\}=\stackrel{~}{D}\gamma _5\stackrel{~}{D}`$. With the $`\gamma _5`$-hermiticity of Dirac operators $`D^{}=\gamma _5D\gamma _5`$, it follows that $`\stackrel{~}{D}`$ is a normal operator,

$$\stackrel{~}{D}^{}\stackrel{~}{D}=\stackrel{~}{D}\stackrel{~}{D}^{}.$$
(6.2)
Normal operators can be written as diagonal matrices in an orthonormal basis of eigenvectors. Thus we can write down the spectral decomposition of a matrix element of the massless normal Dirac operator

$$\stackrel{~}{D}_{ij}=\underset{\lambda }{}\lambda \stackrel{~}{\varphi }_\lambda (i)\stackrel{~}{\varphi }_\lambda ^{}(j),$$
(6.3)
where $`\lambda `$ and $`\stackrel{~}{\varphi }_\lambda `$ are the eigenvalues and the corresponding eigenvectors of $`\stackrel{~}{D}`$, and the indices $`i`$ and $`j`$ contain spin, color and space-time degrees of freedom. Since the inversion of a diagonal matrix is trivial, the spectral decomposition of the quark propagator $`\stackrel{~}{G}_{ij}(\stackrel{~}{D}^1)_{ij}`$ is given by

$$\stackrel{~}{G}_{kl}=\underset{\lambda }{}\frac{1}{\lambda }\stackrel{~}{\varphi }_\lambda (k)\stackrel{~}{\varphi }_\lambda ^{}(l).$$
(6.4)
The orthonormality condition on the eigenvectors reads

$$\underset{l}{}\stackrel{~}{\varphi }_\lambda ^{}(l)\stackrel{~}{\varphi }_\lambda ^{}(l)=\delta _{\lambda \lambda ^{}}.$$
(6.5)
We define the subtracted quark propagator by summing only over the non-zero eigenmodes of $`\stackrel{~}{D}`$,

$$\stackrel{~}{G}_{kl}^{(\mathrm{sub})}\underset{\lambda 0}{}\frac{1}{\lambda }\stackrel{~}{\varphi }_\lambda (k)\stackrel{~}{\varphi }_\lambda ^{}(l).$$
(6.6)
In presence of a finite mass, we will invert a matrix of the form $`a\stackrel{~}{D}+b`$ with $`a,b𝐑`$, and Eq. (6.6) can be written as

$$\stackrel{~}{G}_{kl}^{(\mathrm{sub})}(a,b)=\stackrel{~}{G}_{kl}(a,b)\underset{\lambda =0}{}\frac{1}{b}\stackrel{~}{\varphi }_\lambda (k)\stackrel{~}{\varphi }_\lambda ^{}(l).$$
(6.7)
#### 6.1.2 Basis Transformation

For practical applications, it is not convenient to work with $`\stackrel{~}{D}`$ as defined in Eq. (6.1) due to the appearance of the square root of $`2R`$, whose calculation is a non-trivial numerical problem. A simple basis transformation

$$\widehat{\varphi }S\stackrel{~}{\varphi },$$
(6.8)
with $`S=(2R)^{1/2}`$ helps to get rid of this square root. The eigenvalues $`\lambda `$ remain unchanged under this transformation. The application of the basis transformed Dirac operator

$$\widehat{D}S\stackrel{~}{D}S^1=D2R,$$
(6.9)
on a vector is then reduced to multiplications with both D and R. From the definition (6.6),(6.7) we can read off the subtracted propagator in the new basis,

$$\widehat{G}_{ij}^{(\mathrm{sub})}(a,b)S\stackrel{~}{G}_{ij}^{(\mathrm{sub})}(a,b)S^1=\widehat{G}_{ij}(a,b)\underset{\lambda =0}{}\frac{1}{b}\widehat{\varphi }_\lambda (i)\widehat{\varphi }_\lambda ^{}(k)2R_{kj},$$
(6.10)
where $`\widehat{\varphi }`$ are eigenvectors of $`\widehat{D}`$ and $`\widehat{G}_{kl}(\widehat{D}^1)_{kl}`$. We have to remark that due to the basis transformation, the Dirac operator $`\widehat{D}`$ is no longer normal, and its eigenvectors are not orthogonal, but instead fulfill the generalized orthonormality condition

$$\underset{k,l}{}\widehat{\varphi }_\lambda ^{}(k)2R_{kl}\widehat{\varphi }_\lambda ^{}(l)=\delta _{\lambda \lambda ^{}}.$$
(6.11)
#### 6.1.3 A Cookbook Recipe

The above discussion is valid for general massless Ginsparg-Wilson Dirac operators. We consider in the following the overlap-improved FP Dirac operator (4.1). After basis transformation (6.9), the mass dependence of $`\widehat{D}(m)`$ is the usual one,

$$\widehat{D}(m)\left(1\frac{m}{2}\right)\widehat{D}+m,$$
(6.12)
where the $`2R`$ in the denominator has disappeared.

With these ingredients, we are ready to give a cookbook recipe for the calculation of zero-mode subtracted quark propagators with the massive overlap-improved FP Dirac operator. These are the solutions of the equation

$$D_{kl}^{(\mathrm{sub})}(m)g^{(\mathrm{sub})}(k;m)=b(l),$$
(6.13)
where the subtracted Dirac operator is defined through Eqs. (6.9) and (6.10). The steps to find the solution $`g^{(\mathrm{sub})}(k;m)`$ are the following:

1. Calculate the few smallest eigenmodes of $`\widehat{D}`$,

$$\widehat{D}\widehat{\varphi }=D2R\widehat{\varphi }=\lambda \widehat{\varphi },$$
(6.14)
identify the zero modes, where $`\lambda =0`$, and normalize them according to (6.11).
2. Invert $`\widehat{D}(m)`$ on a source $`b(y)`$ by solving the equation

$$\widehat{D}_{xy}(m)\widehat{g}(x;m)=b(y),$$
(6.15)
for the vector $`\widehat{g}(x;m)`$.
3. Subtract the zero modes from the solution $`\widehat{g}(x;m)`$ as derived in Eq. (6.10):

$$\widehat{g}^{(\mathrm{sub})}(x;m)=\widehat{g}(x;m)\underset{\lambda =0}{}\frac{1}{m}\widehat{\varphi }_\lambda (x)\widehat{\varphi }_\lambda ^{}(z)2R_{zy}b(y).$$
(6.16)
4. To get the inverse of $`D`$, multiply the result by $`2R`$:

$$g^{(\mathrm{sub})}(x;m)=2R\widehat{g}^{(\mathrm{sub})}(x;m).$$
(6.17)

We have to remark that in principle it would be more elegant to perform the basis transformation (6.8) in the inverse direction with $`S=2R`$, because then in the above steps 6.16 and 6.17 the factor of $`2R`$ disappears. The reason we do not follow this apparently simpler path is that in the orthonormality condition (and in hermitean forms in general, which are used in certain parts of our code), a multiplication with the inverse of $`2R`$ would show up, which is of course numerically much more demanding.

### 6.2 Zero Mode Contributions in Meson Propagators

Consider first the case of the normal Dirac operator $`\stackrel{~}{D}`$, where a spectral decomposition into a complete set of eigenstates is possible. Inserting the spectrally decomposed quark propagator (6.4) into the meson propagator (5.11), we get

$$𝒪^\mathrm{M}(x)𝒪^\mathrm{M}(0)=\underset{\lambda ,\lambda ^{}}{}\frac{\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}(x)\mathrm{\Gamma }\gamma _5\stackrel{~}{\varphi }_\lambda ^{}(x)\right]\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}^{}(0)\mathrm{\Gamma }^{}\gamma _5\stackrel{~}{\varphi }_\lambda (0)\right]}{(\lambda +m)(\lambda ^{}+m)},$$
(6.18)
where the color and spin indices have been suppressed and only the spatial index is given explicitly. Splitting the sum into three contributions and summing over $`\stackrel{}{x}`$ to project out zero momentum states leads to

$$\begin{array}{c}\underset{\stackrel{}{x}}{}𝒪^\mathrm{M}(x)𝒪^\mathrm{M}(0)=\underset{\stackrel{}{x}}{}[\underset{\lambda =\lambda ^{}=0}{}\frac{\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}(x)\mathrm{\Gamma }\gamma _5\stackrel{~}{\varphi }_\lambda ^{}(x)\right]\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}^{}(0)\mathrm{\Gamma }^{}\gamma _5\stackrel{~}{\varphi }_\lambda (0)\right]}{m^2}\hfill \\ \hfill +\underset{\lambda =0,\lambda ^{}0}{}\frac{\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}(x)\mathrm{\Gamma }\gamma _5\stackrel{~}{\varphi }_\lambda ^{}(x)\right]\mathrm{Tr}\left[\stackrel{~}{\varphi }_\lambda ^{}^{}(0)\mathrm{\Gamma }^{}\gamma _5\stackrel{~}{\varphi }_\lambda (0)\right]}{m(\lambda ^{}+m)}]+_n\frac{|0|𝒪^\mathrm{M}|n|^2}{2E_n}e^{E_n\tau }.\end{array}$$
(6.19)
The first two summands, which are the pure and mixed zero mode contributions, are divergent in the chiral limit $`m0`$. The zero mode contributions are finite volume artifacts expected to diverge with $`1/\sqrt{V}`$ , and the meson correlator diverges at finite volume with $`𝒪(1/m^2)`$ and $`𝒪(1/m)`$ terms in the chiral limit. If the zero modes in the quark propagator have been subtracted as in Eq. (6.10), these divergences are absent, and also at small quark masses a good signal for the exponential decay of the meson mass can be obtained.

Let us discuss what happens for different meson operators. For the pseudoscalar, we have $`\mathrm{\Gamma }=\mathrm{\Gamma }^{}=\gamma _5`$ with our definition of the Clifford matrices, while for the scalar meson $`\mathrm{\Gamma }=\mathrm{\Gamma }^{}=1`$. The zero modes are chiral eigenstates with $`\gamma _5\stackrel{~}{\varphi }_{\lambda =0}=\pm \stackrel{~}{\varphi }_{\lambda =0}`$, hence it follows that the pseudoscalar and the scalar correlator have the same contributions from zero modes. By taking the difference between the two correlators, the zero mode contributions are therefore exactly cancelled. For the axial vector current correlator $`\mathrm{\Gamma }=\gamma _4\gamma _5`$, the first summand in Eq. (6.19) is zero, since $`\gamma _4`$ couples different chiral sectors and all zero modes on a given gauge configuration have equal chirality. So for this correlator, the zero modes contribute only in the mixed term, and therefore it diverges only with $`1/m`$ in the chiral limit.

If we consider the non-normal Dirac operator $`D`$, we have $`G=\sqrt{2R}\stackrel{~}{G}\sqrt{2R}`$, and in the numerators of Eqs. (6.18) and (6.19) additional factors of $`\sqrt{2R}`$ show up. After subtracting the zero modes of $`\stackrel{~}{D}`$ according to (6.14)–(6.17), again the divergence in the chiral limit is no longer present.

### 6.3 Numerical Results at Small Volume

In order to examine the finite-volume zero mode effects in pion propagators, we perform an exploratory study with the overlap-improved FP Dirac operator on a very small lattice of size $`6^3\times 16`$ at gauge coupling $`\beta =3.0`$, with a set of 90 independent gauge configurations. This lattice volume amounts to a spatial extension of 1 fm, therefore the signal in the pseudoscalar channel is strongly affected by the zero modes. The Dirac operator is constructed with a third order Legendre expansion for the overlap and with exact treatment of the smallest 10 eigenmodes. Two different strategies are examined to remove the zero mode effects from the pion propagator:

The first strategy is to remove the zero modes from the pseudoscalar (P) correlator by building the difference of the pseudoscalar and the scalar (P-S) correlator as described in Section 6.2. We further calculate the zero component axial vector (A) correlator, which also has the quantum numbers of the pion, and for which the zero mode contributions are partially cancelled. We show in Fig. 6.3 the effective masses at the smallest quark mass $`am=0.01`$ and in Fig. 6.4 the squared pion mass as a function of the quark mass. Obviously the different correlators give very different results at small quark mass. While the chiral limit of the pseudoscalar correlator clearly deviates from zero, it is consistent with zero for the P-S correlator. For the axial correlator, the zero mode effects are of the same order as for the pseudoscalar. For comparison, we show in Figs. 6.5 and 6.6 the same plots for a larger volume of $`12^3\times 24`$, with unchanged gauge coupling $`\beta =3.0`$, the same Dirac operator and 100 independent gauge configurations. There the zero mode effects are much smaller due to the larger physical volume. However, the pion mass from the pseudoscalar correlator still clearly deviates from zero in the chiral limit, while for the P-S correlator it nicely goes to zero. From these results we find that removing the zero mode effects in the pion at small quark mass by using the P-S correlator works fine.

The second strategy is to explicitly calculate the zero modes and remove them from the quark propagator as shown in Section 6.1. This is only possible if the Dirac operator allows to identify zero modes unambiguously. To check whether this requirement is fulfilled for the overlap-improved FP Dirac operator, we show in Table 6.1 for a gauge configuration with $`Q=1`$ the position of the zero eigenvalue on the real axis and the first non-zero eigenvalue as a function of the overlap order. We also list the chirality of the corresponding eigenvector $`\varphi ^{(i)}`$, defined by

$$\chi ^{(i)}=\varphi ^{(i)}\gamma _52R\varphi ^{(i)},$$
(6.20)
which is $`\pm 1`$ for zero modes. Obviously, it is easily possible to separate the zero mode and the first non-zero mode, even if the eigenvalue $`\lambda _0`$ is not very close to zero, because the chirality of the zero mode approaches $`|\chi _0|1`$ very quickly with increasing order of the overlap expansion.

We therefore calculate on every gauge configuration the eigenvectors $`\varphi ^{(i)}`$, $`i=1,\mathrm{},10`$ of $`\widehat{D}`$ corresponding to the 10 smallest eigenvalues and computed their chirality, and the eigenvectors with chirality $`|\chi ^{(i)}|0.9`$ are considered to be zero modes. The distribution of the topological charge for the set of 90 gauge configurations is peaked at $`Q=1`$, as can be seen in Table 6.2. All zero modes are subtracted from the quark propagator according to (6.10), and with the subtracted quark propagators, pion propagators are constructed in the usual way. In the pseudoscalar correlator, the zero mode effects should then disappear, while the P-S correlator should remain unchanged.<sup>2</sup><sup>2</sup>2Indeed our P-S correlators change only marginally when built from zero mode subtracted quark propagators. First we examine the correlators at the smallest quark mass $`am=0.01`$ in Fig. 6.7. The subtraction of the zero modes leads both to a strong decrease in the amplitude and the mass, as can be seen from the fitted lines. This effect is illustrated by the Monte Carlo time evolution of the pseudoscalar correlators in Fig. 6.8. At small time $`t=1`$, the correlator fluctuates wildly with Monte Carlo time. The zero mode subtraction removes most of the peaks at topologically non-trivial gauge configurations. The average over Monte Carlo time is therefore much smaller for the subtracted correlator. At large time $`t=7`$, the picture changes and more peaks remain after zero mode subtraction. Some of the most prominent peaks (no. 23 and 42) even are at $`Q=0`$ gauge configurations and therefore do not come from zero modes. The average correlator decreases only moderately after subtraction. Combining the observations at the two times, it follows that the pion mass gets smaller at $`am=0.01`$ when determined from the subtracted propagators.

The situation is different at very large quark mass $`am=0.32`$. The subtracted and unsubtracted pseudoscalar correlators in Fig. 6.9 are almost equivalent at small $`t`$, but clearly disagree at larger $`t`$. Again, the Monte Carlo history in Fig. 6.10 helps illuminating this observation. At $`t=1`$, there is essentially no effect from the zero mode subtraction. At $`t=7`$, after subtraction there appear some large peaks at topologically non-trivial gauge configurations which increase the average considerably, while the full correlator behaves quite smoothly. There are several explanations for this strange behavior of the meson correlator at large mass and time. First of all, the third-order overlap-improved Dirac operator is not exactly chiral, therefore the subtraction of the approximate zero modes leads to small numerical deviations from the exact zero-mode subtracted quark propagators. These deviations become important at large mass and large $`t`$, where the meson correlator is small, and might thus cause the observed flattening in Fig. 6.9. Second, removing the zero modes amounts to a modification of the quenched theory. The meson correlators then do not necessarily have to be a sum of exponential functions. To rule out the first possibility, it would be necessary to repeat this examination with a larger order $`𝒪(10)`$ of the overlap expansion in the Dirac operator. We are however not mainly interested in the large mass behavior of the meson propagator, where a reliable pion mass can easily be extracted from the P correlator. Therefore we do not investigate this effect further.

The results of this small volume study are summarized in Fig. 6.11, where the chiral limit of the squared pseudoscalar meson mass is plotted for the P correlator on $`Q=0`$ configurations, the zero mode subtracted pseudoscalar correlator P<sub>sub</sub> and the P-S correlator. All of them agree within errors at small quark masses and go to zero for $`m0`$.

### 6.4 Conclusion

At small quark masses and fixed volume, the pseudoscalar meson masses measured from P correlators are distorted by topological finite volume effects. Forming the difference P-S, the chiral limit of $`(am_{\mathrm{PS}})^2`$ is consistent with the expectations, confirming that the observed distortion is due to the topological finite size effects. Indeed, in the P-S correlators these effects cancel, up to small chiral symmetry breaking contributions. Furthermore, P-S is a sum over exponentials with physical meson masses, although both the scalar and the pseudoscalar mesons enter. For small quark masses however, the pion dominates.

Unlike P-S, the pseudoscalar correlator P<sub>sub</sub> built from zero mode subtracted quark propagators is a strange quantity and does not need to be a sum of exponential functions. It is therefore better not to use P<sub>sub</sub> in actual calculations.

At intermediate quark masses, the P and P-S correlators agree, as expected. To extract the pseudoscalar meson mass, we therefore use the P-S correlator at small quark masses, where the P correlator would be contaminated by the zero modes, and the P correlator at large quark masses, where it would be difficult to disentangle the contributions of the scalar meson to the P-S correlator. This provides a reliable determination of the pseudoscalar mass over the whole range of quark masses, and as we will also see in Chapter 7, an intersection of $`(am_{\mathrm{PS}})^2`$ with the horizontal axis which is consistent with other determinations.

The best way to avoid any problems with zero mode effects is to work at large enough lattice volumes. As will be shown in Section 7.2, at our largest lattice size $`L_s2.5`$ fm, the zero modes no longer contaminate the pion propagator significantly, and it is possible to get unambiguous answers concerning the chiral limit of pseudoscalar mesons.

## Chapter 7 The Light Hadron Spectrum with Fixed-Point Fermions

Since the first attempts in 1981 , many lattice studies of the light hadron spectrum in quenched QCD have been performed, with quality increasing with time. The first systematic calculation was done in 1993 by the GF11 collaboration , but today’s benchmark is the CP-PACS calculation  from 1998, which included a thorough examination of the chiral and continuum extrapolations and very high statistics. In their study, the most simple choice of actions was taken, namely the Wilson plaquette and fermion actions, and a full year of runs on the dedicated CP-PACS computer with a peak performance of 614 GFLOPS was necessary in order to obtain the quenched particle spectrum in a controlled manner. Because the cut-off effects for the Wilson action are known to be large, the simulation was performed at rather small lattice spacings in the range $`a0.05`$$`0.1`$ fm. To avoid finite volume effects, the physical size was chosen to be 3 fm, which required to run on lattices of sizes up to $`64^3\times 112`$. The Monte Carlo average was taken from up to 800 independent gauge configurations. Their result is a physical particle spectrum with very small statistical errors, which are on the order of 1–2% for mesons and 2–3% for baryons, and with systematic errors from the extrapolations that are estimated to be even smaller. The calculated hadron masses agree qualitatively with the experimentally observed spectrum, but the mass values deviate by up to 11% or $`7\sigma `$, which is thought to be the error introduced by the quenched approximation.

A crucial part in the analysis of the CP-PACS data was the chiral extrapolation. The quark mass was pushed down to a value corresponding to $`m_{\mathrm{PS}}/m_\mathrm{V}0.4`$, which is very small for the Wilson action with its inherent chiral symmetry breaking<sup>1</sup><sup>1</sup>1Typically, simulations with Wilson fermions do not go lower than $`m_{\mathrm{PS}}/m_\mathrm{V}0.5`$, and only with the point at the lowest quark mass it was possible to resolve the non-analytic contributions predicted by quenched chiral perturbation theory (Q$`\chi `$PT) . It is however important to investigate whether including such Q$`\chi `$PT terms leads to significantly different mass values in the chiral limit than when using just low-order polynomial forms, as usually done in earlier works.

Exact or approximately chiral fermion actions allow to go to much smaller quark masses than the Wilson action. It is therefore possible to explore the chiral limit directly, to check the significance of the Q$`\chi `$PT terms and thus to increase the reliability of the hadron mass extrapolation to the physical quark mass. In this Chapter, we report the results of a spectroscopy calculation with the parametrized FP fermion action which probes deep into the chiral limit and also includes investigations of the scaling properties of the hadron masses and their finite-volume dependence. With a total amount of computer resources of about 20 GFLOPS$`\times `$years theoretical peak and an effective amount of about 6 GFLOPS$`\times `$years, this study in the framework of the BGR collaboration  by no means attempts to compete with the above mentioned high-statistics calculations. However, in addition to an independent test for spectroscopy simulations with a new formulation of the lattice QCD action, we get important information from the region of small quark masses, where non-chiral actions do not allow to perform simulations. Furthermore, this is one of the first spectroscopy studies with a chiral symmetric fermion action that examines cut-off effects.

Before we started this simulation, we carried out some tests for spectroscopy with the FP action on a smaller scale, the results of which are published in . The lattice parameters of these tests are listed for completeness in Table 7.1.

### 7.1 Simulation Parameters

For generating the gauge configurations, we use the parametrized isotropic Fixed-Point gluon action from . While this action is relatively expensive compared to the Wilson action<sup>2</sup><sup>2</sup>2A factor of $`60`$ in computer time is estimated in ., generating the gauge configurations was a comparably small effort in terms of computer time in the context of this work, and therefore we could afford using a FP action also for the gluon sector. It was shown in  that the FP gauge action has small scaling violations in gluonic quantities and that it reproduces topological properties well. With this gauge action, we produced a set of configurations at various lattice sizes and gauge couplings as shown in Table 7.2, where we also list the lattice spacing determined from the Sommer scale $`r_00.5`$ fm for the different values of the gauge coupling<sup>3</sup><sup>3</sup>3The detailed analysis of the scale determination for the FP gluon action is given in . This combination of parameters was chosen to allow for a scaling analysis at small physical spatial lattice size $`L_s1.3`$ fm and a finite-volume analysis at gauge coupling $`\beta =3.0`$. The largest lattice at $`\beta =3.0`$ has a physical volume large enough to accommodate hadrons with negligible finite-volume effects. While this value of the gauge coupling is quite far away from the continuum, we expect to get on this lattice precise numbers for hadron masses which can serve as good estimates for the continuum values, because as we will show, the FP action has small scaling violations. We use alternating Metropolis and pseudo-overrelaxation sweeps over the lattice, with 2000 sweeps for thermalization and 500 sweeps to separate between different configurations. The number of separation sweeps is a worst-case estimate based on measurements of autocorrelation times for simple gluonic operators . The configurations are then smeared with the RG inspired two-level hypercubic smearing described in  and fixed to Coulomb gauge<sup>4</sup><sup>4</sup>4The order of smearing and gauge fixing is a matter of choice. In earlier studies we first fixed the gauge and then smeared the links, as we considered the smearing to be part of the definition of the Dirac operator. For this spectroscopy study, the order was reversed. The argument was that because after smearing the configurations are much smoother, the gauge fixing algorithm might have less problems with Gribov copies. with the algorithm presented in Appendix A. We use periodic boundary conditions for both the gauge and fermion fields.

For the fermion action, we take the parametrized FP Dirac operator from Chapter 3, except for one lattice, where the overlap-improved FP Dirac operator from Chapter 4 is used. With the third-order overlap expansion, we decrease the small residual chiral symmetry breaking of the parametrization even further and are able to check what effect the overlap construction has on the mass spectrum. The quark masses cover a very large range, with the smallest value, where $`m_{\mathrm{PS}}/m_\mathrm{V}0.21`$, lying close to the physical point. This provides us with meaningful data for chiral fits and allows for a thorough examination of the chiral limit. To enhance the signal for the hadron correlators, we use Gaussian smeared sources located at the center of timeslice $`t=0`$. The source extension parameter $`\gamma `$ in Eq. (5.26) is chosen to correspond to a source size of $`0.5`$ fm.<sup>5</sup><sup>5</sup>5On the $`12^3\times 24`$ lattice at $`\beta =3.4`$, the source size is $`0.3`$ fm. We use point sinks and project to zero momentum by summing over all spatial sink locations. Quark propagators are calculated with the multi-mass BiCGstab inversion algorithm  (see also Appendix B.3). Source vectors are normalized to 1, and as a stopping criterion, we require the residual to be smaller than some threshold $`|r|\sigma `$.

Due to the fact that a multi-mass solver inverts the Dirac matrix at all quark masses simultaneously, the result of the inversion at larger masses is more accurate than at the smallest one, where the condition number of the Dirac matrix is worst. The value of $`\sigma `$ therefore determines the accuracy of the quark propagators at the smallest quark mass, while the propagators at larger masses are calculated to much higher precision, as shown in Fig. 7.1. In order to demonstrate that the error from the truncation of the inversion algorithm is much smaller than the statistical error from the Monte Carlo estimate of the hadron correlators, we show in Table 7.3 the dependence of pseudoscalar and vector meson correlators on the stopping criterion. Based on this data, we choose a precision of $`\sigma =10^6`$ for the inversion of the Dirac operator, which leads to an error from the truncation of the iterative inversion algorithm that is negligible compared to the statistical error. At our intermediate and larger quark masses, the quark propagator is almost calculated to machine precision. The average number of iterations needed for the inversion on the various lattices is also given in Table 7.2.<sup>6</sup><sup>6</sup>6Interpreting these numbers, one has to keep in mind that the BiCGStab algorithm requires two matrix-vector multiplications with the Dirac operator per iteration.

To check whether the computed hadron propagators are statistically independent, we calculate the statistical error of effective meson masses in dependence of the number of configurations $`N`$ used for the Monte Carlo average. If the gauge configurations are independent, the bootstrap error is proportional to $`1/\sqrt{N}`$. In Fig. 7.2 we plot the error in effective pseudoscalar meson masses at various values of $`t`$ on the $`12^3\times 24`$ lattice at $`\beta =3.4`$. As we can see, the curves nicely agree with the expected $`1/\sqrt{N}`$ behavior. For a second check, we collect bins of $`N_{\mathrm{bin}}`$ successive propagators, take the average of each bin and calculate the statistical error in effective masses from the $`N/N_{\mathrm{bin}}`$ binned propagators. The resulting bootstrap error, plotted in Fig 7.3, turns out to be independent of the bin size $`N_{\mathrm{bin}}`$. These two observations confirm that the number of separation sweeps used in the generation of the gauge configurations is sufficient to ensure statistical independence.

We list our hadron mass results on the different lattices for the pseudoscalar and vector mesons and the octet and decuplet baryons from quarks with degenerate masses, together with the temporal fit range and the resulting value of $`\chi ^2/df`$ for the fit, in Appendix D.1. To account for a possibly biased determination of the fitted masses we apply bootstrap bias correction, as defined in Eq. (5.40), to the fit results.<sup>7</sup><sup>7</sup>7 Bias correction has proven relevant only at the very smallest quark masses, where in some cases the fitted hadron masses showed a noticeable bias. At intermediate and large quark masses, the bias is negligible, as can be seen from the data in Appendix D.1.

Fig. 7.4 gives an overview of the masses of the hadrons in dependence of the quark mass on the various lattices in our simulation. We remind the reader that due to the use of a multi-mass solver, the results at different quark masses are highly correlated. This fact has to be taken into account when interpreting the data presented in this chapter.

The following analysis of our spectroscopy runs has to be considered preliminary. We try to give in this work some first answers to the main questions that arise from the use of such a new chiral action for the determination of the hadron spectrum.

### 7.2 Zero Mode Effects

In Chapter 6, we showed that divergent zero mode contributions appear in pion propagators and examined these quenching artifacts and possible ways to remove them on a small lattice volume. We concluded that it seems sensible to work with different correlators at different quark masses. We return to this issue here and present the zero modes effects seen in the data of our spectroscopy simulation.

In Figs. 7.57.10 the Monte Carlo history of the three correlators P, A, and P-S and the resulting pseudoscalar meson masses are shown on the three different lattice volumes at gauge coupling $`\beta =3.0`$. The values of the correlators in the Monte Carlo time plots are taken from timeslice $`t=8`$ and quark mass $`am_q=0.021`$. As expected, in the case of the smallest lattice of size $`8^3\times 24`$ in Fig. 7.5, there are a few very prominent peaks in the P correlator which dominate the gauge average completely. Also the A correlator is dominated by these peaks, but to a somewhat smaller extent. In contrast to these heavily contaminated correlators, many more configurations contribute to the average of the P-S correlator. Considering the absolute scale, it is obvious that the zero mode contributions, which are present in both the pseudoscalar and scalar correlators, are cancelled to a large extent in the difference. The resulting pseudoscalar meson masses from these three correlators, plotted in Fig. 7.6, differ significantly at the smallest quark masses.

On the lattice with intermediate size $`12^3\times 24`$, the same effects in the Monte Carlo history of the correlators can be observed in Fig 7.7, but they are much less pronounced. Again, the largest peaks in the P correlator are cancelled in the P-S correlator. The masses in Fig. 7.8 agree much better, but still a systematic deviation is visible. On the largest lattice (Fig. 7.9), all three correlators show a fairly smooth behavior. The mass of the pseudoscalar meson (Fig. 7.10) depends only very little on the choice of correlator. We summarize this analysis by listing the relative difference in meson masses from the P and P-S correlators in dependence of the lattice volume and quark mass in Table 7.4, a quantity which serves as a crude estimate for the size of the zero mode effects.

We conclude from these observations that on our largest lattice with spatial extension $`L_s2.5`$ fm, the zero mode contributions are sufficiently small to allow for a reliable determination of the mass in the pseudoscalar channel and the quantitative examination of its chiral limit. Only at the two smallest quark masses, the results from the different correlators do not agree within the (already small) statistical errors. On the smaller lattices, the uncertainty grows rapidly with decreasing quark mass, and therefore it becomes increasingly hard to keep the zero mode effects under control. We will in the following adopt the strategy suggested in Chapter 6 and construct the pseudoscalar meson from the P-S correlator at small and the P correlator at large quark mass. At the smaller lattice volumes, this introduces a systematic uncertainty at the few smallest quark masses, where we can not control how well the subtraction in the meson correlators works. But as our results will show, this strategy turns out to be successful in our simulations.

### 7.3 Chiral Extrapolations and Quenched Chiral Logarithms

With the approximately chiral Fixed-Point fermion action, we are in a position to study the chiral limit of hadron masses in detail. After measuring the residual additive quark mass renormalization, we examine in this section the quark mass dependence of the pseudoscalar meson mass and check for the presence of quenched chiral logarithms. On our largest lattice, we also calculate the resulting spectrum of $`K`$, $`K^{}`$ and $`\varphi `$ mesons and $`N`$, $`\mathrm{\Delta }`$ and $`\mathrm{\Omega }`$ baryons after extrapolating the measured data for the degenerate hadrons to the physical values of the average up and down quark and the strange quark mass, respectively.

#### 7.3.1 Residual Quark Mass

We first determine the residual additive quark mass renormalization introduced by the parametrization of the Fixed-Point Dirac operator. From the axial Ward identity (AWI), the quark mass $`m_q^{\mathrm{AWI}}`$ is given by the large $`t`$ limit of

$$Z_mm_q^{\mathrm{AWI}}(t)=\frac{1}{2}\frac{Z_A_\stackrel{}{x}_4A_4(\stackrel{}{x},t)C(0)}{Z_P_\stackrel{}{x}P(\stackrel{}{x},t)C(0)},$$
(7.1)
where $`C(0)`$ is a source operator with the quantum numbers of the pion and $`A_4(\stackrel{}{x},t)`$ and $`P(\stackrel{}{x},t)`$ are the local fourth component axial vector and pseudoscalar currents. As seen in Fig. 7.11, the ratio of the two correlators in Eq. (7.1) is flat already at small $`t`$ and can easily be fitted to a constant also at very small quark masses. The unrenormalized AWI quark masses $`m_q^{\mathrm{AWI}}`$ extracted from our data are listed in Appendix D.2. Although we do not know the renormalization factors $`Z_m`$, $`Z_S`$ and $`Z_A`$,<sup>8</sup><sup>8</sup>8For exactly chiral actions, $`Z_m^1=Z_S`$, hence $`Z_A`$ can be extracted from Eq. (7.1). the residual quark mass can be determined from the value of the bare quark mass $`m_q`$ where the AWI quark mass vanishes, by linearly extrapolating the measured ratio of correlators to zero.

The resulting values of $`m_{\mathrm{res}}^{(\mathrm{AWI})}=m_q(m_q^{\mathrm{AWI}}=0)`$ from a linear fit to the six smallest masses on each lattice are shown in Table 7.5. At $`\beta =3.0`$, the residual mass for both the parametrized and the overlap-improved FP Dirac operator is very close to zero, whereas at larger $`\beta `$, its value is clearly non-zero. This is not surprising, since we have optimized the parametrization to $`\beta 3.0`$ by choosing a set of gauge configurations with $`2.7<\beta <3.4`$ to fit the FP relation. While the fluctuations of the smallest eigenvalues of the Dirac operator decrease with increasing $`\beta `$, the central value is slightly shifted away from zero, leading to a non-vanishing value of $`m_{\mathrm{res}}^{(\mathrm{AWI})}`$ at larger $`\beta `$. Obviously the fluctuation polynomials in the parametrized FP Dirac operator do not fully absorb this $`\beta `$-dependence of the eigenvalue spectrum.

#### 7.3.2 The Quenched Chiral Log Parameter $`\delta `$

Chiral perturbation theory ($`\chi `$PT)  allows to predict the quark mass dependence of hadrons. In particular the pion, taking the role of the Goldstone boson of spontaneously broken chiral symmetry, should be massless in the chiral limit, and to lowest order in $`\chi `$PT, it depends linearly on the quark mass. Taking into account the quenched approximation, additional logarithmic terms appear: Quenched chiral perturbation theory predicts the following dependence of the pseudoscalar meson mass $`m_{\mathrm{PS}}`$ on the two constituent quark masses $`m_1`$ and $`m_2`$ :

$$\begin{array}{c}m_{\mathrm{PS}}^2=A(m_1+m_2)\left[1\delta \left(\mathrm{ln}\frac{2Am_1}{\mathrm{\Lambda }_\chi ^2}+\frac{m_2}{m_2m_1}\mathrm{ln}\frac{m_2}{m_1}\right)\right]+B(m_1+m_2)^2\hfill \\ \hfill +𝒪(m_q^3),\end{array}$$
(7.2)
where $`\mathrm{\Lambda }_\chi `$ is a Q$`\chi `$PT scale and of order 1 GeV. The term proportional to $`\delta `$ is divergent in the chiral limit and is only present due to quenching. In the case of degenerate quark masses, the divergence can be absorbed into a redefinition of the quark mass by resummation of the leading logarithms . This yields a power form

$$m_{\mathrm{PS}}^2m_q^{\frac{1}{1+\delta }},$$
(7.3)
for the pseudoscalar mass. To leading order in a $`1/N_c`$-expansion, the value of the parameter $`\delta `$ is given by

$$\delta =\frac{m_\eta ^{}^2+m_\eta ^22m_K^2}{48\pi ^2f_\pi ^2}0.18,$$
(7.4)
where the pion decay constant is normalized such that $`f_\pi =93`$ MeV.

CP-PACS has estimated a value of $`\delta =0.10(2)`$ from their analysis of pseudoscalar mesons with non-degenerate quarks and $`\delta =0.09(3)`$ from a fit of Eq. (7.2) to pseudoscalar mesons with degenerate quarks, where the scale was varied in the range $`0.6`$ GeV$`<\mathrm{\Lambda }_\chi <1.4`$ GeV. The drawback of their data is the use of the Wilson action with its explicit breaking of chiral symmetry, which allows to cover only a limited range of quark masses because of the appearance of exceptional configurations. Recently a technique has been proposed for the Wilson action to shift the would-be zero modes, that fluctuate far along the positive real axis, back to zero . This pole-shifting procedure amounts to a modification of the quenched theory and prevents exceptional configurations, thus allowing to go almost down to the physical pion mass even with Wilson fermions. However, one has to assume that the most important effect of explicit chiral symmetry breaking in the Wilson action is the resulting fluctuation of the zero modes on the real axis. The chiral properties are not improved for other parts of the spectrum, and therefore it is not clear how much such a punctual modification helps. In , the chiral log parameter $`\delta `$ has been determined in various ways, amongst others also from the pseudoscalar meson mass, with the pole-shifted Wilson action and similar simulation parameters like we use here. The quoted value averaged over the different determinations is $`\delta =0.065(13)`$, which is a factor of three smaller than the theoretical estimate, but fairly consistent with the CP-PACS result.

It would obviously be interesting to determine the value of $`\delta `$ with a chiral symmetric action, which is free of problems related to exceptional configurations and does not only cure the explicit chiral symmetry breaking punctually. In , investigations of the chiral limit with the Wilson overlap action were performed, and varying values of $`\delta `$ were found, with a tendency towards larger values than what was obtained with non-chiral actions. With our data, we are not only able to quite precisely determine the value of $`\delta `$ from the pseudoscalar meson mass, but we can also compare the results for the approximately chiral parametrized FP and the overlap-improved Dirac operator, providing two at least partially independent results with chiral actions.

First we examine the mass of the pseudoscalar meson with degenerate quarks, calculated from the P and P-S correlators on various lattices at $`\beta =3.0`$. We fit the masses to the polynomial forms

$`(am_{\mathrm{PS}})^2`$ $`=`$ $`2Aa^2(m_q+m_{\mathrm{res}})+4Ba^3(m_q+m_{\mathrm{res}})^2,`$ (7.5)
$`(am_{\mathrm{PS}})^2`$ $`=`$ $`2Aa^2(m_q+m_{\mathrm{res}})+4Ba^3(m_q+m_{\mathrm{res}})^2`$ (7.6)
$`+8Ca^4(m_q+m_{\mathrm{res}})^3,`$
and to the forms inspired by quenched chiral perturbation theory

$`(am_{\mathrm{PS}})^2`$ $`=`$ $`2Aa^2(m_q+m_{\mathrm{res}})\left[1\delta \left(\mathrm{ln}{\displaystyle \frac{2Aa^2(m_q+m_{\mathrm{res}})}{a^2\mathrm{\Lambda }_\chi ^2}}+1\right)\right]`$ (7.7)
$`+4Ba^3(m_q+m_{\mathrm{res}})^2,`$
$`(am_{\mathrm{PS}})^2`$ $`=`$ $`2Aa^2(m_q+m_{\mathrm{res}})^{\frac{1}{1+\delta }},`$ (7.8)
$`(am_{\mathrm{PS}})^2`$ $`=`$ $`2Aa^2(m_q+m_{\mathrm{res}})^{\frac{1}{1+\delta }}+4Ba^3(m_q+m_{\mathrm{res}})^2,`$ (7.9)
allowing for a residual additive quark mass renormalization $`am_{\mathrm{res}}`$. The resulting fit parameters for all three lattices are listed in Table 7.6. The errors are determined by bootstrap resampling, calculating on each of the 500 bootstrap samples the meson masses from a correlated fit and using these for the chiral fit<sup>9</sup><sup>9</sup>9We do not take into account correlations at different quark masses here..

The most reliable results are obtained on the largest lattice. There it is obvious that a quadratic fit with a value of $`\chi ^2/df=11.5`$ misses the clear negative curvature at small quark masses in our data. Also a cubic fit with $`\chi ^2/df=2.9`$ can not account for this curvature well. In contrast, the various Q$`\chi `$PT fits agree perfectly with our data, with $`\chi ^2/df0.5`$. In Fig. 7.12 we compare the logarithmic fit with $`\mathrm{\Lambda }_\chi =1`$ GeV and the quadratic fit, showing the clear discrepancy when the quenched chiral log is neglected. In Table 7.7 we demonstrate that the results for $`\delta `$ on the largest lattice do not differ significantly when choosing one single correlator P, A or P-S instead of extracting the pseudoscalar masses from the P and P-S correlators at different quark mass according to our proposition.

On the $`12\times 24`$ lattice, the results for $`\delta `$ with the parametrized FP Dirac operator agree very well with those at $`16^3\times 32`$, with somewhat larger errors. Also with the overlap-improved operator, we get consistent values, but the errors are then even larger due to the smaller statistics. We compare our values of $`\delta `$ from the degenerate mesons with those from  in Table 7.8. From the results of the fits to the logarithmic form (7.7) on the largest lattice, we estimate a value of $`\delta =0.17(3)`$, where the error mostly comes from the unknown scale $`\mathrm{\Lambda }_\chi `$.

The quenched chiral log parameter $`\delta `$ can also be determined from an analysis of the non-degenerate pseudoscalar meson masses. The two quantities

$`x=2+{\displaystyle \frac{m_1+m_2}{m_1m_2}}\mathrm{ln}\left({\displaystyle \frac{m_2}{m_1}}\right),`$ (7.10)
$`y={\displaystyle \frac{2m_1}{m_1+m_2}}{\displaystyle \frac{m_{\mathrm{PS},12}^2}{m_{\mathrm{PS},11}^2}}{\displaystyle \frac{2m_2}{m_1+m_2}}{\displaystyle \frac{m_{\mathrm{PS},12}^2}{m_{\mathrm{PS},22}^2}},`$ (7.11)
are related by $`y=1+\delta x+𝒪(m^2)`$, and the terms proportional to $`B(m_1+m_2)`$ in Eq. (7.2) cancel . To avoid problems with the residual additive renormalization, we use the axial Ward identity quark masses (7.1) as an input for $`m_1`$ and $`m_2`$. Fig. 7.13 shows our data for the non-degenerate pseudoscalar mesons from the P correlator on the largest lattice. Only the points from mesons with light enough quarks that the linear dependence of $`y(x)`$ remains valid are plotted. For mesons with two heavy quarks, we encountered a systematic deviation towards smaller values of $`y`$ for all $`x`$, therefore they are omitted here. The slope of $`y(x)`$, giving $`\delta =0.170(20)`$, is beautifully consistent with the above value for the equal quark mass case and with the theoretical prediction in Eq. (7.4).

Taking the intersection of the two determinations of the chiral log parameter from degenerate and non-degenerate pseudoscalar mesons, we obtain a value of $`\delta =0.17(2)`$. In Table 7.9, we compare this final result to previous determinations from other groups. While our result is measured at finite lattice spacing $`a0.16`$ fm, the $`a`$-dependence of $`\delta `$ appears to be small .

As a by-product of this analysis, we estimate a value of $`am_{\mathrm{res}}^{(\mathrm{PS})}=0.002(2)`$ for the residual quark mass of the parametrized FP Dirac operator, covering the various quenched fits to the pseudoscalar meson mass at $`\beta =3.0`$. This value agrees well with the residual mass determined from the axial Ward identity in Table 7.5, which is clearly not the case for $`am_{\mathrm{res}}`$ resulting from quadratic or cubic fits to $`(am_{\mathrm{PS}})^2`$. At the other values of the gauge coupling $`\beta =3.4`$ and $`3.7`$, we get larger values of $`am_{\mathrm{res}}^{(\mathrm{PS})}=0.016(4)`$ and $`0.012(7)`$ respectively, in reasonable agreement with the results from the axial Ward identity. This coincidence confirms that measuring the pseudoscalar particle from the P correlator at large and the P-S correlator at small quark mass is a reasonable solution to avoid zero mode effects in the pseudoscalar mass.

#### 7.3.3 Chiral Extrapolations for Vector Mesons and Baryons

For vector mesons and baryons, Q$`\chi `$PT predicts in the continuum limit

$$m(m_{\mathrm{PS}})=m_0+\delta C_{1/2}m_{\mathrm{PS}}+C_1m_{\mathrm{PS}}^2+C_{3/2}m_{\mathrm{PS}}^3+\mathrm{},$$
(7.12)
where the $`C_i`$ are functions of the coefficients in the quenched chiral Lagrangian. Like in the formula for the pseudoscalar channel (7.2), the term proportional to $`\delta `$ appears only in the quenched theory, and according to Q$`\chi `$PT its coefficient $`C_{1/2}`$ is negative.

On our largest lattice, we perform chiral extrapolations for the vector meson and the octet and decuplet baryons both with and without the quenched term. The resulting fits are plotted in Fig 7.14. For the vector meson, there is some evidence for the presence of a term linear in $`m_{\mathrm{PS}}`$ with a negative coefficient, as expected from quenched chiral perturbation theory. The $`\rho `$ meson mass from the extrapolation to the physical quark mass, defined by $`m_{\mathrm{PS}}/m_\mathrm{V}0.18`$, is used to fix the lattice scale. For the fit including the quenched term, we get $`a^1=1104(32)`$ MeV or $`a=0.179(4)`$ fm, while when setting $`C_{1/2}=0`$, we obtain $`a^1=1152(15)`$ MeV or $`a=0.171(2)`$ fm. The two values do not completely coincide, showing that the functional form of the chiral extrapolation indeed can lead to different results. We use in the following the value obtained from the quenched fit, which includes in its error some of the uncertainty in the presence and size of the quenched term.<sup>10</sup><sup>10</sup>10The lattice spacing fixed from the $`\rho `$ meson mass turns out to be somewhat larger than when fixed from $`r_0`$, where we obtain $`a=0.16(1)`$ fm. This uncertainty in the scale determination is a well-known problem in quenched QCD.

For the baryons, neither of the two functional forms is clearly favored by our data. The decuplet mass shows some upward curvature at the smallest quark masses, but the errors are quite large and do not cover the systematic uncertainty in choosing the fit range, which is increasingly difficult at small quark mass. What is however evident is the negative curvature in both baryon masses that gets absorbed by the $`𝒪(m_{\mathrm{PS}}^3)`$ term, which is also present in ordinary chiral perturbation theory.

From a partial analysis of our data, we show in Fig. 7.22 the mass spectrum of hadrons with degenerate light quarks at $`\beta =3.0`$ and lattice size $`16^3\times 32`$. The familiar ambiguity from fixing the strange quark mass either with the $`K`$ or the $`\varphi `$ meson is evident, and for both choices the meson hyperfine splitting turns out to be too small. The chiral extrapolation for the baryons, where the quenched term is included, leads to a nucleon mass which almost agrees with the experimental value, while the $`\mathrm{\Delta }`$ and $`\mathrm{\Omega }`$ baryon masses come out too small. However, since we are at finite lattice spacing, this discrepancy will not necessarily persist after a continuum extrapolation.

### 7.4 Physical Finite Size Effects

When the spatial lattice size is too small to accommodate the wave function of a hadron, its mass suffers from finite-size corrections. Results from quenched simulations show that for lattices larger than 2 fm, the finite volume effects are smaller than 2% . As only our largest lattice with size $`L_s2.5`$ fm fulfills this requirement, we expect to see significant finite volume effects in our data on the smaller lattices. Plotting the ratio of octet baryon to vector meson masses in an Edinburgh plot (Fig. 7.15), we see that indeed on the lattices with spatial size $`L_s1.3`$ fm, this ratio stays more or less constant at $`m_{\mathrm{Oct}}/m_\mathrm{V}1.5`$ over the whole range of quark masses, thus being far too large in the chiral limit. However, already at lattice size $`L_s1.9`$ fm, the results do not differ anymore from those on the largest lattice within statistical errors.

We investigate the finite volume effects more closely in Fig. 7.16 by plotting the masses of all hadrons against the spatial lattice size. For the pseudoscalar meson, only the zero mode effects can be seen, leading to different results for different correlators. Also for the vector meson, there are no obvious finite-volume effects. The situation is different for the baryons: On the smallest lattice, the octet mass increases strongly, independent of the choice of correlator. The same happens with the decuplet mass, where the effect is slightly smaller.

We conclude from our data that at $`\beta =3.0`$, the finite-volume effects become comparable to our statistical errors already at $`N_s=12`$, corresponding to a spatial size of $`L_s1.9`$ fm. At $`L_s1.3`$ fm, the baryons are strongly affected by the small physical volume. The size of the largest lattice $`L_s2.5`$ fm is big enough to provide results which are not affected by finite-volume effects.

### 7.5 Scaling Properties

Let us turn to the investigation of the scaling behavior of the parametrized and overlap-improved FP fermion actions. A standard test for scaling violations of a given lattice action is to plot the vector meson mass, which is known to be particularly sensitive to cut-off effects, against the lattice spacing. Most groups use in this context the string tension to fix the scale. Except at $`\beta =3.4`$, we do not have a direct calculation of the string tension for our gauge action. Therefore we use the determination of $`r_0/a`$ from  instead, where the interpolating formula<sup>11</sup><sup>11</sup>11At the value of $`\beta =3.7`$, an extrapolation is needed, because the formula was fitted to the range $`2.36\beta 3.40`$. We took this into account by applying increasing errors of 0.5%, 0.6% and 1% to the value of $`r_0/a`$ at increasing $`\beta `$.

$$\mathrm{ln}(a/r_0)=1.1539(18)1.0932(68)(\beta 3)+0.132(11)(\beta 3)^2,$$
(7.13)
is given and $`r_0/\sqrt{\sigma }`$ is measured, to set the scale from the string tension. Table 7.10 shows our measurements of the vector meson mass interpolated to<sup>12</sup><sup>12</sup>12The interpolation of the hadron masses is done by fitting the Q$`\chi `$PT formulae (7.7) and (7.12) to the measured masses, with bootstrap resampling to determine the errors. $`m_{\mathrm{PS}}/m_\mathrm{V}=0.7`$ and the resulting value of $`m_V/\sqrt{\sigma }`$.

Our data is plotted in comparison with results from other fermion actions in Fig. 7.17. Wilson fermions have large $`𝒪(a)`$ and unimproved staggered fermions large $`𝒪(a^2)`$ effects, which are clearly seen in the scaling of the vector meson mass. Fat links do not help to improve the situation. For the various types of clover actions shown, only small $`𝒪(a^2)`$ effects remain, as it also seems to be the case for the FP action. While the parametrization of the FP action introduces cut-off effects to all orders, we do not see evidence for $`𝒪(a)`$ effects here.

Because the conversion to the string tension introduces an additional error, it is not clear from this plot whether the scaling violations of the FP action are significant. We therefore determine the hadron masses in units of $`r_0`$ and plot them in Fig. 7.18. For all hadrons under consideration, the point at $`\beta =3.0`$ does not coincide within errors with the data at smaller lattice spacings. One has to take into account that the statistical error in $`a/r_0`$ does not fully cover the systematic uncertainty in determining the lattice scale at this $`\beta `$. Furthermore a scaling study in such a small physical volume encompasses the danger of substantial contributions from small differences in the lattice volume at the various values of the gauge coupling, because the volume dependence of the hadron masses at $`L_s1.3`$ is on the onset of getting strong . Whether our data indicates the presence of small $`𝒪(a^2)`$ cut-off effects in our parametrization of the FP fermion action or results from the determination of the Sommer scale $`r_0`$ and the related uncertainty in the lattice volume therefore needs further investigation.

A striking observation is that the vector meson mass turns out to be significantly larger for the overlap-improved Dirac operator than for the FP operator. We illustrate this by plotting $`am_\mathrm{V}`$ against $`am_{\mathrm{PS}}`$ in Fig. 7.19. Over the whole range of quark masses covered, a clear discrepancy is seen. This is to some extent surprising, since the overlap expansion is only carried out to third order and therefore one might expect that the properties of the input operator are only slightly changed. The discrepancy vanishes in the Edinburgh plot Fig. 7.15, where the results on the $`12^3\times 24`$ lattice agree for the FP and overlap-improved operators. The difference in the cut-off effects thus cancels in these mass ratios.

From our results we therefore find that the overlap-improvement changes the scaling behavior of the fermion action. In Fig. 7.17, the point for the overlap-improved operator at $`\beta =3.0`$ seems to agree with the continuum value obtained by extrapolating all the data to a single point at $`a^2\sigma =0`$, but this coincidence can be misleading due to a possible overall shift from the scale determination. To conclude whether the scaling violations are decreased or increased by the overlap, additional measurements with the overlap-improved FP Dirac operator at different values of $`\beta `$ are needed.

### 7.6 Hadron Dispersion Relations

Another quantity that measures the magnitude of discretization errors of a given action is the energy-momentum dispersion relation for hadrons $`E(\stackrel{}{p})=m^2+\stackrel{}{p}^2`$, or equivalently the squared speed of light

$$c^2(\stackrel{}{p})=\frac{E^2(\stackrel{}{p})m^2}{\stackrel{}{p}^2},$$
(7.14)
which should be $`c^2=1`$ for all momenta according to the continuum dispersion relation. At large lattice spacings, unimproved lattice actions are known to suffer from substantial deviations from the continuum relation. Even for $`𝒪(a)`$-improved clover fermions, $`c^2`$ deviates by 20%–30% at $`a=0.25`$ fm and $`p0.6`$ GeV for pseudoscalar and vector mesons . Our coarsest lattice spacing is $`a=0.16`$ fm, therefore we can not compare our results directly to this data. However, the energy-momentum dispersion relation for our pseudoscalar and vector mesons calculated with the parametrized FP Dirac operator on the lattice of size $`16^3\times 32`$, shown in Fig. 7.20, is consistent with $`c^2=1`$ over the whole range of momenta considered. The parametrization of the Fixed-Point Dirac operator therefore seems to conserve the classically perfect properties very well also for the dispersion relation.

It is interesting to check how the overlap-improved FP Dirac operator performs. We plot the squared speed of light at the smallest non-zero momentum on the $`12^3\times 24`$ lattice with $`\beta =3.0`$ in Fig. 7.21, and compare the results from the parametrized FP and the overlap-improved FP Dirac operator at the three largest quark masses. While the data from the FP operator again agrees with $`c^2=1`$ within errors for both mesons, the results from the overlap-improved operator are too large by 7(3)% for the pseudoscalar and 14(5)% for the vector meson. The overlap construction therefore seems to drive the hadronic dispersion relation of the FP operator away from the continuum form. The situation appears to be analogous to the case of the free Wilson operator, where the dispersion relation is also deteriorated by the overlap .

## Chapter 8 Conclusions and Prospects

After the beautiful properties of classically perfect Fixed-Point actions have been verified in various models, and recently a sophisticated FP $`SU(3)`$ gauge action has been parametrized and successfully tested, we have constructed in this work and in a parallel thesis a FP fermion action for lattice QCD. In this thesis it is applied to calculations of the quenched light hadron spectrum.

##### Properties of Fixed-Point Fermions

We find that the parametrization of FP fermions in QCD is a feasible task. The resulting Dirac operator is rather complex, including the full Clifford structure and hypercubic couplings, and uses two-level hypercubic RG smeared links as a gauge input. We have presented a way to efficiently build the gauge paths that are stored in the matrix representing the Dirac operator. Technically the multiplication of such a Dirac operator with a vector is 36 times more costly than for the Wilson operator, but the overall factor might be smaller in actual simulations due to the faster convergence of the iterative inversion. The computational cost also depends strongly on the architectural details of the computer.

What one gets at this moderately higher price is a lattice fermion action which preserves chiral symmetry to a high level and has largely reduced cut-off effects. The presence of approximate chiral symmetry manifests itself in the eigenvalue spectrum, which is very close to the exact Ginsparg-Wilson case, in the eigenvalue flow of discretized exact instantons, and in the small residual additive quark mass renormalization determined in spectroscopy calculations. Due to the good chiral properties, we are able to perform lattice simulations at very small quark masses, corresponding to $`m_{PS}/m_V<0.3`$.

Like chiral symmetry, reduced cut-off effects are a consequence of the close approximation to the perfect action that is achieved by our parametrization. A high level of improvement compared to standard actions is observed in the scaling of hadron masses, where only small—if any—$`𝒪(a^2)`$ effects seem to remain. Moreover, the hadron dispersion relation shows an impressive agreement with the continuum form.

For applications where very good, but not exact chiral symmetry is needed, the parametrized FP Dirac operator therefore provides a highly competitive alternative to domain wall or Wilson overlap fermions, with the additional advantage of reduced cut-off effects. Furthermore, the parametrized FP fermion action is ultralocal, thus avoiding possible problems with suboptimal localization properties observed for the Wilson overlap.

##### Overlap with Fixed-Point Fermions

If exact chiral symmetry is needed, the overlap construction is the only known solution. We examined the consequences of taking the FP instead of the Wilson Dirac operator as a starting point for the overlap. First of all, since the FP operator is almost chiral, only an expansion to low order is needed for the overlap, illuminated also by the far better behavior of the FP kernel in the exact treatment of the smallest eigenmodes for the calculation of the inverse square root. The need for only a low-order overlap expansion partially compensates the higher computational cost of the FP kernel in the overall cost for the FP overlap operator. In this work, we considered a overlap expansion to third order and examined the effect on hadron spectroscopy. The localization properties of this operator are much better both for the couplings and in reproducing the zero mode of a discretized exact instanton.

The cut-off effects of the overlap-improved FP Dirac operator are another issue. While the overlap construction removes effects to $`𝒪(a)`$, the higher-order effects might even increase. In our measurements, the hadron dispersion relation for instance seems to get worse compared to the FP operator. The scaling violations for the vector meson mass at $`a=0.16`$ fm can not be adequately judged due to uncertainties in the scale determination. We do not have enough data here to make a definite statement about the cut-off effects of the overlap with FP kernel, but the first results indicate that the improvement from the overlap needs not be large and might even be negative for some quantities. This however has to be taken under the premise that the starting operator is already highly improved.

Overall, we find that the overlap with FP kernel yields an operator which has better locality properties than the standard Wilson overlap. Whether in applications where exact chiral symmetry is needed it is better suited than the Wilson overlap remains to be examined. We refer to the thesis of Thomas Jörg for a further discussion of applications and results with the FP overlap.

##### Physical Results

Exploiting the chiral symmetry of FP fermions, we have studied topological finite-volume effects in quenched pion propagators, which are induced by zero modes of the Dirac operator. A recently suggested solution, amounting to the subtraction of the effects in the meson propagators, turns out to be the most practical and efficient way to remove these effects. The explicit calculation and subtraction of the zero modes in the quark propagators does not appear to be competitive. Applying the former method, we find that the intercept of the squared pseudoscalar meson mass with zero is consistent with the determination from the axial Ward identity quark mass.

Having clarified the complication from zero modes, we confirm the presence of the quenched chiral logarithm in the squared pseudoscalar meson mass and measure its coefficient. The resulting value is significantly larger than for previous measurements with non-chiral actions, but consistent with the theoretical expectation.

In a preliminary analysis of our data, we find some hints, but no clear evidence, for quenched terms in the chiral extrapolations of vector mesons and baryons. If however such terms have to be taken into account, the value of the hadrons at the physical mass of the light quarks might be substantially affected at least at fixed lattice spacing.

##### Prospects

The results for hadron spectroscopy with FP fermions are encouraging. Further work leads into various directions. First, our experiences from the parametrization and the simulations give insight into the strengths and weaknesses of the current parametrization, which might help in finding a set of parameters that describes the perfect action even better. Second, the application of FP fermions to physical problems like pion scattering is promising. Third, the construction of chiral currents and FP operators might complement the FP QCD action in the future.

In the context of the BGR collaboration, a comparison of these results to hadron spectroscopy simulations with a different chiral formulation of lattice fermions is under way.

## Appendix A Non-Perturbative Gauge Fixing

In lattice QCD, the functional integrals used to determine physical observables are well-defined due to the finite number of space-time points and the gauge fields being elements of a compact group. In general it is therefore not necessary to fix the gauge. For certain applications however, amongst which are the computation of gauge-dependent quantities like gluon propagators or matrix elements used in non-perturbative renormalization techniques , it is unavoidable to work in a fixed gauge background. Fixing the gauge is also a way to make life easy when using extended sources for calculating quark propagators, as it is then not necessary to ensure gauge invariance by hand. Thus, it is nice to have a fast and reliable algorithm to numerically fix a lattice gauge configuration to a certain gauge. Although lattice gauge fixing is mainly a technical aid for doing calculations, a fair amount of work has been done on this subject by the lattice community, as a recent review of the current status shows .

### A.1 Gauge Fixing and the Lattice

On the lattice, the fundamental variables for the gauge degrees of freedom are not the continuum fields $`𝒜_\mu (x)`$ themselves, but the matrices $`U_\mu (x)`$, which are group elements of $`SU(3)`$ in the fundamental representation and are formally defined as parallel transporters of the color interaction between lattice sites,

$$U_\mu (x)e^{iag𝒜_\mu (x)},$$
(A.1)
so that the $`U_\mu (x)`$ fields live on lattice links. We define the lattice gauge potential

$$A_\mu (x)\frac{1}{2iag}.\left[U_\mu (x)U_\mu ^{}(x)\right]|_{\mathrm{traceless}},$$
(A.2)
which is suggested by the formal relation (A.1) between lattice and continuum gauge fields $`𝒜_\mu (x)`$. Note that $`A_\mu (x)`$ is equivalent to the continuum gauge field $`𝒜_\mu (x)`$ only in the continuum limit $`a0`$. While Eq. (A.2) is a common way to define the lattice gauge potential, it is not unique, and other definitions which differ only by irrelevant terms are perfectly allowed. The choice of one particular definition then also leads to a particular solution of a given gauge fixing condition. However, by comparing Green’s functions, it has been checked that in the continuum limit the continuum gluon field described by different definitions of $`A_\mu (x)`$ on the lattice is unique .

Under a local gauge transformation, the lattice gauge field $`U_\mu (x)`$ transforms like

$$U_\mu (x)U_\mu ^G(x)=G(x)U_\mu (x)G(x+\widehat{\mu }),$$
(A.3)
where $`G(x)`$ are elements of the gauge group SU(3) living on lattice sites. To fix the gauge, a condition $`f(U_\mu ^G(x))=0`$ is introduced, which should pick out one configuration per gauge orbit. In general however, there are multiple solutions to this equation for a given gauge configuration. These solutions belonging to the same gauge orbit are called Gribov copies. The question to what extent this Gribov ambiguity introduces systematic uncertainties in lattice results is not definitely answered. While for certain problems like the calculation of the photon propagator in compact U(1)  or studies of center vortices , Gribov copies are known to distort measurements heavily, the effect seems to be barely distinguishable from the statistical noise for measurements of the axial vector renormalization factor $`Z_A`$  and $`B`$ meson correlators with smeared sources . In the latter study, which was performed on $`10^3\times 20`$ lattices with the Wilson gauge action at $`\beta =6.0`$, the Gribov noise could actually be identified, and it was argued that for larger lattices, its size could become significant compared to the statistical noise. Hence one has to keep in mind that the noise from Gribov copies is a possible source of errors in measurements of hadronic correlators, if we use gauge fixing in combination with extended, gauge-dependent operators.

The gauge fixing condition can be freely chosen. The most common choices in lattice QCD belong to the general class of $`\lambda `$-gauges, which are characterized by a continuum gauge fixing condition

$$\lambda _0𝒜_0+_i𝒜_i=0.$$
(A.4)
For $`\lambda =1`$ one gets Landau and for $`\lambda =0`$ Coulomb gauge. These two conditions are equivalent to finding the extremal value of the lattice functional

$$F_U[G]=\text{Re}\mathrm{Tr}\underset{x}{}\underset{\mu =1}{\overset{l}{}}U_\mu ^{G(x)}(x),$$
(A.5)
where the second sum runs over the spacial indices only ($`l=3`$) for Coulomb gauge and over all space-time indices ($`l=4`$) for Landau gauge. Again, exact equivalence between the lattice and continuum gauge fixing conditions holds only in the continuum limit, and one might consider to reduce discretization error with improved gauge fixing conditions . Since it is numerically impossible to find the global minimum of the functional (A.5), which would be a unique solution up to a global gauge transformation, we specify to take any local minimum as our gauge-fixed configuration. Local minima of (A.5) are numerically found in an iterative procedure. There are several gauge fixing algorithms on the market, we present here the Los Alamos method with improved convergence by stochastic overrelaxation.

### A.2 The Los Alamos Algorithm with Stochastic Overrelaxation

In this method, introduced by De Forcrand and Gupta , the minimizing functional is rewritten using the auxiliary variable

$$w(x)=\underset{\mu =1}{\overset{l}{}}\left(U_\mu (x)+U_\mu ^{}(x\widehat{\mu })\right),$$
(A.6)
so that the sum over the lattice sites in (A.5) is replaced by a sum over only half the lattice sites. If we assign the colors red and black to the lattice sites in a checkerboard manner, the functional reads

$$F_U[G]=\frac{1}{2}\text{Re}\mathrm{Tr}\underset{x\mathrm{red}\mathrm{or}\mathrm{black}}{}w^G(x).$$
(A.7)
The basic idea is now to subsequently transform the gauge fields on the red and black lattice sites separately in a way that the minimizing functional monotonically decreases in every iteration step. The gauge transformation $`G(x)`$ is therefore chosen to be unity on the red(black) lattice sites at even(odd) iteration steps. Under a local gauge transformation, the field $`U_\mu (x)`$ then transforms like

$`U_\mu (x)`$ $``$ $`U_\mu ^G(x)=G(x)U_\mu (x),`$
$`U_\mu ^{}(x)`$ $``$ $`U_\mu ^G(x)=U_\mu ^{}(x)G^{}(x).`$ (A.8)
This gauge transformation amounts to one step in the iterative process. For the variable $`w(x)`$ introduced above, the gauge transformation reads

$$w(x)w^G(x)=G(x)w(x).$$
(A.9)
The transformation $`G(x)`$ is now chosen independently on every other lattice site such that

$$\mathrm{Re}\mathrm{Tr}G(x)w(x)\mathrm{Re}\mathrm{Tr}w(x).$$
(A.10)
We choose $`G(x)`$ to be the projection of $`w(x)`$ onto the $`SU(3)`$ group manifold which maximizes the left hand side of the equation. This is done by iterative maximization of $`SU(2)`$ subgroups . After each step the roles of the red and black sites are interchanged.

To overcome problems with critical slowing down when fixing the gauge on large lattices, several acceleration methods are discussed in the literature, amongst which are Fourier preconditioning , overrelaxation  and multigrid schemes . A very simple method is stochastic overrelaxation , which is based on the idea to overdo the local maximization once in a while. More precisely, with a probability $`0p_{\mathrm{or}}1`$ one applies a local gauge transformation $`G^2(x)`$ instead of $`G(x)`$. For $`p_{\mathrm{or}}=0`$ there is no change to the algorithm, while for $`p_{\mathrm{or}}=1`$ the algorithm does not converge. For intermediate values of $`p_{\mathrm{or}}`$, a dramatic speedup in the convergence can be reached . However, the optimal value for the parameter $`p_{\mathrm{or}}`$ depends quite strongly on the lattice volume and the smoothness of the gauge configuration. It is therefore necessary to optimize $`p_{\mathrm{or}}`$ for every lattice size and lattice spacing separately to get fast convergence.

#### A.2.1 Convergence Criterion

It remains to define a criterion which tells us when the desired accuracy is reached and the algorithm can be stopped. Besides monitoring the value of the functional $`F_U[G]`$ during the process, we determine the gauge fixing accuracy by measuring its first derivative

$$\sigma =\underset{x}{}\mathrm{tr}\left[\mathrm{\Delta }^G(x)\mathrm{\Delta }_{}^{G}{}_{}{}^{}(x)\right],$$
(A.11)
where

$$\mathrm{\Delta }^G(x)=\underset{\mu }{}\left[A_\mu ^G(x)A_\mu ^G(x\widehat{\mu })\right],$$
(A.12)
is a discretized version of the derivative in the continuum gauge fixing condition (A.4), and $`A_\mu ^G(x)`$ is the gauge transformed lattice gauge potential (A.2). At a local minimum of the gauge fixing functional, $`\sigma `$ vanishes. For our purposes, we stop the algorithm if $`\sigma <10^8`$.

#### A.2.2 Tuning of the Overrelaxation Parameter $`p_{\mathrm{or}}`$

It is easily possible to gain a factor of 4 in the number of iterations needed to reach a given gauge fixing accuracy when accelerating the algorithm by stochastic overrelaxation, as the gauge fixing procedure is very sensitive to the value of $`p_{\mathrm{or}}`$. Fig. A.1 shows a plot of the convergence history on one configuration in dependence of the stochastic parameter. The number of iterations needed to fix the gauge is reduced from over 8000 to below 2000 when setting $`p_{\mathrm{or}}=0.75`$. Unfortunately, the optimization of $`p_{\mathrm{or}}`$ is not a very stable procedure, and it can happen that for $`p_{\mathrm{or}}>0`$, the gauge fixing takes much longer on certain gauge configurations than on the others. As seen in Fig. A.2, this is related to the fact that shortly after the beginning, $`\sigma `$ starts to fluctuate wildly for some time, and only then a monotonic decrease is observed. A possible explanation of this behavior is that for a certain time it is not clear which Gribov copy the algorithm is going to choose. In Fig. A.3 we plot the difference of the gauge fixing functional (A.7) from its final value together with $`\sigma `$. Several plateaus in the value of $`F_U`$ show up before the algorithm decides which local minimum to take.

As soon as the region of monotonic decrease of $`\sigma `$ is reached, the convergence is much faster with stochastic overrelaxation than without. It is also obvious from Fig A.2 that small changes in $`p_{\mathrm{or}}`$ can lead to a very different behavior in this fluctuating region, resulting in factors of 2 in the total number of iterations. Still, compared to the $`p_{\mathrm{or}}=0`$ case, a significant reduction is achieved.

In our tests, it appeared that a value of $`p_{\mathrm{or}}`$ slightly below the point where convergence is lost worked well in most cases. As the computer time needed for the gauge fixing was considerably smaller than for the generation of the gauge configurations, a further optimization of the procedure was not necessary.

### A.3 Coulomb vs. Landau Gauge

When one decides to work in a fixed gauge background in order to measure the hadron spectrum with extended operators, one has to choose a specific gauge. There is a theoretical argument why Landau gauge is not well-suited for the measurement of time-dependent correlation functions $`C(t)=𝒪_f(t)𝒪_i(0)`$ which are used to extract masses: Consider a gauge configuration $`U`$, which is fixed to Landau gauge:

$$UU_{\mathrm{gf}}.$$
(A.13)
Suppose we measure a spatially extended, gauge-dependent operator $`𝒪(U_{\mathrm{gf}},t=0)`$ on the time slice $`t=0`$ of the gauge-fixed configuration. Such an operator might for example be the smeared source of a hadron correlator. Now change one single gauge link of the configuration $`U`$ on a time slice $`t0`$, for example by rotating it into the opposite direction:

$$U_\mu ^{}(x)=\{\begin{array}{cc}U_\mu (x)\hfill & \text{if }\mu =1\text{ and }x=(x_0,t_00)\text{ },\hfill \\ U_\mu (x)\hfill & \text{otherwise}.\hfill \end{array}$$
(A.14)
Fixing the resulting gauge configuration $`U^{}`$ to Landau gauge,

$$U^{}U_{\mathrm{gf}}^{},$$
(A.15)
the gauge-fixed configuration $`U_{\mathrm{gf}}^{}`$ differs globally from $`U_{\mathrm{gf}}`$, since in the Landau gauge fixing process, spatial and temporal links enter in the minimizing functional (A.5). It follows that although the operator $`𝒪`$ is thought to be defined only on the time slice $`t=0`$, its value depends on the gauge fields on $`t0`$ time slices, and it is different on the two gauge fixed configurations:

$$𝒪(U_{\mathrm{gf}}^{},t=0)𝒪(U_{\mathrm{gf}},t=0).$$
(A.16)
Landau gauge fixing therefore introduces a non-local interaction in the time direction which might spoil the signal of hadron correlators.

The situation is different with Coulomb gauge: In the minimizing functional the links in the time direction do not enter, hence the gauge fixing is performed independently on each time slice. If a link on a time slice $`t0`$ is changed, the gauge fixed configuration on time slice $`t=0`$ is not affected, and the measured operator will not change.

In a study where hadron correlators with wall sources were compared on configurations fixed to both Coulomb and Landau gauge, a significant difference of order $`3\sigma `$ has been found for the mass of the $`\mathrm{\Delta }`$ meson , which might be an effect of this non-local interaction in the Landau gauge fixing. For our hadron spectroscopy study, we therefore fix the configurations to Coulomb gauge.

## Appendix B QCD on Large Computers

The non-perturbative approach to Quantum Chromodynamics as defined by the lattice regularization yields a theoretical formulation of the strong nuclear force which is ideal for treatment on large computers. The simulation of a realistic problem, like the scattering of two pions or the decay $`K\pi \pi `$, turns out to be a priori numerically very demanding due to the large number of degrees of freedom involved. The following arguments make it clear why lattice QCD calculations are hard: First, the spatial lattice volume has to be large enough, so that the wave functions of all involved particles fit into the box of side length $`L_s=N_sa`$ without getting squeezed. The temporal size $`L_t=N_ta`$ also has to be large, as one is mainly interested in the asymptotic behavior of correlation functions at large Euclidean time (see Chapter 5). Second, to get rid of discretization errors, the continuum limit $`a0`$ has to be taken by performing simulations at several lattice spacings $`a`$ and extrapolating the measured observables to the continuum. This implies that it is necessary to work with lattice spacings $`a`$ small enough to ensure a controlled extrapolation, and therefore $`N_s`$ and $`N_t`$ get large quickly. Third, as the computational effort grows with inverse powers of the quark mass, simulations are mostly performed at values larger than the physical mass. Hence, another extrapolation from results calculated at several higher quark masses is necessary to obtain physical values.

As shown in Section 5.1.1, the numerically demanding part in quenched lattice QCD calculations is the inversion of the Dirac operator, which is a large sparse complex matrix of rank $`4N_cV`$ on a lattice with volume $`V=N_s^3\times N_t`$. This matrix inversion can be done very efficiently on massively parallel computers, because it boils down to many complex multiplication and addition operations, and only minimal communication is necessary between different processors. Large LQCD calculations are done either on commercial machines or on custom-built computers dedicated to lattice QCD. Commercial supercomputers are multi-purpose machines and have to perform well on a very diverse range of problems. Their architecture is therefore highly sophisticated in order to process all kinds of complex code efficiently. Due to the simple computational structure of lattice QCD calculations, building a dedicated computer is in general the much more cost-effective alternative. The lattice groups which have been running the largest simulations in the last few years all use custom-built machines: The Japanese CP-PACS computer in Tsukuba , several machines of the European APE project  and the Anglo-American QCDSP  computers and its successor QCDOC . For our work, we had in contrary access to various multi-purpose machines spanning a large range from workstations to supercomputers. Tables B.1 and B.2 show the different platforms that were used for our simulations and their most relevant properties.<sup>1</sup><sup>1</sup>1Unfortunately, we could not run simulations on the IBM Power4 in Manno yet, as it took up user operations later than expected. The runs on the Hitachi were performed in the framework of the BGR collaboration .

In Appendix B.1, we summarize some technical details of the two computers we mostly worked on. Then we present in Appendix B.2 benchmark measurements from the spectroscopy simulations on the Hitachi and in particular the performance of our code under MPI parallelization. Finally we give in Appendix B.3 a brief introduction to modern matrix inversion algorithms and their extension to shifted linear systems, which is very useful for performing QCD simulations at many quark masses.

### B.1 Specifications of Utilized Supercomputers

It is well-known that what is called a supercomputer at a given time might become inferior to a desktop PC only a few years later, and therefore technical details outdate quickly and are of marginal interest. But as for our simulations the architecture of the involved computers differed quite a lot, and the process of making the simulations run efficiently on the available platforms depends strongly on the technical details of the machines, we present the main properties of the two mostly used computers in some more detail.

#### B.1.1 The NEC SX-5/16

The Swiss Center for Scientific Computing in Manno installed in early 2000 a NEC SX-5 parallel-vector computer with eight processors and 64 GB of shared main memory. In 2001, two more processors were added. At the time it went operational, it was ranked at position 242 in the list of the top 500 supercomputers . This computer differs from the other machines considered here by its vector architecture, which allows to obtain a very high single-processor performance: The basic building block is a 250 MHz CPU with 16 vector pipelines capable of processing two floating-point instructions each per clock cycle, leading to a theoretical peak performance of 8 GFLOPS per CPU. Although it is possible to run parallel processes, we only worked with scalar code due to the small number of installed processors and the long waiting time of the parallel queues.

The SX-5 is a good choice if one needs a lot of memory and does not want to write parallel code, as long as the problem is well-suited for vectorization. The compiler supports automatic vectorization, so it is possible to easily migrate existing scalar code from a workstation to the SX-5. However, to reach good performance, it is essential to optimize the time-critical parts by hand or with compiler directives. Moreover, code which does not vectorize well, like the algorithms used for the Monte Carlo update of the gauge configurations, ends up running extremely slow. The reason is that vector instructions need some time to be initialized, therefore the vectors have to be as long as possible in order to make the overhead irrelevant. In the program code, this manifests itself in the length of loops in which vector instructions appear. If these loops are not long enough, as it happens when doing manipulations on $`SU(3)`$ matrices, it is then even faster to run the program in scalar mode, which means that the machine is in this case inferior to every desktop computer due to its comparably low clock frequency.

#### B.1.2 The Hitachi SR8000-F1

A very different kind of architecture is provided by the Hitachi SR8000-F1 at Leibniz Rechenzentrum in München. This massively parallel scalar computer was also installed in early 2000 and was at that time the world’s fastest computer dedicated to academic research, ranked at 5th position of all supercomputers . At installation time, it comprised of $`112\times 8`$ modified Power3 processors running at a clock frequency of 375 MHz. Each unit is able to process two multiply/add instructions simultaneously. Eight processors are grouped in a node, which can be treated by the programmer like a single CPU, and each node has access to 8 GB of main memory<sup>2</sup><sup>2</sup>2A few special nodes have 16 GB of main memory. The theoretical peak performance per node is $`8\times 1.5=12`$ GFLOPS. In January 2002, the machine was upgraded from 112 to 168 nodes, reaching more than 2 TFLOPS peak performance. The inter-node communication is realized with a multi-dimensional crossbar delivering a bi-directional peak bandwidth of 1 GB/s.

Programs running on the SR8000 have to be parallelized in order make use of multiple nodes. The parallelization within the eight processors of one node is automatically done by the compiler, which also provides a hardware-based ’pseudo-vectorization’ facility that imitates a vector processor. In our applications, it seemed that compared to a real vector computer, the SR8000 was more tolerant when running code that does not vectorize well, resulting in significantly better performance for such programs. However, to make full use of the computer’s capabilities, it is also necessary to tune the programs by optimizing the arrangement of array elements and loops and by placing appropriate compiler directives for parallelization and pseudo-vectorization at the time-critical code segments.

### B.2 Measurements of Parallel Performance

For larger simulations on the Hitachi SR8000, explicitly parallel code is required. The common standard for programming the communication between different nodes of a computer or even different computers is the Message Passing Interface (MPI). In lattice QCD, explicit parallelization of the quark propagator code for $`n`$ nodes is in general a very simple task, since with an iterative solver the inversion of the Dirac operator reduces to matrix-vector multiplications. The Dirac operator can be split into $`n`$ parts containing $`12V/n`$ rows, and it only remains to write a distributed matrix-vector multiplication. However, for our implementation of the FP Dirac operator, parallelization is not completely trivial, because in the process of constructing the matrix it can occur that a gauge path calculated on a certain node needs to be stored on a different node. This is merely a consequence of the way our low-level routines for building the gauge paths in the Dirac operator are designed, and not a problem of the FP Dirac operator itself. Therefore, to make the construction of the Dirac operator work in parallel, some inter-node communication is needed. The Fixed-Point R operator is free of this complication, so parallelizing its construction is trivial. In order to have a completely parallel code, also the vectors of size $`12V`$ appearing in the different algorithms like the matrix inverter or the eigenvalue solver have to be distributed.

Not only in terms of computation time, but also in terms of storage the parallelization of the FP Dirac operator is crucial. The memory needed for storing the Dirac operator on a lattice of volume $`V`$ is $`12\times 81\times 12\times V\times 16`$ bytes<sup>3</sup><sup>3</sup>3We use double precision complex numbers, requiring 16 bytes, in all our code., which exceeds the shared main memory available on most machines even at moderate lattice sizes. Table B.3 lists the storage requirements for various elements of our code at the lattice volumes used in the simulations.

It is clear that the communication between different nodes introduces an overhead, and that this parallelization overhead increases with the number of nodes. Therefore it is important to find a reasonable balance between the gain in wall-clock time and the loss in CPU time due to parallelization. To quantify the parallelization overhead, we list in Tables B.4B.7 for several computational tasks the wall-clock time and the overhead factor $`\omega `$, which is defined as

$$\omega =n\frac{t_n}{t_{\mathrm{ref}}},$$
(B.1)
where $`n`$ is the number of nodes, $`t_n`$ is the wall clock time for the task running on $`n`$ nodes and $`t_{\mathrm{ref}}`$ is the wall-clock time for the smallest $`n`$ on which it was possible to run the task due to memory limitations. As shown in Table B.4, the dot product of two vectors does not really profit from parallelization on the smallest lattice. The situation changes drastically on the larger lattices, where the distribution of the vectors is crucial if one does not want to have severe slowing down in algorithms which perform many vector operations. For the construction of the FP Dirac operator (Table B.5), parallelization works well. The overhead introduced by the above mentioned inter-node communication is small for reasonable ratios of the lattice volume $`V`$ and the number of nodes $`n`$. That the construction of the $`R`$ operator parallelizes trivially can be seen in Table B.6.

The crucial quantity for quenched QCD simulations is the time needed for one iteration of the matrix inversion algorithm, which is listed in Table B.7. This is by far the most time-consuming task, and as one sees the overhead is still considerable. On the $`12^3\times 24`$ lattice, for example, there is only $`50\%`$ gain in time when going from 4 to 8 nodes, so it is advisable to run the simulations on the smallest $`n`$ possible for a given $`V`$. Comparing the results on the different lattice sizes running on 8 nodes shows that the time increases proportionally to the volume. Hence it seems that the code does not profit anymore from longer loops as they appear for larger volumes.

Typical benchmark measurements of the overall performance for the quark propagator inversion in hadron spectroscopy are given in Table B.8. These values depend on the number of iterations needed for the inversion of the Dirac operator, thus they only give an estimate for the performance of simulations in which the smallest quark mass is given by those used in our runs. The bottom line is that our code is reasonably efficient, running at an overall rate of around 30% of peak performance in the production runs. This number noticeably increases when considering only the matrix inversions, since the construction of the $`D`$ and $`R`$ operators and the time for I/O and MPI setup decrease overall performance: The matrix-vector multiplications without communication run at 6.3 GFLOPS per node for the Dirac operator $`D`$ and at 8.6 GFLOPS per node for the $`R`$ operator, which is remarkably fast.

### B.3 Matrix Inversion Techniques

The key element for efficient simulations of quenched QCD is a fast matrix inversion algorithm. As the rank of the Dirac operator matrix is far too large to perform an exact inversion, the methods of choice are iterative procedures, and the most widely used algorithms for QCD are variants of Krylov subspace methods. Consider the linear system of equations

$$Mx=b,$$
(B.2)
where $`M𝐂^{R\times R}`$ is in general a non-hermitean matrix and $`b𝐂^R`$ is the source vector on which the inversion is carried out. Choosing an initial guess $`x_0`$, the initial residual $`r_0=bMx_0`$ is defined. The Krylov subspace $`K^n`$ is then given by

$$K^n=\mathrm{span}\{M^mr_0;m=0,\mathrm{},n\}.$$
(B.3)
The common feature of all variants of Krylov space solvers is that the solution of (B.2) is iteratively approximated using an orthogonal basis of the Krylov subspace $`K^n`$. The iterative solution $`x_n`$ is of the form

$$x_n=x_0+q_{n1}(M)r_0,$$
(B.4)
where $`q_{n1}`$ is a polynomial of maximum degree $`n1`$. The iterative residual $`r_n`$ is therefore

$$r_n=\left(1Mq_{n1}(M)\right)r_0.$$
(B.5)
The best-known Krylov space solver is the Conjugate Gradient (CG) algorithm, which however only works for hermitean matrices. Refined types of algorithms are the Conjugate Gradient Squared (CGS) , Quasi-Minimal Residual (QMR) , Generalized Minimal Residual (GMRES)  or Bi-Conjugate Gradient (BiCG)  algorithms. These more sophisticated methods show either faster convergence or increased stability and are also applicable to non-hermitean matrices. We compared the convergence of the matrix inversion for several algorithms on a toy lattice of size $`4^4`$, with $`M`$ given by the parametrized FP Dirac operator. The results in Fig. B.1 show that the stabilized Bi-Conjugate Gradient (BiCGStab), transpose-free QMR and CGS algorithm perform comparably well in this test. The order $`l`$ in the BiCGStab($`l`$) variants denotes the order of the subspaces which are intermediately orthogonalized in the iteration process. Closer investigations of the properties of the different methods have shown that BiCGStab is in general a good and reliable choice .

As QCD simulations are in generally done at several quark masses, a significant computational gain can be obtained using multi-mass solvers , which exploit the fact that it is possible to get the solution of the shifted linear system

$$(M+\sigma )x=b,$$
(B.6)
for a whole set of values $`\sigma 𝐂`$ at the cost of only one inversion. In QCD this implies that the cost of a multi-mass inversion is equivalent to the cost of a single inversion at the smallest quark mass. For our spectroscopy calculations, we worked with the multi-mass BiCGStab algorithm in <sup>4</sup><sup>4</sup>4A typo in the algorithm had to be corrected.. The drawback of this method is that it is no longer possible to improve the condition number of $`M`$ by using preconditioning techniques, because the starting guess $`x_0`$ for the multi-mass inversion is required to be zero . If one however wants to invert at a large number of quark masses, this disadvantage is more than compensated by the gain from needing only one inversion.

## Appendix C Conditions on the Dirac Operator from Discrete Symmetries

While a lattice transcription of the continuum quark action allows many possible discretizations of the Dirac operator, it is essential that any lattice Dirac operator has the same properties under discrete symmetry transformations as in the continuum. We derive here the transformation properties under reflection of a coordinate axis and charge conjugation from the basic properties of the Dirac spinors and the gauge fields. The representation of the Clifford algebra which is used is given in Appendix E.

### C.1 Reflection of an Axis

A reflection in direction of a coordinate axis $`\eta `$ with $`\eta =1,\mathrm{},4`$ can be described by a unitary operator $`𝒫_\eta =𝒫_\eta ^\mathrm{\Psi }𝒫_\eta ^U`$. The operator $`𝒫_\eta ^\mathrm{\Psi }`$ acts on the fermion fields,

$`𝒫_\eta ^\mathrm{\Psi }\mathrm{\Psi }(n)𝒫_{\eta }^{\mathrm{\Psi }}{}_{}{}^{1}`$ $`=P_\eta \mathrm{\Psi }(\stackrel{~}{n}),`$ (C.1)
$`𝒫_\eta ^\mathrm{\Psi }\overline{\mathrm{\Psi }}(n)𝒫_{\eta }^{\mathrm{\Psi }}{}_{}{}^{1}`$ $`=\overline{\mathrm{\Psi }}(\stackrel{~}{n})P_\eta ^1,`$ (C.2)
where $`\stackrel{~}{n}`$ is the reflected lattice space-time variable,

$$\stackrel{~}{n}_\mu \{\begin{array}{cc}n_\mu \hfill & \text{for }\mu =\eta ,\hfill \\ n_\mu \hfill & \text{for }\mu \eta ,\hfill \end{array}$$
(C.3)
and $`P_\eta `$ is a matrix in Dirac space (in our representation $`P_\eta =\gamma _\eta \gamma _5`$). On the other hand, $`𝒫_\eta ^U`$ acts on the gauge fields,

$$𝒫_\eta ^UU_\mu (n)𝒫_{\eta }^{U}{}_{}{}^{1}=U_\mu ^{𝒫_\eta }(n),$$
(C.4)
where the reflected gauge field $`U_\mu ^{𝒫_\eta }`$ is defined as

$$U_\mu ^{𝒫_\eta }(n)\{\begin{array}{cc}U_\eta ^{}(\stackrel{~}{n}\widehat{\eta })\hfill & \text{for }\mu =\eta ,\hfill \\ U_\mu (\stackrel{~}{n})\hfill & \text{for }\mu \eta .\hfill \end{array}$$
(C.5)
The lattice fermion action

$$S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]=\underset{n,n^{}}{}\overline{\mathrm{\Psi }}(n)D(n,n^{},U_\mu (n))\mathrm{\Psi }(n),$$
(C.6)
has to be symmetric under reflections of the $`\eta `$-axis:

$$S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]=𝒫_\eta S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]𝒫_\eta ^1.$$
(C.7)
Inserting the action (C.6) into (C.7) and using $`𝒫_\eta 𝒫_\eta ^1=\mathrm{𝟏}`$ twice, we get

$$S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]=\underset{n,n^{}}{}𝒫_\eta \overline{\mathrm{\Psi }}(n)𝒫_\eta ^1𝒫_\eta D(n,n^{},U_\mu (n))𝒫_\eta ^1𝒫_\eta \mathrm{\Psi }(n)𝒫_\eta ^1.$$
(C.8)
The fermion and gauge fields transform as specified in Eqs. (C.1)–(C.4):

$$S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]=\underset{n,n^{}}{}\overline{\mathrm{\Psi }}(\stackrel{~}{n})P_\eta ^1D(n,n^{},U_\mu ^{𝒫_\eta }(n))P_\eta \mathrm{\Psi }(\stackrel{~}{n}).$$
(C.9)
It remains to reorder the summation over the $`\eta `$-component of the lattice variables $`n`$ and $`n^{}`$. As it does not matter whether the sum over a variable is performed from above or from below,

$$\underset{n_\eta =\mathrm{}}{\overset{\mathrm{}}{}}F(n)=\underset{n_\eta =\mathrm{}}{\overset{\mathrm{}}{}}F(n)=\underset{n_\eta =\mathrm{}}{\overset{\mathrm{}}{}}F(\stackrel{~}{n}),$$
(C.10)
the argument of the summand can be reflected without changing the value of the sum. Applying this to Eq. (C.9), we get

$$S[\mathrm{\Psi },\overline{\mathrm{\Psi }},U_\mu ]=\underset{n,n^{}}{}\overline{\mathrm{\Psi }}(n)P_\eta ^1D(\stackrel{~}{n},\stackrel{~}{n}^{},U_\mu ^{𝒫_\eta }(\stackrel{~}{n}))P_\eta \mathrm{\Psi }(n),$$
(C.11)
which provides us the condition for the Dirac operator

$$D(n,n^{},U_\mu (n))=P_\eta ^1D(\stackrel{~}{n},\stackrel{~}{n}^{},U_\mu ^{𝒫_\eta }(\stackrel{~}{n}))P_\eta ,$$
(C.12)
when comparing with the original action (C.6).

### C.2 Charge Conjugation

The unitary charge conjugation operator $`𝒞`$ acts on the fermion and gauge fields like

$`𝒞\mathrm{\Psi }𝒞^1`$ $`=C\overline{\mathrm{\Psi }}^T,`$ (C.13)
$`𝒞\overline{\mathrm{\Psi }}𝒞^1`$ $`=\mathrm{\Psi }^TC^1,`$ (C.14)
$`𝒞U_\mu 𝒞^1`$ $`=U_\mu ^{},`$ (C.15)
where the charge conjugation matrix fulfills $`C\gamma _\mu ^TC^1=\gamma _\mu `$ and can be expressed in our representation by $`C=\gamma _2\gamma _4`$. Invariance of the action under $`𝒞`$-transformations means

$$\overline{\mathrm{\Psi }}D(U_\mu )\mathrm{\Psi }=𝒞\overline{\mathrm{\Psi }}D(U_\mu )\mathrm{\Psi }𝒞^1.$$
(C.16)
Inserting $`𝒞^1𝒞`$ twice and using $`C^T=C^1`$, we get

$`\overline{\mathrm{\Psi }}D(U_\mu )\mathrm{\Psi }`$ $`=\mathrm{\Psi }^TC^1D(U_\mu ^{})C\overline{\mathrm{\Psi }}^T`$
$`=\overline{\mathrm{\Psi }}C^1D(U_\mu ^{})^TC\mathrm{\Psi },`$ (C.17)
where the transposition extends over all index spaces. The transformation property of the Dirac operator under charge conjugation is therefore

$$D(U_\mu )=C^1D(U_\mu ^{})^TC,$$
(C.18)
completing the set of conditions from the $`𝒞`$, $`𝒫`$ and (Euclidean) $`𝒯`$ symmetries.

## Appendix D Collection of Data

### D.1 Hadron Masses

For each lattice, we list the bare input quark masses, the bias-corrected masses from correlated fits to the hadron propagators, the value of $`\chi ^2/df`$ for the fit and the fit range. The numbers in brackets and the superscripts denote bootstrap errors and bias (5.40), respectively.

#### D.1.1 Pseudoscalar Mesons

Pseudoscalar mesons are determined from the pseudoscalar (P), fourth component axial vector (A), and pseudoscalar minus scalar (P-S) correlators.

#### D.1.2 Vector Mesons, $`m_{\mathrm{PS}}/m_\mathrm{V}`$ and $`m_{\mathrm{Oct}}/m_\mathrm{V}`$

Together with the bias-corrected fitted masses for the vector meson, we list the mass ratios $`m_{\mathrm{PS}}/m_\mathrm{V}`$ and $`m_{\mathrm{Oct}}/m_\mathrm{V}`$, where the pseudoscalar mass is taken from fits to the P-S correlator at small and the P correlator at large quark mass.

#### D.1.3 Octet Baryons

For each lattice, the bias-corrected fitted masses for the N and N0 correlators, the value of $`\chi ^2/df`$ for the fit and the fit range are given. The correlator which was used for the ratio $`m_{\mathrm{Oct}}/m_\mathrm{V}`$ is marked with a star in each table.

#### D.1.4 Decuplet Baryons

For each lattice, the bias-corrected fitted masses for the D and D0 correlators, the value of $`\chi ^2/df`$ for the fit and the fit range are given.

### D.2 Unrenormalized AWI Quark Masses

The listed values and bootstrap errors of the unrenormalized quark masses from the axial Ward identity are determined by averaging the measurements of $`am_q^{\mathrm{AWI}}(t)`$ in Eq. (7.1) over the time range $`t[t_1,t_2]`$.

## Appendix E Conventions

### E.1 Dirac Algebra in Minkowski Space

In Minkowski space, the Dirac algebra is defined by the anticommutation relation

$$\{\gamma _M^\mu ,\gamma _M^\nu \}=2g^{\mu \nu }\mathrm{𝟏}.$$
(E.1)
From the elements $`\gamma _M^\mu `$ of the Dirac algebra, we construct the tensor

$$\sigma _M^{\mu \nu }\frac{1}{2i}[\gamma _M^\mu ,\gamma _M^\nu ],$$
(E.2)
and the pseudoscalar

$$\gamma _M^5i\gamma _M^0\gamma _M^1\gamma _M^2\gamma _M^3,$$
(E.3)
which satisfies $`\gamma _M^5\gamma _M^5=\mathrm{𝟏}`$ and $`\gamma _{M}^{5}{}_{}{}^{}=\gamma _M^5`$. The set of 16 elements

$$\mathrm{\Gamma }_M\{\mathrm{𝟏},\gamma _M^\mu ,\sigma _M^{\mu \nu },\gamma _M^\mu \gamma _M^5,i\gamma _M^5\},$$
(E.4)
with $`\mu \nu `$ forms a basis of the Dirac algebra and satisfies $`\gamma ^0\mathrm{\Gamma }_M^{}\gamma ^0=\mathrm{\Gamma }_M`$.

The Weyl representation of the Dirac algebra is given by the four-dimensional matrices

$$\gamma _M^0=\left(\begin{array}{cc}0& \mathrm{𝟏}\\ \mathrm{𝟏}& 0\end{array}\right),\gamma _M^i=\left(\begin{array}{cc}0& \sigma ^i\\ \sigma ^i& 0\end{array}\right),$$
(E.5)
where $`\sigma ^i`$ are the Pauli matrices

$$\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(E.6)
With this definition, the basis element $`\gamma _M^5`$ is diagonal:

$$\gamma _M^5=\left(\begin{array}{cc}\mathrm{𝟏}& 0\\ 0& \mathrm{𝟏}\end{array}\right).$$
(E.7)
We conclude this section by listing some transformation properties of the Dirac matrices in our convention. Under hermitean conjugation, we have

$$\gamma _{M}^{0}{}_{}{}^{}=\gamma _M^0;\gamma _{M}^{i}{}_{}{}^{}=\gamma _M^i.$$
(E.8)
Under complex conjugation the Dirac matrices transform as

$$\gamma _{M}^{0}{}_{}{}^{}=\gamma _M^0,\gamma _{M}^{1}{}_{}{}^{}=\gamma _M^1,\gamma _{M}^{2}{}_{}{}^{}=\gamma _M^2,\gamma _{M}^{3}{}_{}{}^{}=\gamma _M^3,$$
(E.9)
and finally, the transposition properties are

$$\gamma _{M}^{0}{}_{}{}^{T}=\gamma _M^0;\gamma _{M}^{1}{}_{}{}^{T}=\gamma _M^1;\gamma _{M}^{2}{}_{}{}^{T}=\gamma _M^2;\gamma _{M}^{3}{}_{}{}^{T}=\gamma _M^3.$$
(E.10)
### E.2 Analytic Continuation to Euclidean Space

Euclidean space-time is reached from Minkowski space-time by analytic continuation, rotating the time direction onto the imaginary axis:

$$x^0ix^4;x^ix^i.$$
(E.11)
In Euclidean space-time, the Dirac (or Clifford) algebra satisfies the anticommutation rule

$$\{\gamma _\mu ,\gamma _\nu \}=2\delta _{\mu \nu }\mathrm{𝟏}.$$
(E.12)
Due to the trivial Euclidean metric, upper and lower indices are the same. We find that the Euclidean matrices

$$\gamma _4\gamma _M^0,\gamma _ii\gamma _M^i,$$
(E.13)
satisfy Eq. (E.12).

From the properties of the Dirac matrices in Minkowski space, it then follows that all the $`\gamma _\mu `$ are hermitean. In the Euclidean version of the Weyl representation, $`\gamma _2`$ and $`\gamma _4`$ are real and symmetric, while $`\gamma _1`$ and $`\gamma _3`$ are purely imaginary and antisymmetric. Like in Minkowski space, we define the tensor element

$$\sigma _{\mu \nu }\frac{1}{2i}[\gamma _\mu ,\gamma _\nu ].$$
(E.14)
The pseudoscalar $`\gamma _5`$ is taken to be the same as in Minkowski space:

$$\gamma _5\gamma _M^5=\gamma _1\gamma _2\gamma _3\gamma _4.$$
(E.15)
The set

$$\mathrm{\Gamma }\{\mathrm{𝟏},\gamma _\mu ,\sigma _{\mu \nu },i\gamma _\mu \gamma _5,\gamma _5\},$$
(E.16)
with $`\mu \nu `$ then forms a hermitean basis $`\mathrm{\Gamma }^{}=\mathrm{\Gamma }`$ with elements

$$\mathrm{\Gamma }=\{S,V,T,A,P\},$$
(E.17)
transforming like scalars (S), vectors (V), tensors (T), axial vectors (A) and pseudoscalars (P).

## Acknowledgements

Many people have played an important role in making this work possible. Peter Hasenfratz introduced me to the subject of perfect actions and gave me the opportunity to work in the fascinating field of lattice QCD. Ferenc Niedermayer provided a wealth of ideas and advice from which I could profit.

During my thesis, I had the chance to collaborate with a number of people: Primarily I have to mention my fellow PhD student Thomas Jörg, with whom I worked closely on the construction of the Fixed-Point fermions. Kieran Holland participated during an essential phase in the project, parametrized the Fixed-Point $`R`$ operator and implemented the multi-mass inverter used in this work. In the first period of my thesis, I could rely on Urs Wenger and Philipp Rüfenacht, who were so kind to share their experience with perfect gauge actions and to answer my naive questions. Thanks to these people, I could also enjoy not only the scientific part of the lattice conferences and workshops.

I would like to thank Tom DeGrand for providing me the possibility of a memorable six-weeks stay in Boulder, for giving me insight into a range of technical problems and for sharing his scaling data. Thanks to Fernando Perez, with whom I shared the office there, I could experience some spectacular rock climbing in the Boulder area. Christof Gattringer initiated the BGR collaboration, consisting of the groups in Regensburg, Graz and Bern, due to which we could run simulations on the Hitachi computer in München, and he contributed gauge configurations and results from the chirally improved operator. I further thank Christian Lang for discussions on the gauge fixing algorithm. I thank the Swiss Center for Scientific Computing in Manno and the Leibniz Rechenzentrum in München for granting access to their supercomputers and offering technical support. Thanks go to Ottilia Hänni for dealing with the administrative hassles. At last, I thank the present and former members of the institute in Bern, who created an enjoyable atmosphere.

Most importantly, I was so fortunate to enjoy the social environment, support and love of my friends and family in all these years. I am deeply indebted to them.

This work has been supported by the Schweizerischer Nationalfonds and by the European Community’s Human Potential Programme under HPRN-CT-2000-00145 Hadrons/Lattice QCD, BBW Nr. 99.0143.