format: text

text:
# Contents

## Chapter 1 Introduction

The task of physics is the determination of all relevant quantities of Nature from basic “principles”. In the case of elementary particle physics the relevant guidelines which have proven to be successful are the symmetry and unification principles. Their application leads to the so-called standard model \[Gla61\], incorporating 3 of the 4 known basic forces, namely the electro-magnetic, weak and strong interactions. Basic degrees of freedom of the standard model are the fundamental fermions (quarks and leptons) and their interaction with each other is described via the exchange of bosonic gauge quanta. An important part is formed by the Higgs mechanism, generating mass for fermions and massive gauge bosons and introducing the scalar Higgs particle into the model. The theory is formulated as a quantised relativistic local gauge theory coupling the symmetry principle with the dynamics of fields \[Zin94, Che84\]. Main building blocks of the standard model are Quantum Electrodynamics (QED) with symmetry group $`U(1)`$ (as part of the electro-weak sector $`U(1)SU(2)`$) and Quantum Chromodynamics (QCD), relying on the symmetry group $`SU(3)`$.

Yet the standard model leaves a lot to be explained. It has as many as 19 free parameters \[Don92\] and no reason is given why a symmetry of $`SU(3)SU(2)U(1)`$ is realized in nature. Thus we would like a more advanced theory with fewer parameters and an even more symmetric structure. Possible ways of extending the standard model include the search for more fundamental constituents of quarks and leptons, string theory including supersymmetric partners to the known elementary particles and grand unified theories. In these higher dimensional symmetry groups are applied at very large energy scales, which for lower energies are broken down to the familiar standard model structure \[Sal90\].

One main research focus is the relativistic $`SU(3)`$ gauge theory including dynamical fermions, Quantum Chromodynamics (QCD) \[Fri73\], which describes hadrons and nuclei under strong interaction with the help of quarks and gluons as fundamental degrees of freedom. The beta-function is known at 3-loop level, leading in the high energy limit (or equivalently at small distances) to asymptotic freedom, i.e. the appearance of almost unbound quarks. In that case, the running coupling constant is small enough to allow perturbative methods. Applying these, there exist in this regime convincing predictions of the theory, which agree well with experiment \[Mut87\].

On the energy scale of the hadrons themselves, i.e. for the long distance part of the strong interaction, quasi-free quarks were never observed. The quark coupling in this regime is too large for perturbation theory around the free solution to be sensible. We therefore need non-perturbative methods to shed light on e.g. the questions of quark confinement and hadron masses.

Among the attempts to reach non-perturbative insight into the structure of QCD approaches using QCD sum rules, large $`N`$ expansions, the Bethe-Salpeter equation or Discretized Light-Cone Quantisation \[Bro91\] can be mentioned besides effective purely phenomenological models. For an overview discussion we refer to the textbook of Narison \[Nar89\].

Compared to these approaches, lattice gauge calculations offer a way to results stemming directly from the use of first principles. Lattice gauge theory uses a path integral representation with a space-time grid regularisation with lattice spacing $`a`$ \[Wil74\]. The theory is set up in such a way that the naive continuum limit $`a0`$ yields the desired continuum Lagrangian to leading order in the lattice spacing. This of course offers a certain freedom in the choice of the lattice theory. In this approach all quantities are measured in units of the lattice spacing or inverses thereof. Results have to be extrapolated to the continuum limit, in which the lattice spacing $`a`$ is set to zero and the dimensionless lattice size $`L`$ to infinity. In order to still have meaningful results, a continuous phase transition point has to be chosen for this procedure. At this point, all correlation functions diverge, so that results in physical units can still be finite. In this limit, we will therefore retain ratios of (e.g. mass) observables corresponding to relations between renormalised quantities. An important point is that certain continuum symmetries (like e.g. Lorentz symmetry) can be broken on a space time grid. The restoration of these symmetries has to be checked in order to justify extrapolation results. For completeness we mention that in order to exclude a contamination of the results by finite size effects, one would like to send also the physical lattice extent $`aL`$ to infinity, a limit called the thermodynamic limit.

One non-perturbative approach to lattice gauge theory obtains results with the help of Monte Carlo methods. Introductions can be found in \[Cre83, Kal90, Mon94\]. To illustrate the success of this approach, we mention that a recent work was able to determine ratios of hadron masses to an error of $`2\%`$ \[Aok97\]. Nevertheless, we are still far away from a complete understanding of the physics of hadronic systems using lattice gauge theory.

A severe problem appearing in simulations is the phenomenon of critical slowing down, i.e. the fact that the CPU cost of a simulations grows more than proportional to the space-time volume of the lattice \[Sok89\]. Especially fermion simulations are very time consuming \[For97a\]. They face the problem that the fermion interaction is inherently non-local, being described by a determinant in the path integral. Thus, up to now most large-scale calculations apply the so-called quenched approximation, replacing this determinant by a constant. Thereby the dynamics of fermionic vacuum fluctuations interacting with a bosonic gauge field are ignored.

The standard algorithm used to include dynamical fermions is the so-called Hybrid Monte Carlo algorithm \[Dua87\]. Here a trajectory of small size update steps is generated introducing a fictitious computer time coordinate, conjugate momenta and the corresponding classical equations of motion. At the end of the trajectory an acceptance step is used to make the algorithm exact. Other approaches, e.g. based on Kramers equation \[Hor91, Jan95\], have not shown significantly better behaviour than Hybrid Monte Carlo.

Recently, an alternative approach was proposed by M. Lüscher \[Lue94\]. This so-called Local Bosonic Algorithm (LBA) applies an $`n`$-th order polynomial approximation to the fermion determinant and thus makes it possible to rewrite the determinant as a set of $`n`$ local and bosonic integrals plus a non-local correction term. As simulations of locally coupled degrees of freedom are much easier, this could potentially reduce the computational task. This algorithm has generated considerable interest \[Jeg95a, For95, Bor96, Jan95a, Jan96a, Wol97\]. Examples for applications of the local bosonic algorithm are Monte Carlo simulations of lattice QCD \[For95, Jan96a\], Supersymmetry \[Mon97b\], the Hubbard model \[Saw97\] and the Schwinger model both with staggered \[Pea94\] and Wilson fermions \[Els96a, Els97\].

Of course, this new algorithm has to be tested and its efficiency compared to the standard Hybrid Monte Carlo algorithm (HMC) \[Bor96a, Jeg95a, Jeg96\]. In this context it may be noted that it is well known how to tune the parameters of Hybrid Monte Carlo routines, whereas this is still under investigation for the local bosonic case.

For the proposed studies of algorithms for dynamical fermions that are rather time-consuming, we want to advocate the massive 2-flavour Schwinger model (2-dimensional QED) \[Sch62, Het95\] as a low-cost laboratory \[Bar97\]. With the appearance of an axial anomaly, confinement, light pseudo-scalar and heavy scalar mesons, it has enough rich physics to be similar to QCD. One of the aims is to set up a program package facilitating further studies.

This work introduces the Schwinger model variants commonly used, and gives some details on the lattice 2-flavour Schwinger model. We will briefly describe the Hybrid Monte Carlo algorithm applied to set the scale for the CPU cost comparison and the four variants of the Hermitean version of the local bosonic algorithm we are testing. We will discuss some problems regarding numerical instabilities and topological metastabilities encountered in our simulations. Finally, we give results of our investigation of the CPU cost of the plaquette and the pion correlator.

Throughout this work, figures and tables are defered to the end of each chapter.

## Chapter 2 The model: Schwinger model – 2D QED

### 2.1 Continuum

Quantum Electrodynamics (QED) in one space and one time dimension ($`d=2`$) with $`N_f`$ flavours of fermions having mass $`m_a`$, $`a=1\mathrm{}N_f`$, is generally defined by the continuum Lagrangian in Euclidean space

$`={\displaystyle \frac{1}{4}}F^2+{\displaystyle \underset{a=1}{\overset{N_f}{}}}\left[\overline{\psi }^a(/+ie/A+m_a)\psi ^a\right],`$ (2.1)
where we use the standard conventions $`\mathrm{}=c=1`$. The fermion fields $`\psi `$ are Grassmann 2-spinors having dimension $`[m]^{\frac{1}{2}}`$ and the electromagnetic tensor

$`F^{\mu \nu }=\left(\begin{array}{cc}0& E\\ E& 0\end{array}\right)`$ (2.2)
in 1+1 dimensions only includes the electric field $`E`$, as because of the missing transverse directions no magnetic field is existing. We point out that the mass and the electric charge are of the same dimension.

For an introduction to Euclidean space we defer the reader to \[Mon94\]. Conventions are generally collected in App. A.

#### 2.1.1 Massless 1-flavour model

2-dimensional 1-flavour QED with massless fermions, i.e. $`N_f=1,m_1=0`$, was presented by Julian Schwinger as an analytically soluble model in 1962 \[Sch62\]. This version is generally called the Schwinger model.

Schwinger was able to show that a theory of a massless vector gauge field does not prevent the existence of a massive particle. The analytic solution allows only uncharged states. The thus formulated theory was therefore regarded as a model for complete charge screening by vacuum polarisation \[Rot79\]. This so-called “quark” trapping became an interesting object of study in view of quark confinement \[Rot79\]. The fact that no electron excitations exist was confirmed in numerous works \[Low71, Kog75\].

An equivalent formulation of the Schwinger model in boson fields using the mass parameter

$`m_B=e/\sqrt{\pi }`$ (2.3)
is given by a Lagrange density of the form \[Col75a, Col75, Col76, Kog75\]

$`={\displaystyle \frac{1}{2}}(_\mu \varphi ^\mu \varphi m_B^2\varphi ^2),`$ (2.4)
thus enabling us to interpret the spectrum as that of free bosons, the one-particle state with mass $`m_B`$ corresponding to a quark-antiquark pair, i.e. a meson. Above that we find a continuum of two meson, three meson $`\mathrm{}`$ states beginning at the minimum energies of $`2m_B`$, $`3m_B,\mathrm{}`$, where the corresponding mesons are at rest relative to each other. A pedagogical introduction can be found in \[Dit86\].

A further feature of the Schwinger model is the appearance of a non-vanishing vacuum condensate $`\overline{\psi }\psi `$. Naively we would expect conservation of electric as well as axial charge. Yet the phenomenon of an axial anomaly is observed \[Bro63, Kog75\].

#### 2.1.2 Massive 1-flavour model

Much interest was concentrated on the variant of the Schwinger model with massive fermions, $`N_f=1,m_1=m0`$, often called the massive Schwinger model. Using quasi-classical approximations, Coleman was able to confirm that also in this model no charged particles are allowed as asymptotic states \[Col75\]. One finds in the massless case decoupled, degenerate vacua upon which for each value of the parameter $`\theta `$, corresponding to a constant electric background field $`E_c=\frac{e}{2\pi }\theta `$, the same spectra are built up. In four dimensions this field would be eliminated due to pair creation. In the case of two dimensions this is not possible; electron-positron pairs are only created until the background field falls below a critical value of $`E_c^{\mathrm{crit}}=\frac{e}{2}`$, so that $`\theta `$ can be chosen to lie in the range $`[\pi ,\pi ]`$. Coleman showed in a semi-classical calculation that also in the massive case the background electric field is decisive for the structure of the spectrum. In the weak coupling limit we find for the number of stable states

$`N={\displaystyle \frac{4m^2}{\pi e^2}}{\displaystyle \frac{1}{1\frac{\theta ^2}{\pi ^2}}}\left(2\sqrt{3}\mathrm{ln}(2+\sqrt{3})\right)+𝒪(1),`$ (2.5)
which at vanishing coupling diverges as expected. In the strong coupling limit the result is one, two and three stable particles for the ranges $`|\theta |>\frac{\pi }{2}`$, $`0<|\theta |<\frac{\pi }{2}`$ and $`\theta =0`$ respectively \[Col75\].

Basic to the expansion around the massless case is the proof that this limit is allowed and leads to the soluble massless Schwinger model \[Kog75\]. The connection of these phenomena and topology was investigated in \[Par93, Man85\].

The massive Schwinger model can equivalently be described by an interacting boson field with the Lagrange density \[Col75\]

$`={\displaystyle \frac{1}{2}}\left[_\mu \varphi ^\mu \varphi m_B^2\varphi ^2+mm_B{\displaystyle \frac{e^\gamma }{\pi }}\mathrm{cos}(2\sqrt{\pi }\varphi )\right],`$ (2.6)
where $`\gamma =0.577\mathrm{}`$ is the Euler constant and $`m_B`$ the mass parameter Eq. 2.3. An expansion in the fermion mass yields through comparison of the coefficients of $`\varphi ^2`$ for the mass of the lowest state \[Kog75, Var94\]

$`M_1^2=m_B^2\left(1+2{\displaystyle \frac{m}{m_B}}e^\gamma \right)+𝒪(m^2).`$ (2.7)
In the framework of lattice gauge theory the massive Schwinger model was investigated numerically using Hamiltonian methods \[Kog75, Cre80\]. Here besides mass eigenvalues of the lowest-lying states also the coefficients of the linear potential in the limit of small mass were obtained. The results are in good agreement of some percent error with the continuum theory \[Cre80\]. Discretized Light-Cone Quantisation \[Bro91\] was able to calculate whole mass spectra in some Fock space approximation \[Ell87, Els95\]. These data were used to study finite temperature quantities and critical exponents of the massive Schwinger model \[Els96\].

#### 2.1.3 Massless $`N`$-flavour model

The $`N`$-flavour model can be solved analytically in the limit $`m_a=0,a=1\mathrm{}N`$ \[Het88\]. Main predictions are the existence of a massive pseudoscalar isosinglet state with mass

$`m=\sqrt{N}{\displaystyle \frac{e}{\sqrt{\pi }}}`$ (2.8)
and $`N^21`$ massless (Goldstone-like) states. We would like to point out that earlier work identifies only $`N1`$ pion states \[Aff86\]. Both the fermion condensate $`\overline{\psi }\psi `$ and the pseudoscalar density $`\overline{\psi }\gamma ^5\psi `$ are predicted to be zero.

#### 2.1.4 Massive 2-flavour model

Most useful for lattice investigations with Wilson fermions is the massive degenerate 2-flavour model $`N_f=2,m_1=m_2=m`$. It is first of all convenient for technical reasons as the effective probability after integration of the fermionic degrees of freedom is manifestly positive. On the other hand, the 2-flavour model has the nice feature of being rather QCD-like, as was observed in several classic papers \[Col75, Het95, Hos96\].

The model describes light pseudoscalar isotriplet states ($`\pi `$, analogous to the Goldstone particles of QCD), a pseudoscalar isosinglet state ($`\eta `$, much like the $`\eta ^{}`$ of QCD) and scalar mesons ($`a_0,f_0`$). For most of these the mass perturbation theory is known to first order \[Het95, Har94\]

$`m_\pi `$ $`=`$ $`2.066m^{\frac{2}{3}}e^{\frac{1}{3}}`$
$`m_{f_0}`$ $`=`$ $`\sqrt{3}m_\pi `$
$`m_\eta `$ $`=`$ $`\sqrt{2}{\displaystyle \frac{e}{\sqrt{\pi }}}.`$ (2.9)
The rich meson physics therefore allows to set the scale via the pion mass and measure physical mass ratios, e.g. $`\frac{m_\pi }{m_\eta }`$. Fermionic densities are generally zero with the exception of the fermion condensate, which is found to be real.

### 2.2 Lattice formulation

The lattice version of the Schwinger model with two flavours of fermions of identical mass is defined using the positive effective distribution for the compact $`U(1)`$ link variables $`U_{x,\mu }`$

$`P_{\mathrm{eff}}[U]`$ $``$ $`\mathrm{det}M^2e^{S_g[U]}`$ (2.10)
one obtains after integration over the fermionic Grassmann variables in the path integral Eq. 2.1 and discretization of the fields \[Wil74\]. We use the standard plaquette action

$`S_g[U]=\beta \mathrm{Re}{\displaystyle \underset{P}{}}(1U_P)`$ (2.11)
with $`\beta =\frac{1}{e^2a^2}`$ the dimensionless lattice gauge coupling and plaquette variables $`U_{Px}`$ defined starting at the lower left corner of each plaquette

$`U_{Px}=U_{x,\mu }^{}U_{x+\widehat{\mu },\overline{\mu }}^{}U_{x+\widehat{\overline{\mu }},\mu }U_{x,\overline{\mu }},`$ (2.12)
so that the lowest order expansion in the lattice spacing matches to the naive continuum limit. Barred Greek symbols signify the orthogonal direction. The Wilson fermion matrix is given by

$`M_{x,y}=\delta _{x,y}\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\widehat{\mu },y}(1+\gamma ^\mu )U_{x\widehat{\mu },\mu }+\delta _{x+\widehat{\mu },y}(1\gamma ^\mu )U_{x,\mu }^{}\right),`$ (2.13)
where $`\kappa =\frac{1}{2(ma+d)}`$ is the dimensionless Wilson parameter. Gamma matrices and details of notation are given in App. A.4. We remark that the lattice spacing $`a`$ is only explicitly shown in exceptional cases to avoid confusion. Generally it is set to 1 in all formulae.

Chiral symmetry is explicitly broken by the Wilson term. We remark that, as $`1`$ is in the centre of the $`U(1)`$ group, periodic or anti-periodic boundary conditions for the links should be irrelevant.

The 2-flavour case ensures that the effective distribution to be simulated is manifestly positive. We work with the Hermitean fermion matrix \[Lue94\]

$`Q`$ $`=`$ $`c_0\gamma ^5M`$
$`=`$ $`c_o\gamma ^5\delta _{x,y}c_o\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\widehat{\mu },y}\gamma ^5(1+\gamma ^\mu )U_{x\widehat{\mu },\mu }+\delta _{x+\widehat{\mu },y}\gamma ^5(1\gamma ^\mu )U_{x,\mu }^{}\right),`$
scaled so that its eigenvalues are in $`[1,1]`$ using the scaling factor

$`c_0`$ $`=`$ $`{\displaystyle \frac{m+d}{m+2d}}{\displaystyle \frac{1}{c_M}}={\displaystyle \frac{1}{1+2d\kappa }}{\displaystyle \frac{1}{c_M}};c_M1.`$ (2.15)
Explicit formulae are given in App. A.4.

The critical Wilson parameter is known in the weak coupling limit

$`\kappa _c(\beta \mathrm{})={\displaystyle \frac{1}{4}}`$ (2.16)
and expected to be larger for lower $`\beta `$. It is still an open question whether simulations above $`\kappa _c`$ are possible. As this was not the focus of this work, we did not investigate this point. Preliminary results (e.g. Sec. 6.1) suggest that reasonable results are difficult but possible in this regime. The second order phase transition point (and thus the continuum limit) is reached for $`\beta \mathrm{}`$ with lines of constant physics fixed via the demand to keep certain (e.g. mass) ratios constant.

Observables
A measurement of the average plaquette $`\frac{1}{LT}_xU_{Px}`$ and of the temporal and spatial Polyakov loops averaged over the orthogonal direction $`P_L=\frac{1}{T}_{x_2=1}^T_{x_1=1}^LU_{x,1}`$ and $`P_T=\frac{1}{L}_{x_1=1}^L_{x_2=1}^TU_{x,2}`$ is directly possible from the definitions.

Averages of local fermion bilinears $`\overline{\psi }\mathrm{\Gamma }\psi `$ can be measured applying a noisy estimator scheme, using random spinors $`\eta (x)`$ as sources for the solver as shown in App. A.5. One expects the pseudoscalar density (also called the chiral condensate) $`\overline{\psi }\gamma ^5\psi `$ and the further densities $`\overline{\psi }\gamma ^\mu \psi `$ to vanish within errors. The fermion condensate $`\overline{\psi }\psi `$ remains non-zero even in the massless limit because Wilson fermions break chiral invariance explicitly \[Wil74\].

For the determination of meson masses operators with various quantum numbers are defined

* flavour triplet - $`(\overline{\psi }\gamma ^5\tau \psi )`$$`\pi `$’, $`(\overline{\psi }\tau \psi )`$$`a_0`$’,
* flavour singlet - $`(\overline{\psi }\gamma ^5\psi )`$$`\eta `$’, $`(\overline{\psi }\psi )`$$`f_0`$’ .

Moreover, insertion of a $`\gamma ^0`$, e.g. $`(\overline{\psi }\gamma ^0\gamma ^5\tau \psi )`$ for the $`\pi `$, leads to alternative operators with the same quantum numbers in the rest frame. This exhausts the Dirac algebra in two dimensions. The result is summed up in Tab. A.1.

The calculation of temporal correlators $`C^{\mathrm{\Gamma }T}(\mathrm{\Delta }t)`$ involves point sources at randomly chosen positions $`(x,t)`$ and summation over spatial displacements $`y`$ in the other time slice $`(y,t+\mathrm{\Delta }t)`$ to project out the zero momentum states. As to the flavour–singlet channels, the disconnected piece is subtracted with the aid of the noisy inversions performed earlier for the measurement of the condensate \[Oka95\]. For more details consult App. A.5.

All measurements are analysed taking into account covariances and auto-correlations of the correlators \[Bra92\]. A binning method is applied to reduce data size.

The scaling idea we are using is based on the following strategy. Imagine that $`m_\pi `$ and the ratio $`\frac{m_\pi }{m_\eta }`$ are fixed, e.g. by experiment. Working at a certain value of $`\beta `$, we will then have to adjust $`\kappa `$ such that the ratio is the physical value, while $`m_\pi `$ sets the scale, determining the lattice spacing in physical units. The size of the lattice has to be chosen so that finite size effects are small.

Repeating this procedure at larger $`\beta `$ results in a smaller pion mass in lattice units, which again sets the scale for the lattice spacing and thus determines the scaling factor achieved.

We remark that motivated by the continuum results we expect the masses to scale like $`m_\eta 1/\sqrt{\beta }`$, so that we have to increase the value of $`\beta `$ enormously if we want to double the lattice size and reduce the lattice spacing by a factor of 2, keeping physical ratios constant by adjusting the $`\kappa `$ value. This unfortunately limits the usefulness of the Schwinger model, as very high $`\beta `$ simulations encounter topological metastabilities \[Joo90, Dil93, Els97a\].

## Chapter 3 The standard way: Hybrid Monte Carlo

We are ultimately interested in deciding whether the new local bosonic algorithm can be faster than the “workhorse” Hybrid Monte Carlo (HMC) \[Dua87\] used for more than 10 years. We thus implemented a HMC code working with the Hermitean matrix $`Q`$ to set the scale. A brief description of the formulae is given in Sec. 3.1.

In order to make a fair comparison we ulilize the optimisation procedures now common in the application of the HMC. Namely, the implementation includes trajectory length set by $`n\mathrm{\Delta }\tau =1`$ and acceptance $`70\%`$. We further use the fact that extending the solver solution in the trajectory via the condition $`\stackrel{\mathrm{}}{x}=0`$ to generate an optimal new start vector for the inverter results in a gain of approximately 20%.

The inverter algorithm applied throughout is a standard Conjugent Gradient algorithm (CG) \[Ral72, Sto90\] used with precision $`10^6`$. The decision not to use advanced inverters like BiCGgamma5 for the observables was taken in order to have a guaranteed convergence for all simulation parameters. In order to avoid trouble with bad pseudo-random numbers, we use Lüscher’s high-quality random number generator (RG) \[Lue93\], which has been shown to have a very long recursion period and efficient decorrelation properties. It is used in a vectorized form applying 257 generators in parallel.

The standard way to get portable results for CPU cost analyses is to count in matrix multiplication operations. These are the most CPU-time intensive parts of any realistic QCD calculation. We describe our explicit counting in Sec. 3.2.

Generally, one uses this information to define the cost of an algorithm by quantifying the cost needed to generate two independent (i.e. decorrelated) configurations on which measurements are executed. Assuming an exponential decorrelation of configurations \[Bra92\] an estimator for the distance of independent configurations is given by twice the integrated auto-correlation time $`\tau _{\mathrm{int}}1`$. Thus this quantity measured in units of molecular dynamics time or respectively matrix multiplications,

$`2\tau _{\mathrm{int}}[\mathrm{Q}\mathrm{ops}]=2\tau _{\mathrm{int}}N_{\mathrm{Q}\mathrm{ops}/\mathrm{update}},`$ (3.1)
can be regarded as a reasonable estimator for the CPU cost of an algorithm. We would like to point out that although the autocorrelation time is dependent on the observable, it is generally of the same order of magnitude characterizing the update algorithm.

### 3.1 Implementation

The idea of HMC is based upon rewriting the partition function given by

$`Z={\displaystyle D[U]\mathrm{det}Q^2e^{S_G(U)}}`$ (3.2)
using complex pseudo-fermion fields (neglecting constant factors)

$`Z={\displaystyle D[U]D[\eta ]D[\eta ^{}]e^{S_G(U)\eta ^{}Q^2\eta }}`$ (3.3)
and introducing real valued momenta $`\pi _{x,\mu }`$ conjugate to the links $`U_{x,\mu }`$ with a kinetic term which is not changing the path integral

$`Z={\displaystyle D[U]D[\eta ]D[\eta ^{}]D[\pi ]e^{S_G(U)\eta ^{}Q^2\eta \frac{1}{2}\pi ^2}}.`$ (3.4)
The thus defined action is called in analogy to classical mechanics the Hamiltonian of the system

$`H=S_G(U)+\eta ^{}Q^2\eta +{\displaystyle \frac{1}{2}}\pi ^2.`$ (3.5)
Applying this distribution, the complete dynamical update of the link configuration is done in a three step procedure

1. * a complex Gaussian update for the field $`\chi =Q^1\eta `$
    * a real Gaussian update for the field $`\pi `$
2. micro-canonical steps via discretised but reversible leapfrog integration of the equations of motion in fictitious computer time $`\tau `$ leaving the Hamiltonian approximately constant

$`\dot{H}={\displaystyle \frac{H}{\tau }}0`$ (3.6)
3. a Metropolis acceptance step correcting for discretisation errors in the leapfrog integration with

$`\mathrm{min}\{1,e^{\mathrm{\Delta }H}\}`$ (3.7)

#### 3.1.1 Micro-canonical update

We need the update equations, also called equations of motion (EOMs), for the fields $`U_{x,\mu }`$ and $`\pi _{x,\mu }`$. The Hamiltonian EOMs can be chosen to fulfil the two requirements

* $`U_{x,\mu }`$ should stay in $`U(1)`$,
* the Hamiltonian is constant in computer time.

yielding a simple update process for the links

$`\dot{U}_{x,\mu }=i\pi _{x,\mu }U_{x,\mu }`$ (3.8)
while the equations for the momenta are fixed via

$`\dot{H}=0.`$ (3.9)
Writing this out, we find

$`\dot{H}=\dot{S}_G(U)+\eta ^{}{\displaystyle \frac{}{\tau }}[Q^2]\eta +\pi \dot{\pi }.`$ (3.10)
Our aim is to express $`\dot{S}_G(U)`$ and $`\eta ^{}\frac{}{\tau }[Q^2]\eta `$ in terms of the derivatives $`\dot{U}_{x,\mu }`$ and $`\dot{\pi }_{x,\mu }`$ to allow numerical integration of the EOMs. We rewrite the second term

$`\eta ^{}{\displaystyle \frac{}{\tau }}\left[Q^2\right]\eta =\eta ^{}Q^2\left[\dot{Q}Q+Q\dot{Q}\right]Q^2\eta `$ (3.11)
and introduce abbreviations for parts independent of the derivatives

$`\omega =Q^2\eta \mathrm{and}\xi =Q\omega `$ (3.12)
to rewrite

$`\dot{H}=\dot{S}_G(U)2\mathrm{R}\mathrm{e}[\xi ^{}\dot{Q}\omega ]+\pi \dot{\pi }.`$ (3.13)
We now calculate the dependency of the first term $`\dot{S}_G`$ of the derivatives

$`\dot{S}_G`$ $`=`$ $`\beta {\displaystyle \underset{x}{}}\mathrm{Re}(\dot{U}_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}+U_{x,1}\dot{U}_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}`$ (3.14)
$`+U_{x,1}U_{x+\widehat{1},2}\dot{U}_{x+\widehat{2},1}^{}U_{x,2}^{}+U_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}\dot{U}_{x,2}^{}).`$
Using the EOMs we can rewrite

$`\dot{S}_G`$ $`=`$ $`\beta {\displaystyle \underset{x}{}}\mathrm{Re}(i\pi _{x,1}U_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}+U_{x,1}i\pi _{x+\widehat{1},2}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}`$ (3.15)
$`U_{x,1}U_{x+\widehat{1},2}i\pi _{x+\widehat{2},1}U_{x+\widehat{2},1}^{}U_{x,2}^{}U_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}i\pi _{x,2}U_{x,2}^{})`$
and reorder

$`\dot{S}_G`$ $`=`$ $`\beta {\displaystyle \underset{x}{}}\pi _{x,1}\mathrm{Re}i(U_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}U_{x\widehat{2},1}U_{x\widehat{2}+\widehat{1},2}U_{x,1}^{}U_{x\widehat{2},2}^{})`$ (3.16)
$`+\pi _{x,2}\mathrm{Re}i(U_{x\widehat{1},1}U_{x,2}U_{x\widehat{1}+\widehat{2},1}^{}U_{x\widehat{1},2}^{}U_{x,1}U_{x+\widehat{1},2}U_{x+\widehat{2},1}^{}U_{x,2}^{}))`$
$`=`$ $`{\displaystyle \underset{x,\mu }{}}\pi _{x,\mu }\mathrm{Re}\mathrm{\Phi }_{x,\mu }`$
defining thus another abbreviation $`\mathrm{\Phi }`$ independent on the derivatives

$`\mathrm{\Phi }_{x,1}=i\beta (U_{Px}U_{Px\widehat{2}})\mathrm{and}\mathrm{\Phi }_{x,2}=i\beta (U_{Px\widehat{1}}U_{Px}).`$ (3.17)
The dependence of $`\dot{Q}`$ on the derivatives is given by

$`\dot{Q}`$ $`=`$ $`c_o\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\mu ,y}\gamma ^5(1+\gamma ^\mu )\dot{U}_{x\mu ,\mu }+\delta _{x+\mu ,y}\gamma ^5(1\gamma ^\mu )\dot{U}_{x,\mu }^{}\right)`$
$`=`$ $`c_o\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\mu ,y}\gamma ^5(1+\gamma ^\mu )i\pi _{x\mu ,\mu }U_{x\mu ,\mu }\delta _{x+\mu ,y}\gamma ^5(1\gamma ^\mu )i\pi _{x,\mu }U_{x,\mu }^{}\right).`$
In the Hamiltonian this is used in the combination

$`2\mathrm{R}\mathrm{e}[\xi ^{}\dot{Q}\omega ]`$ $`=`$ $`2c_0\kappa \mathrm{Re}\left[i{\displaystyle \underset{x,\mu }{}}\pi _{x,\mu }\left(U_{x,\mu }\xi _{x+\widehat{\mu }}^{}\gamma ^5(1+\gamma ^\mu )\omega _xU_{x,\mu }^{}\xi _x^{}\gamma ^5(1\gamma ^\mu )\omega _{x+\widehat{\mu }}\right)\right]`$ (3.18)
$`=`$ $`{\displaystyle \underset{x,\mu }{}}\pi _{x,\mu }\mathrm{Re}\mathrm{\Omega }_{x,\mu }`$
defining the abbreviation $`\mathrm{\Omega }`$. The final Hamiltonian formula is therefore

$`\dot{H}={\displaystyle \underset{x,\mu }{}}\pi _{x,\mu }(\mathrm{Re}\mathrm{\Phi }_{x,\mu }+\mathrm{Re}\mathrm{\Omega }_{x,\mu }+\dot{\pi }_{x,\mu })`$ (3.19)
giving the required EOMs

$`\mathrm{Re}\{\mathrm{\Phi }_{x,\mu }+\mathrm{\Omega }_{x,\mu }\}+\dot{\pi }_{x,\mu }=0.`$ (3.20)
#### 3.1.2 Leapfrog integration

Heuristically one finds optimal behaviour of the algorithm if the EOMs

$`\dot{U}_{x,\mu }`$ $`=`$ $`i\pi _{x,\mu }U_{x,\mu }`$
$`\dot{\pi }_{x,\mu }`$ $`=`$ $`\mathrm{Re}\{\mathrm{\Phi }_{x,\mu }+\mathrm{\Omega }_{x,\mu }\}`$ (3.21)
are integrated for $`N`$ steps of length $`\mathrm{\Delta }\tau `$ under the constraint

$`N\mathrm{\Delta }\tau 1.`$ (3.22)
Integrating the link EOMs, we find

$`U_{x,\mu }(\tau _0+\delta \tau )=U_{x,\mu }(\tau _0)e^{i\pi _{x,\mu }\delta \tau }`$ (3.23)
to explicitly stay in the group manifold and

$`\pi _{x,\mu }(\tau _0+\delta \tau )=\pi _{x,\mu }(\tau _0)(\mathrm{Re}\{\mathrm{\Phi }_{x,\mu }+\mathrm{\Omega }_{x,\mu }\})\delta \tau .`$ (3.24)
Expressing this in steps $`n=1\mathrm{}N`$ in the link case and $`n=1\mathrm{}N+1`$ in the momentum case, we can write the result in a computer-friendly way

$`U_{x,\mu }(n)`$ $`=`$ $`U_{x,\mu }(n1)e^{i\pi _{x,\mu }(n)\mathrm{\Delta }\tau }`$ (3.25)
$`\pi _{x,\mu }(n)`$ $`=`$ $`\pi _{x,\mu }(n1)(\mathrm{Re}\{\mathrm{\Phi }_{x,\mu }(n1)+\mathrm{\Omega }_{x,\mu }(n1)\})\mathrm{\Delta }\tau (n1)`$
with

$`\mathrm{\Delta }\tau (0)`$ $`=`$ $`\mathrm{\Delta }\tau (N)={\displaystyle \frac{\mathrm{\Delta }\tau }{2}}`$
$`\mathrm{\Delta }\tau (n)`$ $`=`$ $`\mathrm{\Delta }\tau n\{1\mathrm{}N1\}.`$ (3.26)
### 3.2 Counting in $`Q`$ operations

We give the number of $`Q`$ matrix multiplications necessary for one update step of our HMC implementation. In the formulae we denote with $`N_{CG}`$ the number of matrix multiplications needed by the Conjugent Gradient routine applying $`Q^2`$.

The Hybrid Monte-Carlo algorithm applies

* a number $`N_{\mathrm{it}}`$ of leapfrog iterations per trajectory
* a CG routine to invert the matrix $`Q^2`$ using $`N_{\mathrm{CG}}`$ iterations to invert

and thus uses

* $`2N_{\mathrm{it}}N_{\mathrm{CG}}`$ Q operations for the main inverter,
* $`2N_{\mathrm{CG}}`$ Q operations for the initial inverter.
* The final inverter can be omitted as the data is known from the trajectory.
* $`2N_{\mathrm{it}}`$ effective Q operations to update the links

so that we finally obtain a total number of matrix multiplications of

$`N_{\mathrm{Q}\mathrm{ops}/\mathrm{update}}=2N_{\mathrm{it}}N_{\mathrm{CG}}+2N_{\mathrm{CG}}+2N_{\mathrm{it}}.`$ (3.27)
## Chapter 4 The new way: Local bosonic algorithm

### 4.1 Basic idea

As an alternative to the Hybrid Monte Carlo algorithm, M. Lüscher proposed a local bosonic formulation \[Lue94\]. The one main variant of the local bosonic algorithm (LBA) we consider uses the fact that we are dealing with a rescaled squared Hermitean fermion matrix, so that the eigenvalues are real and between 0 and 1.

We can then identically rewrite the path integral

$`𝒵={\displaystyle D[U]det[Q^2]P_n(det[Q^2])\left(P_n(det[Q^2])\right)^1e^{S_g[U]}}`$ (4.1)
using for $`P_n(s)`$ any real polynomial of degree $`n`$ which approximates $`1/s`$ in $`(0,1]`$. Then

$`det[Q^2P_n(Q^2)]=det[1R]1`$ (4.2)
is a small correction which can be treated via a correction scheme, provided we can sample the remaining distribution.

We now assume that the approximation polynomial is of the form

$`P_n(Q^2)=N_{\mathrm{norm}}{\displaystyle \underset{k=1}{\overset{n}{}}}(Q^2z_k),`$ (4.3)
where $`N_{\mathrm{norm}}`$ is a normalisation constant and the $`z_k`$ are the roots of the polynomial coming in complex conjugate pairs $`z_k,\overline{z}_k`$ with non-vanishing imaginary part. They define $`\mu _k,\nu _k`$ and $`\overline{\mu }_k,\overline{\nu }_k`$ via

$`\sqrt{z_k}=\mu _k+i\nu _k`$ $``$ $`\overline{\mu }_k=\mu _k\overline{\nu }_k=\nu _k`$ (4.4)
if one defines the branch of the square root to be taken by $`\nu _k>0`$. Using this we can rewrite

$`P(Q^2)=N_{\mathrm{norm}}{\displaystyle \underset{k=1}{\overset{n}{}}}(Q^2z_k)=N_{\mathrm{norm}}{\displaystyle \underset{k=1}{\overset{n}{}}}\left[(Q\mu _k)^2+\nu _k^2\right].`$ (4.5)
As this factorises, sampling this distribution can be done using

$`{\displaystyle D[\varphi ]D[\varphi ^{}]e^{\frac{1}{2}\varphi ^{}A\varphi }}(det[A])^1`$ (4.6)
to transform the distribution to a bosonic integral

$`P_{\mathrm{eff}}(U)`$ $``$ $`det[1R]e^{S_g[U]}{\displaystyle 𝒟\varphi e^{_k[\varphi _k^{}(Q\mu _k)^2\varphi _k+\nu _k^2\varphi _k^{}\varphi _k]}}`$ (4.7)
with $`n`$ complex bosonic Dirac fields $`\varphi _k`$.

The main criteria for the polynomials $`P_n(s)`$ is that they should be fast convergent to $`\frac{1}{s}`$ in $`(0,1]`$ and have roots coming in complex conjugate pairs. We chose as approximation polynomials the Chebyshev-derived polynomials proposed by Bunk et al.\[Bun95a\]

$`P_n(s)={\displaystyle \frac{1+\rho _nT_{n+1}(x)}{s}}={\displaystyle \frac{1R_n(s)}{s}}`$ (4.8)
with

$`x=2{\displaystyle \frac{sϵ}{1ϵ}}1,`$ (4.9)
such that $`s[ϵ,1]`$ is mapped to $`x[1,1]`$. The polynomials $`T_n(z)`$ of degree $`n0`$ are the standard Chebyshev polynomials given by

$`T_n(z)`$ $`=`$ $`\mathrm{cosh}n\varphi \mathrm{with}z=\mathrm{cosh}\varphi .`$ (4.10)
For completeness, we give some of the first polynomials

$`T_0(z)=1,T_1(z)=z,T_2(z)=2z^21,`$
$`T_3(z)=4z^33z,T_4(z)=8z^48z^2+1.`$ (4.11)
An important feature is the Clenshaw recursion relation

$`T_{n+1}(z)+T_{n1}(z)=2zT_n(z)\text{for}n1,`$ (4.12)
known to be numerically stable. For more information we defer the reader to \[Fox68\].

The polynomials $`P_n(s)`$ approximate $`1/s`$ for real $`s[ϵ,1]`$. The rest term $`R_n(s)`$ quantifies the quality of the approximation and the normalisation factor $`\rho _n`$ is defined through the condition that $`P_n(s)`$ exists even at $`s=0`$ via

$`0=sP_n(s)|_{s=0}=1+\rho _nT_{n+1}\left({\displaystyle \frac{1+ϵ}{1ϵ}}\right)\rho _n={\displaystyle \frac{1}{T_{n+1}(\frac{1+ϵ}{1ϵ})}}.`$ (4.13)
The explicit expression for the error term

$`R(s)={\displaystyle \frac{T_{n+1}(\frac{2s}{1ϵ}\frac{1+ϵ}{1ϵ})}{T_{n+1}(\frac{1+ϵ}{1ϵ})}}`$ (4.14)
can be used to deduce limits for the approximation quality

$`|R_n(s)|`$ $``$ $`|\rho _n|=|{\displaystyle \frac{1}{\mathrm{cosh}[(n+1)\chi ]}}|\mathrm{with}\mathrm{cosh}\chi ={\displaystyle \frac{1+ϵ}{1ϵ}}`$ (4.15)
yielding

$`|R_n(s)|`$ $``$ $`2{\displaystyle \frac{1}{\mathrm{cosh}[(n+1)\chi ]+\mathrm{sinh}[(n+1)\chi ]}}=2\left({\displaystyle \frac{1\sqrt{ϵ}}{1+\sqrt{ϵ}}}\right)^{n+1}2e^{2\sqrt{ϵ}n}`$ (4.16)
proving an exponential convergence to $`1/s`$ for $`s[ϵ,1]`$. Explicitly, we use as the approximation quality factor \[Lue94, Bun95a\]

$`\delta =2\left({\displaystyle \frac{1\sqrt{ϵ}}{1+\sqrt{ϵ}}}\right)^{n+1}.`$ (4.17)
Not necessary for the application in the local bosonic algorithm, but convenient is the fact that the roots of this approximation polynomial can be given explicitly observing that

$`P_n(s)=0`$ $``$ $`T_{n+1}(x)=T_{n+1}({\displaystyle \frac{1+ϵ}{1ϵ}})`$ (4.18)
for complex $`s,x`$. This yields

$`z_k={\displaystyle \frac{1}{2}}(1+ϵ){\displaystyle \frac{1}{2}}(1+ϵ)\mathrm{cos}({\displaystyle \frac{2\pi k}{n+1}})i\sqrt{ϵ}\mathrm{sin}({\displaystyle \frac{2\pi k}{n+1}}),`$ (4.19)
i.e. $`n`$ roots on an ellipse in the complex plane around the foci $`ϵ`$ and $`1`$ on the real axis.

We further can deduce the global normalisation constant $`N_{\mathrm{norm}}`$ from the symmetry of the Chebyshev polynomial in $`[ϵ,1]`$

$`T_{n+1}\left({\displaystyle \frac{1+ϵ}{2}}\right)=0{\displaystyle \frac{1+ϵ}{2}}P_n\left({\displaystyle \frac{1+ϵ}{2}}\right)1=0`$ (4.20)
to

$`N_{\mathrm{norm}}^1={\displaystyle \frac{1+ϵ}{2}}{\displaystyle \underset{k=1}{\overset{n}{}}}\left({\displaystyle \frac{1+ϵ}{2}}z_k\right).`$ (4.21)
The updating process consists of exact heatbath sweeps for the $`\varphi `$’s and $`U`$’s \[Bes79\], followed by a number of over-relaxation iterations. The formulae for the updates are given in App. C. We confirmed in preliminary runs that the reflection sweeps for the $`\varphi `$’s and $`U`$’s have to be combined in pairs, as was observed before \[Jeg95, Jan96\]. In total, this introduces the parameters $`n`$, $`ϵ`$ and the number of reflections per heatbath into the algorithm.

### 4.2 Correction step

Finally, the approximation $`det[1R]1`$ has to be controlled. This in principle can either be done via a Metropolis accept/reject step with correction probability

$`P_{U\varphi U^{}\varphi ^{}}^C=\mathrm{min}[1,{\displaystyle \frac{det[1R^{}]}{det[1R]}}]`$ (4.22)
or by reweighting

$`O={\displaystyle \frac{Odet[1R]_L}{det[1R]_L}},`$ (4.23)
where $`\mathrm{}_L`$ denotes the expectation values in the local bosonic ensemble and primed quantities are of the new configuration.

We would like to point out that algorithms including an acceptance step need the reversibility of the update trajectory to guarantee detailed balance. This leads to a symmetrization of the update trajectories in these cases, which has to be included in the CPU cost derivation.

On the other hand, the ensembles generated by the local bosonic distribution with and without correction are distinctly different. In simulations $`det[1R]_L`$ can in fact be far from 1 as will be shown in Tab. 6.5. The integrated autocorrelation time can in these cases be small, while the error from each measurement grows to very large values. We thus can no longer use the autocorrelation time measured in matrix multiplications

$`2\tau _{\mathrm{int}}[\mathrm{Q}\mathrm{ops}]`$ (4.24)
described in Ch. 3 as a good transportable indicator for the CPU cost of the algorithms. We have to regress to the basic effective cost defined by

$`C_{\mathrm{eff}}=N_{\mathrm{total}\mathrm{Q}\mathrm{ops}}{\displaystyle \frac{\sigma _{\mathrm{tot}}^2(A)}{<A>^2}}`$ (4.25)
applicable to all measurements as only the total work and the real result (i.e. the error obtained for a certain observable) are used. The cost factor is normalised by the mean value of the observable to cancel the trivial dependency. We remark that this cost factor is even more observable-dependent than the standard choice $`\tau _{\mathrm{int}}`$, prohibiting even a rough comparison of different observables.

In a first step we approximatively included the correction via the lowest 8 eigenvalues of $`Q^2`$, using them to apply a global Metropolis correction step \[Ale95\]. Recently, stochastic methods which are exact were suggested \[Bor96\]. Some variants of these will be discussed in the next sections.

#### 4.2.1 Reweighting

The correction factor can be treated exactly using a stochastic method. Following a suggestion by Lüscher \[Lue96\], we rewrite

$`O={\displaystyle \frac{Odet[1R]_L}{det[1R]_L}}`$ (4.26)
into

$`O={\displaystyle \frac{[d\eta ]Oe^{\eta ^{}[1R]^1\eta }_L}{[d\eta ]e^{\eta ^{}[1R]^1\eta }_L}}={\displaystyle \frac{[d\eta ]Oe^{\eta ^{}[1R]^1\eta +\eta ^{}\eta \eta ^{}\eta }_L}{[d\eta ]e^{\eta ^{}[1R]^1\eta +\eta ^{}\eta \eta ^{}\eta }_L}}.`$ (4.27)
We then interpret the factor $`e^{\eta ^{}\eta }`$ as a Gaussian weight, resulting in

$`O={\displaystyle \frac{Oe^{\eta ^{}\left(1[1R]^1\right)\eta }_{L,\eta }}{e^{\eta ^{}\left(1[1R]^1\right)\eta }_{L,\eta }}}={\displaystyle \frac{Oe^{\eta ^{}R[1R]^1\eta }_{L,\eta }}{e^{\eta ^{}R[1R]^1\eta }_{L,\eta }}},`$ (4.28)
where $`\mathrm{}_{L,\eta }`$ denotes the ensemble generated by the LBA distribution and the Gaussian noise vectors. Note that rewriting

$`1(1R)^1R(1R)^1`$ (4.29)
is numerically more stable as it avoids subtracting two almost identical objects.

#### 4.2.2 Noisy acceptance step method I

The idea behind this noisy method \[Bor96\] is to generate a new configuration $`(U^{},\varphi ^{})`$ with the standard Lüscher action and then to accept it according to a correction probability that is a function of a noise vector $`\eta `$

$`P_{(U,\varphi )(U^{},\varphi ^{})}^C(\eta )`$ (4.30)
satisfying detailed balance on average over the $`\eta `$ distribution

$`{\displaystyle \frac{[d\eta ]P_{(U,\varphi )(U^{},\varphi ^{})}^C(\eta )}{[d\stackrel{~}{\eta }]P_{(U^{},\varphi ^{})(U,\varphi )}^C(\stackrel{~}{\eta })}}`$ $`=`$ $`{\displaystyle \frac{det[1R^{}]}{det[1R]}}.`$ (4.31)
In order to do this one generates random numbers according to the distribution $`P^{HB}`$ which is constructed to give the required determinant factor. The acceptance of the Metropolis $`P^A`$ is fixed so that the total acceptance

$`P_{(U,\varphi )(U^{},\varphi ^{})}^C(\eta )=P_\eta ^{HB}P_{(U,\varphi )(U^{},\varphi ^{})}^A(\eta )`$ (4.32)
satisfies detailed balance. The almost obvious choice for the heatbath distribution

$`P_\eta ^{HB}={\displaystyle \frac{N_\mathrm{n}}{det[1R]}}e^{\eta ^{}[1R]^1\eta }`$ (4.33)
with $`N_\mathrm{n}`$ a normalisation constant requires an acceptance

$`P_{(U,\varphi )(U^{},\varphi ^{})}^A(\eta )=\mathrm{min}(1,e^{\eta ^{}[1R^{}]^1\eta +\eta ^{}[1R]^1\eta }).`$ (4.34)
The proof of detailed balance can be seen from the symmetrical behaviour of the involved integrals

$`{\displaystyle \frac{[d\eta ]P_{(U,\varphi )(U^{},\varphi ^{})}^C(\eta )}{[d\stackrel{~}{\eta }]P_{(U^{},\varphi ^{})(U,\varphi )}^C(\stackrel{~}{\eta })}}={\displaystyle \frac{[d\eta ]\frac{N_\mathrm{n}}{det[1R]}e^{\eta ^{}[1R]^1\eta }\mathrm{min}(1,e^{\eta ^{}[1R^{}]^1\eta +\eta ^{}[1R]^1\eta })}{[d\stackrel{~}{\eta }]\frac{N_\mathrm{n}}{det[1R^{}]}e^{\stackrel{~}{\eta }^{}[1R^{}]^1\stackrel{~}{\eta }}\mathrm{min}(1,e^{\stackrel{~}{\eta }^{}[1R]^1\stackrel{~}{\eta }+\stackrel{~}{\eta }^{}[1R^{}]^1\stackrel{~}{\eta }})}}`$
$`{\displaystyle \frac{det[1R^{}]}{det[1R]}}{\displaystyle \frac{[d\eta ]\mathrm{min}(e^{\eta ^{}[1R]^1\eta },e^{\eta ^{}[1R^{}]^1\eta })}{[d\stackrel{~}{\eta }]\mathrm{min}(e^{\stackrel{~}{\eta }^{}[1R^{}]^1\stackrel{~}{\eta }},e^{\stackrel{~}{\eta }^{}[1R]^1\stackrel{~}{\eta }})}}={\displaystyle \frac{det[1R^{}]}{det[1R]}}.`$ (4.35)
In order to make the algorithm practically working, we transform the heatbath distribution $`P^{HB}`$ into a simple Gaussian using the substitution

$`\eta =B\chi `$ (4.36)
with $`B`$ and $`B^{}`$ given by the constraint

$`B^{}B=1R.`$ (4.37)
We remark that while this transformation can be done trivially in the non-Hermitean case \[Bor96\], we here need a little trick. The decomposition can be done easily using the product formula

$`1R=N_{\mathrm{norm}}Q^2{\displaystyle \underset{k}{\overset{n/2}{}}}(Q^2z_k)(Q^2\overline{z}_k)`$ (4.38)
taking one root of each complex conjugate pair to yield

$`B=\sqrt{N_{\mathrm{norm}}}Q{\displaystyle \underset{k}{\overset{n/2}{}}}(Q^2z_k).`$ (4.39)
Finally we obtain

$`P_\chi =N_\mathrm{n}e^{\chi ^{}B^{}[BB^{}]^1B\chi }=N_\mathrm{n}e^{\chi ^{}\chi }`$ (4.40)
and

$`P_{(U,\varphi )(U^{},\varphi ^{})}^A`$ $`=`$ $`\mathrm{min}(1,e^{\chi ^{}B^{}[1R^{}]^1B\chi +\chi ^{}\chi }).`$ (4.41)
#### 4.2.3 Noisy acceptance step method I with adapted precision

This method makes use of a possible optimisation for Metropolis-type acceptance step schemes using a restartable solver like the Conjugent Gradient chosen for this study.

As the random number governing the Metropolis decision is known before the solver is applied, it is in principle possible to interrupt the solver iterations and check whether the quality of the solver solution is already good enough to fulfil the requirements of a clear-cut decision, i.e. to distinguish the result from the chosen random number. If not, the solver is restarted and this procedure iterated.

We implemented this idea in a simplified way. We demand a very limited solver precision of $`10^2`$ in the first step and check if this quality is good enough. If not, we in a second step demand full solver precision of $`10^6`$.

#### 4.2.4 Noisy acceptance step method II

Knowing of the idea of using Gegenbauer polynomials \[Bun97\] to solve equation of the type

$`A^{\frac{1}{2}}x=y`$ (4.42)
as explained in App. E, we want to rewrite the formulae of the noisy Metropolis update section, resulting not in the acceptance Eq. 4.41 but in the corresponding

$`P_{(U,\varphi )(U^{},\varphi ^{})}^A`$ $`=`$ $`\mathrm{min}(1,e^{\chi ^{}[1R^{}]^{\frac{1}{2}}(1R)[1R^{}]^{\frac{1}{2}}\chi +\chi ^{}\chi }),`$ (4.43)
where now we effectively do not need the square root of $`Q^2P(Q^2)`$ but the inverse square root.

It is again easy to see the idea. One starts with a heatbath distribution including a determinant factor (but now in the numerator)

$`P_\eta ^{HB}=N_\mathrm{n}det[1R^{}]e^{\eta ^{}(1R^{})\eta }.`$ (4.44)
The acceptance needed in this case is

$`P_{(U,\varphi )(U^{},\varphi ^{})}^A(\eta )=\mathrm{min}(1,e^{\eta ^{}(1R)\eta +\eta ^{}(1R^{})\eta }).`$ (4.45)
We skip the proof of detailed balance as it is identical to the method I case.

We apply the same transformation trick eliminating the Gaussian vector $`\eta `$ via the substitution

$`\eta =[1R^{}]^{\frac{1}{2}}\chi `$ (4.46)
given by the application of the Gegenbauer solver mentioned above.

Thus we again constructed a purely Gaussian distribution for the $`\chi `$ variables

$`P_\chi =N_\mathrm{n}e^{\chi ^{}[1R^{}]^{\frac{1}{2}}(1R^{})[1R^{}]^{\frac{1}{2}}\chi }=N_\mathrm{n}e^{\chi ^{}\chi },`$ (4.47)
yet the acceptance is now

$`P_{(U,\varphi )(U^{},\varphi ^{})}^A`$ $`=`$ $`\mathrm{min}(1,e^{\eta ^{}(RR^{})\eta }).`$ (4.48)
This formula makes it evident that in this scheme the main work is done by the Gegenbauer solver constructing the $`\eta `$ vectors. Thus its usefulness relies on the fact that the convergence rate of the Gegenbauer solver is equivalent to that of the Conjugent Gradient \[Bun97\]. In Ch. 8 we will give results of the first simulations with the Gegenbauer solver. Details of polynomials and solver scheme are given in App. E.

We would like to remark that an optimisation of this method analogous to that of the acceptance step method I with adapted precision is possible. The Gegenbauer inverter is not restartable, but a slight alteration storing some shift vectors in the iteration could be implemented without creating too much overhead. We thus could interrupt the solver iteration, check whether the solver quality is already good enough for the requested Metropolis decision, and continue if this is not yet the case. For reasons of limited computer resources this was not included in this study.

### 4.3 Counting in $`Q`$ operations

Our implementation of the Lüscher local bosonic algorithm applies in the corrected measurement case

* $`N_{\mathrm{refl}}`$ reflections
* $`1`$ heatbath step

for the $`n`$ bosonic fields and the link field, while in the acceptance step cases twice that amount is necessary to symmetrize the trajectory. We use an abbreviation

* $`N_{\mathrm{sym}}=1`$ for corrected measurement and
* $`N_{\mathrm{sym}}=2`$ for acceptance step schemes.

For the boson field correction part, both ways have to invert the approximation polynomial $`Q^2P(Q^2)`$ which takes

$`2N_{\mathrm{CG}}(n+1)`$ (4.49)
$`Q`$ operations. Heatbath and reflections for the boson fields both require

* 1 $`Q`$ operation to initialise the stored auxiliary fields
* 1 $`Q`$ operation in the force calculation routine
* 1.5 $`Q`$ operations in the update routine

so that we obtain

$`3.5nN_{\mathrm{sym}}(N_{\mathrm{refl}}+1)+2N_{\mathrm{CG}}(n+1).`$ (4.50)
The link updates take in 2 dimensions for heatbath and reflection each

* 2 $`Q`$ operations for the staples
* 2 $`Q`$ operations for the bosonic force parts

so that we end up with the total number of matrix multiplications necessary for one update step of

$`N_{\mathrm{Q}\mathrm{ops}/\mathrm{update}}=[3.5n+4]N_{\mathrm{sym}}(N_{\mathrm{refl}}+1)+2(n+1)N_{\mathrm{CG}}.`$ (4.51)
## Chapter 5 The problems I: Instabilities

In the update algorithms of Ch. 4 we encounter the basic numerical problem of evaluating a matrix-valued polynomial of high order. To be specific, we will in the following consider the problem evaluating a polynomial approximation of the inverse determinant

$$\mathrm{det}A\left[\mathrm{det}P_n(A)\right]^1,$$
(5.1)
using the Chebyshev polynomials as defined in Ch. 4.1.

The naive idea would be to use the factorized form Eq. 4.3. However, the numerical construction of a polynomial using the product representation, can – due to rounding errors – easily lead to a loss of precision or even to numerical instabilities. This holds in particular if computers with 32-bit floating point precision are used.

Often, as in the case of Chebyshev polynomials, numerically stable recursion relations are available (viz. Ch. 4.1) \[Fox68\]. However, exact versions of the LBA \[Bor96\] or related approaches like the Polynomial Hybrid Monte Carlo (PHMC) algorithm \[For97b, Fre97a\] often need the factorized form of the polynomial $`P_n`$. Especially a decomposition of the polynomial into two or more (e.g. complex conjugate) parts can in general not be done without recursion to the factorized form. This is the numerical problem we are faced with in the Hermitean LBA variant using a Metropolis acceptance step. It is our intention to investigate in these cases several ordering schemes for the complex roots.

### 5.1 Factorized Chebyshev polynomial

Let us consider the Chebyshev approximation of a function $`f(s)`$ depending on a real variable $`s`$ by a polynomial $`P_n(s)`$ of degree $`n`$. The motivation to initially study a single degree of freedom is that we might think of the matrix $`A`$ as being diagonalized. Then the problem, Eq. 5.1, reduces to finding a polynomial that approximates each $`\lambda ^1(A)`$ separately, where $`\lambda (A)`$ is a real eigenvalue of $`A`$. We therefore expect that studying a single degree of freedom can provide information also about the qualitative behaviour of rounding error effects when the matrix valued polynomial $`P_n(A)`$ is numerically computed.

In principle evaluating the partial product $`P_q(s)`$

$$P_q(s)=N_{\mathrm{norm}}^q\underset{k=1}{\overset{q}{}}(sz_k)=\underset{k=1}{\overset{q}{}}n_k(sz_k),$$
(5.2)
one has to define a normalisation $`N_{\mathrm{norm}}^q`$ for the partial products, effectively distributing the global normalisation constant $`N_{\mathrm{norm}}`$ to $`n`$ $`k`$-dependent normalisation constants $`n_k`$.

Considering precision losses or numerical instabilities can be understood by the following argument: The absolute values of the, in general complex, subsequent partial products $`|P_q(s)|`$ and $`|P_{q+1}(s)|`$ can be different by orders of magnitude if $`sz_q`$. If now $`|P_{q+1}(s)||P_q(s)|`$ then this must have been achieved by subtracting two large numbers, which bears the danger of a significant loss of precision. The problem itself suggests, however, its solution: The monomial factors in Eq. 4.3 or equivalently the roots $`z_k`$ should be ordered, if possible, in such a way that the absolute values of all partial products $`|P^q(s)|`$ have the same order of magnitude. Regarding an application to vectors (where $`s`$ is a priori unknown), this has to be achieved in an $`s`$-independent way.

To investigate the effect of reordering, we propose to use a simple criterion to determine the effects of these rounding errors. The idea is to evaluate in a first step for given $`q`$ the maximal and the minimal value of $`|P_{q,ϵ}(s)|`$ over the spectral interval $`0s1`$. The ratio of the maximum to the minimum value, i.e. $`\stackrel{~}{R}_q=\mathrm{max}_{s[0,1]}|P_{q,ϵ}(s)|/\mathrm{min}_{s[0,1]}|P_{q,ϵ}(s)|`$, is then a measure of how large the fluctuations of the absolute values of the partial products can become. If these fluctuations are very large, the polynomial can not be constructed in a safe way. Building the maximum of $`\stackrel{~}{R}_q`$ with respect to $`q`$, we arrive at the final quantitative measure

$$R_{\mathrm{max}}=\underset{q\{1,\mathrm{},n\}}{\mathrm{max}}\left\{\frac{\underset{s[0,1]}{\mathrm{max}}|P_{q,ϵ}(s)|}{\mathrm{min}_{s[0,1]}|P_{q,ϵ}(s)|}\right\}.$$
(5.3)
It is clear that $`R_{\mathrm{max}}`$ has to be smaller than the inverse relative accuracy on a given computer as a necessary condition for the stability of the evaluation of the full polynomial.

Another quantity of interest is the maximum value of the partial products itself

$$M_{\mathrm{max}}=\underset{s[0,1],q\{1,\mathrm{},n\}}{\mathrm{max}}|P_{q,ϵ}(s)|.$$
(5.4)
This has to be be smaller than the largest representable number in order not to run into overflow. To avoid underflow one should thus also study $`M_{\mathrm{min}}`$ by replacing the maximum in Eq. 5.4 by the corresponding minimum. We, however, restrict ourselves to the maximum quantity.

Note that $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ are computed for $`s[0,1]`$, whereas the Chebyshev polynomial has an exponential convergence only in the interval $`s[ϵ,1]`$. However, as will be explicitly demonstrated below, our results for $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ do not depend very much on the choice of the lower end of the interval.

### 5.2 Ordering schemes

In this section we introduce the different ordering schemes used for the roots $`z_k`$ Eq. 4.19. In principle, one could try to also distribute the normalisation constants $`n_k`$ introduced above in a $`k`$-dependent way to reduce rounding errors. However, we found no improvement over the naive homogeneous distribution

$`n_k=(N_{\mathrm{norm}})^{\frac{1}{n}}`$ (5.5)
which we therefore adopt throughout the rest of this work.

Naive ordering
As naive we regard the ordering given by

$$z_k^{\mathrm{naive}}=z_k,k=1,\mathrm{},n$$
(5.6)
with the $`z_k`$ given in Eq. 4.19. The roots $`z_k`$ form an ellipse in the complex plane and in the naive ordering the roots are selected from this ellipse by starting at the origin and moving around anti-clockwise. This is indicated in Fig. 5.1a, where the roots are shown labelled according to the order in which they are used in the evaluation of the Chebyshev polynomial Eq. 4.3. As we will see later, ordering the roots in this naive way gives rise to substantial rounding error effects, even leading to numerical overflow.

Zürich group scheme
Recently, De Forcrand brought to our attention the scheme used for the simulations of the Zürich group \[For97\], which consists of using the complex pairs $`s_k,s_{(n+1)k}`$ as the $`2k1`$ and $`2k`$-th roots. Because of results generally worse than the naive scheme, we decided to simply state that this scheme is not to be recommended, and not to include it in the plots. We would further like to stress that the studies undertaken by the Zürich group were of a kind not influenced by the quality of the ordering scheme.

Pairing scheme
A first improvement over the naive ordering is to use a simple pairing scheme, reordering the roots

$`z_k^{\mathrm{pair}}=z_{j(k)},k=1,\mathrm{},n.`$ (5.7)
We give the reorder index $`j(k)`$ for the example of $`n/2`$ being a multiple of $`4`$ and $`n^{}=n/8`$. In the lower half plane, $`\mathrm{Im}z_k<0`$, the pairing scheme is achieved by

$`j`$ $`=`$ $`\{1,{\displaystyle \frac{n}{2}},{\displaystyle \frac{n}{4}}+1,{\displaystyle \frac{n}{4}},`$ (5.8)
$`2,{\displaystyle \frac{n}{2}}1,{\displaystyle \frac{n}{4}}+2,{\displaystyle \frac{n}{4}}1,`$
$`\mathrm{}`$
$`n^{},{\displaystyle \frac{n}{2}}n^{}+1,{\displaystyle \frac{n}{4}}n^{},{\displaystyle \frac{n}{4}}+n^{}1\}`$
and for $`\mathrm{Im}z_k>0`$ correspondingly. An illustration of the ordering in the pairing scheme is shown in Fig. 5.1b. The label of the roots indicates in which order they are used in the numerical construction of the polynomial.

In case of $`n/2`$ not divisible by 4, we search for the next integer $`m`$ smaller than $`n/2`$ and divisible by 4. We then repeat the above described procedure on these $`m`$ roots and simply multiply the remaining roots $`z_{m+1}\mathrm{}z_{n/2}`$ at the end. To make our procedure explicit, we include the Fortran code to reorder the first half of the roots in App. D.1. The code for the second half of the roots is constructed analogously. The Fortran code also contains the case of $`n/8`$ being odd, where the above described procedure has to be slightly modified.

Subpolynomial scheme
The problem with rounding errors in computing the polynomial in the product representation arises most severely for a high degree $`n`$ of the polynomial. In order to decrease the effects of rounding errors one may therefore be guided by the following intuitive observation. Let us consider the polynomial $`P_{n,ϵ}(s)`$ with roots $`z_k`$ and $`n1`$. If $`m`$ is an integer divisor of $`n`$, the roots $`z_1,z_{1+m},z_{1+2m},\mathrm{}z_{1+(n/m1)m}`$ turn out to be close to the roots characterising the polynomial $`P_{n^{},ϵ}(s)`$ of degree $`n^{}=n/m`$ (note that we keep the same $`ϵ`$). Moreover, the normalisation constants $`c_k=\left(N_{\mathrm{norm}}(n)\right)^{\frac{1}{n}}`$ and $`n_k^{}=\left(N_{\mathrm{norm}}(n^{})\right)^{\frac{1}{n^{}}}`$ are of the same order (the dependence on $`n`$ of $`n_k`$ turns out to be negligible for large $`n`$). Then the product

$$u=\underset{j=0}{\overset{n/m1}{}}[c_{j+1}(sz_{1+jm})]$$
(5.9)
is a good approximation of $`P_{n^{},ϵ}(s)`$, $`|uP_{n^{},ϵ}(s)|1`$ for all $`ϵ<s1`$, and $`|u|ϵ`$. The same argument may be repeated for the other similar sequences of roots, like $`z_2,z_{2+m},z_{2+2m},\mathrm{}z_{2+(n/m1)m}`$, …, $`z_m,z_{2m},z_{3m},\mathrm{}z_n`$.

This means that the product Eq. 4.3 may be split in a product of $`m`$ subproducts, in such a way that each of them approximates the factorized form of a polynomial $`P_{n^{},ϵ}(s)`$ of lower degree $`n^{}=n/m`$. Because of the lower degree of the polynomial given by the products such as Eq. 5.9, one may expect that much smaller fluctuations occur in the intermediate steps of the evaluation of each of these subproducts.

The reordering of the subpolynomial scheme

$`z_k^{\mathrm{sp}}=z_{j(k)};k=1,\mathrm{},n`$ (5.10)
can be represented by

$`j`$ $`=`$ $`\{1,1+m,2+2m,\mathrm{},1+({\displaystyle \frac{n}{m}}1)m,`$ (5.11)
$`2,2+m,2+2m,\mathrm{},2+({\displaystyle \frac{n}{m}}1)m,`$
$`\mathrm{}`$
$`m,m+m,m+2m,\mathrm{},m+({\displaystyle \frac{n}{m}}1)m\},`$
where $`m`$ is an integer divisor of $`n`$.

We found that $`m`$ has to be chosen $`m\sqrt{n}`$ to avoid severe loss of precision in the construction of the polynomial and to reduce rounding error effects to a tolerable level. We remark that the naive ordering is reproduced by the two extreme choices $`m=1`$ and $`m=n`$. The Fortran code which generates the subpolynomial ordering of the roots may be found in App. D.2.

Bitreversal scheme
The subpolynomial scheme can be generalised, leading to what we will call the bitreversal scheme. To illustrate how this scheme works, let us assume that the degree $`n`$ of the polynomial is a power of two. One now starts with the $`n`$ monomial factors in Eq. 4.3, chooses $`m=n/2`$ and applies the subpolynomial scheme resulting in $`m`$ binomial factors. We then proceed in choosing a $`m^{}=m/2`$ and again applying the subpolynomial scheme to these $`m`$ binomial factors which leaves us with $`m^{}/2`$ subpolynomials each of degree four. The procedure can be iterated until we are left with only one subpolynomial having the degree of the polynomial itself. The above sketched procedure can be realized in practise by first representing the integer label (counting from $`0`$ to $`n1`$) of the roots in the naive order by its bit representation. The desired order is then obtained by simply reversing the bits in this representation. The resulting reordering of the roots is shown in Fig. 5.1c with $`n=16`$ as an example.

For $`n`$ not a power of 2, we pad with dummy roots, chosen to be zero for instance, until the artificial number of roots is a power of 2. The bitreversal procedure can then be applied as described above. Afterwards, the dummy roots have to be eliminated from the sequence.

To make the procedure of reordering the roots explicit, we again give the Fortran code used to generate the bitreversal ordering in App. D.3.

Montvay scheme
Recently, Montvay \[Mon97\] suggested to order the roots according to an optimisation procedure which can be implemented numerically, e.g. using algebraic manipulation programs such as Maple. Let us shortly sketch how Montvay’s ordering scheme works and refer to \[Mon97\] for details. We assume that we have already the optimised order of the roots for the partial product $`P_{q<n}(s)`$, Eq. 5.2. Then the values of $`|sP_q(s)(sz)|`$ are computed for all $`z`$ taken from the set of roots not already used. The values of $`s`$ are taken from a discrete set of points, $`\{s_1,\mathrm{},s_N\}`$ which are all in the interval $`[ϵ,1]`$. Now, the maximal ratio over $`s\{s_1,\mathrm{},s_N\}`$ of all values $`|sP_q(s)(sz)|`$ is computed for each root $`z`$ separately. Finally, that root is taken which gives the lowest of these maximal values. Starting with the trivial polynomial $`P_0(s)=1`$, this procedure obviously defines a scheme with which the roots can be ordered iteratively. We show in Fig. 5.1d the resulting order of the roots using Montvay’s scheme by again labelling the roots in the order as they are used to compute the Chebyshev polynomial Eq. 4.3.

Clenshaw recursion
We repeat the statement that the evaluation of a Chebyshev polynomial normally does, of course, not rely on the product representation, Eq. 4.3. A numerically safe way to construct a Chebyshev polynomial is to use a recursion relation, such as the Clenshaw recursion Eq. 4.12 described in Sec. 4.1 \[Fox68\]. The recursion is known to be numerically very stable and will serve us in the following as a reference procedure for the numerical evaluation of a Chebyshev polynomial.

### 5.3 Scalar numerical tests

In this section we will give our numerical results for the quantities $`R_{\mathrm{max}}`$, Eq. 5.3, and $`M_{\mathrm{max}}`$, Eq. 5.4. When these quantities assume large values, the numerical evaluation of the Chebyshev polynomial Eq. 4.3 is affected by rounding errors and precision losses or numerical instabilities can easily be encountered. All results presented in this section are obtained on computers using 64-bit arithmetic. Montvay’s ordering scheme is constructed applying a Maple program \[Mon97a\] requiring 40 digit precision. We want to emphasise at this point that the quantities $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$, which we are using to test the different ordering schemes, do not directly correspond to the values $`|sP_q(s)(sz)|`$, which serve to optimise the root ordering in the Montvay scheme.

In order to compute the values for $`R_{\mathrm{max}}`$, Eq. 5.3, and $`M_{\mathrm{max}}`$ Eq. 5.4 we take $`5000`$ values of $`s`$, equally spaced in the interval $`[0,1]`$. We check explicitly that the values of $`R_{\mathrm{max}}`$ do not depend very much on the lower end of the interval $`[s_{\mathrm{min}},1]`$ from which $`s`$ is taken. In Fig. 5.2 we show $`R_{\mathrm{max}}`$ as a function of the lower end of the interval $`s_{\mathrm{min}}`$ measured in units of the parameter $`ϵ`$. We present data for three polynomial degrees $`n=30`$, $`n=86`$ and $`n=146`$ at fixed approximation quality $`\delta =0.001`$ using the bitreversal ordering scheme. As Fig. 5.2 shows, the dependence of $`R_{\mathrm{max}}`$ on $`s_{\mathrm{min}}`$ is very weak. For the other root ordering schemes we find a very similar behaviour of $`R_{\mathrm{max}}`$ as a function of $`s_{\mathrm{min}}`$. This justifies the use of $`s_{\mathrm{min}}=0`$ which we have used for the numerical tests described in this section.

We start by comparing the subpolynomial and the bitreversal schemes, as they are closely related to each other. In Fig. 5.3 we show $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ as a function of the degree $`n`$ of the Chebyshev polynomial, keeping the maximal fit accuracy $`\delta =0.1`$ constant by adjusting the parameter $`ϵ`$ accordingly. For the subpolynomial scheme, the divisor $`m`$ is chosen to be $`m\sqrt{n}`$. Fig. 5.3 clearly confirms the expectation that the bitreversal scheme, considered as a generalisation of the subpolynomial ordering scheme, gives smaller values of $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$. For degrees of the Chebyshev polynomial $`n>40`$ rounding error effects are substantially suppressed in the bitreversal scheme compared to the subpolynomial scheme.

In Fig. 5.4 we show the values of $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ for the bitreversal, the naive, the pairing and the Montvay schemes for the case of $`\delta =0.001`$ as a function of $`n`$. We have also performed numerical tests at $`\delta =0.1`$ and $`\delta =0.01`$ where we found the same qualitative behaviour of $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ as a function of $`n`$ for the different ordering schemes. The first striking observation in Fig. 5.4 is that with the naive ordering one obtains already for moderate degrees $`n30`$ of the polynomial large values of $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$, indicating that rounding error effects are becoming severely problematic. Clearly, in the naive ordering scheme rounding errors can lead to very large ratios of particular partial products. Using the naive scheme, especially on machines with only 32-bit precision, a safe evaluation of the Chebyshev polynomial in the product representation can certainly not be guaranteed.

The behaviour of the values of $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ obtained by using the naive scheme demonstrates the necessity of finding better root ordering schemes which are able to reduce rounding error effects in evaluating the factorized Chebyshev polynomial. That such ordering schemes do exist is also demonstrated in Fig. 5.4. For $`n<100`$, the values for $`R_{\mathrm{max}}`$ and $`M_{\mathrm{max}}`$ obtained from the pairing, bitreversal and Montvay schemes are close to each other and many orders of magnitude below the ones of the naive scheme. However, for $`n>120`$, the values of $`R_{\mathrm{max}}`$ from these ordering schemes also start to deviate from each other. Taking $`R_{\mathrm{max}}`$ as a measure of the effects of rounding errors, it seems that the bitreversal scheme can reduce the rounding error effects most efficiently among the ordering schemes investigated here.

### 5.4 Matrix valued numerical tests

We report on a direct test of the ordering schemes described above for matrix-valued polynomials. As these tests were done using an existing program running on the powerful massively parallel Alenia Quadrics (APE) machines, we briefly detour to 4-dimensional lattice QCD \[Mut87\], with details given in App. F.

The task is to construct the matrix valued Chebyshev polynomial of the operator $`\widehat{Q}^2`$

$$P_n(\widehat{Q}^2)=\underset{k=1}{\overset{n}{}}n_k(\widehat{Q}^2z_k),$$
(5.12)
where the roots $`z_k`$ and the normalisation constant $`n_k`$ are given by Eq. 4.19 and Eq. 5.5, respectively. We evaluated the polynomial, Eq. 5.12, using the Clenshaw recurrence as well as using the different ordering of roots described above.

The numerical tests are performed on thermalized configurations on $`8^3\times 16`$ lattices using the APE computers, which have only 32-bit precision. Simulation parameters were chosen to be $`\beta =6.8`$, $`\kappa =0.1343`$ and $`c_{\mathrm{sw}}=1.42511`$. They correspond to realistic parameter values as actually used in simulations to determine values of $`c_{\mathrm{sw}}`$ non-perturbatively \[Jan97\]. We adopt the same (Schödinger functional) boundary conditions as described in \[Lue97\] for the evaluation of $`c_{\mathrm{sw}}`$. For the above choice of parameters and setting $`c_M=0.735`$, the lowest eigenvalue of $`\widehat{Q}^2`$ is $`\lambda _{\mathrm{min}}=0.00114(4)`$ and the largest is $`\lambda _{\mathrm{max}}=0.8721(3)`$. Investigations are performed at values of $`(n,ϵ)`$ $`(16,0.003)`$, $`(32,0.003)`$, $`(64,0.0022)`$ and $`(100,0.0022)`$. At each of these values of $`(n,ϵ)`$ we have $`O(50)`$ configurations.

We apply the matrix $`\widehat{Q}^2P_{n,ϵ}(\widehat{Q}^2)`$, which should be close to the unit matrix for our choices of $`n`$ and $`ϵ`$, to a random Gaussian vector $`G_{alpha,s}(x)`$, which is a complex vector, located at a lattice point $`(x)`$ and carrying colour $`\alpha =1,2,3`$ and spinor $`s=1,\mathrm{},4`$ indices. We then compute the vectors

$$\mathrm{\Phi }_{\mathrm{Clenshaw}}=\widehat{Q}^2P_{n,ϵ}(\widehat{Q}^2)G,$$
(5.13)
where $`P_{n,ϵ}(\widehat{Q}^2)`$ is constructed via Clenshaw’s recurrence relation and

$$\mathrm{\Phi }_{\mathrm{order}}=\widehat{Q}^2P_{n,ϵ}(\widehat{Q}^2)G,$$
(5.14)
where now $`P_{n,ϵ}(\widehat{Q}^2)`$ is evaluated using different root ordering schemes, and order stands for naive, pairing, bitreversal and Montvay. On a given configuration and for given $`G`$ we finally determine

$$\mathrm{\Delta }=\mathrm{\Phi }_{\mathrm{Clenshaw}}\mathrm{\Phi }_{\mathrm{order}}^2.$$
(5.15)
Since the Clenshaw recurrence is the numerically most stable method to evaluate the Chebyshev polynomial, the values of $`\mathrm{\Delta }`$ provide a measure for the effects of rounding errors. The result for $`\mathrm{\Delta }`$ as a function of $`n`$ is shown in Fig. 5.5. Using the naive ordering scheme, we could not run the cases of $`n=64`$ and $`n=100`$ as we hit numerical overflows. For the pairing scheme we find large values for $`\mathrm{\Delta }`$ at $`n=64`$ and $`n=100`$. The bitreversal scheme gives small, but non-negligible values of $`\mathrm{\Delta }`$ for all values of $`n`$ used. Finally, Montvay’s scheme gives $`\mathrm{\Delta }10^6`$ for all values of $`n`$. Surprisingly, using roots ordered by the Montvay scheme, the construction of the Chebyshev polynomial can be done with a stability that is comparable to the one using the Clenshaw recurrence.

## Chapter 6 The test: Physical results

### 6.1 Consistency tests

In general, the results of valid simulation algorithms must be independent of technical simulation parameters, random number generator types, irrelevant boundary conditions and measurement details like source or noise types used for the correlators or fermion densities.

To check whether we can rely on the results of the program package used to generate the CPU cost data of Ch. 8, we present various tests comparing observables from different measurement strategies to analytical results, those from HMC code or other groups.

Critical kappa

The critical kappa was estimated using the peak in the number of Conjugent Gradient solver iterations needed to invert $`Q`$. We perform dynamical simulations using the LBA with acceptance step method I on $`16\times 16`$ lattices with values of $`\beta =2`$, $`\beta =6`$ and $`\beta =10`$. Results shown in Fig. 6.1 to 6.3 lead to estimated $`\kappa _c`$ of $`0.287`$ for $`\beta =2.0`$, $`0.268`$ for $`\beta =6.0`$ and $`0.264`$ for $`\beta =10.0`$ with quite large errors of about $`0.02`$. The estimated values are indicated via a vertical line in the figures. We thus obtain the expected behaviour of $`\kappa _c\frac{1}{4}`$ for high $`\beta `$ and an increase in $`\kappa _c`$ as one decreases the value of $`\beta `$ described in Sec. 2.2.

An independent and usually better way to estimate the critical kappa is the extrapolation to the kappa value at which the pion mass is vanishing. Using dynamical data from $`16\times 40`$ and $`16\times 32`$ lattices respectively, we estimate from Fig. 6.11 a value of $`0.266`$ for $`\beta =6.0`$ and from Fig. 6.6 a value of $`0.262`$ at $`\beta =10.0`$ with errors of about the same order of magnitude as above.

Both independent measurement schemes thus yield the same numbers within errors. Moreover, results also agree with expectations.

Noisy scheme tests

We simulate full dynamical fermions on $`4\times 4`$ lattices and $`\beta =1.0`$ using the LBA with different approximation qualities and approximation regions. The acceptance step method I acceptances or respectively the reweighting factors are calculated both in a noisy and exact way.

Plaquette. We compare in Tab. 6.1 to 6.4 plaquette averages to the analytic values obtained using a hopping parameter expansion and the exact plaquette for pure $`U(1)`$ gauge theory as described in App. B.4 and B.3.

All calculations compare simulations with parameter choices $`n=0,2,10`$ and $`ϵ=0.5,0.01`$ to the leading order hopping parameter expansion result which is believed to be reasonable up to $`\kappa 0.15`$. We generally accumulated statistics of about $`10^5`$ sweeps. Only in the case of the acceptance step method I used with $`n=0`$, where the Metropolis acceptance step is correcting for the full fermion determinant estimated stochastically, this was clearly not enough. We accumulated in this case statistics of up to $`10^7`$ sweeps with integrated autocorrelations of $`O(1000)`$.

The agreement of the simulation results is evident in all cases. Analytical results are reproduced up to $`\kappa =0.15`$. For larger values of $`\kappa `$ systematic deviations from the analytic results are encountered as expected.

Correction factor. Checking the influence of the noisy scheme on the determination of the correction factor we compare in Tab. 6.5 factors calculated from eigenvalues to those estimated. We show results for the correction factor from the eigenvalues of $`Q^2`$ and from 1000 noisy estimation steps. All calculations were done on a $`4\times 4`$ lattice with $`\beta =1.0,\kappa =0.1,ϵ=0.01`$ and different $`n`$.

The results show the compatibility of the exact and the noisy estimation scheme. The $`n=2`$ results nicely demonstrate the importance of importance sampling, i.e. that simulating an almost flat distribution is impossible for dynamical fermions as corrections become too large. They also show that the correction factor can in fact deviate decisively from 1.

Different noises. We compare in Tab. 6.6 $`Z_2`$ and $`Z_4`$ complex Ising variables used in the noisy estimators for fermion densities as described in App. A.5. We calculate the fermion condensate $`\overline{\psi }\psi `$ and pseudo-scalar density $`\overline{\psi }\gamma ^5\psi `$ in the quenched case on $`16\times 16`$ lattices for $`\beta =2.0`$ and $`\kappa =0.26`$.

Results nicely agree with each other. We remark that the imaginary part of the pseudoscalar density was found to be zero to about $`10^{10}`$.

$`U=1`$ tests

For $`U=1`$ the condensate can be calculated analytically as described in App. B.1. We show results for this basic test from runs with links fixed to $`U=1`$ on $`16\times 16`$ lattices for a range of $`\kappa `$ values in Tab. 6.7. They obviously agree with the analytic predictions.

Condensate tests

In order to check the routines, a comparison to completely independent results is desirable.

For the fermion condensate in dynamical simulations the group in Graz obtained data for the 2 flavour Schwinger model with Wilson fermions using a HMC update algorithm on $`16\times 16`$ lattices at $`\beta =2.0`$ \[Lan96\].

In Fig. 6.4 we show the comparison of LBA with acceptance step method I data to these independent results. The agreement is evident.

Random generator tests

As mentioned in Ch. 3, we throughout this study used the high-quality random number generator (RG) by M. Lüscher \[Lue93\] to ensure that results are independent from random numbers. To further verify this, we compare quenched results on $`16\times 16`$ lattices at $`\beta =10`$ obtained using a simple vectorized XOR random number generator (XOR RG) \[Pre92\] and the vectorized Lüscher generator (Lüscher RG). In both cases periodic boundary conditions (BCs) are used.

In Tab 6.8 and Tab 6.9 we compare the condensate and pseudo-scalar density results to each other. These results clearly show that agreement is almost perfect up to the critical kappa $`\kappa _c\frac{1}{4}`$. Both the condensate and pseudoscalar density values above $`\kappa _c`$ are almost meaningless, as errorbars are of the order of the results themselves. Still, even those results are consistent with each other.

We remark that the pseudoscalar density results for large $`\kappa `$ do not agree with zero as was expected. This effect is due to topological effects discussed in Ch. 7. Still, pseudoscalar results are smaller than the fermion condensate by orders of magnitude.

Boundary conditions tests

In a gauge $`U(1)`$ theory we expect no influence of periodic or anti-periodic boundary conditions as $`1`$ is in the centre of the group. To check this, we repeat in Tab. 6.11 and 6.10 the periodic boundary condition data discussed in the paragraph above with anti-periodic boundary conditions.

Agreement of results with those of the periodic boundary conditions runs is evident in all cases.

Pion mass

We check the implementation of the dynamical fermion update comparing mass results with dynamical fermions against an independent calculation. Using a Hybrid Monte Carlo code, Irving \[Irv96\] gives a pseudoscalar vector (pion) mass of $`m_\pi =0.369(3)`$ for a $`32\times 32`$ lattice with $`\beta =2.29`$, $`\kappa =0.26`$. We obtain $`m_\pi =0.377(4)`$ from about 2000 measurements using acceptance step method I, thus agreeing with high precision with the independent value.

Symmetry tests

We test that the alternative meson operators in the various channels give compatible results.

Quenched. Meson masses in the quenched approximation were obtained on $`16\times 32`$ lattices for $`\beta =6.0`$ and $`\kappa =0.20\mathrm{}0.275`$ performing high-statistics runs with about 2000 independent measurements each. Fig. 6.5 shows the clear signal for the pion mass decreasing as $`\kappa `$ is increased, its alternative operator gives masses which agree within errors and is therefore not depicted in the plot. For the $`\eta `$ and its variant, consistency is also verified with larger errors. The higher states unfortunately give noisy results.

Dynamical. Masses were determined on $`16\times 32`$ lattices for $`\beta =10.0`$, $`\kappa =0.20\mathrm{}0.25`$ and are shown in Fig. 6.6. We observe the same symmetry behaviour as in the quenched case, with the alternative operators for the eta and pion masses again agreeing within errorbars. Due to the increased complexity of the computations, the data at large $`\kappa `$ and for higher states is very noisy compared to the quenched case.

### 6.2 Choice of parameters

We want to simulate at relevant physical parameter choices, in optimal cases understanding the scaling towards the continuum and chiral limits. To justify the choice of parameters for the CPU cost studies, we show that finite size effects do not affect the mass results at the chosen lattice sizes, $`\beta `$ and $`\kappa `$ values. In a second step, we demonstrate that the simulation parameters chosen on our larger lattices correspond to the same physical situation with a lattice spacing which is smaller by a factor of about $`1.6`$.

Finite size effects

In order to check the influence of the finite lattice extent, we regard the meson mass spectrum on small $`8\times 20`$ lattices. Naively, one would expect that masses can be obtained in the region $`0.5<m<2`$.

We simulate on $`8\times 20`$ lattices with a very conservative $`\beta =3.0`$ because of the topological metastability problem for high $`\beta `$ values as mentioned in Sec. 2.2 and a run length of generally $`>1000\tau `$. The topological problems are discussed in more detail in Ch. 7.

As can be seen in Fig. 6.7, finite size effects are small up to a $`\kappa `$ value of about 0.24. For larger $`\kappa `$, a deviation from the approximately linear behaviour of the pion mass is detectable.

To justify the linear fits for the pion mass, we included in Fig. 6.7 also a fit assuming the $`m^{\frac{2}{3}}`$ behaviour suggested by perturbation theory. It is evident that this describes the data far worse than the linear fit. We are obviously not yet in the regime where the leading order result is applicable.

This results in a minimal pion mass $`m_\pi (8\times 20)`$ = 0.629 possible on this lattices and a mass ratio $`\frac{m_\pi }{m_\eta }=0.807`$. These parameter choices are used for the CPU cost tests in Ch. 8.

The scaling procedure described in Sec. 2.2 relies on the possibility to determine for a certain $`\beta `$ the appropriate pion mass or $`\kappa `$ value corresponding to a fixed ratio $`\frac{m_\pi }{m_\eta }`$. This is illustrated in Fig. 6.8. The data shows a clear dependency on the pion mass even taking the errorbars into account. A linear fit and an inversion of the functional dependency is thus feasible.

Scaling to larger lattices

To achieve scaling towards the continuum limit, we aim to to reduce the lattice spacing going from $`\beta =3.0`$ and lattices of size $`8\times 20`$ to $`\beta =5.0`$ and appropriate larger lattices as described in Sec. 2.2.

We limit ourselves to these values of $`\beta `$, as for higher $`\beta `$ values topological metastabilities contaminate the data as discussed in Ch. 7. Mass spectra and fits for $`\beta =4,5,6`$ on $`16\times 40`$ lattices are shown in Fig. 6.9 to 6.11. This sequence of $`\beta `$ values is included to show the deterioration of the data going to higher $`\beta `$.

We have to determine the pion mass or respectively the $`\kappa `$ value corresponding to the same mass ratio $`\frac{m_\pi }{m_\eta }`$ as described above. To illustrate this, we include in Fig. 6.12 a plot of $`\frac{m_\pi }{m_\eta }`$ versus the pion mass. The lower quality compared to the $`8\times 20`$ results of Fig. 6.8 is an indication how fast CPU costs become prohibitively large for high-precision studies at larger lattices.

We fix the parameters for the optimal CPU cost search runs at $`\beta =5.0`$ and $`\kappa =0.245`$, yielding a pion mass of $`m_\pi (16\times 40)=0.384`$ and a ratio of $`\frac{m_\pi }{m_\eta }=0.826`$. As this ratio is compatible with the one obtained on the smaller $`8\times 20`$ lattices, we thus achieve a scaling by a factor of $`\frac{m_\pi (8\times 20)}{m_\pi (16\times 40)}=1.64`$. Finite size effects are again small on the lattices we used as the pion correlation length is about 3.

## Chapter 7 The problems II: Topology

Constructing decorrelated configurations in lattice simulations can be difficult if large energy barriers exist between regions of the configuration space to be sampled. Unfortunately, one such possibility are topological sectors for the $`U(1)`$ model \[Gat95\], part of the Schwinger model \[Joo90, Gat97a, Dil93\] we are considering in this study.

While local updates do not work well, we can use in two dimensions a global heatbath update to obtain reasonably high tunnelling rates. This trick obviously works also in the quenched case, while for full dynamical simulations global heatbath updates are not known.

A possible way of dealing with such systems which has been suggested \[Gut97\] is to stay in one fixed topological charge sector and define quantities in that way. In order to study this for relevant observables, we measure pion correlators in fixed (low) topological sectors. In the quenched case we compare local and global update schemes. For dynamical fermions we compare low and high $`\beta `$ results, using both the local bosonic algorithm (LBA) \[Lue94\] and Hybrid Monte Carlo (HMC) \[Dua87\].

### 7.1 Algorithm

Quenched case: local update. We use a local link update consisting of one exact heatbath \[Bes79\] and three over-relaxation steps per trajectory.

Quenched case: global update. Alternatively, it is possible to use a global heatbath for the plaquettes. As the new configuration is constructed without recursion to the old one, there are no problems with metastabilities.

The global update uses the fact that nearly all plaquettes are independent even on a finite lattice. The plaquettes only have to satisfy the constraint

$`{\displaystyle \underset{P}{}}U_P=1.`$ (7.1)
We are therefore able to update $`LT1`$ plaquettes with a heatbath algorithm. The last plaquette is then determined by the condition Eq. 7.1 and the configuration has to be accepted or rejected according to a Metropolis decision to ensure the correct distribution. Finally, this plaquette configuration has to be translated into a valid link configuration. To achieve this, we utilise the freedom to choose a gauge. We use a maximal tree prescription, setting $`LT2`$ links to 1. Then $`LT`$ links can be recursively determined from the plaquettes. The unconstrained two remaining links correspond to the free Polyakov loops $`P_T`$ and $`P_L`$ in 2 dimensions.

We want to state that the problem of slow topological charge fluctuations could also be solved by explicit topological updates \[Dil93\]. Unfortunately, this trick is not applicable in the presence of dynamical fermions. Improvements claimed to increase tunnelling for staggered fermions \[Dil93\] did not work for the Wilson fermions used in this study. Whether a recently suggested reweighting method to enhance topological updates \[For97c\] really reduces the problem is still under discussion.

Full dynamical simulations. In this case we use the Hermitean version of Lüscher’s local bosonic algorithm with noisy acceptance step method I as described in Ch. 4. To have an independent check, we compare to a Hybrid Monte Carlo algorithm implemented as described in Ch. 3.

### 7.2 Topological sectors

The integer-valued topological charge functional

$`{\displaystyle \frac{e}{4\pi }}{\displaystyle d^2xϵ_{\mu \nu }F_{\mu \nu }}`$ (7.2)
can be represented on the lattice by

$`Q={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{P}{}}\varphi _P`$ (7.3)
with plaquette angle $`\varphi _P=\mathrm{Im}\mathrm{ln}(U_P)(\pi ,\pi )`$ \[Lue82\]. We denote by tunnel events all updates resulting in a change of the topological charge and as tunnel probability the number of tunnel events divided by the total number of updates. The topological susceptibility is defined via

$`\chi _{\mathrm{top}}={\displaystyle \frac{1}{N_P}}[<Q^2><Q>^2].`$ (7.4)
We demonstrate the relevance of topological sectors on a $`16\times 32`$ lattice showing in Fig. 7.1 the tunnel rate plotted against $`\beta `$ for local and global updates. The exponential decrease of the tunnel probability for this local update algorithm gives rise to metastabilities in simulations at large $`\beta `$. We therefore consider it worthwhile to investigate observables within fixed topological sectors.

### 7.3 Simulations

Simulation parameters. We simulate on $`8\times 20`$ and $`16\times 40`$ lattices at a beta value of $`\beta =12`$ generating 10000 configurations. We only show data for the larger $`16\times 40`$ lattices. For the fermion part, we choose $`\kappa =0.24`$, where the pion correlation length (in the quenched case) is found to be around 3 and finite size effects can be expected to be small.

Local updates. We monitor the topological charge during the simulations. Due to metastability, no tunnelling of $`Q`$ is observed in the runs with local updates. At $`\beta =12`$ we therefore are able to perform a simulation in a given topological sector by using an initial configuration with this particular charge. This is done generating a classical homogeneous plaquette configuration of the desired charge and converting this to the links as described in Sec. 7.1. The two Polyakov loops $`P_T,P_L`$ are chosen according to a flat random distribution.

Global updates. In the limit of independent plaquettes the topological susceptibility can be calculated analytically as shown in App. B.2. The result in the large $`\beta `$ approximation for $`\beta =12`$ is

$`\chi _{\mathrm{top}}|_{\beta =12}{\displaystyle \frac{1}{4\pi ^2\beta }}|_{\beta =12}=0.0022.`$ (7.5)
To check the global update, we simulate on a $`16\times 32`$ lattice obtaining $`\chi _{\mathrm{top}}=0.0021(1)`$. This shows that the global update works well without metastabilities in the topological charge.

Observables. We generally measure the pion correlator using the prescription detailed in App. A.5.

Quenched results

The results of the quenched runs are shown in Fig. 7.2. Obviously, for $`Q=0,4`$ we are not able to find a plateau in the effective mass plots. We find a valley-like structure in the mass for the low $`Q`$ cases. For high $`Q`$ this turns into a hill-like structure with the peak situated at half of the temporal lattice extent. To show that in the intermediate region the valley and hill structures can approximately cancel and suggest a fake plateau, we include the $`Q=1`$ plot.

In the quenched case, we are able to compare to the results using a global update scheme. For the global update we find a tunnel probability of $`P=0.76`$ and therefore do not expect any influence of topology. The effective mass is shown in the lower right plot of Fig. 7.2. A plateau is clearly more reasonable than in the fixed $`Q`$ cases.

Dynamical fermion results

Effective masses from full dynamical simulations are shown in Fig. 7.3. We do not expect quantitatively the same results as in the quenched case, yet the behaviour is qualitatively similar.

To establish that dynamical results are not influenced by the chosen parameters of the local bosonic algorithm, we repeat the calculation with topological charge $`Q=4`$ using a standard HMC algorithm. This is also shown in Fig. 7.3. The results agree nicely within errorbars.

From these results, we conclude that we need to average over the topological sectors to obtain a plateau in the effective mass.

### 7.4 Projections to topological sectors

To gain further insight, we now use a slightly different approach. In principle, we could also restrict ourselves to definite topological sectors by selecting measurements with fixed topological charge from a simulation, i.e. effectively simulating the path integral given by

$`Z[q]={\displaystyle D[U]D[\overline{\psi }]D[\psi ]\delta _{Q,q}e^S}.`$ (7.6)
To this end we need simulations with a reasonably high fluctuation. Such simulations can be done e.g. with dynamical simulations at low $`\beta `$ or quenched simulations using global updates.

Full dynamical case. We work at low $`\beta =1`$ with a slightly smaller $`\kappa =0.22`$. Effective masses are depicted in Fig. 7.4. We can detect no discrepancy between masses calculated in different topological sectors. This result was also reported by a group working with staggered fermions, which concluded that topological sectors are of no importance to mass estimates in the Schwinger model \[Gut97\].

Quenched case. Here we exploit the opportunity to use the same parameters $`\beta =12,\kappa =0.24`$ as in Sec. 7.3. The results are shown in Fig. 7.5. At this $`\beta `$, we do not find agreement. Rather the effective masses are nearly the same as in the simulations without tunnelling presented in Fig. 7.3. We would like to point out that they do not agree within errorbars. On the other hand, we remark that the statistical sample was very much smaller for the projected data due to the fact that only a part of the generated configurations is projected into the appropriate sectors.

The striking difference between low and high $`\beta `$ results makes a sound understanding of this behaviour highly desirable. As can be seen from the quenched results here, it is evidently not just the averaging over the topological sectors (as found in Sec. 7.3) which is lacking in high $`\beta `$ simulations. There seems to be some more subtle dynamical effect involved.

### 7.5 External gauge configurations

In order to investigate how much averaging is necessary, we plot in Fig. 7.6 effective masses for external configurations with fixed topological charge $`Q=0`$ and $`Q=4`$. These are generated from homogeneous plaquette configurations in the way described in Sec. 7.1. For the fermions we use $`\kappa =0.24`$.

We stress that we use random Polyakov loops $`P_T,P_L`$ in both cases, so that one should not expect free fermion behaviour in the $`Q=0`$ case. For $`Q=4`$ a completely irregular behaviour is observed. It is thus not possible to measure a meaningful pion correlator from one (even very smooth) configuration alone. The translation invariance is manifestly broken even for that smooth configuration.

The result of averaging over 10 values of $`P_T`$ and $`P_L`$ is shown in Fig. 7.7. We clearly regain the qualitatively expected regular valley and hill structure observed in Fig. 7.3 in both the $`Q=0`$ and $`Q=4`$ case.

## Chapter 8 The aim: CPU cost optimisation

We remind the reader of the schemes included in our investigation. Besides
$``$ Hybrid Monte Carlo (HMC) to set the scale,
we compare
$``$ LBA with reweighting,
$``$ LBA with acceptance step method I ,
$``$ LBA with acceptance step method I with adapted precision and
$``$ LBA with acceptance step method II
as described in Ch. 4, investigating the CPU cost behaviour changing $`n,ϵ`$ and the number of reflections for the simulation runs.

We would like to point out again that reweighting and acceptance step algorithms result in two different ensembles, so that one has to use the effective CPU cost $`C_{\mathrm{eff}}=N_{\mathrm{total}\mathrm{Q}\mathrm{ops}}\frac{\sigma _{\mathrm{tot}}^2(A)}{<A>^2}`$ defined in Sec. 4.2.

### 8.1 $`8\times 20`$ lattices

Our search for optimal CPU cost behaviour is conducted on $`8\times 20`$ lattices with parameters $`\beta =3.0`$ and $`\kappa =0.24`$ measuring costs for the average plaquette and the pion correlator at distance $`\mathrm{\Delta }t=3`$.

We generally use anti-periodic boundary conditions and 5 hits for the link heatbath algorithm to provide high acceptance rates. To ensure numerical stability, we apply the bitreversal scheme for ordering of roots as described in Ch. 5.2 whenever a Chebyshev polynomial in the factorized form was necessary. The rescaling factor $`c_M`$ (viz. App. A.4, \[Lue94\]) is set to a conservative $`c_M=1.02`$, so that there is no slowing down from the large eigenvalues of $`Q^2`$. The precision of the inverters demanded is generally $`10^6`$, with $`10^2`$ for the reduced precision of the adapted precision scheme. We accumulated statistics of 10000 calculation sweeps, applying 1000 sweeps thermalization with integrated autocorrelations generally below 50.

To illustrate our search for optimal parameters, we depict in Fig. 8.1 and Fig. 8.2 the CPU cost in the $`nϵ`$ plane, where the number of reflections is optimised for each $`n,ϵ`$ pair. The figure clearly shows that we obtain flat optima.

In Tab. 8.1 we show the CPU cost for the Hybrid Monte Carlo algorithm run with different acceptance probability. This illustrates that the cost is not very dependent on the acceptance. To set the scale, we take the optimal values for an acceptance of $`0.84`$.

We in Tab. 8.2 compare the CPU cost of only the optimal parameter sets of the various LBA variants with the HMC scale.

We find that the LBA is doing better than the HMC by a factor of about 3 for the plaquette and about 2 for the pion correlator at distance $`\mathrm{\Delta }t=3`$. Thus the gain in the plaquette is typically better than that in the pion correlator. The reweighting method is a bit worse than the optimised acceptance step methods, losing for both observables by about a factor of 1.3. The gain using the adapted precision trick is quite considerable for the pion correlator (1.6), while it is even a slight detoriation for the plaquette as the overhead is too large. We remark that the configurations found for the optimised and unoptimised acceptance step method I scheme runs are not identical as one might expect in the case of no numerical errors. Overall, the number of reflections is important to the optimisation. A fixed number of 1 would not have reproduced the real optima. This leaves room for improvement in simulations where this has not yet been utilised.

### 8.2 $`16\times 40`$ lattices

For the larger lattices we choose the parameter values $`\beta =5.0`$ and $`\kappa =0.245`$ on $`16\times 40`$ lattices as detailed in Sec. 6.2. Technical simulation parameters are the same as described in Sec. 8.1.

We again illustrate our search for optimal parameters in Fig. 8.3 and Fig. 8.4, showing the CPU cost in the $`nϵ`$ plane, where the number of reflections is optimised for each $`n,ϵ`$ pair. The figure shows that although optima are fairly flat, we do not find as flat behaviour as for the smaller lattices. We admit, though, that somewhat larger $`n`$ should have been included in the study.

In Tab. 8.3 we give the CPU costs for the Hybrid Monte Carlo algorithm run with different acceptance probability. As we did not study this as exhaustively as in the smaller lattices case, we take the optimal values from all HMC runs to set the scale.

We in table 8.4 again compare CPU cost of the optimal parameter sets of the various LBA variants with the HMC scale.

We find that the LBA is again doing better than the HMC by a factor of about 3 for the plaquette and about 2.5 for the pion correlator at distance $`\mathrm{\Delta }t=3`$. Thus the gain in the plaquette is again better than that in the pion correlator, but not by as large a margin as for the smaller lattices. The reweighting method is again less efficient than the optimised acceptance step methods, this time even more so than for the smaller lattices, losing for both observables by about a factor of 1.8. The gain using the adapted precision trick is negligible for both observables. Overall, the number of reflections is not as important to the optimisation as for the small lattices.

## Chapter 9 Conclusions

The massive two-flavour Schwinger model has physical properties similar to QCD in four dimensions. It is much easier to simulate even with dynamical fermions and allows to determine observables with high precision as the cost will be lower by orders of magnitude as compared to the case of QCD in four dimensions. This makes the Schwinger model a reasonable testing ground for dynamical fermion algorithms. Still, in the case of the local bosonic algorithm, this study is the first application using Wilson fermions.

On the other hand, due to problems with topological sectors (as discussed below) and the scaling properties towards the continuum limit, it is definitely not that easy a toy model as one might expect.

### 9.1 Instabilities

In a class of fermion simulation algorithms relying on the local bosonic algorithm (LBA) a matrix valued Chebyshev polynomial is involved. Recursion relations allow the evaluation of these polynomials in a numerically stable way. Yet in a number of cases, like our implementation of the Acceptance step method I in the Hermitean LBA or the Polynomial Hybrid Monte Carlo (PHMC) algorithm \[Fre97a\], the polynomial is needed in the factorized form. Then rounding errors can easily lead to significant precision losses and even numerical instabilities, especially if simulations are done on machines having only 32-bit floating point arithmetic precision.

We investigate the effects of using various ordering schemes of monomial factors, or equivalently the complex roots, on the numerical construction of the Chebyshev polynomials now commonly used for the LBA. We find that different ordering schemes for the roots can lead to rounding error effects ranging from numerical overflow to retaining a precision comparable to the one numerically stable recurrence relations can provide.

In the case of a Chebyshev polynomial of a single real variable $`s`$ approximating the function $`1/s`$, we find that the bitreversal scheme and a scheme suggested by Montvay can keep rounding error effects to a low level for degrees of the Chebyshev polynomial up to $`n220`$.

Applying these reordering schemes to the evaluation of a matrix valued polynomial, we study numerical simulations of 4-dimensional lattice QCD. We find that Montvay’s ordering scheme of the roots seems to be particularly suited for this problem. The rounding errors could be kept on a level which is comparable to the one that is reached when using the stable Clenshaw recurrence relation.

We conclude that the precision with which the numerical evaluation of the Chebyshev polynomial can be performed depends strongly on the chosen ordering of the roots. We expect also severe consequences for the dynamics of the simulation algorithms where Chebyshev polynomials in the product representation are used, depending on the root ordering scheme employed.

As the most important outcome of this investigation, we consider that there exist orderings of the roots which allow a numerically very stable evaluation of a Chebyshev polynomial, even up to degrees $`n`$ of the polynomial of about 200. Since these values of $`n`$ correspond to degrees of the Chebyshev polynomials commonly used in simulations, we consider our findings as promising for future applications of the local bosonic algorithm.

### 9.2 Topology

During these studies we encountered simulations at high $`\beta `$ ($`8`$) exhibiting metastabilities in the topological charge. For dynamical simulations, no cure to this problem is known.

Motivated by this, we studied simulations at extreme $`\beta `$ values ($`\stackrel{>}{}12`$) which are effectively at fixed charge as no tunnelling occurs within these runs. Results for effective pion masses show that a definition of a pion mass from a plateau is not possible for either quenched or dynamical simulations even for vanishing topological charge.

A comparison with quenched global updates exhibiting no metastabilities demonstrates that a plateau can be found in the correct path integral sample. Furthermore, cross-checks against the Hybrid Monte Carlo algorithm (HMC) indicate that this problem is not an artifact stemming from the fermion algorithm.

To gain insight into the mechanism, we studied effective pion masses calculated from projections to fixed topological sectors. Results from fluctuating ensembles show no dependence on the topological charge. On the other hand, those projected from fluctuating quenched ensembles at high $`\beta `$ show approximately the same problematic behaviour as simulations completely without tunnelling. This suggests a rather subtle dynamical effect we do not understand yet.

We find that external homogeneous plaquette configurations with fixed Polyakov loop values $`P_T`$ and $`P_L`$ exhibit completely irregular behaviour. After averaging over $`P_T`$ and $`P_L`$ we regain the effective mass results characteristic for the topological charge sectors these configurations lie in.

We conclude that there is a need to obtain a better understanding of the interplay between the dynamical mass generation of mesonic states and topological sectors. We would like to point out that higher statistics runs could reveal similar phenomena in other models with non-trivial topological structure.

### 9.3 CPU cost optimisation

We optimised the CPU cost of the local bosonic algorithm varying 3 technical parameters of the algorithm, namely the number of over-relaxation steps in each trajectory, the order of the approximating Chebyshev polynomial $`n`$ and the lower cut of the approximation region $`ϵ`$. We simulated at two different lattice sizes $`8\times 20`$ and $`16\times 40`$ keeping the physical mass ratio $`\frac{m_\pi }{m_\eta }=0.81`$ approximately constant. This resulted in pion masses of about $`m_\pi =0.629`$ and $`m_\pi =0.384`$, or equivalently in a scaling of the lattice spacing of about $`1.6`$.

The tuning of the LBA is demonstrated to be surprisingly easy. It became more difficult for the larger lattices, though. Technically, the number of reflections per heatbath update, usually fixed in other studies of the algorithm, is found to be an important optimisation tool.

The CPU cost can be lower than for HMC, but not by a large factor with present techniques. We find a gain for the plaquette of approximately 3 for both lattice sizes and for the pion propagator of $`2`$ on the smaller and $`2.5`$ on the larger lattices. The gain thus differs and estimates from plaquette-like observables are too optimistic. Still, our main point is that the gain is detectable and consistent for both lattice sizes.

We also demonstrate that using a noisy Metropolis acceptance step scheme to make the LBA exact is also possible for the Hermitean case. The use of the Gegenbauer solver, which avoids instabilities in the evaluation of the polynomial, in the method labelled as acceptance step method II, is shown to be competitive to CG in the first real simulation. Costs are virtually identical to those of the acceptance step method I for both lattice sizes.

In general, all investigated acceptance variants performed similarly. The reweighting method, though, had decisively higher costs. Especially for the larger lattices it performed worse than the acceptance step methods for both observables studied by a factor of 1.8 (1.3 for the smaller lattices). Still, as the difference is not really large, the possible advantages of this method regarding exceptional configurations \[Fre97\] make further study worthwhile.

The optimisation of the acceptance step method I, interrupting the solver iterations and checking at intermediate steps whether the solution quality was already good enough for the requested Metropolis decision, did result in a significant gain of a factor of 1.6 compared to the unoptimised version regarding the cost of the pion correlator on the smaller lattices. On the larger lattices, though, we did not find any improvement. As this is nevertheless an optimisation offering possible gain, future studies should have this trick in mind as the gain could be different for larger inverter trajectories.

We conclude that these results give further evidence that the LBA is a competitive algorithm for the simulation of dynamical fermions.

## Appendix A Conventions

### A.1 General conventions

Dimensions. We use the standard conventions

$`\mathrm{}=c=1.`$ (A.1)
In $`d=2`$ dimensions then follow

$`[p]=[e]`$ $`=`$ $`[m]`$
$`[l]=[t]`$ $`=`$ $`[m]^1.`$ (A.2)
The inverse lattice coupling $`\beta =\frac{1}{e^2a^2}`$ and the Wilson parameter $`\kappa =\frac{1}{2(ma+d)}`$ are dimensionless. Gauge links $`U_{x,\mu }`$ are dimensionless, while Grassmann fermion fields $`\psi `$ have the dimension $`[m]^{\frac{1}{2}}`$. In the electromagnetic tensor

$`F^{\mu \nu }=\left(\begin{array}{cc}0& E\\ E& 0\end{array}\right)`$ (A.3)
in 1+1 dimensions only the electric field $`E`$ appears, as because of the missing transverse directions no magnetic field is existing. The electric field has the same dimension as the charge.
Indices. We denote
flavour indices with Latin lower case letters usually starting with $`a`$,
lattice sites with Latin lower case letters usually starting with $`x`$,
spinor indices with Latin lower case letters usually starting with $`s`$.
Directions we abbreviate with Greek lower case letters $`\mu ,\nu ,1,2,\mathrm{}`$
and use barred letters for the orthogonal direction (only one in 2D) $`\overline{\mu },\overline{\nu },\overline{1},\overline{2},\mathrm{}.`$
Generally, 1 signifies the spatial direction, 2 the temporal.
Unit vectors are denoted by a hat $`\widehat{\mu },\widehat{\overline{\mu }},\widehat{1},\widehat{\overline{1}},\mathrm{}`$.
Trace conventions. The different traces are

$`\mathrm{TR}`$ $`=`$ $`{\displaystyle \underset{x,s,a}{}};\text{i.e. sum over sites, Dirac and flavour indices,}`$
$`\mathrm{Tr}`$ $`=`$ $`{\displaystyle \underset{x,s}{}};\text{i.e. sum over sites and Dirac indices,}`$
$`\mathrm{tr}`$ $`=`$ $`{\displaystyle \underset{s}{}};\text{i.e. sum over Dirac indices,}`$
$`\mathrm{tr}_\mathrm{f}`$ $`=`$ $`{\displaystyle \underset{a}{}};\text{i.e. sum over flavour indices.}.`$ (A.4)
In an update trajectory, new variables and configurations are primed, the old ones unprimed.
Links. We use compact $`U(1)`$ link variables $`U_{x,\mu }`$. Plaquettes are defined starting at the lower left corner going anti-clockwise

$`U_{Px}=U_{x,\mu }^{}U_{x+\widehat{\mu },\overline{\mu }}^{}U_{x+\widehat{\overline{\mu }},\mu }U_{x,\overline{\mu }}.`$ (A.5)
Plaquette angles $`\varphi _x[\pi ,\pi ]`$ are defined via

$`e^{i\varphi _x}=U_{Px}.`$ (A.6)
Polyakov loops $`P_T,P_L`$ are used averaged over the orthogonal direction, i.e.

$`P_L={\displaystyle \frac{1}{T}}{\displaystyle \underset{x_2=1}{\overset{T}{}}}{\displaystyle \underset{x_1=1}{\overset{L}{}}}U_{x,1}\mathrm{and}P_T={\displaystyle \frac{1}{L}}{\displaystyle \underset{x_1=1}{\overset{L}{}}}{\displaystyle \underset{x_2=1}{\overset{T}{}}}U_{x,2}.`$ (A.7)
Fermions. Spinors are denoted with Greek letters $`\psi _{x,s}`$ throughout.

### A.2 Gamma matrices

We use Hermitean gamma matrices

$`\gamma ^0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\gamma ^1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ (A.8)
and a Hermitean $`\gamma ^5`$

$`\gamma ^5=i\gamma ^0\gamma ^1=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)`$ (A.9)
with the relations

$`\{\gamma ^\mu ,\gamma ^\nu \}=2\delta ^{\mu ,\nu },\{\gamma ^5,\gamma ^\mu \}=0,\gamma ^5\gamma ^\mu \gamma ^5=\gamma ^\mu .`$ (A.10)
Note that

$`\gamma ^5(1\pm \gamma ^\mu )\gamma ^5(1\pm \gamma ^\mu )=1\pm \gamma ^5\gamma ^\mu \gamma ^5\pm \gamma ^\mu \pm \gamma ^5\gamma ^\mu \gamma ^5\gamma ^\mu =0,`$ (A.11)
i.e. that the $`Q^2`$ matrix (Eq. A.26) does not connect sites with two straight links in-between. For general purposes we introduce the set

$`\mathrm{\Gamma }`$ $``$ $`\stackrel{~}{\mathrm{\Gamma }}=\{\text{1},(i\gamma ^5),\gamma ^0,(i\gamma ^1)\}`$ (A.12)
chosen so that

$`(\overline{\psi }\mathrm{\Gamma }\psi )^{}=(\overline{\psi }\mathrm{\Gamma }\psi )\mathrm{\Gamma }.`$ (A.13)
### A.3 Flavour matrices

The flavour matrices set is

$`T`$ $``$ $`\stackrel{~}{T}=\{\tau ^+,\tau ^{},\text{1},\tau ^0\},`$ (A.14)
where the tau matrices are derived from the Pauli matrices given by

$`\tau ^0`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\tau ^1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\tau ^2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right)`$ (A.15)
via

$`\tau ^+`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}(\tau ^1+i\tau ^2),\tau ^{}={\displaystyle \frac{1}{\sqrt{2}}}(\tau ^1i\tau ^2).`$ (A.16)
The Hermiticity relations

$`\text{1}^{}=\text{1},\tau ^0=\tau ^0,\tau ^+=\tau ^{}`$ (A.17)
hold. Further useful relations are

$`\mathrm{tr}_\mathrm{f}(\tau ^a)=0,\mathrm{tr}_\mathrm{f}(\tau ^a\tau ^a)=2,\mathrm{tr}_\mathrm{f}(\text{1})=2`$ (A.18)
and

$`{\displaystyle \underset{a}{}}\mathrm{tr}_\mathrm{f}(\tau ^a\tau ^a)=3\mathrm{t}\mathrm{r}_\mathrm{f}(\tau ^b\tau ^b)b.`$ (A.19)
### A.4 Fermion matrix

Starting from the standard Wilson fermion matrix

$`M=\delta _{x,y}\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\widehat{\mu },y}(1+\gamma ^\mu )U_{x\widehat{\mu },\mu }+\delta _{x+\widehat{\mu },y}(1\gamma ^\mu )U_{x,\mu }^{}\right)`$ (A.20)
with the Wilson kappa definition

$`\kappa `$ $`=`$ $`{\displaystyle \frac{1}{2(m+d)}},`$ (A.21)
we define via

$`M=\text{1}\kappa H`$ (A.22)
the hopping matrix $`H`$. The abbreviation $`G`$ for the Greens function $`M^1\delta _{\mathrm{flavour}}`$ is used throughout. We construct a Hermitean matrix $`Q`$ scaled such that eigenvalues lay in \[-1,1\] \[Lue94\] via

$`Q=c_o\gamma ^5\delta _{x,y}c_o\kappa {\displaystyle \underset{\mu }{}}\left(\delta _{x\widehat{\mu },y}\gamma ^5(1+\gamma ^\mu )U_{x\widehat{\mu },\mu }+\delta _{x+\widehat{\mu },y}\gamma ^5(1\gamma ^\mu )U_{x,\mu }^{}\right)`$ (A.23)
with the scaling factor

$`c_0`$ $`=`$ $`{\displaystyle \frac{m+d}{m+2d}}{\displaystyle \frac{1}{c_M}}={\displaystyle \frac{1}{1+2d\kappa }}{\displaystyle \frac{1}{c_M}};c_M1.`$ (A.24)
Inserting explicit gamma matrices, this results in the programmed formulae for $`Q`$

$`\left(Q\varphi \right)_{x,1}=`$ $`+`$ $`ic_0\varphi _{x,2}`$
$``$ $`2ic_0\kappa U_{x,0}^{}\varphi _{x+\widehat{0},2}`$
$``$ $`ic_0\kappa U_{x,1}^{}(\varphi _{x+\widehat{1},2}\varphi _{x+\widehat{1},1})`$
$``$ $`ic_0\kappa U_{x\widehat{1},1}(\varphi _{x\widehat{1},1}+\varphi _{x\widehat{1},2})`$
$`\left(Q\varphi \right)_{x,2}=`$ $``$ $`ic_0\varphi _{x,1}`$ (A.25)
$``$ $`ic_0\kappa U_{x,1}^{}(\varphi _{x+\widehat{1},2}\varphi _{x+\widehat{1},1})`$
$`+`$ $`2ic_0\kappa U_{x\widehat{0},0}\varphi _{x\widehat{0},1}`$
$`+`$ $`ic_0\kappa U_{x\widehat{1},1}(\varphi _{x\widehat{1},1}+\varphi _{x\widehat{1},2})`$
and $`Q^2`$

$`Q_{x,z}^2`$ $`=`$ $`c_0^2\delta _{x,z}`$ (A.26)
$``$ $`2c_0^2\kappa {\displaystyle \underset{\mu }{}}[U_{x\widehat{\mu },\mu }\delta _{x\widehat{\mu },z}+U_{x,\mu }^{}\delta _{x+\widehat{\mu },z}]`$
$`+`$ $`c_0^2\kappa ^2{\displaystyle \underset{\mu ,\mu ^{}}{}}`$
$`(`$ $`\delta _{x\widehat{\mu }\widehat{\mu }^{},z}(1\gamma ^\mu )(1+\gamma ^\mu ^{})U_{x\widehat{\mu },\mu }U_{x\widehat{\mu }\widehat{\mu }^{},\mu ^{}}`$
$`+`$ $`\delta _{x+\widehat{\mu }\widehat{\mu }^{},z}(1+\gamma ^\mu )(1+\gamma ^\mu ^{})U_{x,\mu }^{}U_{x+\widehat{\mu }\widehat{\mu }^{},\mu ^{}}`$
$`+`$ $`\delta _{x\widehat{\mu }+\widehat{\mu }^{},z}(1\gamma ^\mu )(1\gamma ^\mu ^{})U_{x\widehat{\mu },\mu }U_{x\widehat{\mu },\mu ^{}}^{}`$
$`+`$ $`\delta _{x+\widehat{\mu }+\widehat{\mu }^{},z}(1+\gamma ^\mu )(1\gamma ^\mu ^{})U_{x,\mu }^{}U_{x+\widehat{\mu },\mu ^{}}^{}).`$
The diagonal parts are

$`Q_{x,x}`$ $`=`$ $`c_0\gamma ^5=\left(\begin{array}{cc}0& ic_0\\ ic_0& 0\end{array}\right)`$
$`Q_{x,x}^2`$ $`=`$ $`c_0^2(1+4d\kappa ^2)=\left(\begin{array}{cc}c_0^2(1+4d\kappa ^2)& 0\\ 0& c_0^2(1+4d\kappa ^2)\end{array}\right).`$ (A.27)
### A.5 Meson operators, fermion densities and correlators

Operators. As meson operators we introduce fermion bilinears

$`O_x^{\mathrm{\Gamma }T}=\overline{\psi }_x\mathrm{\Gamma }T\psi _x,`$ (A.28)
sandwiching both a gamma matrix $`\mathrm{\Gamma }\stackrel{~}{\mathrm{\Gamma }}`$ and a flavour matrix $`T\stackrel{~}{T}`$ from the sets described above. We show the list of operators in 2 dimensions and the corresponding mesonic states in Tab. (A.1).

The 2D peculiarities

$`\gamma ^5\gamma ^0=i\gamma ^1\mathrm{and}i\gamma ^5\gamma ^1=\gamma ^0`$ (A.29)
reduce the possible gamma structures to four. The trivial action of $`\gamma ^0`$, if used on momentum zero states, leads to still further reduction to 2 relevant operators each for the flavour singlet and triplet. We remark that these symmetries are not used in the program. Correlations from all operators are explicitly calculated and the symmetry properties verified.
Fermion densities. To evaluate densities we need the general expression

$`\psi _{x_1}^{a_1}\overline{\psi }_{y_1}^{b_1}\mathrm{}\psi _{x_n}^{a_n}\overline{\psi }_{y_n}^{b_n}={\displaystyle \underset{x_1^{}\mathrm{}x_n^{}}{}}ϵ_{x_1\mathrm{}x_n}^{x_1^{}\mathrm{}x_n^{}}M_{x_1^{}y_1}^1\delta ^{a_1b_1}\mathrm{}M_{x_n^{}y_n}^1\delta ^{a_nb_n}`$ (A.30)
specialised to the 2 fermion isosinglet case

$`{\displaystyle \frac{1}{Vdn_f}}\psi _{x,r}^a\mathrm{\Gamma }_{rs}\overline{\psi }_{x,s}^a={\displaystyle \frac{1}{Vd}}{\displaystyle \underset{rsx}{}}\mathrm{\Gamma }_{rs}M_{x,s;x,r}^1.`$ (A.31)
Corresponding to the four possible gamma matrices $`\mathrm{\Gamma }\stackrel{~}{\mathrm{\Gamma }}`$ we obtain

* for $`\mathrm{\Gamma }=\text{1}`$ the fermion condensate,
* for $`\mathrm{\Gamma }=\gamma ^5`$ the axial condensate or pseudo-scalar density,
* for $`\mathrm{\Gamma }=\gamma ^1,\gamma ^0`$ the $`\gamma ^0`$ and $`\gamma ^1`$ densities.

Rewriting this in terms of $`Q`$ (Eq. A.23) and normalising by a constant $`\frac{1}{Vd}`$ to suppress trivial volume and dimension behaviour we arrive at the desired expression

$`{\displaystyle \frac{c_0}{Vd}}{\displaystyle \underset{rstx}{}}\gamma _{rt}^5\mathrm{\Gamma }_{ts}Q_{x,s;x,r}^1.`$ (A.32)
Noise. We use noise vectors with the properties

$`<\eta _x^{}\eta _y>=\delta _{x,y},<\eta _x>=0,`$ (A.33)
where the components are either continuous, having a Gaussian distribution, or are discrete $`Z_2`$ ($`\pm 1`$) or $`Z_4`$ ($`\pm 1,\pm I`$) Ising variables. The vectors do not carry spinor indices.
Noisy scheme. Using these noise vectors we apply the noisy estimator of the inverse of $`Q`$

$`Q_{xr,xs}^1=(\eta _x^{}{\displaystyle \underset{z}{}}Q_{xrzs}^1\eta _z)`$ (A.34)
where the inversion has to be done for each spinor component separately to obtain the full information.
Correlators. For the meson 2-point functions (correlators) we specialise Eq. A.30 to the 4 fermion case

$`\widehat{C}_{xy}^{\mathrm{\Gamma }T\mathrm{\Gamma }^{}T^{}}=\overline{\psi }_x\mathrm{\Gamma }T\psi _x\overline{\psi }_y\mathrm{\Gamma }^{}T^{}\psi _y`$ (A.35)
can then be expressed as

$`\widehat{C}_{xy}^{\mathrm{\Gamma }T\mathrm{\Gamma }^{}T^{}}`$ $`=`$ $`{\displaystyle \underset{abcd,\alpha \beta \gamma \delta }{}}\overline{\psi }_{x\alpha }^a\mathrm{\Gamma }_{\alpha \beta }T^{ab}\psi _{x\beta }^b\overline{\psi }_{y\gamma }^c\mathrm{\Gamma }_{\gamma \delta }^{}T^{cd}\psi _{y\delta }^d`$ (A.36)
$`=`$ $`{\displaystyle \underset{abcd,\alpha \beta \gamma \delta }{}}\mathrm{\Gamma }_{\alpha \beta }T^{ab}\mathrm{\Gamma }_{\gamma \delta }^{}T^{cd}\psi _{x\beta }^b\overline{\psi }_{y\gamma }^c\psi _{y\delta }^d\overline{\psi }_{x\alpha }^a`$
$`=`$ $`{\displaystyle \underset{abcd,\alpha \beta \gamma \delta }{}}\mathrm{\Gamma }_{\alpha \beta }T^{ab}\mathrm{\Gamma }_{\gamma \delta }^{}T^{cd}\left(G_{y\delta ,y\gamma }^{dc}G_{x\beta ,x\alpha }^{ba}G_{x\beta ,y\gamma }^{bc}G_{y\delta ,x\alpha }^{da}\right).`$
For operators having an expectation value we have to subtract this value from the correlator to get the connected Greens function. This is not explicitly written out here, nor is it used in the analysis routine - there the constant part is simply fitted.

Using the flavour structure of the Greens function, it is possible to write general formulae for all flavour matrices

$`\widehat{C}_{xy}^{\mathrm{\Gamma }T}`$ $`=`$ $`{\displaystyle \underset{\alpha \beta \gamma \delta }{}}\mathrm{\Gamma }_{\alpha \beta }\mathrm{\Gamma }_{\gamma \delta }\left((\mathrm{tr}_\mathrm{f}T)^2G_{y\delta ,y\gamma }G_{x\beta ,x\alpha }(\mathrm{tr}_\mathrm{f}T^{}T)G_{x\beta ,y\gamma }G_{y\delta ,x\alpha }\right),`$ (A.37)
where we introduce the simplified notation

$`C_{xy}^{\mathrm{\Gamma }T}`$ $`=`$ $`\overline{\psi }_x\mathrm{\Gamma }T\psi _x(\overline{\psi }_y\mathrm{\Gamma }T\psi _y)^{}=\widehat{C}_{xy}^{\mathrm{\Gamma }T;\gamma ^0T^{}\mathrm{\Gamma }^{}\gamma ^0}.`$ (A.38)
Using the connection to the scaled matrix $`Q`$

$`G=c_0Q^1\gamma ^5;Q^{}=Q`$ (A.39)
we get

$`C_{xy}^{\mathrm{\Gamma }T}={\displaystyle \underset{\alpha \beta ϵ\iota }{}}c_0^2(\gamma ^5\mathrm{\Gamma })_{\iota \beta }(\gamma ^5\mathrm{\Gamma })_{ϵ\delta }\left((\mathrm{tr}_\mathrm{f}T)^2Q_{y\delta ,yϵ}^1Q_{x\beta ,x\iota }^1(\mathrm{tr}_\mathrm{f}T^{}T)Q_{x\beta ,yϵ}^1Q_{y\delta ,x\iota }^1\right),`$ (A.40)
where the matrices

$`A_{ϵ\delta }={\displaystyle \underset{\gamma }{}}\gamma _{ϵ\gamma }^5\mathrm{\Gamma }_{\gamma \delta }`$ (A.41)
can be calculated and coded in the program.

We generally use correlators with randomly located point source and summation over the target time slice

$`C^{\mathrm{\Gamma }T}(t)`$ $`=`$ $`{\displaystyle \underset{x}{}}C_{(x+x_0,t+t_0)(x_0,t_0)}^{\mathrm{\Gamma }T},`$ (A.42)
where the index structure should be clear from the context.
Noisy scheme. Using Gaussian or Ising noise vectors with the properties Eq. A.33, we are able to sum over the time slice $`S(x_0t_0)`$ including the initial site $`x_0,t_0`$. Using Hermiticity of $`Q`$ we can derive in the triplet case

$`C_{\mathrm{triplet}}^\mathrm{\Gamma }(t)=6c_0^2{\displaystyle \underset{xabcd}{}}(\gamma ^5\mathrm{\Gamma })_{ab}(\gamma ^5\mathrm{\Gamma })_{cd}Q_{x+x_0t+t_0b;x_0t_0c}^1Q_{x_0t_0d;x+x_0t+t_0a}^1`$
$`={\displaystyle \frac{6c_0^2}{L}}{\displaystyle \underset{xabcd}{}}(\gamma ^5\mathrm{\Gamma })_{ab}(\gamma ^5\mathrm{\Gamma })_{cd}({\displaystyle \underset{z}{}}Q_{x+zt+t_0b;zt_0c}^1\eta _z)({\displaystyle \underset{x_0}{}}(Q_{x+x_0t+t_0a;x_0t_0d}^1\eta _{x_0})^{})`$ (A.43)
and, correspondingly, for the singlet

$`C_{\mathrm{singlet}}^\mathrm{\Gamma }(t)`$ $`=`$ $`{\displaystyle \frac{1}{3}}C_{\mathrm{triplet}}^\mathrm{\Gamma }(t)`$
$`+`$ $`{\displaystyle \underset{x\alpha \beta \gamma \delta }{}}c_0^2(\gamma ^5\mathrm{\Gamma })_{\alpha \beta }(\gamma ^5\mathrm{\Gamma })_{\gamma \delta }\left(4Q_{x_0t_0\delta ,x_0t_0\gamma }^1Q_{x_0+xt_0+t\beta ,x_0+xt_0+t\alpha }^1\right),`$
where one would estimate the disconnected part from the full point source.
Pion. For the pion a special trick is possible. The gamma structure allows to use noise vectors including spinor indices with

$`<\eta _{xs}^{}\eta _{yt}>=\delta _{x,y}\delta _{s,t},<\eta _{xs}>=0.`$ (A.45)
Using these only one inversion is necessary, as can be seen from

$`C^{\mathrm{pion}}(t)`$ $`=`$ $`6c_0^2{\displaystyle \underset{xabcd}{}}\delta _{a,b}\delta _{c,d}Q_{x+x_0t+t_0b;x_0t_0c}^1Q_{x_0t_0d;x+x_0t+t_0a}^1`$ (A.46)
$`=`$ $`{\displaystyle \frac{6c_0^2}{L}}{\displaystyle \underset{ax}{}}({\displaystyle \underset{ze}{}}Q_{x+zt+t_0a;zt_0e}^1\eta _{ze})({\displaystyle \underset{x_0c}{}}(Q_{x+x_0t+t_0a;x_0t_0c}^1\eta _{x_0c})^{}).`$
## Appendix B Analytic results

### B.1 Free fermions

We just give the result for the fermion condensate using free fermions where the momenta are restricted using anti-periodic boundary conditions to

$`k_\mu ={\displaystyle \frac{2\pi }{L_\mu }}(n_\mu +{\displaystyle \frac{1}{2}});n_\mu \{0,\mathrm{}L_\mu 1\}.`$ (B.1)
We obtain the condensate

$`{\displaystyle \frac{1}{n_fVd}}\mathrm{TR}(\overline{\psi }\psi )={\displaystyle \frac{1}{V}}{\displaystyle \underset{k}{}}{\displaystyle \frac{(12\kappa \underset{\mu }{}\mathrm{cos}k_\mu )}{_\mu (2\kappa \mathrm{sin}k_\mu )^2+(12\kappa _\mu \mathrm{cos}k_\mu )^2}}.`$ (B.2)
### B.2 Pure gauge topological susceptibility

In the limit of independent plaquettes the topological susceptibility can be calculated analytically. Introducing the generating functional with periodic BCs for $`N_P`$ independent plaquettes we obtain

$`Z(\beta ,\theta )=\left[{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi }{2\pi }}e^{i\theta Q+\beta \mathrm{cos}\varphi }\right]^{N_P},`$ (B.3)
where $`\varphi `$ denotes the plaquette angle defined in Eq. A.6. From this the topological susceptibility is given by

$`\chi _{\mathrm{top}}={\displaystyle \frac{1}{N_P}}_\theta ^2\mathrm{log}Z(\beta ,\theta )|_{\theta =0}`$ (B.4)
which can be evaluated to

$`\chi _{\mathrm{top}}`$ $`=`$ $`_\theta ^2\mathrm{log}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi }{2\pi }}e^{i\theta Q+\beta \mathrm{cos}\varphi }|_{\theta =0}`$ (B.5)
$`=`$ $`{\displaystyle \frac{_\pi ^\pi \frac{d\varphi }{2\pi }(\frac{\varphi }{2\pi })^2e^{\beta \mathrm{cos}\varphi }}{_\pi ^\pi \frac{d\varphi }{2\pi }e^{\beta \mathrm{cos}\varphi }}}=<({\displaystyle \frac{\varphi }{2\pi }})^2>.`$
Approximations. Two convenient approximations are possible. On the one hand the small $`\beta `$ limit yields

$`\chi _{\mathrm{top}}={\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi }{2\pi }}({\displaystyle \frac{\varphi }{2\pi }})^2={\displaystyle \frac{1}{12}}.`$ (B.6)
On the other hand the integral can be approximated by a Gaussian one in the large $`\beta `$ limit using

$`\beta \mathrm{cos}\varphi =\beta {\displaystyle \frac{\beta }{2}}\varphi ^2+\mathrm{}<\varphi ^2>{\displaystyle \frac{1}{\beta }}`$ (B.7)
yielding a topological susceptibility of

$`\chi _{\mathrm{top}}{\displaystyle \frac{1}{4\pi ^2\beta }}.`$ (B.8)
### B.3 Pure gauge plaquette

For the pure $`U(1)`$ theory in 2 dimensions it is possible to derive an analytical expression for the plaquette even on a finite lattice. We label link angles as $`\varphi _{lt\mu }`$ with $`l`$ and $`t`$ the site indices from 1 to $`L`$ and $`T`$ respectively and express the path integral in these variables

$`Z`$ $`=`$ $`e^{\beta \mathrm{\Omega }}{\displaystyle D[\varphi _{lt\mu }]e^{\beta _{lt}\mathrm{cos}(\varphi _{lt1}+\varphi _{l+1t2}\varphi _{lt+11}\varphi _{lt2})}}`$
$`=`$ $`e^{\beta \mathrm{\Omega }}{\displaystyle \underset{lt}{}}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi _{lt1}}{2\pi }}{\displaystyle \frac{d\varphi _{lt2}}{2\pi }}{\displaystyle \underset{n_{lt}}{}}I_{n_{lt}}(\beta )e^{in_{lt}(\varphi _{lt1}+\varphi _{l+1t2}\varphi _{lt+11}\varphi _{lt2})},`$
where we introduced the Fourier-transform of the exponential

$`e^{\beta \mathrm{cos}\varphi }={\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}c_ne^{in\varphi }`$ (B.10)
with coefficients

$`c_n={\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }e^{\beta \mathrm{cos}\varphi }e^{in\varphi }𝑑\varphi ={\displaystyle \frac{1}{\pi }}{\displaystyle _0^\pi }e^{\beta \mathrm{cos}\varphi }\mathrm{cos}(n\varphi )𝑑\varphi =I_n(\beta ).`$ (B.11)
$`I_n`$ signifies the modified Bessel function \[Pre92\]. As this expression for the path integral factorizes we can rewrite $`Z`$ to

$`e^{\beta \mathrm{\Omega }}{\displaystyle \underset{n_{lt}}{}}[{\displaystyle \underset{lt}{}}\left({\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi _{lt1}}{2\pi }}e^{in_{lt}\varphi _{lt1}in_{lt1}\varphi _{lt1}}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d\varphi _{lt2}}{2\pi }}e^{in_{lt}\varphi _{lt2}+in_{l1t}\varphi _{lt2}}\right){\displaystyle \underset{lt}{}}I_{n_{lt}}(\beta )]`$
$`=e^{\beta \mathrm{\Omega }}{\displaystyle \underset{n_{lt}}{}}\left[{\displaystyle \underset{lt}{}}\left(\delta _{n_{lt},n_{lt1}}\delta _{n_{lt},n_{l1t}}I_{n_{lt}}(\beta )\right)\right].`$ (B.12)
This reduces $`Z`$ to a simple sum

$`Z`$ $`=`$ $`e^{\beta \mathrm{\Omega }}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}I_n(\beta )^\mathrm{\Omega }=e^{\beta \mathrm{\Omega }}\left[I_0(\beta )^\mathrm{\Omega }+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}2I_n(\beta )^\mathrm{\Omega }\right],`$ (B.13)
where we used the symmetry $`I_n(x)=I_n(x)`$ valid for all $`n`$ and $`x`$. To finally get the plaquette we use the relation

$`<P>`$ $`=`$ $`1+{\displaystyle \frac{1}{Z\mathrm{\Omega }}}{\displaystyle \frac{Z}{\beta }}`$ (B.14)
and

$`{\displaystyle \frac{I_n(\beta )}{\beta }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[I_{n1}(\beta )+I_{n+1}(\beta )\right]`$ (B.15)
to get after some algebra to

$`<P>`$ $`=`$ $`{\displaystyle \frac{{\displaystyle \frac{I_1(\beta )}{I_0(\beta )}}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[\left({\displaystyle \frac{I_n(\beta )}{I_0(\beta )}}\right)^{\mathrm{\Omega }1}\left({\displaystyle \frac{I_{n1}(\beta )}{I_0(\beta )}}+{\displaystyle \frac{I_{n+1}(\beta )}{I_0(\beta )}}\right)\right]}{1+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{I_n(\beta )}{I_0(\beta )}}\right)^\mathrm{\Omega }}}.`$ (B.16)
This expression is then numerically evaluated using Fortran routines for Bessel functions \[Pre92\].

### B.4 Hopping parameter expansion

We start with the distribution of the links in a dynamic fermion simulation

$`P_{\mathrm{eff}}(det[\text{1}\kappa H])^2e^{\beta _x\mathrm{Re}U_{Px}}`$ (B.17)
with $`H`$ the hopping terms (Eq. A.22) and rewrite using

$`det[M^2]=e^{2\mathrm{T}\mathrm{r}\mathrm{log}M}.`$ (B.18)
We expand the logarithm

$`\mathrm{log}M=\mathrm{log}(1\kappa H)=\kappa H{\displaystyle \frac{\kappa ^2H^2}{2}}{\displaystyle \frac{\kappa ^3H^3}{3}}{\displaystyle \frac{\kappa ^4H^4}{4}}+\mathrm{}`$ (B.19)
and observe that the first link-dependent term appearing in the trace is the $`H^4`$ term as the odd powers of $`H`$ vanish in the trace $`\mathrm{Tr}H=\mathrm{Tr}H^3=\mathrm{}=0`$ and the second order contribution

$`\mathrm{Tr}H^2=\mathrm{const}.+\mathrm{const}.U_{x\widehat{\mu },\mu }U_{x\widehat{\mu },\mu }^{}\mathrm{}=\mathrm{const}.`$ (B.20)
is an irrelevant constant. So we end up with

$`det[M^2]=e^{\frac{1}{2}\kappa ^4\mathrm{Tr}H^4}`$ (B.21)
to leading order.

Writing out the $`H^4`$ term, we see that only closed loops of links can contribute to the trace, i.e. to leading order only plaquettes are contributing. This results in

$`\mathrm{Tr}H^4`$ $`={\displaystyle \underset{x,\mu }{}}(`$ $`\mathrm{tr}[(1+\gamma ^\mu )(1+\gamma ^{\overline{\mu }})(1\gamma ^\mu )(1\gamma ^{\overline{\mu }})]U_{x\widehat{\mu },\mu }U_{x\widehat{\mu }\widehat{\overline{\mu }},\overline{\mu }}U_{x\widehat{\mu }\widehat{\overline{\mu }},\mu }^{}U_{x\widehat{\overline{\mu }},\overline{\mu }}^{}`$ (B.22)
$`+`$ $`\mathrm{tr}[(1+\gamma ^\mu )(1\gamma ^{\overline{\mu }})(1\gamma ^\mu )(1+\gamma ^{\overline{\mu }})]U_{x\widehat{\mu },\mu }U_{x\widehat{\mu },\overline{\mu }}^{}U_{x\widehat{\mu }\widehat{\overline{\mu }},\mu }^{}U_{x,\overline{\mu }}`$
$`+`$ $`\mathrm{tr}[(1\gamma ^\mu )(1+\gamma ^{\overline{\mu }})(1+\gamma ^\mu )(1\gamma ^{\overline{\mu }})]U_{x,\mu }^{}U_{x+\widehat{\mu }\widehat{\overline{\mu }},\overline{\mu }}U_{x\widehat{\overline{\mu }},\mu }U_{x\widehat{\overline{\mu }},\overline{\mu }}^{}`$
$`+`$ $`\mathrm{tr}[(1\gamma ^\mu )(1\gamma ^{\overline{\mu }})(1+\gamma ^\mu )(1+\gamma ^{\overline{\mu }})]U_{x,\mu }^{}U_{x+\widehat{\mu },\overline{\mu }}^{}U_{x+\widehat{\overline{\mu }},\mu }U_{x,\overline{\mu }}).`$
To evaluate the gamma matrix expressions we observe that

$`\mathrm{tr}\text{1}`$ $`=`$ $`2,\mathrm{tr}\gamma ^i=0i,`$
$`\mathrm{tr}\gamma ^i\gamma ^j`$ $`=`$ $`2\delta _{i,j},\mathrm{tr}\gamma ^i\gamma ^j\gamma ^k=0i,j,k`$ (B.23)
and that therefore all gamma terms give the same contribution

$`\mathrm{tr}[\mathrm{gamma}\mathrm{terms}]=222+\mathrm{tr}[\gamma ^\mu \gamma ^{\overline{\mu }}\gamma ^\mu \gamma ^{\overline{\mu }}]=4.`$ (B.24)
Returning to the full expression we abbreviate the plaquette as above and obtain

$`\mathrm{Tr}H^4=4{\displaystyle \underset{x,\mu }{}}\left(U_{Px\widehat{\mu }\widehat{\overline{\mu }}}+U_{Px\widehat{\mu }}^{}+U_{Px\widehat{\overline{\mu }}}^{}+U_{Px}\right).`$ (B.25)
As we sum over the lattice sites, we can reorder the contributions and find that we have each term twice for each direction. We further combine the complex conjugate pairs so that the final expression for the determinant is

$`det[M^2]=e^{16\kappa ^4_x\mathrm{Re}U_{Px}}.`$ (B.26)
This first order result can be combined with the standard link action term so that the result is an effective pure gauge theory with a shifted beta value of

$`\beta ^{}=\beta +16\kappa ^4.`$ (B.27)
## Appendix C Force

### C.1 Force on bosonic fields

Starting with the effective action

$`S_L`$ $`=`$ $`\beta \mathrm{Re}{\displaystyle \underset{p}{}}(1U_p)+{\displaystyle \underset{k,x,y}{}}\varphi _x^k(Q\mu _k)_{x,y}^2\varphi _y^k+{\displaystyle \underset{k,x}{}}\nu _k^2\varphi _x^k\varphi _x^k`$ (C.1)
we want the part dependent on one specified field spinor $`\varphi _x^k`$. Rewriting yields

$`S_L(U,\varphi )`$ $`=`$ $`\mathrm{const}.+\varphi _x^k\left[{\displaystyle \underset{y|yx}{}}(Q^22Q\mu _k)_{x,y}\varphi _y^k\right]+\left[{\displaystyle \underset{y|yx}{}}\varphi _y^k(Q^22Q\mu _k)_{y,x}\right]\varphi _x^k`$ (C.2)
$`+\varphi _x^k(\mu _k^2+\nu _k^2+(Q^22Q\mu _k)_{x,x})\varphi _x^k`$
$`=`$ $`\mathrm{const}.+\varphi _x^kB_x^k+B_x^k\varphi _x^k+\varphi _x^kK_x^k\varphi _x^k.`$
Mapping this to a Gaussian distribution via

$`e^S`$ $``$ $`e^{\varphi _x^kB_x^kB_x^k\varphi _x^k\varphi _x^kK_x^k\varphi _x^k}e^{R^{}R}`$ (C.3)
we easily see heatbath updates

$`R=(K_x^k)^{\frac{1}{2}}\varphi _x^k+(K_x^k)^{\frac{1}{2}}B_x^k`$
$`\varphi _x^k=(K_x^k)^{\frac{1}{2}}R(K_x^k)^1B_x^k`$ (C.4)
and micro-canonical reflections steps

$`RR`$ $`=`$ $`(K_x^k)^{\frac{1}{2}}\varphi _x^k(K_x^k)^{\frac{1}{2}}B_x^k`$
$`\varphi _x^k`$ $`=`$ $`\varphi _x^k2(K_x^k)^1B_x^k.`$ (C.5)
The necessary coefficients $`B_x^k,(K_x^k)^1,(K_x^k)^{\frac{1}{2}}`$ can be calculated using the formulae in App. A.4 giving explicit expressions for $`K_x^k`$

$`K_x^k=c_0^2(1+4d\kappa ^2)2\mu _kc_0\gamma ^5+\nu _k^2+\mu _k^2=K^k,`$ (C.6)
which is not a function of the lattice site (though a Dirac matrix) and $`B_x^k`$

$`B_x^k`$ $`=`$ $`{\displaystyle \underset{yx}{}}(Q^22Q\mu _k)_{x,y}\varphi _y^k`$ (C.7)
$`=`$ $`{\displaystyle \underset{y}{}}[(Q^22Q\mu _k)_{x,y}\varphi _y^k][c_0^2(1+4d\kappa ^2)2\mu _kc_0\gamma ^5]\varphi _x^k.`$
### C.2 Force on links

The force on links is defined by

$`S_L=\mathrm{const}.+ReF_{x,\mu }U_{x,\mu }`$ (C.8)
with the intention to apply this to the heatbath formulae

$`Pe^{b\mathrm{cos}\varphi }.`$ (C.9)
Observing that the links $`U_{x,\mu }`$ are complex numbers, we define

$`b`$ $`=`$ $`|F_{x,\mu }|`$
$`\mathrm{cos}\varphi `$ $`=`$ $`\mathrm{cos}(\varphi _1+\varphi _2)\mathrm{with}e^{i\varphi _1}={\displaystyle \frac{F_{x,\mu }}{|F_{x,\mu }|}};e^{i\varphi _2}=U_{x,\mu }.`$ (C.10)
Starting with the action one easily sees

$`S_L`$ $`=`$ $`\mathrm{const}.\mathrm{Re}\beta U_p+{\displaystyle \underset{k,x,y}{}}\varphi _x^k(Q_{x,y}^22\mu _kQ_{x,y})\varphi _y^k.`$ (C.11)
Expanding and using $`2\mathrm{R}\mathrm{e}z=z+z^{}`$ and the gamma matrix relations

$`\left(\gamma ^5(1+\gamma ^\mu )\right)^{}=\gamma ^5(1\gamma ^\mu ),\left((1\pm \gamma ^\mu )(1\pm \gamma ^\mu ^{})\right)^{}=(1\pm \gamma ^\mu ^{})(1\pm \gamma ^\mu )`$ (C.12)
to combine the complex conjugate pairs and also complex conjugating one of the staple expressions (allowed as only the real part is taken) one obtains

$`F_{x,\mu }`$ $`=`$ $`\beta {\displaystyle \underset{\mu ^{}\mu }{}}(U_{x,\mu ^{}}^{}U_{x+\widehat{\mu }^{},\mu }^{}U_{x+\widehat{\mu },\mu ^{}}+U_{x+\widehat{\mu }\widehat{\mu }^{},\mu ^{}}^{}U_{x\widehat{\mu }^{},\mu }^{}U_{x\widehat{\mu }^{},\mu ^{}})`$ (C.13)
$`+`$ $`4c_0^2\kappa {\displaystyle \underset{k}{}}\left[\varphi _{x+\widehat{\mu }}^k\varphi _x^k\right]`$
$``$ $`2c_0^2\kappa ^2{\displaystyle \underset{k}{}}\left[{\displaystyle \underset{\mu ^{}}{}}\right(`$
$`\varphi _{x+\widehat{\mu }}^{}(1\gamma ^\mu )(1+\gamma ^\mu ^{})\varphi _{x\widehat{\mu }^{}}U_{x\widehat{\mu }^{},\mu ^{}}`$
$`+\varphi _{x+\widehat{\mu }}^{}(1\gamma ^\mu )(1\gamma ^\mu ^{})\varphi _{x+\widehat{\mu }^{}}U_{x,\mu ^{}}^{}`$
$`+\varphi _{x+\widehat{\mu }^{}+\widehat{\mu }}^{}(1\gamma ^\mu ^{})(1+\gamma ^\mu )\varphi _xU_{x+\widehat{\mu },\mu ^{}}`$
$`+\varphi _{x\widehat{\mu }^{}+\widehat{\mu }}^{}(1+\gamma ^\mu ^{})(1+\gamma ^\mu )\varphi _xU_{x\widehat{\mu }^{}+\widehat{\mu },\mu ^{}}^{})]`$
$``$ $`4c_0\kappa {\displaystyle \underset{k}{}}\left[\mu _k\varphi _{x+\widehat{\mu }}^k\gamma ^5(1+\gamma ^\mu )\varphi _x^k\right].`$
## Appendix D Root Ordering Fortran codes

In this appendix we want to give the Fortran codes which are used for the reordering of the roots Eq. 4.19. In the following $`n`$ will always denote the degree of the Chebyshev polynomial and $`m`$ is an integer divisor of $`n`$ as needed for the subpolynomial scheme.

### D.1 Simple pairing scheme

We show the code for the first half of roots. The second half is constructed analogously.

```
      if (abs(n/4 - one*n/4) .ge. 0.01) then
        med1 = n/4 + 1                             ! n/4 odd
        med2 = n/4 + 1
      else
        med1 = n/4                                 ! n/4 even
        med2 = n/4 + 1
      endif

      i=0
      do k=1,n/8                                   ! take complete
        i = i+1 ; j1(i) = k                        ! groups of 4
        i = i+1 ; j1(i) = n/2-k+1
        i = i+1 ; j1(i) = med1+k                   ! -> hashing table
        i = i+1 ; j1(i) = med2-k
      end do

      r = n/2 - n/2 /4*4
      if (r .ge. 2) then                           ! last pair
        k = n/8 + 1
        i = i+1 ; j1(i) = k
        i = i+1 ; j1(i) = n/2-k+1
        r = r-2
      endif
      if (r .ge. 1) then                           ! last lonely one
        i = i+1 ; j1(i) = med1
      endif

```

### D.2 Subpolynomial scheme

```
      kstep = n / m                                ! \# of subpolynomials
      i =0
      do k=1,m
        do l=1,kstep
          i = i+1 ; j1(i) = (l-1)*m + k            ! reset index
        end do
      end do

```

### D.3 Bitreversal scheme

```
      do b=0,20
        if (2**b .ge. n) then                      ! length of array
          bits=b; goto 333
        end if
      end do
 333  continue

      do k=1,2*n
        j(k) = 0                                   ! init
      end do

      do k=1,n
        do b=0,bits
          bit(b)=0                                 ! init
        end do

        nn = k-1                                   ! shift
        do b=bits-1,0,-1
          if (2**b .le. nn) then
            nn     = nn - 2**b
            bit(b) = 1                             ! extract bit
          end if
        end do

        i = 0
        do b=0,bits-1
          i = i + bit(b) * 2**(bits-1-b)           ! reverse
        end do
        i    = i+1                                 ! reshift
        j(i) = k                                   ! -> hashing table
      end do

      i = 0
      do k=1,2*n
        if (j(k) .ne. 0) then                      ! no dummy root
          i        = i+1
          j1(j(k)) = i                             ! reset index
        end if
      end do

```

### D.4 Montvay scheme

```

read(‘approxima.txt‘);

n   :=  12 ;
eps :=  .1000000044703483 ;
c0  :=  52987.39383322105 ;
Root:=[
6.299918614676031E-02   -.146958373796435 *I,
.2375643902633147   -.2602503551965 *I,
.483704827825326   -.313922116577824 *I,
.7450326909011724   -.295678104222787 *I,
.9616809154023317   -.209697801497745 *I,
1.08401800398973   -7.567825958203078E-02 *I,
1.08401800398973 +  7.567825958203071E-02 *I,
.9616809154023317 +  .2096978014977446 *I,
.7450326909011726 +  .2956781042227866 *I,
.483704827825326 +  .3139221165778238 *I,
.2375643902633154 +  .2602503551965006 *I,
6.299918614676020E-02 +  .146958373796435 *I ];

Optimord(n,Root,c0,1,1.,x,eps,1.,100,40,yes);

```

## Appendix E Gegenbauer Solver

### E.1 Gegenbauer Polynomials

We define polynomials $`C_n^\gamma (z)`$ of the complex variable $`z`$ with degree $`n`$ and a real parameter $`\gamma >0`$ via their generating function

$`(1+t^22tz)^\gamma ={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}t^nC_n^\gamma (z).`$ (E.1)
Alternatively, an integral representation

$`C_n^\gamma (z)`$ $`=`$ $`{\displaystyle \frac{2^{12\gamma }\mathrm{\Gamma }(2\gamma +n)}{n!\mathrm{\Gamma }(\gamma )^2}}{\displaystyle _0^\pi }𝑑\varphi (\mathrm{sin}\varphi )^{2\gamma 1}(z+\sqrt{z^21}\mathrm{cos}\varphi )^n`$ (E.2)
is possible. For the many known features of the thus defined Gegenbauer polynomials like trigonometric representation, coefficients of the highest monomial, parity, bounds and large $`n`$ approximations we defer the reader to the literature \[Bun97, Erd53\]. For completeness, we mention the first few polynomials

$`C_0^\gamma (z)`$ $`=`$ $`1`$
$`C_1^\gamma (z)`$ $`=`$ $`2\gamma z`$
$`C_2^\gamma (z)`$ $`=`$ $`2\gamma (\gamma +1)z^2\gamma `$
$`C_3^\gamma (z)`$ $`=`$ $`{\displaystyle \frac{4}{3}}\gamma (\gamma +1)(\gamma +2)z^32\gamma (\gamma +1)z.`$ (E.3)
The relative error of the partial sums is given by

$`R_n(z)`$ $`=`$ $`1(1+t^22tz)^\gamma {\displaystyle \underset{k=0}{\overset{n}{}}}t^kC_k^\gamma (z)`$ (E.4)
$`=`$ $`(1+t^22tz)^\gamma {\displaystyle \underset{k=n+1}{\overset{\mathrm{}}{}}}t^kC_k^\gamma (z).`$
The main item we want to stress is the existence of an recursion relation

$`(n+1)C_{n+1}^\gamma (z)+(n+2\gamma 1)C_{n1}^\gamma (z)=2(n+\gamma )zC_n^\gamma (z)`$ (E.5)
enabling an efficient calculation of these polynomials.

### E.2 Solver

The generating function suggests the use of this polynomial expansion to build up a solver method. Assume $`M`$ a Hermitean matrix with spectrum $`\text{spec}(M)[\lambda _1,\lambda _2]`$ and $`0<\lambda _1\lambda _2`$. We want to solve for $`x`$ in

$`M^\gamma x=b.`$ (E.6)
To map to the Gegenbauer polynomials defined in Eq. E.1, we transform to a normalised matrix $`A`$ with $`\text{spec}(A)[1,1]`$

$`M=c(1+t^22tA)A={\displaystyle \frac{1}{2tc}}M+{\displaystyle \frac{1+t^2}{2t}},`$ (E.7)
defining parameters $`t`$ and $`c`$

$`t={\displaystyle \frac{\sqrt{\frac{\lambda _2}{\lambda _1}}1}{\sqrt{\frac{\lambda _2}{\lambda _1}}+1}},c={\displaystyle \frac{\lambda _2\lambda _1}{4t}}.`$ (E.8)
We can then write the solution $`x`$ as

$`x`$ $`=`$ $`M^\gamma b=c^\gamma (1+t^22tA)^\gamma b`$ (E.9)
$`=`$ $`c^\gamma {\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}t^nC_n^\gamma (A)b={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}t^ns_n`$
with

$`s_n`$ $`=`$ $`c^\gamma C_n^\gamma (A)b,`$ (E.10)
i.e. as a sum over shift vectors $`s_n`$ with exponentially decreasing factors. Moreover, the shifts have an easy recursion relation

$`(n+1)s_{n+1}+(n+2\gamma 1)s_{n1}=2(n+\gamma )As_n`$ (E.11)
with start shifts $`s_1`$ and $`s_0`$ defined by

$`s_1`$ $`=`$ $`0`$
$`s_0`$ $`=`$ $`c^\gamma b.`$ (E.12)
To obtain a valid solver, we have to find a stopping criterion, i.e we have to calculate (at least bounds for) the relative error Eq. E.4 and the rest vector

$`r_n`$ $`=`$ $`bM^\gamma x_n=R_n(A)b.`$ (E.13)
### E.3 Real solver for $`\gamma =\frac{1}{2}`$

For $`\gamma =1`$ one obtains a standard inverter algorithm. Most interesting are non-standard cases, like e.g. $`\gamma =\frac{1}{4}`$, which is applicable to SUSY models \[Mon97b\]. The general cases for $`\gamma `$ and $`z`$ are discussed in more detail in \[Bun97\]. We for obvious reasons at this place only consider the real case with $`\gamma =\frac{1}{2}`$.

Legendre Polynomials The Gegenbauer polynomials in the $`\gamma =\frac{1}{2}`$ case are the Legendre polynomials

$`C_n^{1/2}(z)P_n(z).`$ (E.14)
In this case the integral representation is given by

$`P_n(\mathrm{cos}\varphi )`$ $`=`$ $`{\displaystyle _0^\pi }{\displaystyle \frac{d\varphi }{\pi }}(z+\sqrt{z^21}\mathrm{cos}\varphi )^n`$ (E.15)
with a normalisation $`P_n(1)=1`$. The recursion is identical to Eq. E.5.

The important point is that in this case the relative error can be estimated for real $`z[1,1]`$ to

$`|R_n(t)||t|^{n+1}.`$ (E.16)
For the proof we defer to \[Bun97\]. More involved is the estimation for $`|z|1`$, so we briefly sketch the idea. We regard the Legendre polynomials with complex argument $`z`$ parametrised by

$`z=\mathrm{cosh}(\tau +i\phi ).`$ (E.17)
Assume that $`z`$ an $`t`$ given by Eq. E.8 are inside the ellipse given by

$`Z=\mathrm{cosh}(\theta +i\varphi )`$ (E.18)
with $`\theta 0,\varphi [0,2\pi ]`$. Then we get a uniform bound for the relative error of the expansion

$`|R_n(z)|(|t|e^\theta )^{n+1}`$ (E.19)
and equivalently for the error of the solver assuming $`A`$ is normal with spectrum inside the ellipse.

A priori, $`\theta `$ is unknown. An estimate can be deducted from the shifts solving

$`{\displaystyle \frac{s_n}{s_0}}={\displaystyle \frac{P_n(A)b}{b}}P_n(\mathrm{cosh}\theta )`$ (E.20)
e.g. by Newton-Raphson iterations.

Solver implementation For the Hermitean local bosonic algorithm case $`M=Q^2P_{n_\varphi }(Q^2)`$ we are able to give an realistic upper bound for the spectrum

$`\lambda _{max}(M)=1+\delta (n,ϵ),`$ (E.21)
but we have to use a guess for the unknown $`\lambda _{min}(M)>0`$

$`\lambda _1`$ $`=`$ $`r_{\mathrm{GB}}ϵP_{n_\varphi }(r_{\mathrm{GB}}ϵ)r_{\mathrm{GB}}(n_\varphi +1)\sqrt{ϵ},`$ (E.22)
assuming that the spectrum starts at $`r_{\mathrm{GB}}ϵ,r_{\mathrm{GB}}1`$.

If $`\lambda _{min}(M)>\lambda _1`$, the convergence factor is $`t`$, but a smaller value would be more efficient. In case of $`\lambda _{min}(M)<\lambda _1`$, $`\mathrm{cosh}\theta >1`$ holds. The convergence factor is $`te^\theta `$. If $`te^\theta 1`$, the series does not converge at all. So there exists an optimal choice $`t_{opt}`$ which maps the extreme value $`A=1`$ exactly to $`\lambda _{min}(M)`$

$`\left({\displaystyle \frac{1t_{opt}}{1+t_{opt}}}\right)^2={\displaystyle \frac{\lambda _{min}(M)}{1+\delta }},`$ (E.23)
where the expansion converges with $`t_{opt}`$. Using the knowledge of $`\theta `$, we can determine $`t_{opt}`$ from

$`\left({\displaystyle \frac{1t_{opt}}{1+t_{opt}}}\right)^2={\displaystyle \frac{1+t^22t\mathrm{cosh}\theta }{(1+t)^2}}.`$ (E.24)
In reality, this is done by adjusting parameters $`t`$ and $`c`$ for the next expansion

$`tt^{}=\{\begin{array}{cc}\hfill t^{1.1}\text{if}& \theta =0\hfill \\ \hfill t_{opt}\text{if}& \theta >0\hfill \end{array}\mathrm{and}cc^{}={\displaystyle \frac{1+\delta }{(1+t^{})^2}}`$ (E.27)
using the information from the last and the fact that we expect similar behaviour for configurations close in an updating sequence.

Computer realization

1. start of basic constants

1. read  t
    2. $`c=\frac{1+\delta (n,ϵ)}{(1+t)^2}`$
2. initialise work scalars and vectors

1. $`s_0=c^{\frac{1}{2}}b`$
    2. $`s_1=0`$
    3. $`x_0=s_0`$
    4. $`N=s_0`$
    5. $`n=0`$
    6. $`\mathrm{mult}=1`$
3. recursion

1. $`\mathrm{Ms}Ms_0`$
    2. $`\mathrm{As}\frac{1+t^2}{2t}s_0\frac{1}{2tc}\mathrm{Ms}`$
    3. $`W\frac{2n+1}{n+1}\mathrm{As}\frac{n}{n+1}s_1`$
    4. $`s_1s_0`$
    5. $`s_0W`$
    6. $`\mathrm{mult}\mathrm{mult}t`$
    7. $`x_0x_0+\mathrm{mult}s_0`$
    8. $`nn+1`$
    9. $`\theta `$ from $`\mathrm{max}(\frac{s_0}{N},1)=P_n(\mathrm{cosh}\theta )`$
4. test

1. if $`(|t|e^\theta )^n<\delta _{\mathrm{rec}}`$ $``$ store $`t`$, exit
    2. else $``$ iterate the recursion

## Appendix F 4D QCD conventions

The theory is established in $`d=4`$ dimensions on a Euclidean space-time lattice with size $`L^3\times T`$. A gauge field $`U_\mu (x)SU(3)`$ is assigned to the link pointing from point $`x`$ to point $`(x+\mu )`$, where $`\mu =0,1,2,3`$ designates the 4 forward directions in space-time. The matrix defining the interaction of the fermions is

$`Q(U)_{xy}`$ $`=`$ $`c_0\gamma ^5[(1+{\displaystyle \underset{\mu \nu }{}}[{\displaystyle \frac{i}{2}}c_{\mathrm{sw}}\kappa \sigma _{\mu \nu }_{\mu \nu }(x)])\delta _{x,y}`$ (F.1)
$``$ $`\kappa {\displaystyle \underset{\mu }{}}\{(1\gamma ^\mu )U_\mu (x)\delta _{x+\mu ,y}+(1+\gamma ^\mu )U_\mu ^{}(x\mu )\delta _{x\mu ,y}\}],`$
where $`\kappa `$ and $`c_{\mathrm{sw}}`$ are parameters that have to be chosen according to the physical problem under consideration, $`c_0=[c_M(1+2d\kappa )]^1`$ and $`c_M`$ is a constant serving to optimise simulation algorithms. Typically $`\kappa 1/8`$ and both $`c_{\mathrm{sw}}`$ and $`c_M`$ are $`\mathrm{O}(1)`$. We refer to the following notation section for the definitions of the matrices $`\gamma ^\mu `$, $`\gamma ^5`$, $`\sigma _{\mu \nu }`$ and the anti-Hermitian tensor $`_{\mu \nu }(x)`$.

In order to speed up the Monte Carlo simulation, not the original matrix $`Q`$ but an even-odd preconditioned matrix $`\widehat{Q}`$ is used. We rewrite the matrix $`Q`$ as

$$Qc_0\gamma ^5\left(\begin{array}{cc}1+T_{ee}& M_{eo}\\ M_{oe}& 1+T_{oo}\end{array}\right),$$
(F.2)
where we introduce the matrix $`T_{ee}`$($`T_{oo}`$) on the even (odd) sites as

$$(T)_{xa\alpha ,yb\beta }=\underset{\mu \nu }{}[\frac{i}{2}c_{sw}\kappa \sigma _{\mu \nu }^{\alpha \beta }_{\mu \nu }^{ab}(x)\delta _{x,y}].$$
(F.3)
The off-diagonal parts $`M_{eo}`$ and $`M_{oe}`$ connect the even with odd and odd with even lattice sites respectively. Preconditioning is now realized by writing the determinant of $`Q`$, apart from an irrelevant constant factor, as

$`det[Q]`$ $``$ $`det[1+T_{ee}]det[\widehat{Q}]`$
$`\widehat{Q}`$ $`=`$ $`\widehat{c}_0\gamma ^5(1+T_{oo}M_{oe}(1+T_{ee})^1M_{eo}).`$ (F.4)
The constant factor $`\widehat{c}_0`$ is given by $`\widehat{c}_0=[c_M(1+64\kappa ^2)]^1`$, and the constant $`c_M`$ is chosen such that the eigenvalues of $`\widehat{Q}`$ are in the interval $`[1,1]`$. Since for the simulation algorithms the eigenvalues have to be positive, we finally work with the matrix $`\widehat{Q}^2`$.

Notations in 4D

Gamma matrices. The matrices $`\sigma _{\mu \nu }`$, $`\mu ,\nu =0,\mathrm{},3`$, are defined via the commutator of $`\gamma `$-matrices

$$\sigma _{\mu \nu }=\frac{i}{2}[\gamma ^\mu ,\gamma ^\nu ].$$
(F.5)
The $`44`$ $`\gamma `$-matrices are given by

$$\gamma ^\mu =\left(\begin{array}{cc}0& e_\mu \\ e_\mu ^{}& 0\end{array}\right),$$
(F.6)
with the $`22`$ matrices

$$e_0=1;e_j=i\sigma _j,j=1,2,3$$
(F.7)
and Pauli-matrices $`\sigma _j`$

$$\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).$$
(F.8)
The matrix $`\gamma ^5=\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$ is thus diagonal

$$\gamma ^5=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(F.9)
$`_{\mu \nu }(x)`$. This antisymmetric and anti-Hermitian tensor is a function of the gauge links and given by

$`_{\mu \nu }(x)`$ $`=`$ $`{\displaystyle \frac{1}{8}}[U_\mu (x)U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)`$ (F.10)
$`+`$ $`U_\nu (x)U_\mu ^{}(x+\widehat{\nu }\widehat{\mu })U_\nu ^{}(x\widehat{\mu })U_\mu (x\widehat{\mu })`$
$`+`$ $`U_\mu ^{}(x\widehat{\mu })U_\nu ^{}(x\widehat{\nu }\widehat{\mu })U_\mu (x\widehat{\nu }\widehat{\mu })U_\nu (x\widehat{\nu })`$
$`+`$ $`U_\nu ^{}(x\widehat{\nu })U_\mu (x\widehat{\nu })U_\nu (x\widehat{\nu }+\widehat{\mu })U_\mu ^{}(x)`$
$``$ $`h.c.].`$
Acknowledgement

I wish to thank

the referees Dr. Karl Jansen, Prof. Dr. Michael Müller-Preußker and Prof. Dr. Gerrit Schierholz for their taking on this extra work,

my colleagues Dr. Wolfgang Bock, Dr. Martin Hasenbusch, Dr. Thomas Kalkreuter, Dr. Klaus Scharnhorst, and Andrea Voigt in the Computational Physics group at the Humboldt Universität,

Dr. Andreas Hoferichter, Dr. Roger Horsley and Dr. Ernst-Michael Ilgenfritz in the Institut für Physik der Humboldt Universität

for the general willingness to collaborate, help and for many fruitful discussions and explanations.

For support and help with all computer problems I am indebted to Konrad-Zuse-Zentrum für Informationstechnik and to Deutsches Elektronensynchrotron – Institut für Hochenergiephysik Zeuthen. I would like to mention Dr. Hinnerk Stueben at ZIB and Dr. Hannelies Novak, Dr. Wolfgang Friebel and Dr. Alf Wachsmann at the IfH.

This work was supported by the Deutsche Forschungsgemeinschaft under grants No. WO 389/3-1 and WO 389/3-2.