Dynamical fermion algorithm for variable actions

A new version of the two-step multi-boson algorithm is developed with different fermion actions in the multi-boson and noisy Metropolis steps.


Introduction
Multi-boson algorithms for dynamical fermions [1,2,3,4] offer promising possibilities for numerical simulations of QCD and other similar theories [5]. The two-step multi-boson (TSMB) algorithm [3] has been succesfully applied, for instance, in supersymmetric Yang-Mills theories [6] and for simulations of QCD with SU(2) colour at high densities [7]. In the present paper I shall consider a TSMB algorithm allowing for different actions in the two steps.
For convenience of the reader let us first shortly summarize the main features of TSMB. Let us consider the case of an arbitrary number of identical fermion flavours N f and assume the existence of a hermitean fermion matrixQ, which has the determinant det(Q) appearing in the effective action for the gauge field after the integration over the fermionic variables in the path integral. Multi-boson algorithms are based on the representation [1] .
Here the polynomial P n satisfies lim n→∞ P n (x) = x −N f /2 (2) in an interval x ∈ [ǫ, λ] covering the spectrum ofQ 2 . For the multi-boson representation of the determinant one uses the roots of the polynomial r j , (j = 1, . . . , n) P n (Q 2 ) = r 0 n j=1 (Q 2 − r j ) .
Assuming that the roots occur in complex conjugate pairs, one can introduce the equivalent forms where r j ≡ (µ j + iν j ) 2 and ρ j ≡ µ j + iν j . With the help of complex boson (pseudofermion) fields Φ jx one can write Since for a finite polynomial of order n the approximation in (2) is not exact, one has to extrapolate the results to n → ∞.
The difficulty for small fermion masses in large physical volumes is that the condition number λ/ǫ becomes very large (10 4 − 10 6 ) and very high orders n = O(10 3 ) are needed for a good approximation. This requires large storage and the autocorrelation of the gauge configurations becomes very bad since the autocorrelation length is a fast increasing function of n. One can achieve substantial improvements on both these problems by introducing a two-step polynomial approximation: The multi-boson representation is only used here for the first polynomial P (1) n 1 which provides a first crude approximation and hence the order n 1 can remain relatively low. The correction factor P (2) n 2 is realized in a stochastic noisy Metropolis correction step with a global accept-reject condition during the updating process. In order to obtain an exact algorithm one has to consider in this case the limit n 2 → ∞.
In the two-step approximation scheme for N f flavours of fermions the absolute value of the determinant is represented as After an update sweep over the gauge field a global accept-reject step is introduced in order to reach the distribution of gauge field variables [U] corresponding to the right hand side of (7). This can be done stochastically by generating a random vector η according to the normalized Gaussian distribution and accepting the change where The Gaussian noise vector η can be obtained from η ′ distributed according to the simple Gaussian distribution by setting it equal to η = P (2) In order to obtain the inverse square root on the right hand side of (12), one can proceed with a polynomial approximation This is a relatively easy approximation because P (2) n 2 (x) − 1 2 is not singular at x = 0, in contrast to the function x −N f /2 . A practical way to obtain P (3) is to use some approximation scheme for the inverse square root. A simple possibility is to use a Newton iteration The second term on the right hand side can be evaluated by a polynomial approximation as for P (2) in (6) with N f = 0 and P (1) → P k P (2) . The iteration in (14) is fast converging and allows for an iterative procedure stopped by a prescribed precision. A starting polynomial P can be obtained, for instance, from P (2) in (6) with N f → − 1 2 N f . The TSMB algorithm becomes exact only in the limit of infinitely high polynomial order: n 2 → ∞ in (6). Instead of investigating the dependence of expectation values on n 2 by performing several simulations, it is better to fix some relatively high order n 2 for the simulation and perform another correction in the "measurement" of expectation values by still finer polynomials. This is done by reweighting the configurations.
The reweighting for general N f is based on a polynomial approximation P (4) n 4 which satisfies lim The interval [ǫ ′ , λ] can be chosen by convenience, for instance, such that ǫ ′ = 0, λ = λ max , where λ max is an absolute upper bound of the eigenvalues ofQ 2 . In practice, instead of ǫ ′ = 0, it is more effective to take ǫ ′ > 0 and determine the eigenvalues below ǫ ′ and the corresponding correction factors explicitly [6]. With properly chosen [ǫ ′ , λ] the limit n 4 → ∞ is exact on an arbitrary set of gauge configurations. For the evaluation of P (4) one can use n 4 -independent recursive relations [8], which can be stopped by observing the required precision of the result. After reweighting, the expectation value of a quantity O is given by where η is a simple Gaussian noise like η ′ in (11). Here . . . U,η denotes an expectation value on the gauge field sequence, which is obtained in the two-step process described above, and on a sequence of independent η's. The necessity and importance of reweighting the gauge field configurations depends on the condition number λ/ǫ, which grows roughly as the squared inverse of the fermion mass in lattice units. For relatively heavy fermions one can choose n 2 high enough, such that the effect of reweighting is negligible compared to the statistical errors. This can be checked on a subsample of statistically independent configurations and then the systematic errors due to the finiteness of n 2 are under control. The reweighting becomes necessary only for very light fermions as, for example, in [7]. In cases if reweighting becomes important the computational work for obtaining the reweighting factors is comparable for the calculation of the inverse ofQ 2 on the vectors η in (16). This is typically less than the off-line calculations performed, for instance, for determining some correlators and their matrix elements. Of course, if the reweighting has some effect, then it has to be taken into account in the process of estimating statistical errors (see [9]).

Variable multi-boson update step
The TSMB algorithm is based on the fact that the change of the fermion determinant in the update sequence can be approximated by a substantially lower order polynomial than would be required for the final precision of the simulation. More generally, one can also allow for the fermion action in the multi-boson step to differ from the true action describing the model to be simulated. (This has been done in a special case in ref. [10] where the update step is performed with the pure gauge action.) Of course, the correction steps become in general more involved because one has to correct for the difference of the two actions, too.
Let us denote the auxiliary hermitean fermion action in the multi-boson step in general byQ 0 . Instead of (7) the fermion determinant is now represented by .
The main difference compared to (7) is that here an additional polynomial (P (2) ) appears which is needed in order to compensate for the use of a different actionQ 0 in the first polynomialP (1) . The polynomials appearing in (17) have to satisfy, with R (..) ≃ 1 and R (..) ≃ 1, the relations Here R (2) (x) ≃ 1 andR (12) (x) ≃ 1 have to be much better approximations thanR (1) The factor (detP (1) ) −1 in (17) is generated by the multi-boson update step using the actionQ 0 . The remaining factor, which has to be reproduced by the noisy Metropolis step, can be written, for instance, as Note that here and in what follows the rôles of P (2) andP (2) can also be exchanged. Using the form in (19) one can write the detailed balance condition for the acceptance probability P A in (9) as (20) Here the gauge field dependence of the polynomials is displayed, which reflects the dependence of the fermion actions on the gauge field. The detailed balance condition can be satisfied similarly to (9) with (21) The required Gaussian can now be obtained from a simple one in (11) by For the evaluation of (21) we need the polynomial approximations and (24) The measurement correction can be performed in an analogous manner as before. For this we need the polynomial approximations P (4) (x) andP (4) (x) corresponding to (18): The final precision is given now by the quality of the approximations R (24) (x) ≃ 1 and R (124) (x) ≃ 1.

An application
Apart from the doubling of the number of polynomials in the noisy Metropolis correction step the present algorithm is entirely analogous to normal TSMB. As a first application let us here consider the case where the two hermitean actionsQ 0 andQ are the same Wilson fermion action, except for the values of the parameters β (the gauge coupling) and κ (the hopping parameter). Let us denote the parameters in the multi-boson step by (β 0 , κ 0 ) and in the noisy Metropolis step by (β, κ). An interesting question is how for fixed parameters (β, κ) the change of (β 0 , κ 0 ) does influence the behaviour of the algorithm. To see this I performed a test run with two flavours of Wilson fermions on an 8 3 · 16 lattice at parameters (β = 5.28, κ = 0.160). (This corresponds to a point on the N t = 4 cross-over line of [11].) The main polynomial parameters were as follows: [ǫ, λ] = [0.00875, 2.8], n 1 = 24, n 2 = 70.
The change of the acceptance rate in the noisy Metropolis step for changing β 0 resp. κ 0 is shown by figure 1. As one can see, the acceptance remains quite good in a relatively broad range of parameters. This is partly due to the fact that the distribution of the exponent in the noisy Metropolis step becomes substantially broader when the multiboson update parameters (β 0 , κ 0 ) differ from (β, κ) (see figure 2).
Changing (β 0 , κ 0 ) can be used to introduce more randomness in the update process. In fact, since the final distribution of gauge configurations does not depend on the parameters of the multi-boson update, during a run one can randomly change (β 0 , κ 0 ). This improves the autocorrelations. As an example, the measured integrated autocorrelation of the plaquette in the original TSMB run with the above polynomial parameters is τ plaq int = 2.9(3) · 10 5 MV M (Matrix-Vector Multiplications by the even-odd preconditioned Wilson-Dirac matrix). This is decreased to τ plaq int = 1.2(3) · 10 5 MV M if (for fixed κ 0 = κ) β 0 is changed randomly during the gauge update in the interval β 0 ∈ [5.16, 5.40]. For counting the number of MVM's in an update cycle the following approximate formula is used: Here N B is the number of boson field updates, N G the gauge field updates and N C the number of noisy Metropolis accept-reject steps. The autocorrelations were determined in independent runs which were at least as long as hundred times the measured autocorrelation. In these test runs with relatively heavy quarks the polynomial order n 2 is high enough and the effect of reweighting in (16) is much smaller than the statistical errors. Therefore, the autocorrelations were determined without taking into account the reweighting factors. The decrease of the plaquette autocorrelation as an effect of updating with a "wrong" action is at the first sight against intuition. This unexpected effect can be understood as due to partly compensating the dominance of the bosonic contributions against the pure gauge contributions in the effective gauge action. As a consequence of the gauge field part in (21), the fluctuation of the gauge fields within the natural band of fluctuations is intensified. At the same time there is some decrease of the acceptance, too, but in total the effect of the amplified fluctuations is more important.
An interesting question is the volume dependence of the effect of action variations in the update. In general, the allowed action changes may be expected to decrease in larger lattice volumes. However, one has to observe that the range with good acceptance in figure 1 is much larger than the range which would be allowed for entirely updating with different parameters and performing the necessary reweighting afterwards on the equilibrium configurations. The main reason is that the distance between two noisy Metropolis steps is typically less than 1% of the plaquette autocorrelation distance. In addition, on larger volumes the number of boson fields n 1 is also larger and the dominance of the bosonic part in the effective gauge action is stronger. Therefore the question of the volume dependence is not easy to answer and can be best done by performing actual test runs on larger volumes [12].
Another interesting possibility is to apply TSMB for variable actions in numerical simulations with improved fermion actions. Using a simplified version of the particular improved action in the multi-boson update step facilitates the implementation of TSMB algorithms and may also improve the autocorrelations. Figure 1: The average acceptance of the noisy Metropolis step at (β = 5.28, κ = 0.16). On left (right): changing β 0 (κ 0 ) in the multi-boson update step.