Reweighting towards the chiral limit
aa r X i v : . [ h e p - l a t ] M a y HU-EP 08/15; SFB/CPP-08-27
Reweighting towards the chiral limit
Anna Hasenfratz ∗ Department of Physics, University of Colorado, Boulder, CO-80309-390
Roland Hoffmann † Bergische Universit¨at Wuppertal, Gaussstr. 20, 42219 Wuppertal, Germany
Stefan Schaefer ‡ Institut f¨ur Physik, Humboldt Universit¨at,Newtonstr. 15, 12489 Berlin, Germany
Abstract
We propose to perform fully dynamical simulations at small quark masses by reweighting inthe quark mass. This approach avoids some of the technical difficulties associated with directsimulations at very small quark masses. We calculate the weight factors stochastically, usingdeterminant breakup and low mode projection to reduce the statistical fluctuations. We find thatthe weight factors fluctuate only moderately on nHYP smeared dynamical Wilson-clover ensembles,and we could successfully reweight 16 , (1 . volume configurations from m q ≈ m q ≈ ǫ − regime. We illustrate the strength of the method bycalculating the low energy constant F from the ǫ − regime pseudo-scalar correlator. ∗ Electronic address: [email protected] † Electronic address: hoff[email protected] ‡ Electronic address: [email protected] . INTRODUCTION The steady progress of simulation techniques over the last decade (see e.g. [1, 2, 3, 4, 5])as well as new insights into the reasons for algorithmic failures [6] have profoundly alteredthe status of lattice QCD: With the latest generation of supercomputers essentially all p -regime points, including the point of physical quark masses, have become accessible to directsimulation. However, in the small quark mass regime the challenges are still considerable: • Large volumes are needed for the stability of the algorithms when Wilson fermions areused [6]. • Auto-correlation times increase dramatically towards the chiral limit [7]. • Statistical fluctuations of fermionic correlators become difficult to estimate since con-figurations with large contributions become rare as small Dirac modes are more andmore suppressed.A possible solution to evade these problems is to avoid generating an ensemble with thefermionic weight of the small target quark mass but instead simulate a heavier quark andreweight to the desired ensemble.This approach solves all of the above mentioned issues: The algorithm is more efficient atthe larger mass, and smaller volumes will be sufficient from the algorithmic point of view. Ata larger quark mass the region of small Dirac eigenvalues is oversampled with respect to thetarget distribution and thus observables that receive large contributions there (e.g. pseudo-scalar correlators) will be better estimated. Previously, the Polynomial Hybrid Monte Carloalgorithm [8] has been used as an alternative way to achieve such an oversampling [9, 10].Since the variance of an observable O is again a field theoretical observable, its statisticalerror on an importance sampled ensemble depends only on the auto-correlation of O . On theother hand when reweighting is employed the error is not given by the quantum mechanicalvariance h O w i − h Ow i , with w being the normalized reweighting factor, but rather thestatistical variance h O w i − h Ow i which depends on the variance of the reweighting factoritself as well as its correlation with the observable of interest. If O and w are (strongly)anti-correlated, this controls and limits the statistical error on the reweighted ensemble.Nevertheless, if the overlap (in configuration space) of the generated and desired ensem-bles is insufficient, reweighting will break down due to the fluctuations of the reweightingfactor. This limits the range of quark mass values that can be bridged by reweighting.In this paper we work with two degenerate flavors of Wilson type fermions, thoughgeneralization to other fermionic actions is straightforward. If we have an ensemble ofconfigurations { U i } generated at bare mass m with the Dirac operator D = D ( U ; m ),we can reweight it to the ensemble that corresponds to bare quark mass m by assigning toeach configuration a weight factor w i ∝ det D † [ U i ] D [ U i ] D † [ U i ] D [ U i ] , (1)2nd calculating expectation values as h O i = P i w i O [ U i ] P i w i . (2)Since the reweighting factors and their fluctuations will increase with the volume, reweightingbecomes inefficient on very large lattices. In practice we find that the fluctuations canbe controlled and quite large volumes can be reweighted. Nevertheless, reweighting is atechnique that is most useful at small quark masses and moderate volumes – like in ǫ -regimecalculations.There are two sources of fluctuations for the weights. One is due to the small eigenmodesof the Dirac operator. These physical infrared eigenmodes contribute to the weight factoras log( w i ) (cid:13)(cid:13)(cid:13) low modes = ( m − m ) X λ λ + m + O (( m − m ) ) . (3)The suppression of configurations due to the small eigenvalues is physical. The weight factorscontrol exceptionally large contributions to quark propagators that arise on configurationswith small eigenmodes. Reweighting in fact reduces the statistical fluctuations of manyobservables when compared to the partial quenched case. The ultraviolet, large eigenvaluemodes are not physical, but due to their large number they can dominate the fluctuations.Some of these fluctuations can be removed by smoothing, and with nHYP smeared Diracoperators [5] reweighting is possible also in bigger volumes. However, even on an nHYPsmeared gauge background, the UV fluctuations are still large. They are also closely corre-lated with the fluctuations of the nHYP smeared plaquette, giving us yet an other optionto reduce the noise by absorbing it into the gauge action: Including a term proportional tothe smeared plaquette in the Lagrangian reduces the UV fluctuations of the weight factors.This latter reduction is not essential, especially at smaller mass shifts, but extends the reachof the method at the expense of introducing a (very) small pure gauge term to the action.Calculating the determinant in Eq. (1) to any reasonable accuracy can be very costly.Fortunately it is not necessary to do that, a stochastic estimator is sufficient. In Sect. IIwe describe the stochastic reweighting, and several of its improvements. Sect. III describesthe numerical tests and efficiency of the reweighting technique, and in Sect. IV we presentphysics results using reweighted configurations. II. STOCHASTIC REWEIGHTING
When the Dirac operator corresponds to Wilson or clover fermions it has the form D [ U ] = 1 − κM [ U ] , (4) here M is a suitable combination of the hopping and clover terms κ is the hopping parameter κ = (2 m + 8) − . Reweighting a configuration from κ to κ requires a weight factor w = det D † D D † D = det − (Ω) , Ω = D − D D † ( D † ) − . (5)The determinant can be calculated as an expectation value w = R D ξ e − ξ † Ω ξ R D ξ e − ξ † ξ = h e − ξ † (Ω − ξ i ξ , but obtaining a reliable estimate for w is expensive, especially when Ω is not close to one.An alternative way is to calculate only a stochastic estimator of the true weight factor anddo the average over the ξ fields together with the configuration average. A similar approachis used in the stochastic global Monte Carlo update [11, 12, 13]. We start by writing theexpectation value of an operator O [ U ] at mass m as h O i = 1 Z Z D U e − S g det( D † D ) O [ U ]= 1 Z Z D U e − S g det( D † D )det − (Ω) O [ U ] (6)= Z Z h O [ U ] e − ξ † (Ω − ξ i ,ξ , where Z Z = R D U e − S g det( D † D ) R D U e − S g det( D † D ) (7)= h e − ξ † (Ω − ξ i ,ξ and the expectation value is with respect to both the U and ξ fields at m . If we considerthe configuration ensemble { U i , ξ i } , with the gauge configurations from the original sequenceand ξ i generated independently with weight e − ξ † ξ , the expectation value becomes h O i = P i O [ U i ] e − ξ † i (Ω[ U i ] − ξ i P i e − ξ † i (Ω[ U i ] − ξ i , (8)i.e. the weight factors are replaced by a single estimator s i = e − ξ † i (Ω[ U i ] − ξ i . (9)The obvious advantage of the stochastic approach is that the averages of the noise sourcesand the gauge fields commute. Without introducing any systematic error, we thereforeneed only one estimator of the weight on each configuration. The disadvantage is that asingle estimator might fluctuate too much and introduce large statistical errors in Eq. (8).Fortunately there are several methods that can reduce the fluctuations of s i to acceptablelevels. 4 . Improving the stochastic estimator In this section we describe two methods, the determinant breakup [14, 15] and the lowmode subtraction that we use in calculating the stochastic weight factors. Both methodswere used in [13] in a different context. It will be useful for the discussion to rewrite Eq. (5)as Ω = ((1 − x )( D † ) − + x )((1 − x ) D − + x ) , (10)where x = κ /κ . We consider reweighting to smaller quark masses so x <
1, thougheverything works for x > − x ) ∝ ( m − m ), i.e. apartfrom a gauge field independent additive constant x , Ω ∝ ( m − m ). The exponent of thestochastic estimator in Eq. (9) now can be written as ξ † (Ω − ξ = (1 − x ) ξ † ( D † ) − D − ξ + x (1 − x ) ξ † (( D † ) − + D − ) ξ + ( x − ξ † ξ , (11)which requires one inversion of the Dirac operator D on ξ .
1. Determinant Breakup
The determinant breakup is based on Ref. [15] and has been used extensively in dynamicalsimulations. The idea is to break up the interval ( κ → κ ) to N sections κ → κ + ∆ κ → .... → κ + N ∆ κ = κ , and write the determinant as the product of N terms, each on a ∆ κ interval. Thus the stochastic estimator takes the formexp (cid:16) − N X n =1 ξ † n (Ω n − ξ n (cid:17) , (12)where Ω n is the analogue of Eq. (5) on the κ n → κ n +1 interval. Since the operators Ω n aremuch closer to one then the original Ω, the estimator in Eq. (12) has reduced fluctuations.Calculating it on N intervals requires N applications of D − , so it is expensive. On theother hand it predicts a stochastic weight factor for the configuration at any intermediate κ n value, so the same calculation can be used to reweight to many different mass values.
2. Separating the low eigenmodes
Eqs. (3) and (11) clearly show that the low eigenmodes of the Dirac operator not onlygive large contributions to the reweighting factor but they can dominate the stochasticfluctuations as well. This latter problem can be reduced by explicitly calculating the loweigenmodes and removing them from the stochastic estimator. It is not necessary to workwith the exact eigenmodes of Ω, subtracting approximate eigenmodes works as well.Assume that P is an arbitrary but Hermitian projection operator, P = P , P † = P , and¯ P = 1 − P is the complementary projector. Any operator can be decomposed asΩ = P Ω P P
Ω ¯ P ¯ P Ω P ¯ P Ω ¯ P ! T ! P Ω P Q ! R ! with T = ¯ P Ω P ( P Ω P ) − , R = ( P Ω P ) − P Ω ¯ P and Q = ¯ P Ω ¯ P − ¯ P Ω P ( P Ω P ) − P Ω ¯
P . (13)Now the determinant can be written asdet Ω = det ( P Ω P ) det Q , (14)so the sectors projected by P and ¯ P are separated. If the projection operator P is builtform the eigenvectors of Ω, the second term in Eq. (13) vanishes and Q = ¯ P Ω ¯ P . If P isbuilt only from approximate eigenmodes, both terms in Q are present but the second termgives only a small correction. We will refer to it in the following as correction term. Separating Dirac eigenmodes:
Due to the γ Hermiticity of the Wilson Dirac operator theeigenvalues of D come in complex conjugate pairs and the eigenvectors are γ orthogonal.The eigenvectors of the massless operator D | v λ i = λ | v λ i are the eigenvectors of the massiveone as well, with eigenvalues λ + m . The γ orthogonality implies that one can normalizethe eigenvectors such that h v λ ′ | γ | v λ ∗ i = δ λ ′ λ and the operator P λ = | v λ ih v λ ∗ | γ is a projector. Since P λ is not Hermitian, the discussion of the previous paragraph does notapply. Nevertheless if P λ is built form the eigenvectors of Dirac operator, Eq. (14) is validwith vanishing correction term. If P projects to a few low energy modes of D , the first term,det( P Ω P ), can be calculated explicitly, and the second one, det ¯( P Ω ¯ P ), can be estimatedstochastically using Eq. (11). Since the eigenvectors of the massless Dirac operator can beused at all mass values, the overhead of the low mode subtraction is a one time calculation ofthe eigenmodes. If those modes are subtracted during the inversion, some of the additionalcost can be recouped by the improved convergence as well. Separating Hermitian eigenmodes : An alternative approach is to construct the projectionoperator from the eigenmodes of the Hermitian operator H = γ D .
Eqs. (5) and (10) in terms of the Hermitian operator becomeΩ = H − H H − = (1 − x ) H − + x (1 − x )( γ H − + H − γ ) + x . (15)Now the projector is Hermitian, but it is not constructed from the eigenmodes of Ω so thecorrection term in Eq. (13) is necessary.Working with n normalized eigenvectors of H | w i i = η i | w i i the projector is P = P ni =1 | w i ih w i | , and P Ω P = ω ij | w i ih w j | igure 1: The topological charge of the original configurations. is an n × n matrix with coefficients ω ij = (1 − x ) δ ij η i + x (1 − x ) h w i | γ | w j i ( η − i + η − j ) + x δ ij . (16)Both the determinant and the inverse of P Ω P is easily calculable. The leading term of thestochastic estimator, ¯ P Ω ¯ P requires the evaluation of H − ¯ P | ξ i . Since the inversion is doneon the subspace that is orthogonal to the low eigenmodes, this can be considerably fasterthan calculating H − | ξ i . The correction term does not require a new inversion if one usesthe identity H − | w i i = 1 /η i | w i i .The eigenmodes of H have to be recalculated for every κ interval. While the initialcalculation of the eigenmodes is usually expensive, changing κ by a small amount does noteffect the eigenmodes much, and starting from nearly correct eigenmodes the calculationconverges fast. An alternative is to use the eigenmodes of H of the first interval throughout,but than one needs to calculate H − | w i i on every interval. The inversion on the low modesis expensive and we found recalculating the eigenmodes on each interval is a better option. III. NUMERICAL TESTS
We have tested the reweighting on a set of 16 configurations generated with a 2-flavornHYP clover action. At the original parameter values, β = 7 . κ = 0 . a = 0 . m P CAC ≈
20 MeV. We have 180thermalized configurations separated by 5 trajectories. The topological charge, as measuredby an overlap operator based on nHYP Wilson fermions ( R = 0 . igure 2: The gap of the Hermitian Dirac operator of the original 16 , κ = 0 . β = 7 . about am q = 0 . κ = 0 . m P CAC ≈ κ = 0 . κ = 0 . ǫ -regime with Wilson fermions. Alreadythe configurations at κ = 0 . × m π L ≈ . lattices. In the ǫ − regimethe eigenmodes of both the Dirac and Hermitian Dirac operators are pushed away fromzero, while the finite volume has little effect on the PCAC quark mass. The ratio of themedian of the gap, ¯ µ , and m P CAC increases as one approaches the ǫ − regime. In infinitevolume the ratio is the renormalization factor Z A , which we expect to be near 1 with nHYPfermions. On large but finite volume the ratio might not have any physical meaning, butthere is indication that in practice it remains close to Z A [6]. Our large volume runs give¯ µ/m P CAC ≈ .
82, while the ratio here is 1 .
1, showing the effect of the finite volume. (In [5]we quoted ¯ µ/m
P CAC ≈ .
91. Those runs on 12 ×
24 lattices with La = 1 . m q = 70 MeV quark mass.) A. Reweighting with Hermitian eigenmodes
With Hermitian subtraction one has to balance the cost of calculating the low energyeigenmodes with the improved convergence of the conjugate gradient iteration. The costof the latter is proportional to the number of intervals between the starting and ending κ values, which in turn determines the statistical fluctuations of the stochastic estimator dueto the ultraviolet modes.We have found that removing more than 6 eigenmodes did not substantially increase theconvergence of the inversion, nor did it decrease the fluctuations of the estimator. Sub-tracting only 3 eigenmodes resulted in a significantly more expensive inversion that quicklyovertook the cost of calculating the eigenmodes. Therefore we have settled on subtracting 6Hermitian eigenmodes. Next we considered the optimal number of steps in the determinantbreakup. This should be such that the stochastic fluctuations of the weight factors are smallcompared to the fluctuations between the weight factors of the different configurations. Wefound that 99 intervals between κ = 0 . κ = 0 . B. Reweighting with Dirac eigenmodes
We have also tested reweighting with Dirac eigenmodes. Using the ARPACK package wecalculated 20 eigenmodes and separated ≈
16 real or complex conjugate pairs. While theseeigenmodes work on every κ interval, the conjugate gradient inversion is still expensive(we project on D modes but invert the operator D † D ) and we found this approach moreexpensive than removing Hermitian eigenmodes. Our tests were done with single precisionDirac eigenmodes and it is quite possible that with better eigenmodes the Dirac eigenmodeseparation becomes competitive. C. Reweighting factors
In this section we present results obtained using Hermitian eigenmode separation.Figure 3 shows the reweighting factors at three different κ values starting form the original κ = 0 . κ = 0 . igure 3: Reweighting factors at κ = 0 . . . configurations. Thereweighting factors are normalized such that their average is h s i = 1. There is another source for the large fluctuations, the ultraviolet noise. The fluctuationof the nHYP plaquette is a good representative of the UV noise, and it correlates closelywith the UV part of the weight factor (top panel of Figure 4 ). We define the UV part as theweight factor without the explicitly calculated low eigenmodes, i.e. the value determinedby the stochastic process. While this definition is not unique, it captures the correlationwith the nHYP plaquette and suggests that at least some of the UV noise can be removedby introducing an nHYP plaquette term in the new action by reweighting from S = S g − ln det( D † D ) to S = S g − ln det( D † D ) + ˜ β X p (3 − Re Tr U p, nHYP ) . (17)The correlation shown in Figure 4 can be captured by an nHYP plaquette term with ˜ β = − . κ = 0 . β = − . κ = 0 . κ = 0 . igure 4: The logarithm of the UV part of the reweighting factor versus the smeared nHYPplaquette (top panel), and the reweighting factors after the observed correlation is removed by annHYP plaquette term in the action (lower panel). The dashed red line in the lower panel is theoriginal reweighting factor. The data is for κ = 0 . here, both for reweighting to κ = 0 . t = n t / C ππ ( n t /
2) and s i C ππ ( n t / s i corresponding to a new action S both with and without the nHYP plaquette term(see Eq. (17)). Reweighting removes the very large spike (corresponding to a configurationwith a nearly zero eigenmode), and reweighting with the nHYP plaquette term reduces thefluctuations by an additional factor of two (observe the scale difference). The reweighteddata has considerably smaller statistical fluctuations than the partially quenched one.11 PSfrag replacements t C ( t ) β =3 . β =3 . β =3 . β =3 . Figure 5: The scalar correlator at κ = 0 . t = n t / κ = 0 . a m P CAC F (MeV)0.1278 0.0119(5) 79(3)(4)0.1279 0.0090(3) 79(4)(4)0.1280 0.0062(3) 81(8)(3)0.1281 0.0027(5) 78(10)(1)Table I: The PCAC quark mass and the low energy constant F on the original and reweightedensembles. The first error of F is statistical, the second is systematic, due to the uncertainty ofthe parameter m Σ V . IV. PHYSICAL RESULTS
Our goal in this paper is to illustrate the effectiveness of the reweighting method. Thephysical results we present in this section are preliminary, they merely serve to illustrate thepower of the reweighting technique.As we have mentioned in Sect. III, the original ensemble consists of 180 16 configurationsgenerated with 2 degenerate flavors of nHYP clover fermions at coupling β = 7 . κ = 0 . a = 0 . r /a = 4 . r = 0 .
49 fm, and the PCAC quark mass is am P CAC = 0 . m P CAC ≈
20 MeV. From the quadratic time dependence of the axial correlator [23], onour volume of V = (1 .
87 fm) we estimate m Σ V ≈ .
1, in or close to the ǫ − regime.Reweighting has no observable effect on the lattice spacing, both with standard andwith nHYP reweighting r /a = 4 .
25 at every κ value, though the error increases to 0 .
15 at κ = 0 . m P CAC shows a linear dependence onthe bare quark mass with κ cr = 0 . ǫ − regime NLO chiral perturbation theory predicts a quadratic form for the mesoncorrelators. For the pseudoscalar meson the prediction to O ( ǫ ) is [24, 25, 26] G P ( t ) = 1 L s Z d x h P ( x ) P (0) i = a p + L t L s b p h ( tL t ) + O ( ǫ ) , (18)where a p = Σ ρ I (2 m Σ V ρ ) (19) b p = Σ F (1 − I (2 m Σ V ρ ))are related to the chiral low energy constants Σ and F , while ρ = 1 + 3 β F √ V igure 7: The m P CAC quark mass as a function of the bare mass, m q = (1 /κ − /κ cr ) / κ cr = 0 . is the shape factor with β = 0 . I can be expressedin terms of Bessel functions, I ( u ) = 8 Y ′ ( u ) / ( uY ( u )). The function h ( τ ) = 12 [( τ −
12 ) −
112 ] (20)describes the quadratic time dependence. Our data follows this expected functional form inthe region t ∈ [4 , a p with 6-8% error, predicting Σ / with 1% accuracy,while the ratio a p /b p has 8-30% errors, predicting F at the 4-15% level. The large errors arenot surprising as the quadratic term measures only finite size effects, and our lattices arenot small. The errors are particularly large at κ = 0 . a p /b p as it is free of renormalization factors. To predict the low energy constantΣ we would need the renormalization factor of the pseudoscalar density. This calculationis in progress and we will report the results in a forthcoming publication [23]. In Table Iwe list the predictions for F as obtained from the pseudoscalar correlator. The first error isstatistical, the second is systematic from the uncertainty of the quantity m Σ V . The centralvalues correspond to m Σ V = 2 . κ = 0 . κ values we rescaledthis according to the PCAC quark mass. It is satisfying that the values we obtain areconsistent with each other, suggesting that all four data sets are governed by the ǫ -regimepredictions. The value is also consistent with recent p -regime overlap action calculations,though somewhat smaller than overlap ǫ -regime calculations [27, 28, 29]. It is possible todetermine F from the eigenvalue distribution of the Dirac operator at imaginary chemical14 igure 8: The pseudoscalar correlator on the original and nHYP reweighted data sets. potential, giving consistent results, though with large finite volume corrections [30]. Othermeson correlators can be used in similar way to predict the low energy constants of thechiral Lagrangian.The advantage of using Wilson-clover fermions is that it is relatively inexpensive to createeven large volume configurations. With the reweighting technique it is possible to probe awhole range of mass values and approach the ǫ − regime without independent simulations. Atpresent we are running simulations on 24 volumes at the same lattice spacing ( L = 2 . V. CONCLUSION
Dynamical simulations with light quarks are still computationally expensive and haveto overcome technical difficulties due to large auto-correlation time, algorithmic instabilityand statistical fluctuations. In this paper we presented an alternative to direct simulations,suggesting that reweighting in the quark mass to reach the desired light mass value might be abetter alternative. We described stochastic reweighting and presented several improvements.Our numerical tests show that with reweighting one can easily approach the ǫ − regimewith Wilson-clover quarks and we presented preliminary data for the low energy chiralconstant F . Reweighting can also provide an alternative to partial quenched studies inthe p − regime, since the statistical fluctuations of the reweighted dynamical lattices can15e considerably smaller than the partially quenched data. As an example we showed howthe scalar correlator, known to become negative in partial quenched simulations, becomespositive and follows the expected theoretical form on the reweighted (and therefore fullydynamical) ensembles.Reweighting might not be efficient in large volumes, or might be limited to small massdifferences. Our experience indicates that this is not a problem on (1 .
87 fm) volumesbetween 20 MeV and 5 MeV quarks or (2 . volumes between 8 and 3 MeV quarks withnHYP smeared Wilson-clover fermions. VI. ACKNOWLEDGMENT
We thank T. Degrand for useful discussions. Some of the numerical calculations werecarried out on the KAON cluster at Fermilab and we acknowledge the allocation fromUSQCD/SciDac. This research was partially supported by the US Dept. of Energy and theDeutsche Forschungsgemeinschaft in the SFB/TR 09. [1] C. Morningstar and M. J. Peardon, Phys. Rev.
D69 , 054501 (2004), hep-lat/0311018.[2] M. Luscher, Comput. Phys. Commun. , 199 (2005), hep-lat/0409106.[3] C. Urbach, K. Jansen, A. Shindler, and U. Wenger, Comput. Phys. Commun. , 87 (2006),hep-lat/0506011.[4] T. Takaishi and P. de Forcrand, Phys. Rev.
E73 , 036706 (2006), hep-lat/0505020.[5] A. Hasenfratz, R. Hoffmann, and S. Schaefer, JHEP , 029 (2007), hep-lat/0702028.[6] L. Del Debbio, L. Giusti, M. Luscher, R. Petronzio, and N. Tantalo, JHEP , 011 (2006),hep-lat/0512021.[7] S. Aoki et al. (JLQCD) (2008), 0803.3197.[8] R. Frezzotti and K. Jansen, Phys. Lett. B402 , 328 (1997), hep-lat/9702016.[9] M. Della Morte, R. Hoffmann, F. Knechtli, and U. Wolff (ALPHA), Comput. Phys. Commun. , 49 (2005), hep-lat/0405017.[10] K. Jansen, A. Nube, A. Shindler, C. Urbach, and U. Wenger, PoS
LATTICE2007 , 084(2006), 0711.1871.[11] A. Hasenfratz and F. Knechtli, Comput. Phys. Commun. , 81 (2002), hep-lat/0203010.[12] A. Alexandru and A. Hasenfratz, Phys. Rev.
D66 , 094502 (2002), hep-lat/0207014.[13] A. Hasenfratz, P. Hasenfratz, and F. Niedermayer, Phys. Rev.
D72 , 114508 (2005), hep-lat/0506024.[14] M. Hasenbusch, Phys. Rev.
D59 , 054505 (1999), hep-lat/9807031.[15] M. Hasenbusch, Phys. Lett.
B519 , 177 (2001), hep-lat/0107019.[16] A. Hasenfratz, R. Hoffmann, and S. Schaefer, JHEP , 071 (2007), 0709.0932.[17] S. Schaefer, A. Hasenfratz, and R. Hoffmann, PoS LATTICE2007 , 132 (2006), 0709.4130.[18] R. Hoffmann, A. Hasenfratz, and S. Schaefer, PoS
LAT2007 , 104 (2007), 0710.0471.
19] A. Stathopoulos, SIAM J. Sci. Comput. , 481 (2007).[20] A. Stathopoulos and J. R. McCombs, SIAM J. Sci. Comput. , 2162 (2007).[21] R. Sommer, Nucl. Phys. B411 , 839 (1994), hep-lat/9310022.[22] A. Hasenfratz, R. Hoffmann, and F. Knechtli, Nucl. Phys. Proc. Suppl. , 418 (2002),hep-lat/0110168.[23] A. Hasenfratz, R. Hoffmann, and S. Schaefer (2008), in preparation.[24] P. Hasenfratz and H. Leutwyler, Nucl. Phys.
B343 , 241 (1990).[25] F. C. Hansen, Nucl. Phys.
B345 , 685 (1990).[26] F. C. Hansen and H. Leutwyler, Nucl. Phys.
B350 , 201 (1991).[27] J. Noaki et al. (JLQCD), PoS
LAT2007 , 126 (2007), 0710.0929.[28] H. Fukaya et al. (JLQCD) (2007), 0711.4965.[29] S. Necco, PoS
LATTICE2007 , 021 (2006), 0710.2444.[30] T. DeGrand and S. Schaefer, Phys. Rev.
D76 , 094509 (2007), 0708.1731., 094509 (2007), 0708.1731.