{'fraction_non_alphanumeric': 0.028701594533029614, 'fraction_numerical': 0.01249050873196659, 'mean_word_length': 4.320339325388811, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Impact of large cutoff-effects on algorithms for improved Wilson fermions

As a feasibility study for a scaling test we investigate the behavior of algorithms for dynamical fermions in the N_f=2 Schroedinger functional at an intermediate volume of 1 fm^4. Simulations were performed using HMC with two pseudo-fermions and PHMC at lattice spacings of approximately 0.1 and 0.07 fm. We show that some algorithmic problems are due to large cutoff-effects in the spectrum of the improved Wilson-Dirac operator and disappear at the smaller lattice spacing. The problems discussed here are not expected to be specific to the Schroedinger functional.


Motivation
For improved Wilson fermions it has long been established that in the quenched approximation cutoff-effects at a lattice spacing of a ≃ 0.1 fm are tolerable and a continuum extrapolation can be started there. Recently more and more evidence has been accumulated that for dynamical improved Wilson fermions in a similar physical condition the cutoff-effects are much larger than expected. As an extreme example, for three flavors the existence of a phase transition in the β-κ-plane has been numerically conjectured and is interpreted as a lattice artifact [1]. A summary of large scaling violations in the two-flavor-theory is given in ref. [2].
In order to quantify those we are preparing a scaling test similar to what was done for the quenched case in ref. [3]. This will also serve as a benchmark for new actions. In the course of the scaling study we plan to calculate the axial current normalization constant Z A (g 2 0 ) and the axial current improvement constant c A (g 2 0 ) using the methods described in refs. [4] and [5], respectively.
On our coarser lattices we encountered algorithmic difficulties in both the molecular dynamics integration of the Hybrid Monte Carlo (HMC) and the efficient simulation of the canonical ensemble. We thus found it advantageous to deviate from importance sampling. Here we will discuss these problems and their link to cutoff-effects at the infrared end of the spectrum of the Dirac operator.
This paper is organized as follows: In the remainder of Section 1 we will briefly describe our setup and give a summary of the parameters of the simulations, which we will quote in the following sections. In Section 2 we describe in more detail the problems encountered at large coupling and also discuss methods to address these. In this context we study the behavior of the Polynomial Hybrid Monte Carlo (PHMC) algorithm [6,7] in this situation and find it a very useful tool for a detailed investigation of the properties of the small eigenvalues of the Dirac operator.
Section 3 is devoted to a comparison of the spectrum with the quenched case at matched physical parameters and in Section 4 we present results from an exploratory run at a smaller lattice spacing. We close with a brief summary of our findings.

Setup
All our simulations were performed in the Schrödinger functional (SF) setup [8,9]. We use non-perturbatively improved Wilson fermions [11,10,12,13] and the plaquette gauge action. For N f = 2 the clover coefficient c SW at β = 5.2 has been set to the value suggested in [12] and recently confirmed by the JLQCD Collaboration [13]. For the additional boundary-improvement coefficients needed in the SF we used the perturbative values for c t (2-loop) [14] andc t (1-loop) [15]. The axial current improvement constant c A is also set to its 1-loop value [16]. The first algorithm used is the HMC with two pseudo-fermion fields as proposed in ref. [17]. We want to note that the physical situation here is quite different from the one where this algorithm was previously applied and its performance tested by the ALPHA Collaboration [18]. In our planned scaling study we are interested in intermediate size physical volumes and lattice spacings between 0.1 and 0.05 fm. As mentioned above the second algorithm we employed is the PHMC, which we will discuss in some detail in Section 2. Apart from global sums all our calculations are carried out in single-precision arithmetics.

Simulation parameters
In Table 1 we list the lattice sizes and bare parameters of our simulations. In all cases we have T = 9/4 L for the time extension T . The bare quark mass m is defined in the appendix of ref. [19]. In the algorithm column 'H 2 ' refers to HMC with two pseudo-fermion fields and 'P n ' stands for PHMC with a polynomial of degree n. The trajectory length is always equal to one and the molecular dynamics integration step-size is denoted by δτ . For each simulation we ran 16 independent replica to gain more statistics. Concerning the SF parameters we work with zero background field and periodic spatial boundary conditions (θ = 0).  2 Sampling problems on coarse lattices

Instabilities in the molecular dynamics integration
Algorithms making use of molecular dynamics (MD) require a numerical integration of the equations of motion. Along a trajectory the Hamiltonian is then only conserved up to powers of the step-size δτ employed in the integration. Apart from these small deviations, under certain conditions the currently used integration schemes can become unstable and produce very large Hamiltonian violations ∆H. For a more detailed discussion see ref. [20], where a connection between these instabilities and large driving forces in the MD is proposed in analogy to a harmonic oscillator model. In this model the integrator becomes unstable when the product of the force and the integration step-size exceeds a certain value. The reversibility of the numerical integration is needed to prove detailed balance for these algorithms, which in turn implies that e −∆H = 1. Here one should note that the average is taken over all proposed configurations (see ref. [18]). Therefore this quantity is also sensitive to those, which were rejected in the Metropolis step following the MD integration, i.e. trajectories resulting in a large value of ∆H. In a histogram of e −∆H these contribute to bins close to zero while the distribution is peaked around one. They can also lead to an unusual autocorrelation of this quantity, making the Monte Carlo error estimate difficult. 1 In particular this applies also to the integrated autocorrelation time of e −∆H itself. This is due to the long periods of rejection in the Metropolis step, which sometimes follow large ∆H values. Fig. 1 shows a histogram of e −∆H and also its integrated autocorrelation time from one of our simulations. In this data set there are several series of large ∆H values, during which the proposed configurations were rejected. In the distribution of e −∆H these lead to an additional peak close to zero. One also sees from the right-hand plot that e −∆H is noticeably autocorrelated only when a large number of proposals were rejected in a row. As argued above in these cases the error of τ int could be underestimated. These two effects might cause some concerns when using e −∆H − 1 as an indicator for the absence of reversibility violations [18].
Spikes in ∆H have been observed by several collaborations using (improved) Wilson fermions in various setups (e.g. different gauge actions and volumes) at relatively large lattice spacings [20,12,22,23]. There these spikes have been traced back to large values of the driving force in the MD evolution and also their dependence on the quark mass has been investigated. Here we want to clarify a point, which is essentially implied by the previous observations [22,24], namely the strong correlation between spikes in ∆H and small eigenvalues of the Dirac operator. 2 In this way we hope to be able to separate physical effects from cutoff-effects, i.e. the occurrence of unphysically small eigenvalues. In Fig. 2 we clearly see a long period of rejection (corresponding to the rightmost data point in Fig. 1) caused by the presence of a very small eigenvalue. Although we did not measure them, this is expected to produce large fermionic contributions to the driving forces since they involve an inverse power of the Dirac operator.
We found the observed average λ min to be close to its tree-level estimate with Schrödinger functional boundary conditions [9]. However, the smallest λ min is an order of magnitude below that and we therefore consider these eigenvalues unphysical and will later establish their nature as cutoff-effects.
Finally, following the procedure of ref. [18], the absence of global reversibility violations is explicitly verified even for trajectories resulting in large values of ∆H. Nevertheless our experience shows that the increased cost of using a smaller δτ such that no long periods of rejection occur is more than compensated by the reduction in autocorrelation time of all observables. The reason is that already a small decrease of the integration step-size greatly reduces the Hamiltonian violations. For example, repeating run I with a step-size of 1/20 instead of 1/16, the longest period of rejection was 4 (instead of 47) consecutive trajectories.

MC estimates of fermionic observables
We concluded in the previous section that unphysically small eigenvalues ofQ 2 produce algorithmic problems only on a practical and not on a theoretical level. But apart from slowing down the algorithm these small eigenvalues also cause problems in the MC evaluation of fermionic Green's functions.
Consider the Schrödinger functional correlation function f 1 as defined in ref. [19]. It is the correlation between pseudo-scalar composite fields at the first and last time-slice, respectively. We will denote its value on a given gauge field configuration by [φ 1 ] F . Fig. 3 shows the MC history of the normalized [φ 1 ] F for the 16 replica of run II. Here t MD refers to the molecular dynamics time for each replicum. While on this scale the bulk of the data are below one and hence not visible there are several peaks, which have a big influence on the mean value. These spikes also affect the error estimate σ(f 1 ) through both the variance and the integrated autocorrelation time [21]. For statistically accessible quantities the error should approach a 1/ √ t MD behavior in the limit t MD → ∞. In this respect we found f 1 and all other fermionic correlation functions we considered to be very hard to measure. Even when using 16 replica, this asymptotic behavior does not set in after t MD ≃ 500.
The reason is the rare occurrence of very large values of [φ 1 ] F , which appear to be correlated with small eigenvalues ofQ 2 . However, this effect is washed out by using several replica. We therefore show in Fig. 4 the MC history of [φ 1 ] F , λ min and our error estimate for f 1 for one replicum of run II with such a spike in [φ 1 ] F . Indeed, for each spike in [φ 1 ] F the smallest eigenvalue drops below its average. That the converse is not true could be ascribed to a lack of overlap of the eigenvector corresponding to λ min with the source needed to compute the quark propagator. Quantitatively, for the correlation between [φ 1 ] F and λ min we measure a value of C [φ 1 ] F ,λ min = −0.33(4) if we use all replica and −0.46(6) from the replicum shown in Fig. 4 alone. Here we used as a definition of the correlation C A,B between two observables A and B Even though in the limit of infinite statistics configurations carrying very small eigenvalues are given the correct weight, depending on the algorithm this might be badly approximated for a typical ensemble size. Similar arguments referring in particular to the HMC algorithm motivated the introduction of the Polynomial Hybrid Monte Carlo (PHMC) algorithm in refs. [7].
Hence the difficulty in measuring fermionic correlation functions might be an efficiency problem related to the choice of the algorithm. To check this conjecture we employ a second algorithm and compare ensembles generated by HMC (with two pseudo-fermion fields) with PHMC ensembles. Indeed, PHMC can be tuned in such a way that it enhances the occurrence of configurations carrying small eigenvalues, thus resulting in a better sampling of this region of configuration space. A reweighting step is introduced to render the algorithm exact. As a preparation for the following discussions we want to recall some properties and introduce the notations concerning the PHMC.

The PHMC algorithm
One of the main ideas of the PHMC algorithm is to deliberately move away from importance sampling by using an approximation to the fermionic part of the lattice QCD action. More precisely, in an HMC algorithm the inverse ofQ 2 is replaced by a polynomial P n,ǫ (Q 2 ) of degree n. Here P n,ǫ (x) approximates 1/x in the range ǫ ≤ x ≤ 1. As a consequence this algorithm stochastically implements the weight dU det P −1 n,ǫ (Q 2 )e −Sg , whereas standard HMC generates ensembles according to dU detQ 2 e −Sg with S g being the gauge part of the action and U the gauge link configuration. Denoting averages over the PHMC ensemble by . . . P , the correct sample average of an observable O can then be written as and we introduce the reweighting factor W as a (partially) 3 stochastic estimate of det{Q 2 P n,ǫ (Q 2 )}. When using Chebyshev polynomials the relative approximation error for ǫ ≤ x ≤ 1 is bounded by δ ≃ 2 exp(−2 √ ǫn). To give an impression of the rôle of ǫ and δ we plot in Fig. 5 a set of polynomials P n,ǫ (x) for typical (in our simulations) values of these parameters and compare them with 1/x in the region of small x. Depending on the smallest eigenvalue of Q 2 the parameters ǫ and n have to be tuned such that the reweighting factor does not fluctuate too much. The authors of ref. [7] suggested to take ǫ of the same order as λ min and in practice used ǫ ≃ 2 λ min and δ 0.01.
Recalling that PHMC replaces detQ 2 in the HMC weight with det P −1 n,ǫ (Q 2 ) and observing from Fig. 5 that P n,ǫ (x) is smaller than 1/x for x ≤ ǫ, the aforementioned property of enhancing the occurrence of small eigenvalues is evident. At this point we would like to note that the fermionic contribution to the driving force in the PHMC is bounded from above since P n,ǫ (x) is finite even at x = 0. In this way the polynomial provides a regularized inversion ofQ 2 , thus also addressing the problems mentioned in Section 2.1.

HMC vs. PHMC
Coming back to the comparison of samples from HMC and PHMC, we repeated run II with PHMC using a polynomial of degree 140 and ǫ = 6 ·10 −4 , resulting in δ ≃ 0.002. The ratio ǫ/ λ min turned out to be around 2.7. In Fig. 6 where ω is the normalized reweighting factor. Our error estimate of f 1 shows the expected scaling behavior as soon as the run is long enough for a reliable extraction of τ int .
We first observe that apart from removing the largest spikes the inclusion of the reweighting factor does not seem to significantly change the relative fluctuations. This means that the parameters of the polynomial have been chosen properly. Events where [φ 1 ] F assumes a value O(10) times larger than f 1 are no longer isolated as in Fig. 3 but happen frequently, which means that the PHMC algorithm can more easily explore the associated regions in configuration space. This is what allows a reliable error estimate as shown in the upper part of Fig. 6, i.e. with 16 replica the asymptotic behavior of the error sets in after t MD ≃ 100.
The advantage of using PHMC instead of HMC can be clearly seen by considering the spread of σ(f 1 ) √ t MD among the replica. We analyzed this quantity in extensions of runs II and III. The result is shown in figure Fig. 7, where the shaded areas represent the range of values covered by the 16 replica as a function of the MD time. In the limit of infinite statistics all replica should converge to the same value, which need not be the same for the two algorithms because of reweighting and different autocorrelation times. We see that the spread for the HMC data is more than twice as large as for PHMC, i.e. the error on f 1 is significantly harder to estimate with HMC. What we are suggesting here is that the algorithm should be chosen depending on the type of observables and the parameter values. From our experience we conclude that PHMC sampling might just be For the two algorithms we show the ranges covered by the 16 replica.
more effective than HMC when computing fermionic quantities that are sensitive to small eigenvalues. To gain some more insight into the difference in sampling we consider the distribution of λ min since this is where we expect the largest effect. The distributions are analyzed by treating Λ bin = χ bin (λ min ) as a primary observable. Here χ bin denotes the characteristic function of each given bin in the histogram. We then perform our normal error analysis for Λ bin , where eq. (2.2) has to be used if it is a PHMC sample. For comparison Λ bin P is also analyzed in this case.
The histograms in the upper part of Fig. 8 compare the results from 200 independent measurements produced by HMC and PHMC (runs II and III, respectively). As expected the distributions agree within errors. For the PHMC run we also plot the unreweighted histogram, i.e. Λ bin P . Here we again confirm that with the parameters we chose for the polynomial the PHMC produces more configurations with small eigenvalues than HMC. As a consequence of the reweighting the errors at the infrared end of the spectrum should be smaller for the PHMC data. This is explicitly verified in the lower part of the plot where we show the ratio of the errors on Λ bin from the two algorithms. The three symbols refer to different bin sizes. The advantage in using PHMC to sample this part of the spectrum is significant and we will make use of this in the following discussion. For run III we also show Λ bin P . Lower panel (from the same data): ratio of the error on Λ bin from HMC to that from PHMC.

Comparison to the quenched case
In the previous section we studied various problems related to the occurrence of small eigenvalues. All the data presented there were produced at bare parameter values, which correspond to relatively large quark masses and small volumes. These small eigenvalues might therefore have a different nature from the "physical" ones expected to show up in large volumes and/or close to the chiral limit. Here and in the next section we will establish them as cutoff-effects.
To this end we made an additional simulation at the parameters of run II and calculated the ten lowest-lying eigenvalues λ i , i = 1 . . . 10. In Fig. 9 the smallest eigenvalue, λ 1 , is denoted by an open symbol. It seems that while λ 2 through λ 10 form a rather compact band, the lowest eigenvalue fluctuates to very small values quite independently of the others. It is expected and has been shown numerically [25] that the spectrum of the Dirac operator depends quite strongly on the bare gauge coupling. A well-defined lower bound should be recovered close to the continuum limit only. Therefore we take the strong fluctuations of λ min as an indication for the presence of large cutoff-effects. Here we should point out that the eigenvalues of the Dirac operator are not on-shell quantities. Hence the Symanzik improvement programme does not necessarily reduce cutoff-effects here. Quenched experience even suggests that the opposite might be true [26].
The occurrence of small eigenvalues at these bare parameters poses a some- what unexpected problem in dynamical simulations. Comparing the quenched situation to the N f = 2 dynamical case, the naïve expectation is that at fixed bare parameters the probability of finding configurations with small eigenvalues should be reduced by the determinant. To us the more relevant question seems to be whether small eigenvalues are suppressed in a situation where the physical parameters (e.g. volume and pseudo-scalar mass) are kept constant. Using the quenched data from ref. [27] and the dynamical data from refs. [28] and [29] (where an estimate of r 0 /a = 5.21(6) for β = 5.2 can be found) we chose the parameters of the quenched run V such that the lattice spacing and the (large volume) pseudo-scalar mass are matched to run IV. This was found to occur at almost equal bare current quark mass (see Lm in Table 1). In Fig. 10 we compare the distributions of λ min for these two runs. Two comments are in order here: • For the dynamical run the mean value is shifted up from 1.44(1) · 10 −4 to 1.72(5) · 10 −4 . This agrees with the naïve expectation but in a physically matched comparison it is a non-trivial observation.
• The distribution itself is significantly broader compared to the quenched case and in particular it is falling off more slowly towards zero. This means that even though λ min is larger for N f = 2 the probability of finding very small eigenvalues is enhanced.
The second point, i.e. that the lower bound of λ min is less well-defined, seems to imply that at a lattice spacing of a ≈ 0.1 fm the cutoff-effects are much larger in the N f = 2 case. To substantiate this we will compare the distribution of λ min to that from a run at finer lattice spacing and matched physical parameters.

Finer lattices
Apart from cutoff-effects, in the massless theory the Schrödinger functional couplingḡ 2 is a function of the box size L only [8,9]. We measured it on a small lattice of extension L/a = 4 at β = 5.2, obtaining a value ofḡ 2 = 3.7(1). We then extrapolated to this value the L/a = 6 data used in ref. [30] as a function of β. Our result from the matching is that for the two-flavor theory a bare gauge coupling of β = 5.5 roughly corresponds to a lattice spacing, which is 1.5 times smaller than at β = 5.2.
Hoping that the algorithmic difficulties arising from cutoff-effects would be much smaller in this situation, we simulated a 12 3 ×27 lattice at this value of β (run VI) using the HMC algorithm. With the κ we chose (and ignoring the change in renormalization factors) the bare quark mass Lm is roughly matched to the heavier runs at β = 5.2. We therefore compare run VI with run III.
Normally, a constant acceptance requires a decrease of the MD integration step-size if ones goes to finer lattices at fixed physical conditions. This argument is based on the scaling of the small eigenvalues, which influence the MD driving force. We found that λ min in run VI is a factor two smaller than in run III. Nevertheless, at β = 5.5 the step-size necessary for a certain (≃ 90%) acceptance is roughly the same as at β = 5.2. This indicates that the value of δτ we had to 2) and VI (β = 5.5). While the data from the coarse lattice stretch almost to zero, the β = 5.5 data seem to have a more well-defined lower bound.
use in the HMC runs at β = 5.2 was dictated by the occurrence of extremely small eigenvalues rather than by the average smallest eigenvalue. In addition, where in run I at the same average acceptance a maximum of 47 proposals were rejected in a row, the maximum for run VI is 4 trajectories. For this reason e −∆H shows no autocorrelation and its distribution is well separated from zero. Concerning fermionic observables, we have not observed spikes and hence expect the error to scale properly. However, for an accurate estimate of the error on e.g. f 1 our present statistics is not yet sufficient. 4 The reason for these effects is the change in the distribution of λ min . To compensate for the different lattice spacing, Fig. 11 compares λ min / λ min from runs III and VI. One can clearly see that at the finer lattice spacing the probability of finding a smallest eigenvalue less than half its average is greatly reduced compared to β = 5.2. The width of the distribution is smaller in this case and in particular the spectrum is now clearly separated from zero. Quantitatively, the normalized variance of λ min is reduced from 0.18(1) to 0.13 (2).
This comparison explicitly shows that the long tail of the eigenvalue distribution we observed at a ≃ 0.1 fm, and which caused the problems we have discussed, is a cutoff-effect. Matching also run VI to a quenched simulation (run VII), we again found an upward shift of λ min for the dynamical case. In addition, at this finer lattice spacing, the tails of the distributions of λ min look already very similar to each other as shown in Fig. 12.

Conclusions
At a lattice cutoff of approximately 2 GeV we have studied the behavior and performance of HMC-type algorithms in an intermediate size volume of 1 fm 4 . We discussed problems related to the occurrence of small eigenvalues in two-flavor dynamical simulations with improved Wilson fermions. We found these small eigenvalues to be responsible for large Hamiltonian violations in the molecular dynamics. Even for integration step-sizes such that the acceptance is 80 ∼ 90% those can still cause long periods of rejection, thus degrading algorithmic performance. However, in spite of employing only single-precision arithmetics we never observed reversibility violations. In addition, those eigenvalues make the estimate of fermionic quantities very difficult. The naïve intuition is that the fermionic determinant should suppress small eigenvalues compared to the quenched case. Through a direct comparison at matched physical parameters we indeed verified that λ min is larger with two dynamical flavors. On the other hand there is no obvious expectation for the tail of the distribution and we observed that it extends further towards zero than in the quenched case. Given the infrared cutoff induced by the Schrödinger functional boundary conditions and the quark mass we interpret this as a lattice artifact. We were able to confirm this picture with a simulation at finer lattice spacing, where the spectrum turned out to have a much sharper lower bound.
In our study we found that the PHMC algorithm is more efficient than HMC (with two pseudo-fermions) in incorporating the contribution to the path integral of configurations carrying small eigenvalues. In other words, the distortion of the spectrum by cutoff-effects actually makes it advantageous to deviate from importance sampling. Also without such special problems we found PHMC at least comparable in performance to HMC (in our implementations).
We want to emphasize that the problems discussed here do not occur only in the Schrödinger functional setup. Without this infrared regulator they are expected to show up already at larger quark masses.