An Algorithm for Gluinos on the Lattice

TL;DR

This paper presents a novel Monte Carlo algorithm for simulating supersymmetric Yang-Mills theories with fermions on the lattice, combining L"uscher's local bosonic approach with polynomial correction steps.

Contribution

It introduces a new algorithm that enhances lattice simulations of fermionic theories, potentially outperforming traditional molecular dynamics methods.

Findings

The algorithm is promising and competitive with existing methods.

Application of polynomial approximation schemes to hopping parameter expansions.

Discussion on optimizing the polynomial approximation for fermionic simulations.

Abstract

L\"uscher's local bosonic algorithm for Monte Carlo simulations of quantum field theories with fermions is applied to the simulation of a possibly supersymmetric Yang-Mills theory with a Majorana fermion in the adjoint representation. Combined with a correction step in a two-step polynomial approximation scheme, the obtained algorithm seems to be promising and could be competitive with more conventional algorithms based on discretized classical (``molecular dynamics'') equations of motion. The application of the considered polynomial approximation scheme to optimized hopping parameter expansions is also discussed.

Tables

Table 1. Table 1: The values of some simple characteristic quantities on 4 3 ⋅ 8 ⋅ superscript 4 3 8 4^{3}\cdot 8 lattice at β = 2.0 , K = 0.150 formulae-sequence 𝛽 2.0 𝐾 0.150 \beta=2.0,\;K=0.150 . The quantities are defined in the text. In the sixth column the number of combined sweeps is given in thousands. For the two-step algorithms the seventh column shows the updating probability of gauge links, and the last column contains the acceptance rate of the noisy correction step.

$(\overline{n},)n$	$\square$	$\left|P_l\right|$	$-\left({\overline{\Psi }}_x{\Psi }_x\right)$	${\chi }_2$	sweep	$p_u$	acc.
2	0.50068(19)	0.05088(23)	9.4887(24)	0.54831(93)	28
4	0.51117(16)	0.05095(23)	11.1486(45)	0.5668(10)	40
8	0.50857(38)	0.05191(51)	11.818(12)	0.5698(24)	20
16	0.50675(8)	0.0488(8)	11.908(24)	0.5657(39)	20
24	0.5075(14)	0.0524(16)	11.897(45)	0.5619(72)	15
4,24	0.50611(29)	0.0506(3)	11.160(4)	0.5561(13)	60	0.1	0.59
6,24	0.50578(32)	0.0514(3)	11.639(8)	0.5563(14)	25	1.0	0.54
8,24	0.50600(40)	0.0508(6)	11.833(12)	0.5556(21)	20	1.0	0.74

Table 2. Table 2: The integrated autocorrelations of different quantities expressed in numbers of combined sweeps on 4 3 ⋅ 8 ⋅ superscript 4 3 8 4^{3}\cdot 8 lattice at β = 2.0 , K = 0.150 formulae-sequence 𝛽 2.0 𝐾 0.150 \beta=2.0,\;K=0.150 . The quantities are defined in the text.

$(\overline{n},)n$	$\square$	$\left|P_l\right|$	$-\left({\overline{\Psi }}_x{\Psi }_x\right)$	${\chi }_2$
2	6.4(3)	1.0(1)	1.0(1)	2.8(1)
4	13.9(3)	1.4(1)	1.5(1)	5.2(1)
8	36(2)	3.3(3)	1.6(1)	15(1)
16	117(14)	8.3(3)	1.9(1)	35(4)
24	165(35)	19(2)	2.3(2)	43(5)
4,24	29(3)	4.0(3)	1.6(3)	12.8(6)
6,24	18(2)	1.7(4)	1.8(3)	7.8(7)
8,24	23(4)	2.7(5)	1.7(4)	9.1(6)

Equations

$$S_f=\sum _x\text{mbox}\{{\overline{\psi }}_x^r{\psi }_x^r-K\sum _{\mu =1}^4\left[{\overline{\psi }}_{x+\hat{\mu }}^r{V}_{rs,x\mu }(1+{\gamma }_{\mu }){\psi }_x^s+{\overline{\psi }}_x^r{V}_{rs,x\mu }^T(1-{\gamma }_{\mu }){\psi }_{x+\hat{\mu }}^s\right]\text{mbox}\}.$$

$$V_{rs,x\mu }\equiv V_{rs,x\mu }[U]\equiv 2\mathrm{Tr}(U_{x\mu }^{†}{T}_r{U}_{x\mu }{T}_s)=V_{rs,x\mu }^{\ast }=V_{rs,x\mu }^{-1T}.$$

$$S=S_g+S_f.$$

$$S_g=\beta \sum _{pl}\left(1-\frac{1}{N_c}\mathrm{Re}\mathrm{Tr}{U}_{pl}\right),$$

$${\Psi }^{(1)}\equiv \frac{1}{\sqrt{2}}(\psi +C{\overline{\psi }}^T),{\Psi }^{(2)}\equiv \frac{i}{\sqrt{2}}(-\psi +C{\overline{\psi }}^T)$$

$${\overline{\Psi }}^{(j)}={\Psi }^{(j)T}C(j=1,2).$$

$$\psi =\frac{1}{\sqrt{2}}({\Psi }^{(1)}+i{\Psi }^{(2)}),{\psi }_c\equiv C{\overline{\psi }}^T=\frac{1}{\sqrt{2}}({\Psi }^{(1)}-i{\Psi }^{(2)}).$$

$$S_f=\frac{1}{2}\sum _x\sum _{j=1}^2\text{mbox}\{{\overline{\Psi }}_x^{(j)r}{\Psi }_x^{(j)r}-K\sum _{\mu =1}^4\left[{\overline{\Psi }}_{x+\hat{\mu }}^{(j)r}{V}_{rs,x\mu }(1+{\gamma }_{\mu }){\Psi }_x^{(j)s}+{\overline{\Psi }}_x^{(j)r}{V}_{rs,x\mu }^T(1-{\gamma }_{\mu }){\Psi }_{x+\hat{\mu }}^{(j)s}\right]\text{mbox}\}.$$

$$Q_{yv,xu}\equiv Q_{yv,xu}[U]\equiv {\delta }_{yx}{\delta }_{vu}-K\sum _{\mu =1}^4\left[{\delta }_{y,x+\hat{\mu }}(1+{\gamma }_{\mu })V_{vu,x\mu }+{\delta }_{y+\hat{\mu },x}(1-{\gamma }_{\mu })V_{vu,y\mu }^T\right].$$

$$S_f=\sum _{xu,yv}{\overline{\psi }}_y^v{Q}_{yv,xu}{\psi }_x^u=\frac{1}{2}\sum _{j=1}^2\sum _{xu,yv}{\overline{\Psi }}_y^{(j)v}{Q}_{yv,xu}{\Psi }_x^{(j)u},$$

$$\int [d\overline{\psi }d\psi ]e^{-S_f}=\int [d\overline{\psi }d\psi ]e^{-\overline{\psi }Q\psi }=\det Q=\prod _{j=1}^2\int [d{\Psi }^{(j)}]e^{-\frac{1}{2}{\overline{\Psi }}^{(j)}Q{\Psi }^{(j)}}.$$

$$\int [d\Psi ]e^{-\frac{1}{2}\overline{\Psi }Q\Psi }=\pm \sqrt{\det Q}.$$

$$\mathrm{pf}(M)\equiv \int [d\varphi ]e^{-{\varphi }_{\alpha }{M}_{\alpha \beta }{\varphi }_{\beta }}=\frac{1}{N!}{ϵ}_{{\alpha }_1{\beta }_1\dots {\alpha }_N{\beta }_N}{M}_{{\alpha }_1{\beta }_1}\dots M_{{\alpha }_N{\beta }_N}.$$

$${\left[\mathrm{pf}(\frac{1}{2}M)\right]}^2=\det M.$$

$$S_f\equiv \frac{1}{2}\overline{\Psi }Q\Psi \equiv \frac{1}{2}\sum _x\text{mbox}\{{\overline{\Psi }}_x^r{\Psi }_x^r-K\sum _{\mu =1}^4\left[{\overline{\Psi }}_{x+\hat{\mu }}^r{V}_{rs,x\mu }(1+{\gamma }_{\mu }){\Psi }_x^s+{\overline{\Psi }}_x^r{V}_{rs,x\mu }^T(1-{\gamma }_{\mu }){\Psi }_{x+\hat{\mu }}^s\right]\text{mbox}\}.$$

$$⟨{\psi }_{y_1}{\overline{\psi }}_{x_1}{\psi }_{y_2}{\overline{\psi }}_{x_2}\cdots {\psi }_{y_n}{\overline{\psi }}_{x_n}⟩=Z^{-1}\int [dU]e^{-S_g[U]}\det Q[U]$$

$$\cdot \sum _{z_1\cdots z_n}{ϵ}_{y_1{y}_2\cdots y_n}^{z_1{z}_2\cdots z_n}Q[U{]}_{z_1{x}_1}^{-1}Q[U{]}_{z_2{x}_2}^{-1}\cdots Q[U{]}_{z_n{x}_n}^{-1},$$

$$Z\equiv \int [dU]e^{-S_g[U]}\det Q[U].$$

$$⟨{\Psi }_{y_1}{\overline{\Psi }}_{x_1}{\Psi }_{y_2}{\overline{\Psi }}_{x_2}\cdots {\Psi }_{y_n}{\overline{\Psi }}_{x_n}⟩=Z_M^{-1}\int [dU]e^{-S_g[U]}\sqrt{\det Q[U]}$$

$$\cdot (\det Q[U]{)}^{-1}\int [d\overline{\psi }d\psi ]e^{-\overline{\psi }Q\psi }({\overline{\psi }}_{x_1}+C{\psi }_{x_1}^T)({\overline{\psi }}_{x_2}+C{\psi }_{x_2}^T)\cdots ({\overline{\psi }}_{x_n}+C{\psi }_{x_n}^T)2^{-n},$$

$$Z_M\equiv \int [dU]e^{-S_g[U]}\sqrt{\det Q[U]}.$$

$$⟨{\Psi }_y{\overline{\Psi }}_x⟩=Z_M^{-1}\int [dU]e^{-S_g[U]}\sqrt{\det Q[U]}\frac{1}{2}\left\{Q_{yx}^{-1}+C^{-1}{Q}_{xy}^{-1}C\right\}.$$

$$⟨{\Psi }_{y_2}{\overline{\Psi }}_{x_2}{\Psi }_{y_1}{\overline{\Psi }}_{x_1}⟩=Z_M^{-1}\int [dU]e^{-S_g[U]}\sqrt{\det Q[U]}$$

$$\cdot \frac{1}{4}\sum _{z_1{z}_2}\{{ϵ}_{y_1{y}_2}^{z_1{z}_2}{Q}_{z_1{x}_1}^{-1}{Q}_{z_2{x}_2}^{-1}+C_{x_1}^{-1}{ϵ}_{x_1{y}_2}^{z_1{z}_2}{Q}_{z_1{y}_1}^{-1}{Q}_{z_2{x}_2}^{-1}{C}_{y_1}+C_{x_2}^{-1}{ϵ}_{y_1{x}_2}^{z_1{z}_2}{Q}_{z_1{x}_1}^{-1}{Q}_{z_2{y}_2}^{-1}{C}_{y_2}$$

$$+C_{x_1}^{-1}{C}_{x_2}^{-1}{ϵ}_{x_1{x}_2}^{z_1{z}_2}{Q}_{z_1{y}_1}^{-1}{Q}_{z_2{y}_2}^{-1}{C}_{y_1}{C}_{y_2}-C_{x_1}^{-1}{ϵ}_{y_1{x}_1}^{z_1{z}_2}{Q}_{z_1{x}_2}^{-1}{Q}_{z_2{y}_2}^{-1}{C}_{y_2}-C_{x_2}^{-1}{ϵ}_{y_2{x}_2}^{z_1{z}_2}{Q}_{z_1{x}_1}^{-1}{Q}_{z_2{y}_1}^{-1}{C}_{y_1}\}.$$

$$Q^{†}={\gamma }_{5}Q{\gamma }_5,$$

$$|\sqrt{\det (Q)}|=\{\det (Q^{†}Q){\}}^{1/4}=\{\det ({\tilde{Q}}^2){\}}^{1/4}.$$

$$\tilde{Q}\equiv {\gamma }_{5}Q={\tilde{Q}}^{†}.$$

$${\int }_{ϵ}^{\lambda }dx{x}^w{\left[x^{-\alpha }-P(x)\right]}^2.$$

$$\delta \equiv {\left\{{\int }_{ϵ}^{\lambda }dx{\left[1-x^{\alpha }P(x)\right]}^2\right\}}^{\frac{1}{2}}.$$

Full text

1 Introduction

Among quantum field theories the supersymmetric ones play a very distinguished rôle: they have much less free parameters than a general renormalizable quantum field theory with the same set of fields, and show remarkable non-renormalization and finiteness properties (for general references see, for instance, [1]). Even more special are supersymmetric Yang Mills (SYM) theories. The Yang Mills theory with a Majorana fermion in the adjoint representation is automatically $N=1$ supersymmetric in the massless case, at least according to perturbation theory. The supersymmetry constraints in $N=2$ SYM theory can be exploited to determine, under some reasonable assumptions, its exact solution in the low energy limit [2].

The non-perturbative properties of supesymmetric quantum gauge field theories enjoyed particular interest already in the 80’s (see the review [3] and references therein). After the recent beautiful exact results of Seiberg and Witten [2, 4] there is a revived interest because one can hope to understand better on these examples some general non-perturbative phenomana as confinement and chiral symmetry breaking.

It is obvious that numerical simulations for studying non-perturbative properties of supersymmetric quantum field theories would be both interesting and desirable. In this way the assumptions beyond the exact solutions could, in principle, be checked and the results perhaps extended. Nevertheless, up to now it turned out to be impossible to find a lattice formulation with exact supersymmetry. The only possibility for reconciling lattice regularization with supersymmetry seems to be to formulate a more general non-supersymmetric model and recover supersymmetry in the continuum limit as a result of fine-tuning the parameters. This approach was pioneered by Curci and Veneziano [5] in the simple case of $N=1$ SYM, where there is only one parameter to tune (namely, the Majorana fermion mass). More general cases can be expected to be dealt with similarly, as shown, for instance, in the case of $N=2$ SYM in ref. [6]. This way of “embedding” supersymmetric theories in more general non-supersymmetric ones is natural since in Nature supersymmetry is also broken. Besides, as argued recently [7], many of the “exotic” dynamical features of supersymmetric theories survive, if small symmetry breaking terms are added.

Since supersymmetry connects bosons and fermions, a numerical simulation also involves the numerical approximation of fermionic Grassmann integrals. This is a notoriously difficult problem, which makes numerical simulations of e. g. QCD with dynamical fermions quite hard. Nevertheless, there are several known fermion algorithms which do the job for QCD, for instance, the popular Hybrid Monte Carlo algorithm [8]. This is not directly applicable in supersymmetric models where Majorana fermions are important. Nevertheless, some related algorithms based on discretized classical “molecular dynamics” equations of motions, as for instance the one in ref. [9], can work. However, recently Lüscher proposed a conceptually rather different local bosonic algorithm [10], which is an alternative for QCD [11], and can also be extended to supersymmetric cases. As recently pointed out by Boriçi and de Forcrand [12], this algorithm gains on attractivity, if it is combined with a “noisy correction” step as proposed in the early days of fermion algorithm developments [13]. This local bosonic formulation and its combination with a noisy correction step is the basis of the fermion algorithm investigated in the present paper. The formulation and tests were done in $N=1$ SYM with SU(2) gauge group, but the methods used are obviously more general: they can be applied in many supersymmetric models with $N=1$ and $N>1$ supersymmetry and, of course, also in other models containing fermions and not related to supersymmetry.

The plan of this paper is the following: In the next section the lattice formulation of models with Majorana fermions is discussed on the example of $N=1$ SYM. In section 3 the optimized polynomial approximation scheme is introduced, which is then used in the local bosonic algorithms defined and tested in section 4. Section 5 is devoted to the discussion of the applications of the optimized polynomial approximations to hopping parameter expansions. A possible fermion algorithm based on optimized numerical hopping parameter expansion is also briefly discussed there. Finally, section 6 contains some concluding remarks.

2 Majorana fermions on the lattice

In order to define the path integral for a Yang Mills theory with Majorana fermions in the adjoint representation, let us first consider the same theory with Dirac fermions. The lattice action in case of Wilson’s lattice fermion formulation has been considered in [6]. If the Grassmanian fermion fields in the adjoint representation are denoted by $\psi^{r}_{x}$ and $\overline{\psi}^{r}_{x}$, with $r$ being the adjoint representation index ($r=1,..,N_{c}^{2}-1$ for SU($N_{c}$) ), then the fermionic part of the lattice action is:

[TABLE]

Here $K$ is the hopping parameter, the irrelevant Wilson parameter removing the fermion doublers in the continuum limit is fixed to $r=1$, and the matrix for the gauge-field link in the adjoint representation is defined as

[TABLE]

The generators $T_{r}\equiv\frac{1}{2}\lambda_{r}$ satisfy the usual normalization ${\rm Tr\,}(\lambda_{r}\lambda_{s})=\frac{1}{2}$. In case of SU(2) ($N_{c}=2$) we have, of course, $T_{r}\equiv\frac{1}{2}\tau_{r}$ with the isospin Pauli-matrices $\tau_{r}$. The normalization of the fermion fields in (1) is the usual one for numerical simulations. The full lattice action is the sum of the pure gauge part and fermionic part:

[TABLE]

The standard Wilson action for the SU($N_{c}$) gauge field $S_{g}$ is a sum over the plaquettes

[TABLE]

with the bare gauge coupling given by $\beta\equiv 2N_{c}/g^{2}$.

In order to obtain the lattice formulation of a theory with Majorana fermions, let us introduce the Majorana field components

[TABLE]

with the charge conjugation matrix $C$. These satisfy the Majorana condition

[TABLE]

The inverse relation of (5) is

[TABLE]

In terms of the two Majorana fields the fermion action $S_{f}$ in eq. (1) can be written as

[TABLE]

For later purposes it is convenient to introduce the fermion matrix

[TABLE]

Here, as usual, $\hat{\mu}$ denotes the unit vector in direction $\mu$. In terms of $Q$ we have

[TABLE]

and the fermionic path integral can be written as

[TABLE]

Note that for Majorana fields the path integral involves only $[d\Psi^{(j)}]$ because of the Majorana condition in (6). This implies for $\Psi\equiv\Psi^{(1)}$ or $\Psi\equiv\Psi^{(2)}$

[TABLE]

For a given gauge field the sign can be taken by definition to be positive, but for different gauge fields one has to care about continuity. The sign convention can be fixed, for instance, at the trivial gauge field $U_{x\mu}\equiv 1$ to be positive. Then the positive sign stays by continuity until one reaches gauge fields with $\det Q[U]=0$.

The squre root of the determinant in eq. (12) is a Pfaffian. This can be defined for a general complex antisymmetric matrix $M_{\alpha\beta}=-M_{\beta\alpha}$ with an even number of dimensions ($1\leq\alpha,\beta\leq 2N$) by a Grassmann integral as

[TABLE]

Here, of course, $[d\phi]\equiv d\phi_{2N}\ldots d\phi_{1}$, and $\epsilon$ is the totally antisymmetric unit tensor. Using the above trick with doubling the number of Grassmann variables one can easily show that

[TABLE]

This establishes the connection between Pfaffians and determinants. According to eq. (6), in case of the above fermion action the antisymmetric matrix is $\bar{Q}\equiv CQ$.

It is now clear that the fermion action for a Majorana fermion in the adjoint representation $\Psi^{r}_{x}$ can be defined by

[TABLE]

The path integral over $\Psi$ is defined by the Pfaffian ${\rm pf}(\frac{1}{2}CQ)=\sqrt{\det(CQ)}=\sqrt{\det(Q)}$. In this way the Majorana nature of $\Psi$ implies that one has to take the square root of the usual fermion determinant [5]. Note that the Pfaffian defined by (13) assigns a unique sign to the path integral for Majorana fermions, therefore there is no sign ambiguity. For positive determinant the Pfaffian is real, but if the determinant is negative then the Pfaffian has to be pure imaginary.

Expectation values of Majorana fermion fields can be calculated as follows. For the Dirac fermion fields $\psi,\overline{\psi}$ we have, as is well known,

[TABLE]

[TABLE]

with the antisymmetrizing unit tensor $\epsilon$ and

[TABLE]

In order to express the expectation value of Majorana fermion fields by the matrix elements of the propagator $Q[U]^{-1}$, one can use again the doubling trick: one can identify, for instance, $\Psi\equiv\Psi^{(1)}$ and introduce another Majorana field $\Psi^{(2)}$, in order to obtain a Dirac field. Then from eqs. (5) and (16) we obtain

[TABLE]

[TABLE]

where now

[TABLE]

In particular, for $n=1$ we have with $Q\equiv Q[U]$

[TABLE]

The important case $n=2$ can be expressed by six terms:

[TABLE]

[TABLE]

[TABLE]

The indices on the charge conjugation matrix $C$ show how the Dirac indices have to be contracted.

Since the path integral over Majorana fermion fields is defined by a Pfaffian, the sign in eq. (12) is uniquely determined. It is easy to show that the fermion matrix in (9) satisfies

[TABLE]

therefore $\det(Q)$ is always real. In principle, in numerical simulations one can take into account the phase of the Pfaffian by performing the Monte Carlo integration with the positive square root $+\sqrt{|\det Q|}$, and correct for the phase by putting it into the expectation values. Although this solves the phase problem in principle, in practice there is still a problem because the numerical evaluation of the Pfaffian from the formula (13) is rather cumbersome. The only practical possibility in a Monte Carlo process is to monitor the lowest eigenvalues of $Q$ and/or explicitly calculate $\det Q$, in order to make it probable that $\det Q$ does not cross zero, where the question of the sign change in eq. (12) and the problem of the non-trivial phase of the Pfaffian emerge. In what follows it will be assumed that one can restrict the Monte Carlo simulation to gauge configurations where ${\rm pf}(\frac{1}{2}CQ)=+\sqrt{|\det(Q)|}$.

3 Optimized polynomial approximations

An important ingredient of Lüscher’s fermion algorithm [10] is the polynomial approximation of the function $1/x$ in some positive interval $x\in[\epsilon,\lambda]$ $(0<\epsilon<\lambda<\infty)$. This can be achieved by Chebyshev polynomials, which minimize the maximum relative error (“infinity-norm”). Here we shall follow another approach by minimizing the deviation in an average sense. The problem can then be reduced to a simple quadratic minimization. The advantage of this approach is its flexibility in choosing weight factors and regions of minimization. In addition, it is easily applicable to more general cases including broad classes of functions and also two-step approximations needed in the next section. In fact, the function $1/x$ refers to $N_{f}=2$ degenerate quark flavours. In the general case one has to approximate the function $x^{-N_{f}/2}$. Since in the path integrals for one species of Majorana fermions the square root of the fermion determinant appears (see previous section), a Majorana fermion is obtained for $N_{f}=1/2$, hence we have to approximate $x^{-1/4}$. This corresponds to writing

[TABLE]

Here $\tilde{Q}$ denotes the Hermitean fermion matrix defined by

[TABLE]

(See eq. (22).)

In what follows, we shall always assume that in the numerical simulations the spectrum of $Q^{\dagger}Q$ is bounded from below, and is contained in the interval $[\epsilon,\lambda]$. In principle, this may induce some restrictions on the bare parameters $\beta,K$, but in practical numerical simulations, e. g. in QCD, this condition seems always to be satisfied. Before discussing in section 4 the ways of implementing the fermion determinant (23) in the path integral, let us introduce some classes of the necessary optimized polynomial approximations.

3.1 Single step approximation

In order to define an optimized polynomial approximation $P(x)$ of the function $x^{-\alpha}\equiv x^{-N_{f}/2}$ in the above interval, the first thing to do is to define a positive norm $\delta$ characterizing the deviation of two functions from each other. The approximation is optimal if $\delta$ is minimal. A simple possibility would be to minimize, for instance,

[TABLE]

Here $x^{w}$ is some arbitrary power giving different weights to different parts of the interval.

In what follows, we shall only consider the special case $w=2\alpha$, which corresponds to considering the relative deviation. This gives the deviation norm

[TABLE]

The square root in the definition of the norm is a matter of convention. One can, as well, minimize

[TABLE]

[TABLE]

The coefficients of the polynomial $c_{\nu}$ are assumed to be real.

The second form in (27) shows the advantage of choosing norms of the type (25)-(26): since $\delta^{2}$ is a quadratic form of the coefficients, the minimum can be easily determined to be at

[TABLE]

where from (27) we have

[TABLE]

[TABLE]

The optimized polynomial approximation is then

[TABLE]

The roots of this polynomial will be needed in the form

[TABLE]

where $\nu_{j}>0$ will be assumed. Since the coefficients $c_{\nu}$ are real, the roots come in pairs with equal $\nu_{j}$ and opposite $\mu_{j}$. One can then write

[TABLE]

For a good approximation usually quite high orders $n\gg 1$ are necessary. In this case the numerical determination of the coefficients $c_{n\nu}(\alpha;\epsilon,\lambda)$ and roots $r_{nj}(\alpha;\epsilon,\lambda)$ can become cumbersome. In particular, very high numerical precision is necessary, which is not available e. g. in Fortran. However, algebraic programs usually provide the possibility to calculate with arbitrarily high precision. For instance, using this feature of Maple, one can obtain the required precision, and use the results as input data for Fortran programs. For the calculation of the inverse matrix in (28) one can use LU-decomposition with pivoting [14]. The roots of the polynomial can be numerically determined by Laguerre iteration [15].

For gluinos on the lattice we need the power $\alpha=1/4$. In this case good approximations can be achieved by polynomials of order $n=32$ or $n=40$. An example for the spectral interval $[\epsilon=0.03,\lambda=4.0]$, which is typical in the test runs discussed in section 4.3, is illustrated in figure 1.

3.2 Two-step approximations

In the next section two-step approximations will be used for the local bosonic description of the fermion determinant. Let us assume that a first polynomial approximation $\bar{P}(x)$ of the function $x^{-\alpha}$ is already known. This can be obtained, for instance, by the method described in the previous subsection for some relatively low order $\bar{n}$. The task is now to find the best approximation in the product form $\bar{P}(x)P(x)$, with a polynomial $P(x)$ of order $n>\bar{n}$.

This can be done, in fact, quite similarly as in subsection 3.1. The equation replacing eq. (27) is now

[TABLE]

Here $\bar{c}_{\bar{\nu}}$ denotes the coefficients in $\bar{P}$. Everything remains the same as in subsection 3.1, only (29) has to be replaced by

[TABLE]

[TABLE]

The coefficients of the optimized polynomial $P(x)\equiv P_{\bar{n},n}(\alpha;\epsilon,\lambda;x)$ are given by

[TABLE]

A somewhat different type of approximation will also be needed when, for the given polynomial $\bar{P}(x)$, the optimized $P(x)$ is searched which is minimizing

[TABLE]

In this case, instead of (34), we have

[TABLE]

[TABLE]

3.3 Approximations in the complex plane

The above polynomial approximation scheme can also be extended to the complex plane. Denoting the complex variable by $(x+iy)$, it is possible to approximate the function $(x+iy)^{-\alpha}$ in some region of the complex plane. This problem occurs, for instance, in the non-hermitian variants of local bosonic algorithms [12]. The approximation regions have to cover the complex spectrum of the fermion matrix $Q$. Since, according to the relation (22), the eigenvalues of $Q$ come in complex conjugate pairs, the region of approximation is symmetric with respect to a reflection on the real axis. Natural shapes of the complex region are circle or ellipse, but a simple symmetric rectangle is also appropriate.

In this latter case, for instance, one can minimize

[TABLE]

Therefore, for $w=2\alpha$, eq. (27) is replaced by

[TABLE]

In the interesting cases, when $\alpha=N_{f}=\frac{1}{2},1,2,\ldots$, the integrals can be performed analytically and give rise to relatively simple expressions in $\lambda,\epsilon,\gamma$. The coefficients $c_{n\nu}(\alpha;\epsilon,\lambda,\gamma)$ of the optimal polynomial approximation $P_{n}(\alpha;\epsilon,\lambda,\gamma;x+iy)$ are given again by a similar expression as eq. (28).

An illustration for the positions of the roots of the polynomial $P_{n}$ for $\alpha=1$, $(\epsilon=0.1,\;\lambda=2.0,\;\gamma=1.0)$ is figure 2. The spectral region is typical for the situation in the tests discussed in section 4.3, similarly to fig. 1. As one can see, the roots wind around the rectangular region in the complex plane. The value of the deviation norm is $\delta=0.0132..$. This is relatively large for this high order, showing that the rectangular region is not optimal for this kind of approximations. Elliptical shapes are more natural from the mathematical point of view. A rectangle is, however, more convenient from the point of view of monitoring the eigenvalue ranges. For instance, one can easily determine the extremal eigenvalues of the Hermitean matrices $(Q+Q^{\dagger})$ and $i(Q^{\dagger}-Q)$, and infer from them the necessary rectangle. (The extremal eigenvalues can be determined, for instance, by gradient methods [16].) In addition, as a closer inspection shows, most of the contribution to $\delta$ comes from the areas of the rectangle near the corners, where there are typically no eigenvalues anyway. In fact, in the important central regions the approximation is several orders of magnitude better.

4 Local bosonic algorithms

Let us now turn to the question, how to implement the fermion determinant factor in eq. (23) in the path integrals appearing, for instance, in the expectation values (18)-(21). Following ref. [10] and taking some polynomial approximations, as the Chebyshev polynomials or the polynomials introduced in section 3.1, we have

[TABLE]

Here $\;\stackrel{{\scriptstyle\delta}}{{\longrightarrow}}\;$ means a polynomial approximation with deviation norm equal to $\delta$. The determinant factors can be written with the help of complex pseudofermion fields $\phi_{jx},\;(j=1,..,n)$ as Gaussian integrals:

[TABLE]

This is how the fermion determinant is represented in path integrals by a local bosonic action.

4.1 Two-step approximation with noisy correction

In order to obtain a very high precision approximation one has to take a large number of pseudofermion fields. This increases the required storage space and might also cause problems with long autocorrelations. The solution is to introduce a two-step approximation scheme, as discussed in section 3.2.

In the first step we make an approximation by the polynomial $\bar{P}(x)$ with a deviation norm $\bar{\delta}$:

[TABLE]

This is improved to a smaller deviation norm $\delta$ in the second step by multiplying with the polynomial $P(x)$:

[TABLE]

The first approximation is realized by pseudofermion fields as in eq. (41), whereas the correction factor $\det[P(\tilde{Q}^{2})]$ is taken into account by a Metropolis correction step, similarly to ref. [13]. Since this latter is realized by Gaussian i. e. “noisy” unbiased estimators, we shall call it here briefly noisy correction.

By the Monte Carlo updating process one has to reproduce the canonical gauge field distribution

[TABLE]

where $S_{g}[U]$ is the pure gauge field part of the action in (4). Using (41) to represent $\det\bar{P}$, the canonical distribution in terms of the gauge field $U$ and pseudofermion field $\phi$ can be written in the form

[TABLE]

The updating of the $\phi$ field can be straightforwardly done by heatbath and/or overrelaxation algorithms.

In order to prove detailed balance for the updating of the gauge field, we have to show that the transition probability $P([U^{\prime}]\leftarrow[U])$ of the Markov process satisfies

[TABLE]

The transition probability is given by

[TABLE]

where $P_{M}([U^{\prime}]\leftarrow[U])$ is the transition probability of the updating of the gauge field with the action $S_{g\bar{P}}[U]$, and $P_{A}([U^{\prime}]\leftarrow[U])$ is the acceptance probability of the correction step. $P_{M}([U^{\prime}]\leftarrow[U])$ is realized by standard Metropolis sweeps. It is assumed to satisfy detailed balance for the action $S_{g\bar{P}}[U]$, that is

[TABLE]

This can be achieved in different ways, for instance, by a randomly chosen order of links or by updating always one half of the links in a checkboard decomposition or by doing the sweeps pairwise with opposite orders of link updates. Comparing eqs. (46)-(48) one can see that the condition for the acceptance probability is

[TABLE]

This could be satisfied by the the standard Metropolis acceptance probability

[TABLE]

However, this would be too expensive to implement numerically.

The idea of the noisy correction [13] is to generate a random vector $\eta$ according to the normalized Gaussian distribution

[TABLE]

and to accept the change $[U^{\prime}]\leftarrow[U]$ with probability

[TABLE]

where

[TABLE]

This means that the acceptance probability is

[TABLE]

Due to the relation

[TABLE]

the detailed balance condition in (49) is satisfied.

In this form the algorithm would be still too expensive, because of the necessity to generate the distribution in eq. (51). One can, however, easily generate the simple Gaussian distribution

[TABLE]

and then obtain the above $\eta$ from

[TABLE]

The exact calculation is still too difficult, but one can obtain the result to a required approximation by the two-step polynomial approximation technique introduced at the end of section 3.2. An approximation with the same deviation norm $\delta$ as in (43) is sufficient.

According to eq. (43) the polynomial approximation $R(x)$ for $P(x)^{-1/2}$ has to obey

[TABLE]

Here we defined the approximation

[TABLE]

which can be obained by a single-step polynomial approximation defined in section 3.1 for $\alpha=-\frac{1}{2}$ in the interval $[\lambda^{-1/4},\epsilon^{-1/4}]$. After this, $R(x)$ in (58) can be determined by minimizing the deviation norm $\delta$ defined in eq. (36) for $\alpha=1/8$ and $\bar{P}(x)$ replaced by $S(\bar{P}(x))$. The weight factor is somewhat arbitrary, but a good choice is, for instance, $w=0$. This means that in the interval $[\epsilon,\lambda]$ the absolute deviation is minimized.

This algorithm based on two-step approximations is not “exact” because for any finite $\delta$ there is still some deviation from the functions required in eqs. (43) and (58-59). This is in contrast to the exact algorithms with noisy correction proposed for QCD in refs. [12, 13]. It is, however, not clear how to extend the exact algorithms to the case with square roots considered in this paper. Nevertheless the two-step approximations can easily be made very precise. Namely, provided that the second approximation is already good, the acceptance rate in the noisy correction step mainly depends on the first polynomial in (42), and only very little on the precision of the approximation achieved by the second one. The application of the polynomials in eq. (53) and (58) is not very time consuming, therefore the decrease of $\delta$ does not lead to a dangerous slowing down of the simulation. This applies both to the gluino model considered in this paper and to QCD with quarks in the fundamental representation. Therefore the two-step algorithm proposed here is also a potential candidate for simulating QCD.

4.2 Noisy estimators

The pseudofermion fields can also be used to obtain the matrix elements of the fermion propagators. The required expectation values are easily obtained during the updating process, and in this way the fermion matrix elements can be determined. The expressions in terms of the pseudofermion fields are usually called noisy estimators.

The starting relation is, for a general matrix $M$,

[TABLE]

In our case we have from this

[TABLE]

Combining this formula with eqs. (40)-(41) we obtain

[TABLE]

Here the notation

[TABLE]

is used. For instance, for the gluino condensate eqs. (20) and (62) give

[TABLE]

(Note that this is the unrenormalized condensate, which has to be subtracted and renormalized, in order to obtain a physical quantity.)

The expectation values of several pairs of fermion variables can also be obtained in a similar way. For instance, for two pairs one can use

[TABLE]

and in our case

[TABLE]

As in (62), it is straightforward to transform this in an expectation value in terms of the pseudofermion field.

4.3 Numerical simulation tests

The first tests of the algorithm described in the previous subsections were performed in the SU(2) Yang Mills model with gluinos. The bare parameters were $\beta=2.0$ and $0.125\leq K\leq 0.175$. The lattice sizes were in the range $4^{3}\cdot 8$ to $8^{3}\cdot 16$.

In order to observe the $n$-dependence of some simple global quantities and the corresponding autocorrelations, a series of tests were done on $4^{3}\cdot 8$ lattice at $(\beta=2.0,K=0.150)$. This is in a reasonable range of parameters for numerical simulations, as the obtained estimates of the value of the string tension in lattice units show (see below). For the spectral interval $[\epsilon=0.03,\;\lambda=4.0]$ was taken, which contains the eigenvalues of $\tilde{Q}^{2}$ confortably: it turned out that in equilibrium configurations the eigenvalues were always larger than 0.05 and smaller than 3.7 (they were monitored at every sweep by a gradient algorithm [16]). The chosen quantities were: the plaquette $\Box\equiv W_{1,1}\equiv\frac{1}{2}{\rm Tr\,}U_{pl}$, the absolute value of the Polyakov line in the time ($L_{t}=8$) direction $|P_{l}|$, the unrenormalized gluino condensate from the noisy estimator in eq. (64), and the $k\otimes k$ Creutz ratio with $k=2$ calculated from the Wilson loops as

[TABLE]

Here $a^{2}\sigma$ is the string tension in lattice units.

The results with $2\leq n\leq 24$ for the single-step algorithm, and with $\bar{n}=4,6,8;\;n=24$ for the two-step algorithm discussed in subsection 4.1 are shown in table 1. The corresponding integrated autocorrelations are given in table 2.

No particular effort was made to optimize the mixture of the different updating steps. This should certainly be an important part of the final optimization [11, 12], which should be done in larger scale applications. Here only the $n$-dependence was investigated for a fixed reasonably looking mixture: 1 heatbath and 6 overrelaxation sweeps for the pseudofermions followed by 2 Metropolis sweeps with 8 hits per link for the gauge field. In case of the two-step algorithm before the accept-reject correction step 2 Metropolis sweeps were done with opposite orders of links. This was repeated 5-times. In order to reach a good acceptance, at every link it was first decided whether to update it or not. By choosing the updating probability $p_{u}\leq 1$, the acceptance rate could be tuned appropriately. In table 1 the numbers of “combined sweeps” always mean the numbers of these sequences of sweeps.

Tables 1 and 2 show that for the given statistical errors $n=24$ is enough, and from the point of view of autocorrelations the algorithm with two-step approximation and noisy correction is better. In addition, compared to a single-step algorithm with the same approximation order $n$, the combined sweeps of the two-step algorithm also need less computer time. In table 2 it is amazing to observe that the physically more interesting Creutz ratio $\chi_{2}$ has definitely smaller integrated autocorrelations than the plaquette $\Box$. The noisy estimators for $-(\overline{\Psi}_{x}\Psi_{x})$ always have very short autocorrelations.

In the two-step algorithm it is an interesting question, how the acceptance probability of the noisy correction is behaving for increasing lattice sizes. In order to see this, the simulation with $\bar{n}=4,\;n=24$ was scaled up to $6^{3}\cdot 12$ and $8^{3}\cdot 16$, always with $n=24$ in the correction step. It turned out that acceptances similar to $\bar{n}=4$ on $4^{3}\cdot 8$ could be reached with $\bar{n}=6$ and $\bar{n}=8$ on $6^{3}\cdot 12$ and $8^{3}\cdot 16$, respectively. With the same mixture of sweeps and link updating probability $p_{u}=0.1$, as in table 1, the acceptance rates were 55% and 65%, respectively. This is good news, which tells that the number of pseudofermion fields for the first approximation has to grow at most linearly with the lattice extension.

Another important question is the behaviour of these local bosonic algorithms for decreasing lattice spacing and/or fermion mass. A detailed investigation of this could not be carried out here. From the data in table 1 follows that the square root of the string tension in lattice units is at $(\beta=2.0,K=0.150)$ roughly $a\sqrt{\sigma}\simeq 0.75$. Defining the physical scale by $\sqrt{\sigma}\simeq 0.45\;GeV$, we obtain $a\simeq 1.7\;GeV^{-1}$. In order to see the behaviour of the algorithm and of $a\sqrt{\sigma}$ as a function of the fermion mass ($K$), numerical simulations on $6^{3}\cdot 12$ lattices were performed at $\beta=2.0$ and $K=0.125,\;0.150,\;0.175$ with the single-step algorithm. The approximations were done by the polynomials $P_{8}(1/4;\epsilon=0.07,\lambda=3.5;x)$, $P_{8}(1/4;\epsilon=0.03,\lambda=4.0;x)$ and $P_{8}(1/4;\epsilon=0.005,\lambda=5.0;x)$, respectively. In all three cases the chosen spectral intervals comfortably covered the observed extrema of the spectrum of $\tilde{Q}^{2}$. About 1000 combined sweeps were done after equilibration. The plaquette autocorrelations were about 10, 30 and 50 complete sweeps, respectively. At $K=0.125$ and $K=0.150$ the square root of the string tension in lattice units were roughly the same as the above value. At $K=0.175$ a small decrease to a value $a\sqrt{\sigma}\simeq 0.7$ was observed.

Another question is the size of the effects of dynamical gluinos on the expectation values. In the considered range of bare parameters this is still small. A short simulation in pure gauge theory gives, for instance, at $\beta=2.0$ on $8^{3}\cdot 16$ lattice $\Box=0.50113(6)$ and $\chi_{2}=0.5490(2)$. These are only slightly below the values in table 1. Larger effects are expected at larger $\beta$ in the critical region of the hopping parameter $K$.

The experience with the local bosonic algoritms on these relatively small lattices at relatively low lattice cut-offs was quite positive. As remarked above, from the point of view of autocorrelations and computational speed the two-step algorithm with noisy correction was clearly better. Of course, detailed tests on larger lattices and at smaller lattice spacings are necessary. A direct comparison with the conventional algorithms based on discretized classical equations of motion [9] would also be interesting.

5 Other applications of the optimized polynomials

5.1 Optimized solvers for matrix inversion

The optimized polynomial approximations defined in section 3 can also be used for matrix inversion instead of the Chebyshev polynomials in Chebyshev iteration schemes [17]. In numerical simulations of fermionic field theories an important task is to calculate the fermion propagator matrix elements in a given gauge field background $Q[U]^{-1}\equiv Q^{-1}$. An approximate solution of the equation

[TABLE]

is the vector

[TABLE]

Here we used the first form of the optimized polynomial $P_{n}$ in (32), but other forms are also possible. Since we are actually considering the non-Hermitean fermion matrix $Q$, the approximation has to be optimized in the region of the complex plane $x\in[\epsilon,\lambda],\;y\in[-\gamma,\gamma]$, which covers the spectrum of $Q$. The deviation norm $\delta$, which gives the precision of the approximation in (69), can be characterized by an expression similar to (39) weighted by the spectral density of the matrix $Q$.

For illustrative purposes figure 3 shows the squared length of the residue vector

[TABLE]

obtained from (69) on a typical $4^{3}\cdot 8$ configuration at $(\beta=2.0,K=0.150)$. For comparison, the same quantity is also shown as a function of the number of iterations in the popular conjugate gradient (CG) algorithm. As the figure shows, an application of the optimized solver $P_{32}$ reduces the residuum roughly as much as 32 CG iterations.

An iterative scheme is obtained if the optimized polynomials are successively applied to the residuum vectors. The resulting cyclic iteration process [17] can be represented by the following simple iteration equations: start by defining $p_{n}^{(1)}\equiv p_{n}$ and $r_{n}^{(0)}\equiv r_{n}$ and then set for $k=1,2,\ldots$

[TABLE]

[TABLE]

$p_{n}^{(k)}$ tends to the solution $p$ of (68) in the limit $k\to\infty$ because

[TABLE]

Therefore the iteration is realizing the simple identity

[TABLE]

For an illustration of this cyclic optimized solver (COS) iteration scheme see figure 4, where instead of $Q$ the positive matrix $Q^{\dagger}Q$ is inverted. In the figure the square of the residuum vector is shown as a function of the number of matrix multiplications (in this case a CG iteration step requires only one matrix multiplication). As the figure shows, the COS iteration with a high enough order polynomial has a similar performance as CG iteration.

The application of the optimized solver $P_{n}$ requires $n$ matrix multiplications and no vector norm calculations. Therefore the amount of arithmetics in COS iterations compares favourably with the CG iteration. Nevertheless there are also known methods which are better than CG iteration [18], therefore the best choice of the inversion method at given bare parameters and lattice sizes has to be decided on a case by case basis.

The optimized solver can also be considered as an optimized hopping parameter expansion. Let us write the fermion matrix in (9) as

[TABLE]

with the hopping matrix $M$. Then we have

[TABLE]

This can be compared to the $n$’th order hopping parameter expansion of the inverse matrix

[TABLE]

The only change in (75) with respect to this is the different (optimized) choice of the roots of the $n$’th order polynomial. From the point of view of the numerical iterative hopping parameter expansion [19, 20] this is practically no change in the arithmetics, but a considerable gain in precision. In addition, there is no convergence problem because for $n\to\infty$ we always have $\delta\to 0$.

5.2 An algorithm with optimized hopping parameter expansion

The optimized hopping parameter expansion of the previous subsection can also be used in the old fermion algorithm based on the direct evaluation of the change of the fermion determinant [21, 20]. Let us remember that according to (9) the change of the fermion matrix is

[TABLE]

if the link $z\rho$ is changed according to

[TABLE]

The ratio of the two fermion determinants is

[TABLE]

The $24\otimes 24$ matrix needed for the calculation of the determinant on the right hand side can be approximately calculated by the iterative numerical hopping parameter expansion, or by the improved variant discussed above using optimized solvers for the inverse.

6 Discussion

The local bosonic algorithms for gluinos on the lattice investigated in section 4 seem suitable for numerical experiments aiming at an understanding of Yang Mills theories with massive gluinos. The study of the limit when the gluino mass goes to zero should be possible and could reveal the character of the $N=1$ supersymmetric limit by comparing to the expectations based on previous analytical work [3, 5]. In particular, the combination of Lüscher’s local bosonic approach [10] with the noisy correction method [13] turned out to be rather effective in the performed tests. Since this combined algorithm has several algorithmic parameters, the final optimization for given lattice sizes and bare parameter values is a non-trivial task, which has to be carried out in the particular applications.

The required polynomial approximations are defined in a scheme which is different from the usual approach leading to Chebyshev polynomials. The deviation norm is defined as an integral which is quadratic in the polynomial coefficients (see section 3). This allows a great flexibility in choosing weight factors and approximation regions. The optimal approximation in this average sense is well suited for the present purposes, in particular, in the two-step approximation scheme applied in section 4.

Besides the local bosonic algorithms, other applications of the optimized polynomial approximation technique of section 3 have also been briefly discussed. In particular, the use of optimized solvers for matrix inversion and the application of the optimized hopping parameter expansion for a fermion algorithm without pseudofermion fields have been considered in section 5. The test of this latter algorithm goes beyond the scope of the present paper.

An important question in these algorithms is the rôle of the deviation norm $\delta$, which is characterizing the quality of the polynomial approximations. Since exact results are obtained only in the limit $\delta\to 0$, this could necessitate an extrapolation of the expectation values to $\delta=0$. Nevertheless, in the two-step approximation scheme a relative precision with $\delta\simeq 10^{-5}-10^{-6}$ seems possible with reasonable effort, therefore the practical need for an extrapolation is not clear. For instance, such a precision comes already close to the machine precision on 32-bit computers, which is generally considered to be enough in present numerical simulations.

Although the methods used here have been formulated and tested in case of the SU(2) gauge theory with a Majorana fermion in the adjoint representation, they are also applicable in other quantum field theories with fermion fields, as QCD or Higgs-Yukawa models including also scalar fields.

**Acknowledgements **

I thank Martin Lüscher and Rainer Sommer for discussions and helpful suggestions contributing in an essential way to the content of this paper.