Low-mode deflation for twisted-mass and RHMC reweighting in lattice QCD

We propose improved estimators to compute the reweighting factors which are needed for lattice QCD calculations that rely on twisted-mass reweighting for the light quark contribution and the Rational Hybrid Monte Carlo (RHMC) algorithm for non-degenerate quark masses. This is the case for a number of modern large-scale simulations based on O(a) improved Wilson fermions. We find a significant reduction of uncertainties for the reweighting factors at similar computational cost compared to the conventional estimation. This leads to a significant increase in precision for phenomenologically relevant observables with high correlation to the low eigenmodes of the Wilson-Dirac operator in the presence of exceptionally small eigenvalues. Supplementary details regarding the spectral gap of the light quark Dirac operator on the 2+1 flavor large-volume ensembles explored in this study can be found in an accompanying appendix.


Introduction
Lattice QCD is the only known framework for the calculation of ab initio predictions in the low-energy regime of Quantum Chromodynamics, the theory of the strong interactions in the Standard Model of particle physics.The field has reached the stage where high precision lattice QCD calculations can have a significant impact on phenomenologically relevant observables.Prime examples are the strong coupling constant α s and the hadronic contributions to the anomalous magnetic moment of the muon [1,2].For a recent review of lattice results on flavor physics observables see [3].
To control all sources of systematic uncertainties, simulations with resolutions finer than 0.05 fm on lattices with spatial extents larger than 4 fm and physical quark masses are desirable and by now state of the art.The generation of gauge ensembles in lattice QCD using the Wilson formulation and the Hybrid Monte Carlo (HMC) algorithm [4] is known to suffer from instabilities in large volumes and towards physical light quark masses [5][6][7].This is caused by accidental near-zero modes of the light quark lattice Dirac operator.Twisted-mass reweighting of the light quark determinant as introduced in [8,9] greatly reduces instabilities and sampling inefficiencies and thus allows for the simulation of QCD with Wilson fermions at the physical quark mass and in large volumes.
Twisted-mass reweighting is used in a number of modern large-scale simulations with O(a) improved Wilson fermions [10][11][12][13].Among them, the CLS consortium has generated a large set of large-volume, 2 + 1 flavor gauge ensembles across six lattice spacings below 0.1 fm and a range of quark masses down to the physical values of the masses of the up and down quarks [10,14,15] using O(a) improved Wilson fermions [16,17] and the tree-level improved Lüscher-Weisz gauge action [18]. 1 Difficulties in accurately estimating the twisted-mass reweighting factor have been encountered on a number of gauge configurations, as discussed in Appendix G of [19], see also [20].As pointed out in [10], for observables that have a strong correlation with the reweighting factor, i.e., observables that are sensitive to the low eigenmodes of the light quark Dirac operator, insufficient precision in the determination of the reweighting factor can lead to large fluctuations in the Monte Carlo history and correspondingly to significantly increased statistical errors.This behavior is particularly pronounced in the case of an anticorrelation of reweighting factor and observable, since a large value of the observable and a small reweighting factor must cancel to high precision on configurations that are affected by strong fluctuations.
Traditionally, the calculation of the twisted-mass reweighting factor has been performed via stochastic estimation.It was already anticipated in [8,9] that improved estimators might be needed for high-precision studies at near-physical light quark masses.A viable solution to improve the estimation of reweighting factors with high correlation with the low eigenmodes of the Dirac operator has been introduced in [21] in the context of quark mass reweighting: It combines the exact computation of the reweighting factor in the subspace of the lowest eigenmodes of the Dirac operator with a stochastic estimation in the corresponding orthogonal subspace.Such an estimator improves the computation because the low modes contribute significantly to the reweighting factor and its stochastic fluctuations.In this work, we develop and apply low-mode deflation to the computation of twisted-mass reweighting factors.We find that the improved setup effectively solves the accuracy problems in the computation with similar or even less cost than the purely stochastic setup.
Another reweighting factor has to be considered when the Rational Hybrid Monte Carlo (RHMC) algorithm [22][23][24] is employed to simulate fermions with non-degenerate masses, as it is the case in state-of-the-art simulations, see for example [11][12][13][25][26][27][28].It corrects for the approximation of the inverse square root of the Dirac operator in the RHMC.Whereas this reweighting factor is not expected to fluctuate strongly for appropriately chosen parameters, exceptionally small eigenvalues of the strange quark Dirac operator can hinder the precise computation.We find that a low-mode deflated setup may be used in these cases to resolve all precision issues.This work is structured as follows: We discuss the computation of twisted-mass reweighting factors in section 2, where we first recapitulate the underlying regularization of the quark determinant and the stochastic estimation of the reweighting factor before incorporating exact low-mode deflation into the computation.We then turn to the case of the RHMC reweighting factor in section 3, where we again give an overview of the stochastic estimation before introducing the application of low-mode deflation.The newly developed estimators are presented in sections 2.3 and 3.3, respectively.In section 4, we present the effect of low-mode deflation on the computation of reweighting factors for the case of 2 + 1 flavor CLS ensembles.We summarize the expected impact of our findings on ongoing and future lattice calculations in section 5. Appendix A gathers some information on the behavior of the spectral gap of the light quark Dirac operator for O(a) improved Wilson fermions obtained in the course of this work.

Twisted-mass reweighting
Twisted-mass reweighting is employed to stabilize HMC simulations of lattice QCD with very light Wilson quarks.The fermion determinant is split up into two factors and only one of them is taken into account in the HMC, whereas the factor that includes the contribution of the low eigenmodes of the Dirac operator is treated as reweighting factor.Before we turn to the application of low-mode deflation to the determination of the twistedmass reweighting factor, we will briefly recapitulate the definition of the regulator and the stochastic estimation of the reweighting factor.

Twisted-mass determinant reweighting
There exists some freedom in the exact definition of the twisted-mass regulator.For definiteness, we here choose to work with the setup that is employed on the 2 + 1 flavor CLS ensembles [10], where the regularization is applied to the Schur complement of the asymmetric even-odd preconditioning of the hermitian Dirac operator [29,30], with the light quark Wilson-Dirac operator D l .The second of the two variants in [8] (see Table 1 of that reference) is employed for the regularization of the light quark determinant, as it is less affected by fluctuations from the ultraviolet part of the spectrum of the Dirac operator [9].This amounts to the replacement where µ > 0 is the infrared regulator.The reweighting factor, needed to restore the target distribution, is defined via and is computed on the gauge configurations that have been generated with the regularized determinant.When computing expectation values of primary observables O that have been computed on ensembles with the modified action, the reweighting factor W l is taken into account by the redefinition of the gauge expectation value via Frequency splitting, similar to that introduced by Hasenbusch [31,32], may be applied to the computation of W l [21].In this case, the interval between µ = 0 and µ = µ 0 is split into N sp smaller intervals with µ 0 = μ0 > μ1 > • • • > μNsp = 0 and W l can be factorized, . (2.5) The determinants entering this product can be computed stochastically from the estimator where η i are complex-valued Gaussian distributed random fields of unit variance, defined on the even lattice sites.The determinant of R is related to this stochastic estimator via up to an irrelevant constant.The variance of R can be computed from and the relative variance increases significantly when R deviates strongly from unity.In the context of [10], only one set of stochastic estimators, i.e., N sp = 1 is used in the computation of the reweighting factor.The use of frequency splitting for the computation of twisted-mass reweighting factors has been investigated in [33] on several configurations, where the relative uncertainty of the stochastic estimator was found to be above 100%.A reduction of this uncertainty to the O(10%) level could be achieved by the use of frequency splitting and a significantly increased computational effort.

Deflated determinant reweighting
Low-mode deflation may help to reduce significant statistical fluctuations of the stochastic estimator R in the presence of exceptionally low eigenvalues of the Dirac operator.Following [21], where low-mode deflation has been applied in the context of quark mass reweighting, we define a general reweighting factor W based on an operator Ω via In the following, we will only consider the case where eigenmodes of Ω, determined to high precision, are employed for the deflation.Inexact deflation has not been implemented in this work.For the computation of det(R) this would result in a reduced effort for the computation of eigenmodes but an increased number of inversions of the Dirac operator in the computation of the mixed terms.

Deflated twisted-mass determinant reweighting
In the case of twisted-mass reweighting, as introduced in section 2.1, we work with or, if we consider frequency splitting, as operator that defines the reweighting factor in eq.(2.9).Eigenmodes of the evenodd preconditioned, hermitian light quark Dirac operator Ql are therefore eigenmodes of Ω and may be used to construct the projector P and thereby deflate Ω without mixed contributions of P , Ω and P .We can compute the first term on the right hand side of eq.(2.15) from N L eigenvalues λ i of Ql via The second factor det( P Ω P ) in eq.(2.15) may be estimated stochastically, similar to the standard stochastic setup of eq.(2.7).To compute the determinant in the subspace which is orthogonal to the space of the low modes, defined via P , it is enough to project out the low modes.The only change compared to eq. (2.7) is the use of deflated random sources P η i in the inversion. 2The stochastic estimators for det( P Ω P ) are then defined by and we can compute where, in the second line, we have taken the limit of one single factor in the frequency splitting.The variance of the stochastic estimator R , according to eq. (2.8), is now solely based on the high-mode contribution to the reweighting factor.

RHMC reweighting
The use of the rational approximation when including quarks in the simulation that do not come in mass degenerate pairs requires the inclusion of an additional reweighting factor to correct for the small deviation from the target action.In the context of 2 + 1 and 2 + 1 + 1 flavor simulations with Wilson fermions, this concerns the strange and charm quark components.If non-degenerate up and down quark masses are simulated, the RHMC algorithm may be employed for all quark flavors [27,28].The idea of deflated reweighting factors, as discussed in section 2.2, may also be applied to the RHMC reweighting factor.We begin this section with a brief summary of RHMC reweighting [22][23][24], based on the notation used in [27,34], before including low-mode deflation.

The RHMC quark determinant
We consider the quark determinant for flavor q in its even-odd preconditioned form, for the application in the RHMC algorithm.The odd-odd part of the Dirac operator, D q,oo , can be directly included in the molecular-dynamics evolution.Therefore, we only need to handle the even-odd preconditioned operator, similar to the twisted-mass case.We write where R q is an operator that approximates ( D † q Dq ) −1/2 and the reweighting factor is defined from To approximate the inverse square root function 1/ √ y in the range ≤ y ≤ 1, the Zolotarev function [35] with known coefficients A, a i , of degree [n, n] may be used.In simulations based on openQCD [34], as the ones described in [10], the representation of eq.(3.4) is employed to approximate the operator R q of eq.(3.2) via The scalar factors r a , r b define the approximation region [r a , r b ], which has to cover the spectral range of the operator ( D † q Dq ) 1/2 = γ 5 Dq in order to achieve a good approximation, resulting in a small deviation from the target distribution.Due to the rather large mass of the strange quark, Qs was assumed to have a safe spectral gap, although chiral symmetry is broken for Wilson fermions at finite lattice spacing.In [36] it was shown, that this assumption is problematic at coarse lattice spacing, where configurations with negative strange quark determinants can be found.Another reweighting factor has been introduced to correct for the neglected sign and to restore the proper Boltzmann weight in the path integral.
The Zolotarev rational function can be split into several factors to achieve a frequency splitting of the quark determinant.We do not discuss this modification in the following, since the application of exact low-mode deflation can be trivially extended to the modified estimator.

RHMC reweighting
The reweighting factor W q can be determined by stochastic estimation.This can be done by writing [34], with the newly introduced operator Z q .The stochastic estimation is then performed by computing using N r random quark fields η j (x) with normal distribution, defined on the even lattice sites.The use of Z q is profitable here, because one can use the series to evaluate the square root in the stochastic estimator via provided that the series expansion converges rapidly (i.e., the approximation is good).In [34] it is shown that this is the case, if [r a , r b ] covers the spectral range of Qq .Furthermore, arguments are given that in this case the fluctuations derived from the random sources can be assumed to be subleading with respect to the fluctuations from the gauge field and therefore N r = 1 can be chosen.Simulations of QCD with 2 + 1 flavors of Wilson quarks using the RHMC algorithm have been studied in [36].It was found that regions in the configuration space where the strange quark Wilson-Dirac operator has negative eigenvalues are not negligible, contrary to standard assumptions.Since only the modulus of the fermion determinant is taken into account in the above procedure, it is not affected by negative signs.However, these signs have to be taken into account as an additional reweighting factor to restore the correct probability density in the path integral.Strategies for computing this sign are outlined in [28,36]. 3s explained in [36], changing the number of negative eigenvalues of the strange quark Dirac operator by one unit can only be achieved by passing through configurations with a zero eigenvalue of D s , which have zero probability in the path integral.The use of the rational approximation in the RHMC leads to a regularization of the infinite barrier between positive and negative eigenvalues.In [36] it is proposed to use this feature actively to allow the algorithm to tunnel between the two sectors and thus avoid possible ergodicity problems.
The appearance of configurations with exceptionally small eigenvalues of Qs , which are much smaller than r a , cannot be excluded in this case.The associated reweighting factor is intended to correct for such fluctuations.However, if the modulus of a nearzero eigenvalue is too small, the computation of W s by means of the Taylor expansion in eq.(3.9) fails because the series does not converge within a tolerable number of terms.

Deflated RHMC reweighting
Low-mode deflation of the RHMC reweighting factor can resolve the issue described above.We start with the factorization of the determinant of an operator Ω in eq.(2.13).In the case that the projector P is built out of eigenmodes of Ω we may use eq.(2.15).In the computation of W q , we employ and therefore eigenmodes of Qq = γ 5 Dq are eigenmodes of Ω.Using where v i are the eigenmodes of the hermitian Dirac operator Qq , we can compute from N L real-valued eigenvalues λ i of Qq and the known coefficients A and a j of the Zolotarev function.The computation of det( P Ω P ) can be done by stochastic estimation in the subspace that is orthogonal to the space of the low modes via where the issues related to the convergence of the series expansion are eliminated since all eigenvalues of P Qq are greater than r a .The only change in existing implementations of the stochastic estimation is the projection of the random sources using P and all other refinements, such as the frequency splitting, remain unchanged.

Numerical results
We have implemented the computation of deflated reweighting factors, as described in the previous sections, by modifying existing routines of the openQCD package.The computation of a small number of low eigenmodes of Qq is performed using the PRIMME package [38,39].In this section, we present the properties of the corresponding reweighting factors on a number of 2 + 1 flavor CLS ensembles.For an overview of all gauge ensembles considered in this work see table 1.They cover a wide range of lattice spacings and quark masses, including the physical ones.

Application of deflated twisted-mass reweighting
Stochastic evaluation of twisted-mass reweighting factors can become problematic on configurations where the lowest eigenmode has an exceptionally small eigenvalue, compared to the ensemble average and compared to the parameter µ 0 that has been used in the course of the simulation.This is precisely the situation where we expect low-mode deflation to be most effective.To illustrate this assumption, we examine the relative standard deviation of the stochastic estimators, configuration by configuration, and compare the pure stochastic setup with a deflated computation.
We perform the comparison on the three ensembles N451, D450 and D452 with approximate pion masses 289 MeV, 219 MeV and 156 MeV at a bare inverse gauge coupling of β = 3.46, corresponding to a relatively coarse lattice spacing of ≈ 0.075 fm.In this regime, especially towards physical quark masses, we observe an increased number of configurations with an exceptionally small spectral gap of the light quark Dirac operator, leading to small twisted-mass reweighting factors.The Dirac operator does not have a solid spectral gap here, since chiral symmetry is broken for Wilson fermions at finite lattice spacing.In the stochastic setup, 24 random sources were used to determine W l on N451 and D450 and 36 sources on D452, respectively.For the computation of the deflated reweighting factors, we employ four low modes of Ql and reduce the number of random sources by a factor of two, compared to the stochastic setup.
Fig. 1 shows the relative standard deviation of the stochastic estimator W stoch l,i and of the stochastic contribution to the deflated estimator W dfl l,i , depending on the central value of the estimator on each configuration i of the three Monte Carlo chains.We also show the prediction of this dependence for W stoch l,i based on eq.(2.8).It is clearly visible that the data follow the prediction, such that a drastic increase in the relative standard deviation is observed as the central value of the reweighting factor decreases.The relative standard deviation of the deflated estimator on the other hand is independent of the value of the reweighting factor, since only a small deviation from unity has to be estimated stochastically.Thus, by using low-mode deflation, we obtain a reduction of the relative uncertainty by several orders of magnitude in the regime of small reweighting factors.The gain is therefore maximized on those configurations where we need to determine the reweighting factor precisely because it is significantly different from one.At the same time, we also observe an overall reduction of the variance, as depicted in Figure 2, where we show the Monte Carlo history of the variance of the two estimators on the three ensembles under investigation.Furthermore, the spread of the distribution decreases and no large spikes in the variance are present, when deflation is applied.
In Figure 3 we directly compare twisted-mass reweighting factors determined in the stochastic setup with 24 sources on ensemble E250 (m π = 132 MeV, a ≈ 0.064 fm) [15] with deflated ones, using four low modes and 12 stochastic sources.For better visibility we plot the deviation of the reweighting factors from unity and show only the factors on configurations i where W dfl l,i < 0.93.Despite the approximate overall agreement, we can observe significant deviations between the two sets of estimators for a number of configurations.Note that these deviations do not affect the correctness of the reweighting procedure but that they may lead to increased fluctuations in reweighted observables.
The computational cost of calculating the deflated reweighting factor in the chosen setup is of the same order as in the purely stochastic computation.However, this cost can be significantly reduced, considering that we have never observed a case where two lowest eigenmodes have an exceptionally small eigenvalue.Therefore, any significant deviation of W l from one is due to the lowest eigenmode and deflation of a single mode is sufficient to achieve high precision.At the same time, the number of random sources may be further reduced, both compared to the stochastic setup and also compared to the setup that was chosen here, to achieve sufficient precision.
In Fig. 4 we show the dependence of the variance of the deflated stochastic estima- tor on the number of low modes N L , starting with N L = 0, which is equivalent to the traditional stochastic estimator.32 random sources are used for the stochastic estimation and for the estimation of the variance.We illustrate the behavior on two configurations of the D452 ensemble.W l on configuration D452r002n4 is close to unity.Accordingly, the estimator is well determined without deflation.However, increasing the number of low modes significantly reduces the error of the stochastic part, about an order of magnitude when increasing N L from one to ten.For D452r002n134, where the reweighting factor is significantly different from one, deflation of the lowest mode has a significant impact on the relative error.For a larger number of modes, the behavior is similar to D452r002n4.The reduction in relative error below the per mil level that can be achieved by increasing the number of low modes is not necessary in practical applications.However, to achieve a target precision with minimal cost, it may be beneficial to use more low modes and fewer stochastic sources.
The impact of the improved estimator for the reweighting factor on the expectation value of an observable depends strongly on the observable in question.If its correlation with the low modes of the Dirac operator is large, an observable will be strongly influenced by exceptionally low eigenvalues.One particularly affected observable is the pion decay constant f π which is often used to set the scale in lattice QCD computations, e.g. in [40,41], and may enter the determination of the hadronic vacuum polarization contribution to the anomalous magnetic moment of the muon via f π -rescaling [42,43].
A significant number of exceptionally small eigenvalues of the light Dirac operator can be found on the ensemble D150 (m π = 132 MeV, a ≈ 0.085 fm).The bare decay constant af π , computed as outlined in [43,44] from stochastically estimated correlation functions and reweighting factors (using four stochastic time slice sources per configuration for the correlation functions and 32 sources for the twisted-mass reweighting factor) is determined to be 0.06808(155) on this ensemble.When a deflated twisted-mass reweighting factor is used, the decay constant is evaluated to be 0.06837(74).Thus, we find a reduction of the error by a factor of two, simply by using a more precise reweighting factor.If a robust estimator for the standard deviation is used, e.g., the median absolute deviation [45], we find that the resulting estimate for the error, 5.5 • 10 −4 , is invariant under the exchange of the reweighting factors within the error of the error.In this case, the reduction in uncertainty from using deflation in the reweighting factor is due to the removal of outliers in the Monte Carlo history.
The interplay of observable and reweighting factor also depends on the estimator that is used for the observable.It has been found that correlation functions constructed from stochastic time slice sources, as used in this exemplary case, have significantly fewer outliers in the MC history than those built from point sources, such as those used in hadron structure calculations.This can be explained by the fact that potentially large local fluctuations of the eigenmodes cannot be completely eliminated by the reweighting factor, which is based on the full four-volume of the box.Correlation functions determined by low-mode averaging (LMA) [46][47][48], on the other hand, include a full volume average of the contribution of the lowest modes of the Dirac operator.Therefore, we expect that large fluctuations due to the lowest mode are exactly canceled in combination with deflated reweighting factors.
We illustrate the removal of outliers in the Monte Carlo history by using improved estimators for the case of the light quark two-point function based on pseudoscalar and axial densities (defined as C A in [43]) in figure 5.This correlation function enters the computation of af π .We observe a significant reduction in the size and the number of spikes when low-mode deflation is used for the reweighting factor and a further reduction when LMA is used to compute the correlation function.Note that a larger scale is used to depict the fluctuations for the combination of the stochastically estimated reweighting factor and correlation function in the top panel.The variance derived from the mean absolute deviation is the same in all four cases considered here.

Application of deflated RHMC reweighting
Smoothly fluctuating RHMC reweighting factors for the strange quark determinant can be observed on the vast majority of the CLS ensembles, indicating a proper rational approximation during the simulation.On some of the ensembles, deviations towards smaller values of W s can be observed on a few configurations, indicating a lowest eigenvalue λs 1 below the parameter r a that was chosen for the RHMC.On a handful of configurations with an exceptionally low strange eigenvalue, the convergence of the Taylor series in eq.(3.9) is particularly slow, so that O(100) elements are not sufficient to achieve the required precision.An investigation of the lowest eigenvalue of | Qs | together with the information about its sign from [36] in the vicinity of a problematic configuration helps to understand the behavior of the simulation.In Fig. 6 we show the Monte Carlo history of the lowest eigenvalue of Qs on the D452r002 ensemble around configuration number 190, where convergence issues are found.A departure in the region of negative eigenvalues can be observed for configurations 191 and 192.When switching from positive to negative eigenvalues, the exceptionally small eigenvalue 2.6 • 10 −4 , well below r D452 a = 5 • 10 −3 , is encountered on configuration 190. 4  This is exactly the situation that has been the subject of discussion in section 3.2.[36].The assumed spectral gap r a = 0.005 that was chosen for the simulation is depicted by dashed black lines.Such a small eigenvalue poses two problems for the computation of the reweighting factor W s : A single random source may not be sufficient to estimate W s with precision, similarly to the case of the twisted-mass reweighting factor described above.Secondly the estimation along the lines of section 3.2 may fail because the Taylor series does not converge at all.
We have found that both issues are resolved when low-mode deflation using N L = 1 is applied in the computation of W s , because any significant deviation from the ensemble average is accounted for by the exactly evaluated contribution of the lowest eigenmode.This is the case for all investigated configurations where W s could not be determined at all or not be determined accurately from stochastic estimation.

Conclusions
In the course of simulating lattice QCD with O(a) improved Wilson fermions and twistedmass reweighting it has been found that small but imprecisely determined twisted-mass reweighting factors can significantly affect the ensemble average of observables that are sensitive to the low eigenmodes of the light Dirac operator. 5We have shown that this issue can be resolved by exact deflation of a small number of low modes in the computation of the corresponding reweighting factor.The procedure can be implemented without much effort and may even reduce the overall cost of the computation.
Two other improved estimators for the twisted-mass reweighting factor have been introduced in the past.One is based on a frequency splitting of the quark determinant, as described in section 2.1, and requires a much larger amount of computer time on all configurations of an ensemble to more accurately estimate a small number of reweighting factors.The resulting precision at tolerable cost is not competitive with the observed precision of deflated reweighting factors.The strategy of [19] relies on a larger number of stochastic samples on configurations with small reweighting factors and a more robust estimator for the stochastic estimator.Low-mode deflation, on the other hand, allows to use a small number of noise sources for the estimation of the stochastic remainder on all configurations.
We have applied low-mode deflation to the computation of the twisted-mass reweighting factors on 15 CLS ensembles with pion masses between 360 MeV and 132 MeV at lattice spacings a ≥ 0.05 fm.The newly computed set of reweighting factors will be used in future work.The impact on the overall precision of phenomenologically relevant observables at the physical point depends on the observable in question and its sensitivity to the low modes of the Dirac operator.We have discussed the impact on the computation of af π on a particularly affected ensemble in section 4.1.This directly influences the attainable precision in scale setting and the computation of the hadronic vacuum polarization contribution to the muon g − 2. Indications of stabilizing behavior on the computation have also been observed in the study of I = 1 π-π scattering at the physical point [50].
As pointed out in the introduction, the precise estimation of the reweighting factor alone may not be sufficient to guarantee the cancellation of fluctuations in observables and reweighting factors to the amount that is required for high-precision calculations.Examples of such fluctuations are shown for the nucleon sigma term in [51], where the reweighting factors of this work have been employed, and for nucleon three-point functions in [52].We note in passing that nucleon three-point functions are particularly affected by exceptionally small eigenvalues because, when expressed in terms of eigenmodes, each propagator carries a factor λ−1 1 which may enhance fluctuations of the local magnitude of eigenmodes as the ones described in [53].Improved estimators for the correlation functions may be needed to remove these outliers completely.Low-mode averaging [46][47][48] is expected to have a stabilizing effect because the contributions of the lowest eigenmodes of Q can be exactly taken into account, as it is the case for the deflated reweighting factor.
While the strange quark Wilson-Dirac operator was long time assumed to have a solid spectral gap, it was found that sign changes of the quark determinant have a nonvanishing probability at finite lattice spacing [36].In the vicinity of some of these sign flips, problems in the computation of the strange quark RHMC reweighting factor occur when the lowest eigenvalue of the Dirac operator is exceptionally small.This leads to imprecise estimates or the failure of the computation.We have found that low-mode deflation can be applied here to overcome these problems by exactly computing the contribution of the lowest eigenmode to the reweighting factor.
With the results of this work, in combination with the findings of [36], the reweighting factors that are needed to restore the target distribution for modern simulations of Wilson quarks based on twisted-mass reweighting and the RHMC algorithm can be computed with high precision.This is the basis for precise predictions in the low energy regime of QCD.
gauge averages in this work has been performed using the Γ-method in the implementation of the pyerrors package [55][56][57].Plots have been generated with matplotlib [58].

A The spectral gap of the light quark Wilson-Dirac operator
The distribution of the spectral gap of the Dirac-Wilson operator has been investigated in some detail in the two-flavor theory [5,6], where the first of the two references focused on the theory without O(a) improvement.Information on the behavior of exponentiatedclover Wilson fermions in the 2 + 1 flavor theory is provided in [11]. 6In this appendix, we will provide some insight in the distribution of the lowest eigenvalue λ1 of the even-odd preconditioned hermitian light-quark Wilson-clover Dirac operator Ql on the 2 + 1 flavor CLS ensembles that have been included in this study.
We define spectral gap of Ql as as the absolute value of its smallest eigenvalue λ1 .Because the tail of the distribution of the spectral gap might be poorly sampled along a Monte Carlo chain, we do not rely on the standard deviation to estimate the width of the distribution.Instead, the median and width σ of the reweighted eigenvalue distributions are computed as proposed in appendix D of [11].V /a, versus Φ 2 (from [43]) as proxy for the light quark mass.The uncertainties of the widths have been estimated using a bootstrap procedure.The vertical line denotes physical light quark masses.
In Figure 7 we plot the product of the width σ and the square root of the four volume V versus Φ 2 = 8t 0 m 2 π .The latter may be used as proxy for the light quark mass to leading order in chiral perturbation theory [40].In the two-flavor theory without O(a) improvement, the product σ √ V /a was found to be approximately constant, where σ was the width of the distribution of the spectral gap of Q.In the improved two-flavor theory [6], a dependence of σ √ V /a on the quark mass was observed.This is also the case in our study, where σ√ V /a seems to decrease towards smaller light quark masses.Figure 8 shows the dependence of the mean of the distribution of λ1 , renormalized with 1/Z P from [59], on the light quark proxy Φ 2 .In formulations of lattice QCD that preserve chiral symmetry, the gap of Q l is bounded from below by the bare current-quark mass m 12 [5].While we cannot expect this relation to hold for Wilson quarks, the data in Fig. 8 seems to be approximately proportional to Φ 2 .Overview of ensembles used in this work.The lattice spacing a has been computed in [19].Pion and kaon masses have been determined in the context of [43].L denotes the spatial extent of the box.

Figure 1 :
Figure 1: Dependence of the relative standard deviation of the stochastic estimator W stoch l and the deflated estimator W dfl l on the value of the estimator on the ensembles N451, D450 and D452 at β = 3.46.The dotted line shows the prediction for the standard deviation of W stoch l based on eq.(2.8).Four low modes have been used in the deflated estimator.

Figure 2 :
Figure 2: Monte Carlo history of the stochastic variances of the estimator W stoch l and of the deflated estimator W dfl l on the ensembles N451, D450 and D452 at β = 3.46.

Figure 3 :
Figure 3: Deviation from one for the two sets of twisted-mass reweighting factors on ensemble E250r001 at close-to-physical pion mass.For clarity we show only factors smaller than 0.93.The error bars are derived from the distribution of the stochastic estimates on each configuration.The uncertainties on W dfl l are smaller than the symbol size.

Figure 4 :
Figure 4:Relative standard deviation of the stochastic estimator for W l depending on the number of low modes N L that is projected out, starting with N L = 0, i.e., no deflation, for D452r002n4 and D452r002n134.Lines are drawn to guide the eye.

Figure 5 :
Figure 5: Monte Carlo history of the light quark two-point correlation function based on pseudoscalar and axial density at source-sink separation T /4 on ensemble D150 for four combinations of reweighting factor and correlation function.Note the different scale in the first panel.

1 D452r002Figure 6 :
Figure 6: Part of the Monte Carlo history of the lowest eigenvalue of Qs with the sign according to[36].The assumed spectral gap r a = 0.005 that was chosen for the simulation is depicted by dashed black lines.

Figure 7 :
Figure 7: Width of the reweighted gap distributions, scaled by √V /a, versus Φ 2 (from[43]) as proxy for the light quark mass.The uncertainties of the widths have been estimated using a bootstrap procedure.The vertical line denotes physical light quark masses.

Figure 8 :
Figure 8: Reweighted spectral gap λ1 /Z P versus φ 2 .Note that the symbols of the two ensembles D150 and E250, both approximately at physical light quark mass, overlap.For the majority of data points, the error bars are smaller than the symbol sizes. id